pls_uncnst_tv.m

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

  function [uh] = pls_uncnst_tv(uh,bh,k_hat,alpha,beta, ...
      newt_params,pcg_flag,cg_params)

%
%

  global tvplotfigsize;
  [nx,ny] = size(uh);
  N = nx*ny;
  hx = 1/nx;
  hy = 1/ny;
  
  gradnormvec = [];
  newt_stepnormvec = [];
  newt_iter = 0;
  stop_flag = 0;
  newt_maxiter = newt_params(1);
  newt_steptol = newt_params(2);
  newt_gradtol = newt_params(3);
  newt_tab_flag= newt_params(4);
  newt_fig     = newt_params(5);
  
%  Perform Quasi-Newton iterations. 
  
  while stop_flag == 0
    newt_iter = newt_iter + 1;
    
    %  Compute regularization operator Lh.

    [dcoef_xh,dcoef_yh] = get_dcoef(uh,beta);
    Lh = get_L(dcoef_xh,dcoef_yh);
    
    %  Compute gradient of objective functional.
    
    Khuh = integral_op(uh,k_hat);
    KstarKuh = integral_op(Khuh,conj(k_hat));
    Lhuh = reshape(Lh*uh(:),nx,ny);
    gh = KstarKuh + alpha * Lhuh - bh;

    %  Solve linear system Ah*duh = -gh using PCG iteration.
    
    if (pcg_flag == 0 | pcg_flag == 1)

      [duh,cg_residnormvec,cg_stepnormvec] = pcg_os(...
        -gh,k_hat,dcoef_xh,dcoef_yh,alpha,cg_params);

    elseif (pcg_flag == 2 | pcg_flag == 3)
      
      [duh,cg_residnormvec,cg_stepnormvec] = pcg_ml(...
        -gh,k_hat,dcoef_xh,dcoef_yh,alpha,cg_params);
    
    else
      
      disp(' *** Wrong pcg_flag value in pls_uncnst_tv.m. ***');
      return
    end

    uh = uh + duh;
    
    stepnorm = norm(duh(:))/norm(uh(:));
    gradnorm = norm(gh(:));
    newt_stepnormvec = [newt_stepnormvec; stepnorm];
    gradnormvec = [gradnormvec; gradnorm;];
    cg_conv_rate = [cg_conv_rate; -mean(diff(log10(cg_residnormvec)))];

    %  Check stopping criteria.
    
    if newt_iter >= newt_maxiter
      stop_flag = 1;
    elseif stepnorm < newt_steptol 
      stop_flag = 2;
    elseif gradnorm/gradnormvec(1) < newt_gradtol
      stop_flag = 3;
    end
    
    %  Output results of current iteration.
 
    if newt_tab_flag == 1
      fprintf('\n QN iter%3.0f, ||grad||=%6.4e, CG conv rate=%6.4e\n', ...
            newt_iter, gradnormvec(newt_iter), cg_conv_rate(newt_iter));
    end
    
    if newt_fig > 0
 
      figure(newt_fig)
      set(newt_fig,...
        'units','pixels',...
        'pos',tvplotfigsize ,...
        'numbertitle','off',...
        'name','"Total Variation Non-Linear Solver"');

      colormap(jet)
      subplot(221)
        imagesc(uh), colorbar
        title('TV Reconstr.')
      subplot(222)
        semilogy(newt_stepnormvec,'o')
	xlabel('Newton iter')
        title('||Newton step||')
      subplot(223)
        semilogy(gradnormvec,'o')
	xlabel('Newton iter')
        title('||gradient||')
      subplot(224)
        plot(cg_conv_rate,'o')
	xlabel('Newton iter')
        title('CG Conv. Fact.')
      drawnow
    end
  end

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