经典的变分法图像去噪的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 (doublearray gray level1-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]

%%(defaultvalues arein[])

%%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;



fori=1:iter,%%doiterations

%estimate derivatives

I_x=(I(:,[2:nx nx])-I(:,[11:nx-1]))/2;

I_y=(I([2:ny ny],:)-I([11:ny-1],:))/2;

I_xx=I(:,[2:nx nx])+I(:,[11:nx-1])-2I;

I_yy=I([2:ny ny],:)+I([11:ny-1],:)-2I;

Dp=I([2:ny ny],[2:nx nx])+I([11:ny-1],[11:nx-1]);

Dm=I([11:ny-1],[2:nx nx])+I([2:ny ny],[11:nx-1]);

I_xy=(Dp-Dm)/4;

%compute flow

Num=I_xx.(ep2+I_y.^2)-2I_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+dtI_t;%%evolve image by dt

end%fori



%%returnimage

%J=IImean/mean(mean(I));%normalize to original mean

J=I;



figure(1); imshow(uint8(I0)); title('Noisy image');

%denoise image byusingtvforsome iterations

figure(2); imshow(uint8(J)); title('Denoised image');


另外我在我的图像处理框架程序里实现了这个最经典版本的TV去噪算法，核心代码如下：


//TV去噪函数

boolMyCxImage::TVDenoising(intiter/= 80/)

{

if(my_image==NULL)returnfalse;

if(!my_image->IsValid())returnfalse;

//算法目前不支持彩色图像，所以对于彩图，先要转换成灰度图。

if(!my_image->IsGrayScale())

{

my_image->GrayScale();

//return false;

}



//基本参数，这里由于设置矩阵C为0矩阵，不参与运算，所以就忽略之

intep=1, nx=width, ny=height;

doubledt=(double)ep/5.0f, lam=0.0;

intep2=epep;



doubleimage=newDoubleMatrix(nx, ny);

doubleimage0=newDoubleMatrix(nx, ny);

//注意一点是CxImage里面图像存储的坐标原点是左下角，Matlab里面图像时左上角原点

for(inti=0; i<ny; i++)

{

for(intj=0; j<nx; j++)

{

image0[i][j]=image[i][j]=my_image->GetPixelIndex(j, ny-i-1);

}

}



doubleimage_x=newDoubleMatrix(nx, ny);//I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2;

doubleimage_xx=newDoubleMatrix(nx, ny);//I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2I;

doubleimage_y=newDoubleMatrix(nx, ny);//I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;

doubleimage_yy=newDoubleMatrix(nx, ny);//I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2I;

doubleimage_tmp1=newDoubleMatrix(nx, ny);

doubleimage_tmp2=newDoubleMatrix(nx, ny);



doubleimage_dp=newDoubleMatrix(nx, ny);//Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1



doubleimage_dm=newDoubleMatrix(nx, ny);//Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);



doubleimage_xy=newDoubleMatrix(nx, ny);//I_xy = (Dp-Dm)/4;



doubleimage_num=newDoubleMatrix(nx, ny);//Num = I_xx.(ep2+I_y.^2)-2I_x.I_y.I_xy+I_yy.(ep2+I_x.^2);

doubleimage_den=newDoubleMatrix(nx, ny);//Den = (ep2+I_x.^2+I_y.^2).^(3/2);



//////////////////////////////////////////////////////////////////////////

//对image进行迭代iter次

iter=80;

for(intt=1; t<=iter; t++)

{

//进度条

my_image->SetProgress((long)100t/iter);

if(my_image->GetEscape())

break;

//////////////////////////////////////////////////////////////////////////

//计算I(:,[2:nx nx])和I(:,[1 1:nx-1])

//公共部分2到nx-1列

for(inti=0; i<ny; i++)

{

for(intj=0; j<nx-1; j++)

{

image_tmp1[i][j]=image[i][j+1];

image_tmp2[i][j+1]=image[i][j];

}

}

for(inti=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])-2I;

for(inti=0; i<ny; i++)

{

for(intj=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])-2image[i][j];

}

}



//////////////////////////////////////////////////////////////////////////

//计算I([2:ny ny],:)和I([1 1:ny-1],:)

//公共部分2到ny-1行

for(inti=0; i<ny-1; i++)

{

for(intj=0; j<nx; j++)

{

image_tmp1[i][j]=image[i+1][j];

image_tmp2[i+1][j]=image[i][j];

}

}

for(intj=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],:)-2I;

for(inti=0; i<ny; i++)

{

for(intj=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])-2image[i][j];

}

}



//////////////////////////////////////////////////////////////////////////

//计算I([2:ny ny],[2:nx nx])和I([1 1:ny-1],[1 1:nx-1])

//公共部分分别是矩阵右下角，左上角的ny-1行和nx-1列

for(inti=0; i<ny-1; i++)

{

for(intj=0; j<nx-1; j++)

{

image_tmp1[i][j]=image[i+1][j+1];

image_tmp2[i+1][j+1]=image[i][j];

}

}

for(inti=0; i<ny-1; i++)

{

image_tmp1[i][nx-1]=image[i+1][nx-1];

image_tmp2[i+1][0]=image[i][0];

}

for(intj=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(inti=0; i<ny; i++)

{

for(intj=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(inti=0; i<ny-1; i++)

{

for(intj=0; j<nx-1; j++)

{

image_tmp1[i+1][j]=image[i][j+1];

image_tmp2[i][j+1]=image[i+1][j];

}

}

for(inti=0; i<ny-1; i++)

{

image_tmp1[i+1][nx-1]=image[i][nx-1];

image_tmp2[i][0]=image[i+1][0];

}

for(intj=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(inti=0; i<ny; i++)

{

for(intj=0; j<nx; j++)

{

image_dm[i][j]=image_tmp1[i][j]+image_tmp2[i][j];

}

}



//////////////////////////////////////////////////////////////////////////

//计算I_xy = (Dp-Dm)/4;

for(inti=0; i<ny; i++)

{

for(intj=0; j<nx; j++)

{

image_xy[i][j]=(image_dp[i][j]-image_dm[i][j])/4;

}

}



//////////////////////////////////////////////////////////////////////////

//计算过程：



//计算Num = I_xx.(ep2+I_y.^2)-2I_x.I_y.I_xy+I_yy.(ep2+I_x.^2) 和 Den = (ep2+I_x.^2+I_y.^2).^(3/2);

for(inti=0; i<ny; i++)

{

for(intj=0; j<nx; j++)

{

image_num[i][j]=image_xx[i][j](image_y[i][j]image_y[i][j]+ep2)

-2image_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+dtI_t; %% evolve image by dt

for(inti=0; i<ny; i++)

{

for(intj=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(inti=0; i<ny; i++)

{

for(intj=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);



returntrue;

}

//////////////////////////////////////////////////////////////////////////

//开辟二维数组函数

doubleMyCxImage::newDoubleMatrix(intnx,intny)

{

doublematrix=newdouble*[ny];



for(inti=0; i<ny; i++)

{

matrix[i]=newdouble[nx];

}

if(!matrix)

returnNULL;

return

matrix;

}

//清除二维数组内存函数

boolMyCxImage::deleteDoubleMatrix(doublematrix,intnx,intny)

{

if(!matrix)

{

returntrue;

}

for(inti=0; i<ny; i++)

{

if(matrix[i])

{

delete[] matrix[i];

}

}

delete[] matrix;



returntrue;

}

//////////////////////////////////////////////////////////////////////////


这个代码单独显然是无法运行的，因为还要涉及底层的图像处理的类库，图像的读取显示我用了CxIamge类，而程序界面我是用的MFC的框架。不过代码基本一直都是在做矩阵运算，如果要是能有一个比较好的矩阵运算类库的话，代码会简介许多，效率也会高一些。总体上C++代码还是要比Matlab效率高许多的。

关于变分法的算法原理和基本思想，我这两天再读一些论文在做总结。。

Email：lichao@icst.pku.edu.cn


原文链接: https://www.cnblogs.com/CCBB/archive/2010/12/29/1920884.html

欢迎关注

微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍