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]);