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