# 【图像处理】全分发TV图像去噪

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;



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