pcg_ml.m

采用matlab编写的数字图像恢复程序

  function [uh,residnormvec,stepnormvec] = pcg_ml(...
      bh,k_hat,dcoef_xh,dcoef_yh,alpha,cg_params)

%  [uh,residnormvec,stepnormvec] = pcg_ml(...
%        bh,k_hat,dcoef_xh,dcoef_yh,alpha,cg_params);
%
%  Invert the linear system
%    (Kh'*Kh + alpha*Lh)*uh = bh.
%  Here alpha is a positive parameter, Kh is an integral operator 
%  of convolution type, and Lh is a discretization of the 
%  diffusion operator 
%        L(dcoef)u = -div (dcoef grad u).
%  Inversion is carried out using a preconditioned conjugate gradient 
%  iteration with a preconditioner based on multilevel splitting.

  global cgplotfigsize;
  [nx,ny] = size(bh);
  p = log2(nx);
  N = nx * ny;
  hx = 1/nx;
  hy = 1/ny;
  maxiter = cg_params(1);
  steptol = cg_params(2);
  residtol= cg_params(3);
  cg_tab_flag = cg_params(4);
  cg_fig  = cg_params(5);
  p_type  = cg_params(6);
  p0      = cg_params(7);
  Q_iter  = cg_params(8);
  nu      = cg_params(9);

%  Compute fine grid regularization operator Lh.

  Lh = get_L(dcoef_xh,dcoef_yh);

%  Compute coarse grid regularization operator LH.

  [dcoef_xH,dcoef_yH] = restrict_dcoef(dcoef_xh,dcoef_yh,p-p0);
  LH = get_L(dcoef_xH,dcoef_yH);
  
%  Compute coarse grid convolution integral operator KstarKH.

  kstark_hat = abs(k_hat).^2 * (hx*hy);
  kernel_KstarK = fftshift(real(ifft2(kstark_hat)));
  kerKstarKH = restrict_ker(kernel_KstarK,p-p0);
  kstarkH_hat = fft2(fftshift(kerKstarKH));
  KstarKH = get_Kmat(kstarkH_hat);
  
  AH = KstarKH + alpha*LH;

%  PCG initialization.

  uh = zeros(nx,ny);
  resid = bh;
  residnormvec = norm(resid(:));
  stepnormvec = [];
  pcgiter = 0;
  stop_flag = 0;

%  PCG iterations.
  
  while stop_flag == 0
    pcgiter = pcgiter + 1;

    if p_type == 1
      
    %  Apply Jacobi preconditioner to solve Ah*d = resid.
    
      [dh,Qresidnormvec] = ...
             precond_jac(resid,AH,p0,dcoef_xh,dcoef_yh,alpha,Q_iter,nu);
    elseif p_type == 2
      
    %  Symmetric Gauss Seidel preconditioner.
    
      [dh,Qresidnormvec] = ...
        precond_sgs(resid,AH,p0,kstark_hat,dcoef_xh,dcoef_yh,alpha,Q_iter,nu);
    else
      disp(' *** Error in pcgml.m.  Illegal p_type value. ***');
      return
    end

    %  Compute conjugate direction ph.

    rd = resid(:)'*dh(:);
    if pcgiter == 1,
       ph = dh; 
     else
       betak = rd / rdlast;
       ph = dh + betak * ph;
    end

    %  Form product Ah*ph.
    
    Khph = integral_op(ph,k_hat);
    KstarKhph = integral_op(Khph,conj(k_hat));
    Lhph = reshape(Lh*ph(:),nx,ny);
    Ahph = KstarKhph + alpha*Lhph;
 
    %  Update uh and residual.
    
    alphak = rd / (ph(:)'*Ahph(:));
    uh = uh + alphak*ph;
    resid = resid - alphak*Ahph;
    rdlast = rd;

    residnorm = norm(resid(:));
    stepnorm = abs(alphak)*norm(ph(:))/norm(uh(:));
    residnormvec = [residnormvec; residnorm];
    stepnormvec = [stepnormvec; stepnorm];
    
    %  Check stopping criteria.
    
    if pcgiter >= maxiter
      stop_flag = 1;
    elseif stepnorm < steptol
      stop_flag = 2;
    elseif residnorm / residnormvec(1) < residtol
      stop_flag = 3;
    end
    
    %  Output convergence information.

    if cg_tab_flag == 1
      fprintf('  CG iter%3.0f, ||resid||=%6.4e, ||step||=%6.4e \n', ... 
         pcgiter, residnormvec(pcgiter), stepnormvec(pcgiter));
    end
    if cg_fig > 0
      figure(cg_fig)
      set(cg_fig,...
        'units','pixels',...
        'pos',cgplotfigsize ,...
        'numbertitle','off',...
        'name','"Total Variation - CG Performance"');
        
        subplot(221)
          semilogy(residnormvec,'o')
	  xlabel('CG iteration')
	  title('CG residual norm')
	subplot(222)
	  semilogy(stepnormvec,'o')
	  xlabel('CG iteration')
	  title('CG relative step norm')
        subplot(223)
          semilogy(Qresidnormvec/Qresidnormvec(1),'o')
          title('PCG/MG performance for QLQ*Qu=Qb')
      drawnow
    end

  end

global figure_list;
uicontrol(cg_fig, ...
        'style','push', ...
        'callback','close(gcf),figure_list(3) = -1;', ...
        'string','Close');