function [scrambled,inputImg] = caoPhaseScramble(inputImg,vmSigma,tuningParams,cmplxFlag,graphFlag)

% Usage: scrambled = caoPhaseScramble(inputImg,[vmSigma],[tuningParams],[cmplxFlag],[graphFlag])
%      where  inputImg is the path to the image you want to scramble *or*
%               a variable containing the image.  For now, the code centers
%               the biggest 2^m x 2^n square it can in the image, and works
%               with that, since the transform/inverse transform process
%               seems to move power to the imaginary part of the image if
%               it's not dimensioned in powers of 2.  If imgFile is a path,
%               the code then writes the 2^m x 2^n subset scrambled image
%               in the same format, after matching the pixel intensity
%               distribution.
%             vmSigma is the width of the vonMises distribution used to
%               perturb the phases. If not provided, 10 is used for full
%               scrambling
%             tuningParams is an optional vector with 4 parameters:
%                 [centerFreq freqSigma orientation orientBW]
%                 centerFreq:  for pass band, if only some of image
%                   spectrum is to be scrambled.
%                 freqSigma:  width of frequency passband ( = Inf if
%                   scrambling is to be broadband in frequency)
%                 orientation:  of pass band
%                 orientBW:  set to Inf if scrambling isn't to be tuned
%                   in orientation
%             cmplxFlag:  set to 0 or leave out if you want a real-valued
%               image returned; set to 1 if you want the complex image so
%               you can check if the code messed up your contrast
%             graphFlag: set to 1 to see output; default 0
% Example:  scrambled = caoPhaseScramble('twoGates.jpg',10,[10 10 30 Inf],1);
% 
% cao 11/06

% Get the input image, make sure it's grayscale, and window it down to a 
% 2^m x 2^n size
if isstr(inputImg)
    [FPATH FNAME ext] = fileparts(inputImg);
    outputFile1 = fullfile(FPATH,[FNAME '_scr' ext]);
    outputFile2 = fullfile(FPATH,[FNAME '_cropped' ext]);
    inputImg = double(imread(inputImg));
end
if size(inputImg,3) > 1; inputImg = mean(inputImg,3); end
dims = size(inputImg);
finalDims = 2.^(floor(log2(dims)));
imgCenter = floor(dims/2); imgEdge = imgCenter - finalDims/2;
inputImg = inputImg(imgEdge(1)+(1:finalDims(1)),imgEdge(2)+(1:finalDims(2)));

% Take the Fourier transform and separate out the magnitude and phase
% components
f = fftshift(fft2(fftshift(inputImg))); 
magF = abs(f); phF = angle(f);

% If tuningParams were specified, then build pass filter for scrambling
if exist('tuningParams','var') & ~isempty(tuningParams)
    centerFreq = tuningParams(1); % cycles per image
    freqSigma = tuningParams(2); % Set this to Inf if you want all frequencies
    orientation = tuningParams(3); % degrees
    orientBW = tuningParams(4); % Set this to Inf if you want all orientations
    % make generic matrices, where r represents frequency and theta,
    % orientation
    [x y] = meshgrid(-size(f,2)/2:size(f,2)/2-1,-size(f,1)/2:size(f,1)/2-1);
    r = sqrt(x.^2 + y.^2); % this has inherent units of cycles per image
    % There must be a better way to do this - this is my historic cluge
    y(find(y == 0)) = .01;
    theta = atan(x./y); 
    theta(size(f,2)/2+1:end,:) = theta(size(f,2)/2+1:end,:) - pi;
    theta = theta + 3*pi/2;
    % Build filter
    if freqSigma == Inf;  freqBand = ones(size(r));
    else
        freqBand = exp(-((r-centerFreq).^2/2/freqSigma.^2));
    end
    if orientBW == Inf; passBand = freqBand;
    else
        orientBand = exp(-((theta-pi*(orientation+90)/180).^2/2/(pi*orientBW/180).^2)) ...
            + exp(-(((angle(exp(i*(theta)))+pi)-pi*(orientation+90)/180).^2/2/(pi*orientBW/180).^2));
        passBand = freqBand.*orientBand;
    end
else
    passBand = ones(size(f));
end

% Now calculate phase perturbations.  prsVonMisesRand and
% caoReflectPhaseSpectrum are imbedded in this one file for convenience,
% but also exist as stand-alone functions
if ~exist('vmSigma','var') | isempty(vmSigma)
    vmSigma = 10;
end
randPh = vonMisesrand(prod(size(phF)),0,vmSigma);
phF_vonMises = angle(exp(i*(phF + passBand.*reshape(randPh,size(phF)))));

% Reconstruct filtered image
scrambled = real(fftshift(ifft2(fftshift(magF.*exp(i*caoReflectPhaseSpectrum(phF_vonMises))))));
% and show results
if ~exist('graphFlag','var'); graphFlag = 0; end
if graphFlag
    % show original (windowed) image and passband
    subplot(2,3,1); imagesc(inputImg); axis image off; title('input image')
    subplot(2,3,4); imagesc(passBand); axis image off; title('pass band for scrambling')
    subplot(2,3,2); imagesc(real(scrambled)); axis image off; colorbar; title('real')
    subplot(2,3,3); imagesc(imag(scrambled)); axis image off; colorbar; title('imaginary')
    subplot(2,3,5); imagesc(abs(scrambled)); axis image off; colorbar; title('magnitude')
    subplot(2,3,6); imagesc(angle(scrambled)); axis image off; colorbar; title('phase')
    set(gcf,'Color',[1 1 1]); colormap(gray)
end
if ~exist('cmplxFlag','var') | (cmplxFlag == 0) 
    scrambled = real(scrambled);
end

outDist = sort(inputImg(:));
[temp order] = sort(scrambled(:));
scrambled(order) = outDist;
    
if exist('outputFile1','var')
    imwrite(uint8(scrambled),outputFile1); disp('writting cropped file.')
    imwrite(uint8(inputImg),outputFile2); disp('writting cropped file.')
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function reflectim = caoReflectPhaseSpectrum(im)

dims = size(im);
reflectim = im;
% fix the center column
 reflectim(dims(1):-1:(dims(1)/2+2),dims(2)/2+1) = -reflectim(2:dims(1)/2,dims(2)/2+1);
% reflect the 2nd half of the phase spectrum
reflectim(dims(1):-1:2,dims(2):-1:(dims(2)/2+2)) = -im(2:dims(1),2:dims(2)/2);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vonrand = vonMisesrand(nrand,mu,sigma)
% Usage:   vonrand = vonMisesrand(nrand,mu,sigma)
%  Returns nrand samples from circular normal (vonMises) distribution,
%  centered at mu with standard deviation sigma.
%
% inverse cumulative method, executed by table lookup with
% linear interpolation.  Written by Paul Schrater 2001
%
% C. Olman: 2006 - imbedded vonMisespdf so only have to keep track of one
% file.

if sigma > 0.25,
	% build sampled cdf
	x = (-pi:2*pi/(2e3):pi);
	pofx = vonMisespdf(x,0,sigma);
	cofx = cumsum(pofx/sum(pofx));
	%
	u = rand(1,nrand);
	vonrand = interp1(cofx,x,u)+mu;
	vonrand = min(max(vonrand,-pi+eps),pi-eps);
else
	vonrand = randn(1,nrand)*sigma +mu;
end	

return

function pofx = vonMisespdf(x,mu,sigma)
% For -pi <= x <= pi
% force x-mu within -pi to pi
y = angle(exp(i*(x-mu)));
kappa = 1/(sigma)^2;
pofx = exp(kappa*cos(y))/(2*pi*besseli(0,kappa));
return
