pls_uncnst_fft.m

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

  function u = pls_uncnst_fft(Kstarz,k_hat,alpha,q,fft_fig)
  
%  u = pls_uncnst_fft(Kstarz,k_hat,alpha)
%  u = pls_uncnst_fft(Kstarz,k_hat,alpha,q)
%
%  Uses an FFT-based implementation of Tikhonov order q Regularization 
%  (i.e., Penalized Least Squares with Sobolev H^q penalty) 
%  for 2-D image deblurring.
%
%  Let khat and zhat denote the 2-D FFT's of the PSF and the blurred 
%  data, respectively. Compute the minimizer uhat for
%     ||khat .* uhat - zhat||^2 + alpha * ||(2*pi*i)^q.*uhat||^2.
%  The deblurred image is then the 2-D inverse fft of uhat.
%
%  q = 0 corresponds to squared L^2 penalty, yielding Tikhonov 
%  Regularization with the identity.
%
%  This implementation is very fast, but the deblurred image may
%  contain artifacts due to the periodicity of the Fourier Transform.

  [m,n] = size(Kstarz);
  [m2,n2] = size(k_hat);
  if m ~= n
    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*n
    disp(' *** Error in pls_uncnst_fft. dim(k_hat) must be twice dim(Kstarz)');
    return
  end
  if nargin == 3
    q = 0;
  elseif q < 0
    disp(' *** Error in pls_uncnst_fft. Reg. order q must be nonnegative.');
    return
  end
  
  h = 1/n;
  hsq = h^2;
  
  if q == 0
    u = conv2d(Kstarz, 1./(hsq^2 * abs(k_hat).^2 + alpha));
  else
    hd2 = h/2;
    x = [0:hd2:1-hd2]';
    indx = [0:2*n-1]';
    dvec = ((2*pi) * abs(indx - n)).^(2*q);
    dvec = fftshift(dvec);
    dmat = dvec*ones(1,2*n) + ones(2*n,1)*dvec';
    u = conv2d(Kstarz, 1./(hsq^2 * abs(k_hat).^2 + alpha*hsq*dmat));
  end
figure(fft_fig);
set(fft_fig,...
        'units','pixels',...
        'numbertitle','off',...
        'name','"FFT without PCG"');

        subplot(221)
          imagesc(u)
          colormap(jet), colorbar
          title('Approx. Solution')     


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