function J=imwhiten(I,alpha,cutoff)
% IMWHITEN  Whiten image
%  J = IMWHITEN(I,ALPHA,CUTOFF) approximatively whitens the power
%  spectrum of the natural image I. The algorithm assumes that the
%  modulus of the spectrum decays as 1/f^ALPHA (f is the frequency).
%
%  IMWHITEN(I) uses ALPHA=1 (typical value for natural images).
%
%  IMWHITEN(I,ALPHA,CUTOFF) also applies a low-pass filter with cutoff
%  frequency equal to CUTOFF x FN, where FN is the Nyquist
%  frequency (half of the sampling frequency).
%

if ~exist('alpha'),  alpha = 1 ; end
if ~exist('cutoff'), cutoff = [] ; end

[M,N]=size(I) ;

% Frequency domain
fn = 0.5 ; % Nyquist freq (=1/2T, T=1)
fx_range=linspace(-fn, fn, N) ;
fy_range=linspace(-fn, fn, M) ;
[fx fy]=meshgrid(fx_range, fy_range) ;

% Whitening filter
rho=sqrt(fx.*fx+fy.*fy);
filt=rho ;

% Low-pass filter
if ~isempty(cutoff)
  fcut = cutoff * fn ;
  filt = filt .* exp(-(rho/fcut).^4);
end

% Apply filter
J = real(ifft2(fft2(I).*fftshift(filt))) ;