%Written by: Sukanta Basu (Last updated: March 12, 2007)

%**************************************************************************
%ITERATIVE AMPLITUDE ADJUSTED FOURIER TRANSFORM SURROGATES (2D)
%Function to create 2D Surrogates that preserve PDF (exact) and
%the Spectra (near-accurate/approximate)
%**************************************************************************

function J = MakeSurrogate(varargin)

I = varargin{1};
if nargin == 2
    maxit = varargin{2};
else
    maxit = 100;
end

[NI NJ]	= size(I); N = NI*NJ;
frac 	= 0.95;
amps 	= abs(fft2(I));
sur     = aaft2D(I);

for k=1:maxit
    temp 			= fft2(sur);
    ampsnew 		= abs(temp);
    g               = find(ampsnew == 0);
    if ~isempty(g), ampsnew(g) = ampsnew(g) + 1e-15*randn(1); end;
    ampsshouldbe 	= frac*amps + (1-frac)*ampsnew;
    temp 			= temp.*ampsshouldbe./ampsnew;
    sur 			= real(ifft2(temp));
    sur_col			= reshape(sur,N,1);

    I_col			= reshape(I,N,1);
    S 				= histmatch(sur_col,I_col); sur = reshape(S,NI,NJ);

    if k == 1; PastErr = Inf; Err_min = Inf; end
    Err = mean(mean(abs(amps-ampsnew)));

    if Err < Err_min
        sur_min = sur;
        Err_min = Err;
    end
    if abs((Err-PastErr)/PastErr)*100 <= 1 | k == maxit
        J = sur_min;
        break;
    end
    PastErr = Err;
end
return
%**************************************************************************

%**************************************************************************
function J      = aaft2D(I)
[NI NJ]         = size(I); N = NI*NJ;
X 				= reshape(I,N,1);
Z				= histmatch(X);
Z				= reshape(Z,NI,NJ);
J				= surrfft2(Z);
J				= reshape(J,N,1);
J				= histmatch(J,X);
J				= reshape(J,NI,NJ);
return
%**************************************************************************

%**************************************************************************
%Function to do histogram matching: Must be a COLUMN vector
function Z  = histmatch(X,Y)
if nargin == 1
    Y 			= randn(size(X));
end
[XS,Xi] 		= sort(X);
[XS,Xii]		= sort(Xi);
YS				= sort(Y);
Z				= YS(Xii);
return
%**************************************************************************

%**************************************************************************
%Function to create simple 2D Surrogates: DOES NOT PRESERVE PDF
function J  = surrfft2(I)
[NI NJ] 	= size(I);
F 			= fft2(I);
PI			= 2*pi*rand(ceil(NI/2)-1,1);
PJ			= 2*pi*rand(1,ceil(NJ/2)-1);
PIJ		    = 2*pi*rand(NI-1,ceil(NJ/2)-1);

if rem(NJ,2)==0
    QIJ	= 2*pi*rand(ceil(NI/2)-1,1);		
    if rem(NI,2)==0
        QIJ	= [QIJ;0;-flipud(QIJ)];
    else
        QIJ	= [QIJ;-flipud(QIJ)];
    end
    PHIJ		= [PIJ QIJ -fliplr(flipud(PIJ))];
else
    PHIJ		= [PIJ -fliplr(flipud(PIJ))];
end

if rem(NI,2) == 0
    PHI = [0;PI;0;-flipud(PI)];
else
    PHI = [0;PI;-flipud(PI)];
end

if rem(NJ,2) == 0
    PHJ = [0 PJ 0 -fliplr(PJ)];
else
    PHJ = [0 PJ -fliplr(PJ)];
end
PH		 		= zeros(NI,NJ);
PH(2:NI,2:NJ)   = PHIJ;
PH(1,:)			= PHJ;
PH(:,1)			= PHI;
K               = F .* exp(j*PH);
J               = real(ifft2(K));
return
%**************************************************************************