www.pudn.com > ibd_blindequalizer.rar > ibd.m


% Function for iterative Blind Deconvolution
%
% blurred   : input blurred image 
% fo        : initial estimate of the image
% mxiters   : maximum number of iterations  
% noise     : threshold level for Ftfm(blurred), try between 0-1
% beta      : averaging parameter, usually ~ 0.9
%
%
% reference : G. R. Ayers and J. C. Dainty, "Iterative blind deconvolution
%             method and its application," Optics Letters, vol 13(7), 
%             pp.547-549, July 1988
%
% Alper Corlu 12/11/2002
% corlua@dept.physics.upenn.edu

function [fimg,h] = ibd (blurred, fo, mxiters,noise,beta)


% the size of the blurred image
  [nrow ncol] = size (blurred);

%
% start timing
%
tic

% take the F-tfm of the input image
Cblur = fft2(blurred);

% start the itreation loop

iter = 0;
Fnew = zeros (nrow,ncol);

fimg = fo;

while iter < mxiters
     
     % progress bar
       fprintf(1,'.');
       
     % get the fourier transform of the image
       Fold = fft2(fimg);
     
     % form the estimate for  H
       H = Cblur ./ Fold;
       
     % take the inverse tfm of H
       h = real(ifft2 (H));
      
     % impose blur constraints
       h (find (h < 0.0) ) = 0;
       
     % take the fourier tfm of h
       H = fft2(h);
       
     % form the new estimate from Cblur and H
     % take care of the regions of low value
    
       index = find (abs(H) >= abs(Cblur));
       Fnew(index) = (1-beta).*Fold(index) + beta.*(Cblur(index)./H(index));
       
       index = find (abs(H) < abs(Cblur));
       Fnew(index) = ( (1-beta)./Fold(index)  + beta.*H(index)./Cblur(index) ).^-1;
       
       index = find (abs(Cblur) < noise );
       Fnew(index) = Fold(index);
       
     % get the inverse-fft of Fnew
       fimg = real(ifft2(Fnew));
     
     % impose image constraints 
       fimg (find (fimg < 0.0 )) = 0;
     
      
     iter = iter + 1;
    
end
% stop timing
toc

figure, imagesc (fimg); colormap gray;
daspect([1 1 1]);
figure, imagesc (h); colormap gray;
daspect([1 1 1]);