function noiseSample = caoLandyTexture(imgSize,tuningParams,graphics)

% Usage: noiseSample = caoLandyTexture(imgFile,[vmSigma],[tuningParams],[cmplxFlag])
%      where  imgSize is the size of the output noiseSample (e.g. 64 or [64 128]).
%             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
%             display:  pass this flag if you want to see the filter used
% Example:  noiseSample = caoLandyTexture(128,[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
inputImg = randn(imgSize);

dims = size(inputImg);

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

% 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
% Reconstruct filtered image
noiseSample = real(fftshift(ifft2(fftshift(f.*passBand))));
noiseSample = (noiseSample - mean(noiseSample(:)))./std(noiseSample(:));

if exist('graphics','var') & ~isempty(graphics)
    % show original (windowed) image and passband
    subplot(2,2,1); imagesc(inputImg); axis image off; title('input image')
    subplot(2,2,3); imagesc(passBand); axis image off; title('pass band for scrambling')

    % and show results
    subplot(2,2,2); imagesc(noiseSample); axis image off; colorbar; title('noise sample')
    subplot(2,2,4); imagesc(abs(fftshift(fft2(fftshift(noiseSample))))); axis image off; colorbar; title('FT of noise')
end

set(gcf,'Color',[1 1 1]); colormap(gray)

