{
image_tmp1[ny-1][j] = image[ny-1][j+1]; image_tmp2[0][j+1] = image[0][j]; }
image_tmp1[ny-1][nx-1] = image[ny-1][nx-1]; image_tmp2[0][0] = image[0][0];
//计算Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]); for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++) {
image_dp[i][j] = image_tmp1[i][j] + image_tmp2[i][j]; } }
//////////////////////////////////////////////////////////////////////////
//计算I([1 1:ny-1],[2:nx nx])和I([2:ny ny],[1 1:nx-1]) //公共部分分别是矩阵左下角，右上角的ny-1行和nx-1列 for (int i = 0; i < ny-1; i++) {
for (int j = 0; j < nx-1; j++) {
image_tmp1[i+1][j] = image[i][j+1]; image_tmp2[i][j+1] = image[i+1][j]; } }
for (int i = 0; i < ny-1; i++) {
image_tmp1[i+1][nx-1] = image[i][nx-1]; image_tmp2[i][0] = image[i+1][0]; }
for (int j = 0; j < nx-1; j++) {
image_tmp1[0][j] = image[0][j+1]; image_tmp2[ny-1][j+1] = image[ny-1][j]; }
image_tmp1[0][nx-1] = image[0][nx-1]; image_tmp2[ny-1][0] = image[ny-1][0];
//计算Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]); for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++)
{
image_dm[i][j] = image_tmp1[i][j] + image_tmp2[i][j]; } }
//////////////////////////////////////////////////////////////////////////
//计算I_xy = (Dp-Dm)/4; for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++) {
image_xy[i][j] = (image_dp[i][j] - image_dm[i][j])/4; } }
//////////////////////////////////////////////////////////////////////////
//计算过程：
//计算
Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2) 和 Den = (ep2+I_x.^2+I_y.^2).^(3/2);
for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++) {
image_num[i][j] = image_xx[i][j]*(image_y[i][j]*image_y[i][j] + ep2)
- 2*image_x[i][j]*image_y[i][j]*image_xy[i][j] + image_yy[i][j]*(image_x[i][j]*image_x[i][j] + ep2);
image_den[i][j] = pow((image_x[i][j]*image_x[i][j] + image_y[i][j]*image_y[i][j] + ep2), 1.5); } }
//计算
I: I_t = Num./Den + lam.*(I0-I+C); I=I+dt*I_t; %% evolve image by dt for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++) {
image[i][j] += dt*(image_num[i][j]/image_den[i][j] + lam*(im
age0[i][j] - image[i][j])); } } }
//迭代结束
//////////////////////////////////////////////////////////////////////////
//赋值图像 BYTE tmp;
for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++) {
tmp = (BYTE)image[i][j]; tmp = max(0, min(tmp, 255));
my_image->SetPixelIndex(j, ny-i-1, tmp); } }
//////////////////////////////////////////////////////////////////////////
//删除内存
deleteDoubleMatrix(image_x, nx, ny); deleteDoubleMatrix(image_y, nx, ny); deleteDoubleMatrix(image_xx, nx, ny); deleteDoubleMatrix(image_yy, nx, ny); deleteDoubleMatrix(image_tmp1, nx, ny); deleteDoubleMatrix(image_tmp2, nx, ny); deleteDoubleMatrix(image_dp, nx, ny); deleteDoubleMatrix(image_dm, nx, ny); deleteDoubleMatrix(image_xy, nx, ny); deleteDoubleMatrix(image_num, nx, ny); deleteDoubleMatrix(image_den, nx, ny); deleteDoubleMatrix(image0, nx, ny); deleteDoubleMatrix(image, nx, ny);
return true; }
////////////////////////////////////////////////////////////////////////// //开辟二维数组函数
double** MyCxImage::newDoubleMatrix(int nx, int ny) {
double** matrix = new double*[ny];
for(int i = 0; i < ny; i++) {
matrix[i] = new double[nx]; }
if(!matrix) return NULL; return matrix; }
//清除二维数组内存函数
bool MyCxImage::deleteDoubleMatrix(double** matrix, int nx, int ny) {
if (!matrix) {
return true; }
for (int i = 0; i < ny; i++) {
if (matrix[i]) {
delete[] matrix[i]; } }
delete[] matrix;
return true; }
//////////////////////////////////////////////////////////////////////////
这个代码单独显然是无法运行的，因为还要涉及底层的图像处理的类库，图像的读取显示我用了CxIamge类，而程序界面我是用的MFC的框架。不过代码基本一直都是在做矩阵运算，如果要是能有一个比较好的矩阵运算类库的话，代码会简介许多，效率也会高一些。总体上C++代码还是要比Matlab效率高许多的。
关于变分法的算法原理和基本思想，我这两天再读一些论文在做总结。。