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

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未经允许请勿用于商业用途）