setup_data.m

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

%  setup_data.m  
%
%  Load user defined data.  The following variables are expected:
%	Data - z, 
%	Point Source - kernel,
%	Exact Solution - u_exact  (optional) 
%  on cell-centered mesh, where we assume the model
%
%      z(x,y) = \int_Omega k(x-x',y-y') u_exact(x',y') dx' dy'
%                         + error.
%
%  Omega is the unit square in R^2.
%
%  Also, extend the kernel by zero (to reduce edge effects), compute its 
%  Fourier transform khat, compute the Fourier transform abs(khat).^2 
%  of and K^*K, and evaluate K^* z.

  if n_reduce==0
    n = max(size(z));
  else
 
    for kk = 1:n_reduce
 
    %  Sum over 4 neighbors.
  
      n = max(size(z));
      kx = kernel(1:2:n-1,:) + kernel(2:2:n,:);
      kernel = kx(:,1:2:n-1) + kx(:,2:2:n);
      kernel = kx(:,1:2:n-1) + kx(:,2:2:n);
      zx = z(1:2:n-1,:) + z(2:2:n,:);
      z = zx(:,1:2:n-1) + zx(:,2:2:n);
      clear kx zx
      if ( exist('u_exact') ),
	ux = u_exact(1:2:n-1,:) + u_exact(2:2:n,:);
	u_exact = ux(:,1:2:n-1) + ux(:,2:2:n);
	clear ux
      end
    end
    n = n/2;
  end
  
  nx = n;
  ny = n;
  hx = 1/nx;
  hy = 1/ny;

%  Plot point spread function and blurred, noisy data
if ( figure_list(2) > 0 )
	datafig = figure(figure_list(2));
else
	datafig = figure(length(findobj('type', 'figure'))+1);
	figure_list(2) = datafig;
end
  set(datafig,...
        'units','pixels',...
        'pos',dataplotfigsize ,...
        'numbertitle','off',...
        'name',datasetupflag);    

  colormap(jet)
  subplot(221)
    imagesc(kernel), title('Point Spread Function'), colorbar
  subplot(222)
    kspec = fftshift(abs(fft2(fftshift(kernel))));
    imagesc(log(kspec+sqrt(eps)*max(max(kspec))))
    title('Log Power Spectrum of PSF'), colorbar
  subplot(223)
    imagesc(z), title('Noisy, Blurred Data'), colorbar
  subplot(224)
     zspec = fftshift(abs(fft2(fftshift(z))));
    imagesc(log(zspec+sqrt(eps)*max(max(zspec))))
    title('Log Power Spectrum of Data'), colorbar
  
%  Extend kernel and compute its 2-d Fourier transform. Then use this to 
%  compute K'*z and kstark_hat, the 2-d Fourier transform for K'*K.
 
  m2 = 2*n;
  nd2 = n/2;
  kernele = zeros(m2, m2);
  kernele(nd2+1:n+nd2,nd2+1:n+nd2) = kernel;
  k_hat = fft2(fftshift(kernele));
  clear kernele
  kstark_hat = abs(k_hat).^2 * (hx*hy);
  Kstarz = integral_op(z,conj(k_hat));

%  Compute kernal for integral operator K'*K.
  kernal_KstarK = fftshift(real(ifft2(kstark_hat)));
  return