TV去噪主页：http://visl.technion.ac.il/~gilboa/PDE-filt/tv_denoising.html
可以下载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 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; 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;
另外有博主谢了C++的代码：经典的变分法图像去噪的C++实现 ，不过个人感觉写的不太简洁。稍稍修改了一下，新版代码出炉：
//TV去噪函数
Mat TVDenoising(Mat img, int iter)
{
int ep = 1;
int nx=img.cols;
int ny = img.rows;
double dt = 0.25f;
double lam = 0.0;
int ep2 = ep*ep;
double** image = newDoubleMatrix(nx, ny);
double** image0 = newDoubleMatrix(nx, ny);
for(int i=0;i<ny;i++){
uchar* p=img.ptr<uchar>(i);
for(int j=0;j<nx;j++){
image0[i][j]=image[i][j]=(double)p[j];
}
}
//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_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++){
//for (int i = 0; i < ny; i++){
// for (int j = 0; j < nx; j++){
// //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;
// int tmp_i1=(i+1)<ny ? (i+1) :(ny-1);
// int tmp_j1=(j+1)<nx ? (j+1): (nx-1);
// int tmp_i2=(i-1) > -1 ? (i-1) : 0;
// int tmp_j2=(j-1) > -1 ? (j-1) : 0;
// image_x[i][j] = (image[i][tmp_j1] - image[i][tmp_j2])/2;
// image_y[i][j]= (image[tmp_i1][j]-image[tmp_i2][j])/2;
// image_xx[i][j] = image[i][tmp_j1] + image[i][tmp_j2]- image[i][j]*2;
// image_yy[i][j]= image[tmp_i1][j]+image[tmp_i2][j] - image[i][j]*2;
// image_dp[i][j]=image[tmp_i1][tmp_j1]+image[tmp_i2][tmp_j2];
// image_dm[i][j]=image[tmp_i2][tmp_j1]+image[tmp_i1][tmp_j2];
// image_xy[i][j] = (image_dp[i][j] - image_dm[i][j])/4;
// 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);
// image[i][j] += dt*(image_num[i][j]/image_den[i][j] + lam*(image0[i][j] - image[i][j]));
// }
//}
for (int i = 0; i < ny; i++){
for (int j = 0; j < nx; j++){
int tmp_i1=(i+1)<ny ? (i+1) :(ny-1);
int tmp_j1=(j+1)<nx ? (j+1): (nx-1);
int tmp_i2=(i-1) > -1 ? (i-1) : 0;
int tmp_j2=(j-1) > -1 ? (j-1) : 0;
double tmp_x = (image[i][tmp_j1] - image[i][tmp_j2])/2; //I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2;
double tmp_y= (image[tmp_i1][j]-image[tmp_i2][j])/2; //I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
double tmp_xx = image[i][tmp_j1] + image[i][tmp_j2]- image[i][j]*2; //I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
double tmp_yy= image[tmp_i1][j]+image[tmp_i2][j] - image[i][j]*2; //I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
double tmp_dp=image[tmp_i1][tmp_j1]+image[tmp_i2][tmp_j2]; //Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
double tmp_dm=image[tmp_i2][tmp_j1]+image[tmp_i1][tmp_j2]; //Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
double tmp_xy = (tmp_dp - tmp_dm)/4; //I_xy = (Dp-Dm)/4;
double tmp_num = tmp_xx*(tmp_y*tmp_y + ep2)
- 2*tmp_x*tmp_y*tmp_xy +tmp_yy*(tmp_x*tmp_x + ep2); //Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
double tmp_den= pow((tmp_x*tmp_x + tmp_y*tmp_y + ep2), 1.5); //Den = (ep2+I_x.^2+I_y.^2).^(3/2);
image[i][j] += dt*(tmp_num/tmp_den+ lam*(image0[i][j] - image[i][j]));
}
}
}
Mat new_img;
img.copyTo(new_img);
for(int i=0;i<img.rows;i++){
uchar* p=img.ptr<uchar>(i);
uchar* np=new_img.ptr<uchar>(i);
for(int j=0;j<img.cols;j++){
int tmp=(int)image[i][j];
tmp=max(0,min(tmp,255));
np[j]=(uchar)(tmp);
}
}
//////////////////////////////////////////////////////////////////////////
//删除内存
//deleteDoubleMatrix(image_x, nx, ny);
//deleteDoubleMatrix(image_y, nx, ny);
//deleteDoubleMatrix(image_xx, nx, ny);
//deleteDoubleMatrix(image_yy, 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);
//imshow("Image",img);
//imshow("Denosing",new_img);
return new_img;
}
（转载请注明作者和出处：http://blog.csdn.net/xiaowei_cqu 未经允许请勿用于商业用途）