经典的变分法图像去噪的C++实现
由于这学期的图像处理课程的大作业需要写一个图像处理程序，不能使用古典的线性滤波，或者基于频域（小波）或者基于统计之类的方法。只能用老师讲过的一些方法，诸如变分，PDE，微分几何等。。感觉上简单的变分法稍微要好实现一些，就打算基于最早的TV图像去噪模型，做一个VC的实现。但是找遍了网上也没有TV去噪的C++源码，与之只好自己动手写了。
关于变分法和泛函分析的一些基础原理今天就先不多说了，TV图像去噪经典论文：《Nonlinear Total Variation based noise removal algorithms》Google上可以搜得到。
关于Matlab的程序实现，有一个经典的主页： http://visl.technion.ac.il/~gilboa/PDE-filt/tv_denoising.html
下面是一个Matlab代码实现：复制到记事本用matlab打开就可以运行，要注意图像的名称和路径要对应。如果只是想学学算法思路或者看看处理效果的话，只需要Matlab的代码就行了。
function J=tv(I,iter,dt,ep,lam,I0,C)
%% Private function: tv (by Guy Gilboa).
%% Total Variation denoising.
%% Example: J=tv(I,iter,dt,ep,lam,I0)
%% Input: I - image (double array gray level 1-256),
%% iter - num of iterations,
%% dt - time step [0.2],
%% ep - epsilon (of gradient regularization) [1],
%% lam - fidelity term lambda [0],
%% I0 - input (noisy) image [I0=I]
%% (default values are in [])
%% Output: evolved image
clc
clear
I=imread('grids.bmp'); % load image
I = double(I);
if ~exist('ep')
ep=1;
end
if ~exist('dt')
dt=ep/5; % dt below the CFL bound
end
if ~exist('lam')
lam=0;
end
if ~exist('I0')
I0=I;
end
if ~exist('C')
C=0;
end
[ny,nx]=size(I); ep2=ep^2;
% params
iter=80;
for i=1:iter, %% do iterations
% estimate derivatives
I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2;
I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
I_xy = (Dp-Dm)/4;
% compute flow
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);
I_t = Num./Den + lam.*(I0-I+C);
I=I+dt*I_t; %% evolve image by dt
end % for i
%% return image
%J=I*Imean/mean(mean(I)); % normalize to original mean
J=I;
figure(1); imshow(uint8(I0)); title('Noisy image');
% denoise image by using tv for some iterations
figure(2); imshow(uint8(J)); title('Denoised image');
%% Private function: tv (by Guy Gilboa).
%% Total Variation denoising.
%% Example: J=tv(I,iter,dt,ep,lam,I0)
%% Input: I - image (double array gray level 1-256),
%% iter - num of iterations,
%% dt - time step [0.2],
%% ep - epsilon (of gradient regularization) [1],
%% lam - fidelity term lambda [0],
%% I0 - input (noisy) image [I0=I]
%% (default values are in [])
%% Output: evolved image
clc
clear
I=imread('grids.bmp'); % load image
I = double(I);
if ~exist('ep')
ep=1;
end
if ~exist('dt')
dt=ep/5; % dt below the CFL bound
end
if ~exist('lam')
lam=0;
end
if ~exist('I0')
I0=I;
end
if ~exist('C')
C=0;
end
[ny,nx]=size(I); ep2=ep^2;
% params
iter=80;
for i=1:iter, %% do iterations
% estimate derivatives
I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2;
I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
I_xy = (Dp-Dm)/4;
% compute flow
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);
I_t = Num./Den + lam.*(I0-I+C);
I=I+dt*I_t; %% evolve image by dt
end % for i
%% return image
%J=I*Imean/mean(mean(I)); % normalize to original mean
J=I;
figure(1); imshow(uint8(I0)); title('Noisy image');
% denoise image by using tv for some iterations
figure(2); imshow(uint8(J)); title('Denoised image');
另外我在我的图像处理框架程序里实现了这个最经典版本的TV去噪算法，核心代码如下：
//TV去噪函数
bool MyCxImage::TVDenoising(int iter /* = 80 */)
{
if(my_image == NULL) return false;
if(!my_image->IsValid()) return false;
//算法目前不支持彩色图像，所以对于彩图，先要转换成灰度图。
if(!my_image->IsGrayScale())
{
my_image->GrayScale();
//return false;
}
//基本参数，这里由于设置矩阵C为0矩阵，不参与运算，所以就忽略之
int ep = 1, nx = width, ny = height;
double dt = (double)ep/5.0f, lam = 0.0;
int ep2 = ep*ep;
double** image = newDoubleMatrix(nx, ny);
double** image0 = newDoubleMatrix(nx, ny);
//注意一点是CxImage里面图像存储的坐标原点是左下角，Matlab里面图像时左上角原点
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx; j++)
{
image0[i][j] = image[i][j] = my_image->GetPixelIndex(j, ny-i-1);
}
}
double** image_x = newDoubleMatrix(nx, ny); //I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2;
double** image_xx = newDoubleMatrix(nx, ny); //I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
double** image_y = newDoubleMatrix(nx, ny); //I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
double** image_yy = newDoubleMatrix(nx, ny); //I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
double** image_tmp1 = newDoubleMatrix(nx, ny);
double** image_tmp2 = newDoubleMatrix(nx, ny);
double** image_dp = newDoubleMatrix(nx, ny); //Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1
double** image_dm = newDoubleMatrix(nx, ny); //Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
double** image_xy = newDoubleMatrix(nx, ny); //I_xy = (Dp-Dm)/4;
double** image_num = newDoubleMatrix(nx, ny); //Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
double** image_den = newDoubleMatrix(nx, ny); //Den = (ep2+I_x.^2+I_y.^2).^(3/2);
//////////////////////////////////////////////////////////////////////////
//对image进行迭代iter次
iter = 80;
for (int t = 1; t <= iter; t++)
{
//进度条
my_image->SetProgress((long)100*t/iter);
if (my_image->GetEscape())
break;
//////////////////////////////////////////////////////////////////////////
//计算I(:,[2:nx nx])和I(:,[1 1:nx-1])
//公共部分2到nx-1列
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx-1; j++)
{
image_tmp1[i][j] = image[i][j+1];
image_tmp2[i][j+1] = image[i][j];
}
}
for (int i = 0; i < ny; i++)
{
image_tmp1[i][nx-1] = image[i][nx-1];
image_tmp2[i][0] = image[i][0];
}
//计算I_x, I_xx
// I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2
//I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx; j++)
{
image_x[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/2;
image_xx[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - 2*image[i][j];
}
}
//////////////////////////////////////////////////////////////////////////
//计算I([2:ny ny],:)和I([1 1:ny-1],:)
//公共部分2到ny-1行
for (int i = 0; i < ny-1; i++)
{
for (int j = 0; j < nx; j++)
{
image_tmp1[i][j] = image[i+1][j];
image_tmp2[i+1][j] = image[i][j];
}
}
for (int j = 0; j < nx; j++)
{
image_tmp1[ny-1][j] = image[ny-1][j];
image_tmp2[0][j] = image[0][j];
}
//计算I_xx, I_yy
// I_y = I([2:ny ny],:)-I([1 1:ny-1],:)
//I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx; j++)
{
image_y[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/2;
image_yy[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - 2*image[i][j];
}
}
//////////////////////////////////////////////////////////////////////////
//计算I([2:ny ny],[2:nx nx])和I([1 1:ny-1],[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][j] = image[i+1][j+1];
image_tmp2[i+1][j+1] = image[i][j];
}
}
for (int i = 0; i < ny-1; i++)
{
image_tmp1[i][nx-1] = image[i+1][nx-1];
image_tmp2[i+1][0] = image[i][0];
}
for (int j = 0; j < nx-1; j++)
{
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*(image0[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;
}
//////////////////////////////////////////////////////////////////////////
bool MyCxImage::TVDenoising(int iter /* = 80 */)
{
if(my_image == NULL) return false;
if(!my_image->IsValid()) return false;
//算法目前不支持彩色图像，所以对于彩图，先要转换成灰度图。
if(!my_image->IsGrayScale())
{
my_image->GrayScale();
//return false;
}
//基本参数，这里由于设置矩阵C为0矩阵，不参与运算，所以就忽略之
int ep = 1, nx = width, ny = height;
double dt = (double)ep/5.0f, lam = 0.0;
int ep2 = ep*ep;
double** image = newDoubleMatrix(nx, ny);
double** image0 = newDoubleMatrix(nx, ny);
//注意一点是CxImage里面图像存储的坐标原点是左下角，Matlab里面图像时左上角原点
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx; j++)
{
image0[i][j] = image[i][j] = my_image->GetPixelIndex(j, ny-i-1);
}
}
double** image_x = newDoubleMatrix(nx, ny); //I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2;
double** image_xx = newDoubleMatrix(nx, ny); //I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
double** image_y = newDoubleMatrix(nx, ny); //I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
double** image_yy = newDoubleMatrix(nx, ny); //I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
double** image_tmp1 = newDoubleMatrix(nx, ny);
double** image_tmp2 = newDoubleMatrix(nx, ny);
double** image_dp = newDoubleMatrix(nx, ny); //Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1
double** image_dm = newDoubleMatrix(nx, ny); //Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
double** image_xy = newDoubleMatrix(nx, ny); //I_xy = (Dp-Dm)/4;
double** image_num = newDoubleMatrix(nx, ny); //Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
double** image_den = newDoubleMatrix(nx, ny); //Den = (ep2+I_x.^2+I_y.^2).^(3/2);
//////////////////////////////////////////////////////////////////////////
//对image进行迭代iter次
iter = 80;
for (int t = 1; t <= iter; t++)
{
//进度条
my_image->SetProgress((long)100*t/iter);
if (my_image->GetEscape())
break;
//////////////////////////////////////////////////////////////////////////
//计算I(:,[2:nx nx])和I(:,[1 1:nx-1])
//公共部分2到nx-1列
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx-1; j++)
{
image_tmp1[i][j] = image[i][j+1];
image_tmp2[i][j+1] = image[i][j];
}
}
for (int i = 0; i < ny; i++)
{
image_tmp1[i][nx-1] = image[i][nx-1];
image_tmp2[i][0] = image[i][0];
}
//计算I_x, I_xx
// I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2
//I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx; j++)
{
image_x[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/2;
image_xx[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - 2*image[i][j];
}
}
//////////////////////////////////////////////////////////////////////////
//计算I([2:ny ny],:)和I([1 1:ny-1],:)
//公共部分2到ny-1行
for (int i = 0; i < ny-1; i++)
{
for (int j = 0; j < nx; j++)
{
image_tmp1[i][j] = image[i+1][j];
image_tmp2[i+1][j] = image[i][j];
}
}
for (int j = 0; j < nx; j++)
{
image_tmp1[ny-1][j] = image[ny-1][j];
image_tmp2[0][j] = image[0][j];
}
//计算I_xx, I_yy
// I_y = I([2:ny ny],:)-I([1 1:ny-1],:)
//I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
for (int i = 0; i < ny; i++)
{
for (int j = 0; j < nx; j++)
{
image_y[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/2;
image_yy[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - 2*image[i][j];
}
}
//////////////////////////////////////////////////////////////////////////
//计算I([2:ny ny],[2:nx nx])和I([1 1:ny-1],[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][j] = image[i+1][j+1];
image_tmp2[i+1][j+1] = image[i][j];
}
}
for (int i = 0; i < ny-1; i++)
{
image_tmp1[i][nx-1] = image[i+1][nx-1];
image_tmp2[i+1][0] = image[i][0];
}
for (int j = 0; j < nx-1; j++)
{
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*(image0[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效率高许多的。
关于变分法的算法原理和基本思想，我这两天再读一些论文在做总结。。
Email：lichao@icst.pku.edu.cn
