%	to correct for image nonuniformity
%	first fills in holes, then generates a low frequency 
%	background image and makes correction factor
%	input image is stored in img, output is Iout

%	rev 0	4/8/00		original


foo = size(img);
npix = foo(1);
np2 = npix/2;

thres = input('gimme threshold % [5] =');
if(isempty(thres))
  thres = 5;
end;
thres = thres*.01;

Imax = max(max(img));
Ithr = Imax*thres;

%  fill in holes

mask = (img>Ithr);
Imask = img.*mask;
Iave = sum(sum(Imask))/sum(sum(mask));
Ifill = Imask + (1-mask).*Iave;

%  make low freq image

z = fftshift(fft2(Ifill));

%ndef = 3;
ndef = np2/8;
n = input(sprintf('half size of filter [%d] =', ndef));
if(isempty(n))
  n = ndef;
end;
%i1 = np2 - n + 1;
%i2 = np2 + n;
%zl = zeros(npix, npix);
%zl(i1:i2,i1:i2) = z(i1:i2,i1:i2);

y2 = (1:np2).^2;  % use 15 as default win for guassian
a2 = 1/n;
for x=1:np2
  r2 = y2 +x*x;
  win1(x,:) = exp(-a2*r2);
end
win(np2+1:npix,np2+1:npix) = win1;
win(np2+1:npix,np2:-1:1) = win1;
win(1:np2,:) = win(npix:-1:np2+1,:);
zl = win.*z;

Ilow = abs(ifft2(fftshift(zl)));

%  corrected image

Ilave = sum(sum(Ilow.*mask))/sum(sum(mask));
Icorf = (Ilave./Ilow).*mask;  % correction factor
Iout = img.*Icorf;


