testImg1 = imread('Mable.jpg');

imagesc(testImg1); axis image off % look at image in color
testImg1 = mean(testImg1,3); % throw out color info
imagesc(testImg1); axis image off; colormap(gray) % look at image, grayscale
f = fftshift(fft2(fftshift(testImg1))); 
imagesc(abs(f),[0 1e5]); axis image off % look at Fourier transform (mag)

% Check to see that transforming and re-transforming gives us back the
% original image
testImg = fftshift(ifft2(fftshift(f)));
subplot(2,2,1); imagesc(real(testImg)); colorbar; axis image off; title('real')
subplot(2,2,2); imagesc(imag(testImg)); axis image off; colorbar; title('imaginary')
subplot(2,2,3); imagesc(abs(testImg)); axis image off; colorbar; title('magnitude')
subplot(2,2,4); imagesc(angle(testImg)); axis image off; colorbar; title('phase')

% It doesn't.
% I think the phase roll b/w real and imaginary comes from the
% non-squareness and/or the odd dimensions of the image ...
testImg = testImg1(41:296,51:306); % get a square image

f = fftshift(fft2(fftshift(testImg)));
subplot(1,2,1); imagesc(testImg); axis image off
subplot(1,2,2); imagesc(abs(f),[0 1e5]); axis image off
% Now make sure that everything's fine - all the power in the reconstructed
% image should be in the real part, not the imaginary
testImg = fftshift(ifft2(fftshift(f)));
subplot(2,2,1); imagesc(real(testImg)); colorbar; axis image off; title('real')
subplot(2,2,2); imagesc(imag(testImg)); axis image off; colorbar; title('imaginary')
subplot(2,2,3); imagesc(abs(testImg)); axis image off; colorbar; title('magnitude')
subplot(2,2,4); imagesc(angle(testImg)); axis image off; colorbar; title('phase')

% Now let's phase-scramble the square test image
f = fftshift(fft2(fftshift(testImg - mean(testImg(:)))));
magF = abs(f);
phF = angle(f);
% check to see that we can put magnitude and phase back together perfectly
testImg = fftshift(ifft2(fftshift(magF.*exp(i*phF))));
subplot(2,2,1); imagesc(real(testImg)); colorbar; axis image off; title('real')
subplot(2,2,2); imagesc(imag(testImg)); axis image off; colorbar; title('imaginary')
subplot(2,2,3); imagesc(abs(testImg)); axis image off; colorbar; title('magnitude')
subplot(2,2,4); imagesc(angle(testImg)); axis image off; colorbar; title('phase')

% vonMisesRand was written by Paul Schrater in 2001 to allow us to control
% the distribution of random phases that we added to the phase spectrum of
% a natural image
sigma = 1;
randPh = prsVonMisesRand(prod(size(phF)),0,sigma);
clf; hist(randPh,50); title(['Von Mises distribution, sigma = ' num2str(sigma,2) ])

phF_vonMises = phF + reshape(randPh,size(phF));

% look at the original and perturbed mag and phase parts (mag isn't
% perturbed ...)
subplot(2,2,1); imagesc(magF,[0 1e5]); axis image; title('original mag spectrum')
subplot(2,2,2); imagesc(phF); axis image; title('original phase spectrum')
subplot(2,2,3); imagesc(magF,[0 1e5]); axis image; title('same mag spectrum')
subplot(2,2,4); imagesc(phF_vonMises); axis image; title('perturbed ph spectrum')

% Adding perturbations to the phase spectrum introduces values outside of
% the acceptable -pi to pi range
clf
subplot(1,3,1); hist(phF(:),50); title('original phase spectrum')
subplot(1,3,2); hist(phF_vonMises(:),50); title('perturbed phase spectrum')
% Here's a quick one-liner for wrapping phases back into the acceptable
% range
phF_vonMises = angle(exp(i*phF_vonMises));
subplot(1,3,3); hist(phF_vonMises(:),50); title('perturbed, fixed phase spectrum')

% Now put Mable back together again ...
Mable = fftshift(ifft2(fftshift(magF.*exp(i*phF_vonMises))));

subplot(2,2,1); imagesc(real(Mable)); axis image off; colorbar; title('real Mable')
subplot(2,2,2); imagesc(imag(Mable)); axis image off; colorbar;title('imaginary Mable')
subplot(2,2,3); imagesc(abs(Mable)); axis image off; colorbar;title('magnitude Mable')
subplot(2,2,4); imagesc(angle(Mable)); axis image off; colorbar;title('phase Mable')

% Oops ... by failing to respect the symmetry in the Fourier domain, we
% dumped a bunch of power into the imaginary part of the inverse Fourier
% transform ...
% Put Mable back together again with the right symmetry...
Mable = fftshift(ifft2(fftshift(magF.*exp(i*caoReflectPhaseSpectrum(phF_vonMises)))));

subplot(2,2,1); imagesc(real(Mable)); axis image off; colorbar; title('real Mable')
subplot(2,2,2); imagesc(imag(Mable)); axis image off; colorbar; title('imaginary Mable')
subplot(2,2,3); imagesc(abs(Mable)); axis image off; colorbar; title('magnitude Mable')
subplot(2,2,4); imagesc(angle(Mable)); axis image off; colorbar; title('phase Mable')

% We need to think about contrast energy in the original and perturbed images.
% Although this is an insufficient definition of
% natural image contrast, we'll use RMS contrast to describe the contrast of
% the input image.
testImg1 = mean(imread('Mable.jpg'),3);
testImg = testImg1(41:296,51:306);
RMS = std(testImg(:));
RMScontrast = std(testImg(:))/mean(testImg(:));
subplot(1,2,1); imagesc(testImg); axis image off; colorbar; title(['RMS contrast:  ' num2str(100*RMScontrast,2) '%'])
subplot(1,2,2); hist(testImg(:),50); colormap(gray)
title(['RMS:  ' num2str(RMS,2)])
 
% Normalize pixel intensity distribution
targetDist = randn([1 256*256]); % generate the pixel values we want in the image
targetDist = sort(targetDist); % sort them
[sorted order] = sort(testImg(:)); % sort the input image
testImg(order) = targetDist;
clf; imagesc(testImg); title('Normalized pixel intensity distribution'); axis image off
% Mable looks ugly with a normal intensity distribution; let's pick someone
% else
% testImg2 = imread('Cottonwood.jpg');
% testImg2 = testImg2(121:2:1144,51:2:1074,:);
testImg2 = imread('CumbresToltec.jpg');
testImg2 = testImg2(21:2:1044,51:2:1074,:);

subplot(1,2,1); imagesc(testImg2); axis image off
testImg = mean(testImg2,3); % throw out color info
imagesc(testImg); axis image off

targetDist = sort(randn([1 prod(size(testImg))])); % generate the pixel values we want in the image
[sorted order] = sort(testImg(:)); % sort the input image
testImg(order) = targetDist;
subplot(1,2,2); imagesc(testImg); axis image off

% Back to where we were - now that we have an image with a normal
% distribution and unit standard deviation, RMS contrast is more meaningful
example = testImg*42 + 127; example = min(max(example,0),255);
RMS = std(example(:)); RMScontrast = std(example(:))/mean(example(:));
subplot(3,2,1); imagesc(example); axis image off; colorbar; title(['RMS contrast:  ' num2str(100*RMScontrast,2) '%'])
subplot(3,2,2); hist(example(:),50); axis([0 255 0 2e4]); title(['RMS:  ' num2str(RMS,2)])

example = testImg*30 + 127; example = min(max(example,0),255);
RMS = std(example(:)); RMScontrast = std(example(:))/mean(example(:));
subplot(3,2,3); imagesc(example); axis image off; colorbar; title(['RMS contrast:  ' num2str(100*RMScontrast,2) '%'])
subplot(3,2,4); hist(example(:),50); axis([0 255 0 2e4]); title(['RMS:  ' num2str(RMS,2)])

example = testImg*10 + 127; example = min(max(example,0),255);
RMS = std(example(:)); RMScontrast = std(example(:))/mean(example(:));
subplot(3,2,5); imagesc(example); axis image off; colorbar; title(['RMS contrast:  ' num2str(100*RMScontrast,2) '%'])
subplot(3,2,6); hist(example(:),50); axis([0 255 0 2e4]); title(['RMS:  ' num2str(RMS,2)])

% Do phase perturbation, demonstrating that contrast isn't messed up ...
f = fftshift(fft2(fftshift(testImg - mean(testImg(:))))); magF = abs(f); phF = angle(f);
sigmas = [0 .1 .5 1 5];
for vonMisesSigma = sigmas;
  randPh = prsVonMisesRand(prod(size(phF)),0,vonMisesSigma);
  phF_vonMises = angle(exp(i*(phF + reshape(randPh,size(phF)))));
  scrambled = real(fftshift(ifft2(fftshift(magF.*exp(i*caoReflectPhaseSpectrum(phF_vonMises))))));
  RMS = std(scrambled(:));
  subplot(length(sigmas),1,find(sigmas == vonMisesSigma))
  imagesc(scrambled); axis image off; colorbar; title(['vonMises sigma: ' num2str(vonMisesSigma,2) ', RMS: ' num2str(RMS,2)])
end

for n = 1:length(sigmas)
    vonMisesSigma = sigmas(n);
    randPh = prsVonMisesRand(prod(size(phF)),0,vonMisesSigma);
    phF_vonMises = angle(exp(i*(phF + reshape(randPh,size(phF)))));
    scrambled(:,:,n) = real(fftshift(ifft2(fftshift(magF.*exp(i*caoReflectPhaseSpectrum(phF_vonMises))))));
end 
clf; imagesc(caoMakeImageMontage(scrambled,1,size(scrambled,3))); axis image off

% What about phase scrambling in particular frequency or orientation bands?
testImg1 = mean(double(imread('twoGates.JPG')),3);
testImg = testImg1(1:2:512,1:2:512);
f = fftshift(fft2(fftshift(testImg - mean(testImg(:))))); magF = abs(f); phF = angle(f);
subplot(1,2,1); imagesc(abs(f)); axis image

centerFreq = 30; % cycles per image
freqSigma = 10; % Set this to Inf if you want all frequencies
orientation = 60; % degrees
orientBW = 20; % Set this to Inf if you want all orientations

[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
if freqSigma == Inf;  freqBand = ones(size(r));
else
    freqBand = exp(-((r-centerFreq).^2/2/freqSigma.^2));
end

if orientBW == Inf; passBand = freqBand;
else
    theta = atan(x./y); theta(size(f,2)/2+1,size(f,1)/2+1) = 0;
    theta(size(f,2)/2+1:end,:) = theta(size(f,2)/2+1:end,:) - pi;
    theta = theta + 3*pi/2;
    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
subplot(1,2,2); imagesc(passBand); axis image 

% First try just filtering the image with this to be sure it works ...
filtered = fftshift(ifft2(fftshift(f.*passBand)));
subplot(1,3,1); imagesc(testImg); axis image off 
subplot(2,3,2); imagesc(real(filtered)); axis image off; colorbar; title('real')
subplot(2,3,3); imagesc(imag(filtered)); axis image off; colorbar; title('imaginary')
subplot(2,3,5); imagesc(abs(filtered)); axis image off; colorbar; title('magnitude')
subplot(2,3,6); imagesc(angle(filtered)); axis image off; colorbar; title('phase')
% I don't like the striping in the imaginary part of the filtered image -
% it means that something's not perfect about the symmetry of the pass band
% But it's small, so I'm leaving it for now ...

% Do phase perturbation
%testImg1 = mean(double(imread('MABLE.JPG')),3);
%testImg = testImg1(41:296,51:306);
f = fftshift(fft2(fftshift(testImg - mean(testImg(:))))); 
magF = abs(f); phF = angle(f);

centerFreq = 10; % cycles per image
freqSigma = 10; % Set this to Inf if you want all frequencies
orientation = 30; % degrees
orientBW = Inf; % Set this to Inf if you want all orientations

[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
if freqSigma == Inf;  freqBand = ones(size(r));
else
    freqBand = exp(-((r-centerFreq).^2/2/freqSigma.^2));
end

if orientBW == Inf; passBand = freqBand;
else
    theta = atan(x./y); theta(size(f,2)/2+1,size(f,1)/2+1) = 0;
    theta(size(f,2)/2+1:end,:) = theta(size(f,2)/2+1:end,:) - pi;
    theta = theta + 3*pi/2;
    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
subplot(2,3,1); imagesc(testImg); axis image off; title('input image')
subplot(2,3,4); imagesc(passBand); axis image off; title('pass band for scrambling')

 vonMisesSigma = 10;
 randPh = prsVonMisesRand(prod(size(phF)),0,vonMisesSigma);
 phF_vonMises = angle(exp(i*(phF + passBand.*reshape(randPh,size(phF)))));
 filtered = real(fftshift(ifft2(fftshift(magF.*exp(i*caoReflectPhaseSpectrum(phF_vonMises))))));
subplot(2,3,2); imagesc(real(filtered)); axis image off; colorbar; title('real')
subplot(2,3,3); imagesc(imag(filtered)); axis image off; colorbar; title('imaginary')
subplot(2,3,5); imagesc(abs(filtered)); axis image off; colorbar; title('magnitude')
subplot(2,3,6); imagesc(angle(filtered)); axis image off; colorbar; title('phase')
