pcg_fft.m

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

  function uh = pcg_fft(bh,k_hat,alpha,q,cg_params)

%
%  Use Tikhonov regularization (i.e., Penalized Least 
%  Squares with L^2 penalty or H^1) for 2-D image deblurring.
%
%  Let K denote the blurring operator and z the blurred data. The
%  deblurred image u minimizes
%      ||K*u - z||^2 + alpha * <L*u,u>,
%  where L is either the identity or the negative Laplacian operator.
%  u satisfies the linear system 
%      Ah * u = bh,
%  where Ah = K'*K + alpha*L and bh = K'*z.
%  This system is solved using Preconditined Conjugate Gradient iteration. 
%  The preconditioning matrix is a circulant approximation to 
%  Ah, and is implemented using 2-D FFT's.

%  Check array sizes and parameter values.

  [nx,ny] = size(bh);
  [m2,n2] = size(k_hat);
  if nx ~= ny
    disp(' *** Error in pls_uncnst_fft. Kstarz must be square. ***');
    return
  elseif m2 ~= n2
    disp(' *** Error in pls_uncnst_fft. k_hat must be square. ***');
    return
  elseif n2 ~= 2*nx
    disp(' *** Error in pls_uncnst_fft. dim(k_hat) must be twice dim(Kstarz)');
    return
  end
  if (q ~= 0 & q ~= 1)
    disp(' *** Error in pcg_fft. Reg. order q must be 0 or 1.');
    return
  end

%  Get CG iteration parameters.
  
  maxiter = cg_params(1);
  steptol = cg_params(2);
  residtol= cg_params(3);
  cg_tab_flag = cg_params(4);
  cg_fig = cg_params(5);

%  Set up preconditioner.
  
  h = 1/nx;
  hsq = h^2;
  if q == 0
    KpLinv = 1 ./ (hsq^2 * abs(k_hat).^2 + alpha);
  else
    hd2 = h/2;
    x = [0:hd2:1-hd2]';
    indx = [0:2*nx-1]';
    dvec = ((2*pi) * abs(indx - nx)).^(2*q);
    dvec = fftshift(dvec);
    dmat = dvec*ones(1,n2) + ones(n2,1)*dvec';
    KpLinv = 1 ./ (hsq^2 * abs(k_hat).^2 + alpha*hsq*dmat);
    Lh = laplacian(nx,ny);
  end

%  CG initialization.
  
  uh = zeros(nx,ny);
  resid = bh;
  residnormvec = norm(resid(:));
  stepnormvec = [];
  cgiter = 0;
  stop_flag = 0;

%  CG iteration.

  while stop_flag == 0
    cgiter = cgiter + 1;

    %  Apply preconditioner.
    
    dh = conv2d(resid,KpLinv);

    %  Compute conjugate direction p.
 
    rd = resid(:)'*dh(:);
    if cgiter == 1,
       ph = dh; 
     else
       betak = rd / rdlast;
       ph = dh + betak * ph;
    end

    %  Form product Ah*ph.

    if q == 0
      alpha_Lhph = alpha * ph;
    else
      alpha_Lhph = alpha * reshape(Lh*ph(:),nx,ny);
    end
    Ahph = integral_op(integral_op(ph,k_hat),conj(k_hat)) + 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',...
        'numbertitle','off',...
        'name','"FFT with PCG"');



        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')
	subplot(223)
	  imagesc(uh)
	  colormap(jet), colorbar
	  title('PCG solution')
      drawnow
    end
 
  end

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