fixed_pt.m

matlab语言开发的图像去噪的固定点迭代法

%  Fixed_pt.m
%
%  Use "lagged diffusivity" fixed point iteration to minimize
%      T(u) = ||K*u - d||^2/2 + alpha*J(u),
%  where K is a discretized integral operator, d is discrete data, 
%  ||.|| denotes the l^2 norm, alpha is a positive regularization 
%  parameter, and J is a smooth approximation to the 
%  Total Variation functional.
%      J(u) = sum_i psi(|[D*u]_i|^2,beta) * Delta_x,
%  where D is a discretization of the first derivative operator and
%  beta is a positive smoothing parameter.
%
%  At each iteration, replace u by u + Delta_u, where Delta_u solves
%    (K'*K + alpha*L(u)) * Delta_u = K'(K*u-d) + alpha*L(u)*u,
%  where 
%    L(u)*v = D'* diag(psi'(|[D*u]_i|^2,beta) * D * Delta_x. 
%  A combination of upwind and downwind differencing is used to approximate
%  the gradient.

  max_fp_iter = 5% input(' Max. no. of fixed point iterations = ');
  max_cg_iter =5 % input(' Max. no. of CG iterations = ');
  cg_steptol = 1e-5;
  cg_residtol = 1e-5;
  cg_out_flag = 0;  %  If flag = 1, output CG convergence info.
  upwind_flag = 1;
  reset_flag = 1 %input(' Enter 1 to reset; else enter 0: ');
  if exist('f_alpha','var')
    e_fp = [];
  end
  
 if reset_flag == 1
  alpha = 40  %input(' Regularization parameter alpha = ');
  beta = 0.01 %input(' TV smoothing parameter beta = ');
  
  %  Set up discretization of first derivative operators.
  
  n = nfx;
  nsq = n^2;
  Delta_x = 1 / n;
  Delta_y = Delta_x;
  D = spdiags([-ones(n-1,1) ones(n-1,1)], [0 1], n-1,n) / Delta_x;
  I_trunc1 = spdiags(ones(n-1,1), 0, n-1,n);
  I_trunc2 = spdiags(ones(n-1,1), 1, n-1,n);
  Dx1 = kron(D,I_trunc1);   %  Forward (upwind) differencing in x
  Dx2 = kron(D,I_trunc2);   %  Backward (downwind) differencing in x
  Dy1 = kron(I_trunc1,D);   %  Forward (upwind) differencing in y
  Dy2 = kron(I_trunc2,D);   %  Backward (downwind) differencing in y

  %  Initialization.

  k_hat_sq = abs(k_hat).^2;
  Kstar_d = integral_op(dat,conj(k_hat),n,n);   %  Compute K'*dat.
  f_fp = zeros(n,n);
 end
 
  fp_gradnorm = [];
  snorm_vec = [];
  
  for fp_iter = 1:max_fp_iter
      
    %  Set up regularization operator L.

    fvec = f_fp(:);
    psi_prime1 = psi_prime((Dx1*fvec).^2 + (Dy1*fvec).^2, beta);
    Dpsi_prime1 = spdiags(psi_prime1, 0, (n-1)^2,(n-1)^2);
    L1 = Dx1' * Dpsi_prime1 * Dx1 + Dy1' * Dpsi_prime1 * Dy1;
    if upwind_flag == 1  %  Upwind differencing
      L = L1 * Delta_x * Delta_y;
    else  %  Combine upwind and downwind differencing.
      psi_prime2 = psi_prime((Dx2*fvec).^2 + (Dy2*fvec).^2, beta);
      Dpsi_prime2 = spdiags(psi_prime2, 0, (n-1)^2,(n-1)^2);
      L2 = Dx2' * Dpsi_prime2 * Dx2 + Dy2' * Dpsi_prime2 * Dy2;
      L = (L1 + L2) * Delta_x * Delta_y / 2;
    end
    KstarKf = integral_op(f_fp,k_hat_sq,n,n);    
    G = KstarKf(:) - Kstar_d(:) + alpha*(L*f_fp(:));
    gradnorm = norm(G);
    fp_gradnorm = [fp_gradnorm; gradnorm];
    
    %  Use PCG iteration to solve linear system
    % (K'*K + alpha*L)*Delta_f = r

    fprintf(' ... solving linear system using cg iteration ... \n');
    [Delf,residnormvec,stepnormvec,cgiter] = ...
      cg(k_hat_sq,L,alpha,-G,max_cg_iter,cg_steptol,cg_residtol);
    Delta_f = reshape(Delf,n,n);
      
    %  Update TV regularized reconstruction.

    f_fp = f_fp + Delta_f;
    snorm = norm(Delf);
    snorm_vec = [snorm_vec; snorm];
    if exist('f_alpha','var')
      e_fp = [e_fp; norm(f_fp - f_alpha,'fro')/norm(f_alpha,'fro')];
    end
    
    %   Output fixed point convergence information.
   
    fprintf(' FP iter=%3.0f, ||grad||=%6.4e, ||step||=%6.4e, nCG=%3.0f\n', ...
       fp_iter, gradnorm, snorm, cgiter);
   
    figure(2)
      subplot(221)
        semilogy(residnormvec/residnormvec(1),'o')
        xlabel('CG iteration')
        title('CG Relative Residual Norm')
      subplot(222)
        semilogy(stepnormvec,'o')
        xlabel('CG iteration')
        title('CG Relative Step Norm')
      subplot(223)
        semilogy([1:fp_iter],fp_gradnorm,'o-')
        xlabel('Fixed Point Iteration')
        title('Norm of FP Gradient')
      subplot(224)
        semilogy([1:fp_iter],snorm_vec,'o-')
        xlabel('Fixed Point Iteration')
        title('Norm of FP Step')
    figure(3)
      imagesc(f_fp), colorbar
      title('Reconstruction')
    figure(4)
      plot([1:nfx]',f_fp(ceil(nfx/2),:), [1:nfx]',f_true(ceil(nfx/2),:))
      title('Cross Section of Reconstruction')
    drawnow
       
  end %for fp

 %% rel_soln_error = norm(f_fp(:)-f_true(:))/norm(f_true(:))