pcg_os.m

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

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

%  [uh,residnormvec,stepnormvec] = pcg_os(...
%      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 conjugate gradient iteration.
%  iteration with a preconditioner based on operator splitting.

  global cgplotfigsize;
  [nx,ny] = size(bh);
  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);
  if max(size(cg_params)) == 5
    gam = 0;
  elseif max(size(cg_params)) == 6
    gam = cg_params(6);
  else
    disp(' *** In pcg_os.m, cg_params vector has wrong size. ***');
    return
  end
  
%  CG initialization.

  if gam > 0      %  Set up preconditioner.
    sqrtIpKinv = 1 ./ sqrt((hx*hy)^2 * abs(k_hat).^2 + gam);
  end
  
  Lh = get_L(dcoef_xh,dcoef_yh);  
  uh = zeros(nx,ny);
  resid = bh;
  residnormvec = norm(resid(:));
  stepnormvec = [];
  cgiter = 0;
  stop_flag = 0;

%  CG iterations.
  
  while stop_flag == 0

    cgiter = cgiter + 1;
    if gam <= 0               % No preconditioning
      dh = resid;
    else                      % Apply preconditioning
      
    %  Solve C*d = resid, where preconditioning matrix 
    %  C = 1/gamma * (gamma + K-star*K)^(1/2) 
    %              * (gamma + alpha*L) * (gamma + K-star*K)^(1/2).

      dh = conv2d(resid*gam,sqrtIpKinv);
      dh = (gam*speye(N) + alpha*Lh) \ resid(:);
      dh = conv2d(reshape(resid,nx,ny),sqrtIpKinv);
    end
    
    rd = resid(:)'*dh(:);
    
    %  Compute conjugate direction ph.

    if cgiter == 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 cgiter >= 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', ... 
         cgiter, residnormvec(cgiter), stepnormvec(cgiter));
    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/residnormvec(1),'o')
	  xlabel('CG iteration')
	  title('CG residual norm')
	subplot(222)
	  semilogy(stepnormvec,'o')
	  xlabel('CG iteration')
	  title('CG relative step norm')
      drawnow
    end

  end

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