cgm_2d2.m

一个简单的图象去噪的算法

function [u,iter,res] = cgm_2d2(z, alpha, beta, nit, tol)
%[u,w,it,res] = cgm_2d(z, alpha, beta, nit, tol, iout);
%Solves the unconstrained formulation of TV denoising
%z: noisy image,                u: initial guess
%nit: #nwt Newton iterations    tol: CG tolerance
%beta: beta for |nabla u|^2+beta
%alpha: for  J = alpha*TV(u) + (u-z)^2/2
n=length(z); u=[];m=1;
N=n^2;
shuf=reshape([1:N ; N+1:2*N], 2*N,1);%矩阵按列拉直成为一个列向量
ix=1:2:2*N;
iy=2:2:2*N;

e1=ones(n, 1);
A1=spdiags([e1 -e1], [-1 0], n, n);
% A1(n, n)=0;
A=[ kron(speye(n), A1) kron(A1, speye(n))]; %A:N*2N;[Dx Dy]
A=A(:, shuf);  %A an zhao shuf de yuan su chong xin pai lie,why?  %A:div  A':nabla

I2N=speye(2*N);
IN=speye(N);

eta=0.9;

z=reshape(z, N, 1 );%矩阵按列拉直成为一个列向量
u=z;  %initial value
%w=zeros(2*N, 1);
w=ones(2*N,1)/sqrt(2*N);  % initial value %w=nabla u/sqrt(|nabla u|^2+beita):xiang liang
normrh_old=1;
fprintf('  %-9s %-9s %-12s %-10s %-8s\n',...
    '|rh|', 'rel|rh|', '  cgres', '1-ww', 's_w');
res=0; du=ones(n,n);
for iter=1:nit
  grad_u=A'*u;
  norm_grad_u=kron(sqrt(grad_u(ix).^2+grad_u(iy).^2+beta), [1;1]);%2N*1 sqrt(|nabla u|^2+beita) %grad_u(ix):Dx*u
%   rh=-alpha*A*(grad_u./norm_grad_u)+u-z;
    rh=alpha*A*(grad_u./norm_grad_u)+u-z; % gradient of J:rh=-alpha*A*(grad_u./norm_grad_u)+u-z;?
  normg=norm(rh);
  if ( iter == 1 )
    normg0=normg;
  end
  
  Grad_u=[spdiags(grad_u(ix), 0, N, N);
      spdiags(grad_u(iy), 0, N, N)];  %2N*N
  Grad_u=Grad_u(shuf, :);  %true gradient of u:[Dx u;Dy u]
  
  W=[spdiags(w(ix), 0, N, N);
      spdiags(w(iy), 0, N, N)];
  W=W(shuf, :);
  D=spdiags(1./norm_grad_u, 0, 2*N, 2*N);
  F=I2N-0.5*D*(Grad_u*W'+W*Grad_u');
%   M=-alpha*A*D*F*A' + IN;
    M=alpha*A*D*F*A'+IN;
  %R=cholinc(M,10^(-4));
  %du=pcg(M,-rh,min(0.1,.9*(norm(rh)/normrh_old)),500,R',R);
  du=-M\rh;
  dw=grad_u./norm_grad_u - w + D*F*A'*du;
	jj = j2(u,z,alpha);
        res(iter)=jj;
    if normg < tol*normg0
        fprintf('% .2e % .2e\n', normg, normg/normrh_old);
        break;
    end

  ww=w(ix).^2+w(iy).^2;
  wdw=w(ix).*dw(ix)+w(iy).*dw(iy);
  dwdw=dw(ix).^2+dw(iy).^2;
  idx_dw= dwdw > 0; % elements are 0,1
  ww=ww(idx_dw);  % dimention(ww)=the number of 1 in idx_dw 
  wdw=wdw(idx_dw);
  dwdw=dwdw(idx_dw);
  s_w=min((-wdw+sqrt(wdw.^2+dwdw.*(1-ww)))./dwdw);
  s_w=eta*min(s_w, 1);
  
  w = w + s_w*dw;
  u = u + du;
%  while(m<3) 
%   grad_u=A'*u;
%   norm_grad_u=kron(sqrt(grad_u(ix).^2+grad_u(iy).^2+beta), [1;1]);%2N*1 sqrt(|nabla u|^2+beita) %grad_u(ix):Dx*u
% %   rh=-alpha*A*(grad_u./norm_grad_u)+u-z;
%     rh=alpha*A*(grad_u./norm_grad_u)+u-z; % gradient of J:rh=-alpha*A*(grad_u./norm_grad_u)+u-z;?
%   normg=norm(rh);
%   if ( iter == 1 )
%     normg0=normg;
%   end
%   du=-M\rh;
%   dw=grad_u./norm_grad_u - w + D*F*A'*du;
% 	jj = j2(u,z,alpha);
%         res(iter)=jj;
%     if normg < tol*normg0 
%         fprintf('% .2e % .2e\n', normg, normg/normrh_old);
%         break;
%     end
% 
%   ww=w(ix).^2+w(iy).^2;
%   wdw=w(ix).*dw(ix)+w(iy).*dw(iy);
%   dwdw=dw(ix).^2+dw(iy).^2;
%   idx_dw= dwdw > 0; % elements are 0,1
%   ww=ww(idx_dw);  % dimention(ww)=the number of 1 in idx_dw 
%   wdw=wdw(idx_dw);
%   dwdw=dwdw(idx_dw);
%   s_w=min((-wdw+sqrt(wdw.^2+dwdw.*(1-ww)))./dwdw);
%   s_w=eta*min(s_w, 1);
%   w = w + s_w*dw;
%   u = u + du;
%   m=m+1;
%  end
%  m=1;
%   
 fprintf('%.2e % .2e % .2e % .2e % .2e\n',...
      normg, normg/normrh_old, res(iter), min(1-ww), s_w);
 normrh_old=normg;
end
res=j2(u,z,alpha);
return