function height = GenFractalTexPublic(n, beta, randSeed)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Generates a fractal surface of size n x n and returns an array of values
% from [0 : 1] which can be saved as a greyscale heightmap image.
%
%
% v1.01 Mike Chantler 22/5/2006
% v1.20 Khemraj Emrith 21/09/2006 
% v1.50 Lin Qi 17/11/2008
% v2.00 Thomas Methven & Lin Qi 30/04/2012 **CURRENT VERSION**
%
% Created by the Texture Lab, Heriot-Watt University, original source:
% http://www.macs.hw.ac.uk/texturelab/resources/scripts/
% 
%   n = size of the image (n x n)
%   beta = spetral rolloff factor ('roughness')
%   randSeed = is for seeding the random generator. Leave blank for
%              nondeterministic randomness.
%
%   EXAMPLE USES: height = GenFractalTexPublic(512, 2.2);
%                 (Nondeterministic, relatively smooth surface, 512 x 512
%                 pixels)
%
%                 height = GenFractalTexPublic(1024, 1.5, 12345);
%                 (Deterministic (i.e. will give the same result every 
%                 time), relatively rough surface, 1024 x 1024 pixels)
%
%   N.B. This function has been updated to use the newer random seeding
%   technique as recommended by MathWorks. If the 'rng' command doesn't
%   work on your version of Matlab, swap it with the commented line below.
%   (Remember to comment out the function check just below here too!)
%   ***PLEASE NOTE THAT THE DIFFERENT RANDOM GENERATOR VERSIONS WILL GIVE
%   DIFFERENT FRACTAL SURFACES EVEN WHEN USING THE SAME SEED!***
%   Explaination here:
%   http://www.mathworks.co.uk/help/techdoc/math/bsn94u0-1.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(exist('rng', 'file') == 0)
    disp(' ');
    disp('ERROR: The function ''rng'' does not exist in your version of Matlab!');
    disp('ERROR: This code''s random number generator has been updated to the latest version,');
    disp('ERROR: as recommended by MathWorks. Check N.B. in source code introduction for solution.');
    disp(' ');
    return;
end    

%%%%%%%%%%%%%%%%%%%%%
% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%
if (nargin == 0)
    disp(' ');
    disp('ERROR: MISSING AT LEAST 2 VARIABLE INPUTS. The format is ''GenFractalTexPublic(n, beta, [randSeed]);''');
    disp('ERROR: ''n'' is the size of the image (n x n)');
    disp('ERROR: ''beta'' is the spetral rolloff factor (''roughness'')');
    disp('ERROR: ''randSeed'' is for seeding the random generator. Leave blank for nondeterministic randomness');
    disp(' ');
    return;
elseif (nargin == 1)
    disp(' ');
    disp('ERROR: MISSING AT LEAST 1 VARIABLE INPUT. The format is ''GenFractalTexPublic(n, beta, [randSeed]);''');
    disp('ERROR: ''n'' is the size of the image (n x n)');
    disp('ERROR: ''beta'' is the spetral rolloff factor (''roughness'')');
    disp('ERROR: ''randSeed'' is for seeding the random generator. Leave blank for nondeterministic randomness');
    disp(' ');
    return;
elseif (nargin == 3)
    % Generate random phase with a known seed - Deterministic.
    rng(randSeed);
    %rand('state', randSeed);
    
    theta = rand(n, n) * 2 * pi; 
else
    % Generate random phase - Nondeterministic.
    rng('shuffle');
    %rand('seed', sum(100*clock));
    
    theta = rand(n, n) * 2 * pi; 
end

%%%%%%%%%%%%%%%%%%%%%
% GENERATE FRACTAL ARRAY
%%%%%%%%%%%%%%%%%%%%%

% Generate indices
nq = n / 2;
V = -repmat((-nq:nq-1)', 1,2*nq);
U = repmat((-nq:nq-1), 2*nq,1);
f = sqrt(U.*U + V.*V);

%Generate magnitude
mag = power(f, -beta); 

%Covert to spatial domain
mag = ifftshift(mag);

%Shift & zero d.c.
mag(1,1) = 0; 

%Convert to cartesian
[x,y] = pol2cart(theta,mag);
F = x+1i*y;

%Inverse fft and take real part
f = ifft2(F, 'symmetric');
r = real(f); 

%Check to make sure there is no imaginary component as we are throwing it
%away
i = imag(f);
i = sum(sum(i));

if(i ~= 0)
    disp(' ');
    disp('ERROR: IMAGINARY COMPONENT IS NON-ZERO, CANNOT DISCARD!');
    disp(' ');
    return;
end

%Not really required as its a periodic signal, but helps visualisation when
%debugging.
ht = fftshift(r); 

%Normalise to [0:1] and return array of values (relating to pixels)
im = ht-min(min(ht));
im = im./max(max(im));

height = im;

end
