% Create a randomized height field, and display a surface fit to it
% Calls: FillWithWater, Odd

% method: The method used to create the randomized height field
% sz: The height field is sz x sz
% iters: Number of iterations the height field generator runs for
% randm: Amount of randomness added to the CA height field algorithm
% peak: Amount of jaggedness of the higher peaks
% decay: Amount of decay in the weight of successive iterations of the CA
% jag: Jaggedness of the landscape
% water: Amount of water in the landscape
% customBuild: Whether or not the current landscape is custom built
function [sz,iters,dim] = Build(method,sz,iters,randm,peak,decay,dim,jag,water,customBuild);

clc;
cla;

% Available methods: simple, fractl, subdiv, spline, spplus, random,
% randsp, rndpls, custom

% Default values
if ~customBuild | (method == 'subdiv')
jag = 1;
water = 0;
end

% All methods but fractl use a CA to calculate the height field
if (method ~= 'fractl') | (method == 'normal')

% Set all edge vertices to 0
% Randomly set all inner vertices to either 0 or 1
CA = randn(sz-2);
CA = [zeros(1,sz-2);CA;zeros(1,sz-2)];
CA = [zeros(sz,1) CA zeros(sz,1)];
for i = 2:sz-1
for j = 2:sz-1
CA(i,j) = (CA(i,j) >= 0);
end
end

% Run the CA and save the sum in Height
Height = CA;
Id = eye(sz);
up = [zeros(sz,1) Id(:,1:sz-1)];
down = up';
for k = 1:iters
upCA = up * CA;
downCA = down * CA;
NewCA = CA + upCA + downCA + CA*up + CA*down + upCA*up + upCA*down + downCA*up + downCA*down;
NewCA = (NewCA > 5) | (NewCA == 4);
CA = xor(NewCA,(rand(sz) < (randm/100)));
Height = Height + CA / (1+decay*k/5);
end

% Calculate an estimate of the fractal dimension of the landscape
c = 0;
dim(1) = 0;
maxk = floor(log2(sz+1));
perturbs = zeros(2^maxk);
for k = 1:maxk
size = 2^k;
oldSize = 1 / size;
half = (sz+1) / size;
for i = 1:2:size
for j = 1:2:size
midi = i * half;
midj = j * half;
perturbs(i,j) = ((Height(ceil(midi),ceil(midj))+Height(ceil(midi),floor(midj))+Height(floor(midi),ceil(midj))+Height(floor(midi),floor(midj))) / 4) / (size * 2);
newSize = sqrt(oldSize^2 + perturbs(i,j)^2);
newDim = log(4) / log(2*oldSize/newSize);
if (newDim > 0) & (newDim < 10)
dim(1) = dim(1) + log(4) / log(2*oldSize/newSize);
c = c + 1;
end
end
end
end
dim(1) = dim(1) / c;
dim(1) = (dim(1) + 2.5 + decay/10) / 2;
if length(dim) > 1
set(dim(3),'String',num2str(dim(1)));
end
end

if method == 'subdiv'
% Square the number of vertices by subdividing the matrices
X = zeros(sz*(sz-1)+1);
Y = zeros(sz*(sz-1)+1);
for k = 1:sz*(sz-1)+1
pos = (k-1) / (sz*(sz-1));
X(:,k) = pos;
Y(k,:) = pos;
end
Z = zeros(sz*(sz-1)+1);

for i = 1:sz
% Copy original height of first component of each row
Z((i-1)*sz+1,1) = Height(i,1);
for j = 2:sz
% Copy heights of rest of original components
Z((i-1)*sz+1,(j-1)*sz+1) = Height(i,j);
for k = 1:sz-1
% Calculate heights of inner components of each row by interpolation
Z((i-1)*sz+1,(j-2)*sz+1+k) = Height(i,j-1) + k/sz*(Height(i,j)-Height(i,j-1));
end
end
end

% Calculate heights of inner components of each column by interpolation
for j = 1:sz*(sz-1)+1
for i = 2:sz
for k = 1:sz-1
Z((i-2)*sz+1+k,j) = Z((i-2)*sz+1,j) + k/sz*(Z((i-1)*sz+1,j)-Z((i-2)*sz+1,j));
end
end
end

% sHeight is a miniaturized version of the original height field
sHeight = Height / sz * jag;

% Duplicate sHeight (sz-1) x (sz-1) times, and store the result in aug
augtemp = sHeight(1,:);
for i = 1:sz-1
augtemp = [sHeight;augtemp];
end
aug = augtemp(:,1);
for j = 1:sz-1
aug = [augtemp aug];
end

% Artificially fractalize the height field by superimposing the miniaturized
% versions of the height field on top of the original version
Z = Z + aug;

elseif (method == 'spline') | (method == 'spplus') | (method == 'random') | (method == 'rndpls') | (method == 'randsp')

% divs is the number of subdivisions to make between each pair of original vertices
if (method == 'spline')
divs = ceil(50/sz);
else
divs = sz;
end

% Increase the number of vertices by subdividing the matrices
size = (sz-1)*divs+1;
X = zeros(size);
Y = zeros(size);
for k = 1:size
pos = (k-1) / (size-1);
X(:,k) = pos;
Y(k,:) = pos;
end
Z = zeros(size);

% Fit splines to the original height field
oldbase = linspace(0,1,sz);
newbase = linspace(0,1,size);
for i = 1:sz
Z((i-1)*divs+1,:) = spline(oldbase,Height(i,:),newbase);
end
for j = 1:size
Z(:,j) = spline(oldbase',Z(1:divs:size,j),newbase');
end

if (method == 'random') | (method == 'rndpls')
for i = 1:sz
% Copy original height of first component of each row
Z2((i-1)*sz+1,1) = Height(i,1);
for j = 2:sz
% Copy heights of rest of original components
Z2((i-1)*sz+1,(j-1)*sz+1) = Height(i,j);
for k = 1:sz-1
% Calculate heights of inner components of each row by interpolation
Z2((i-1)*sz+1,(j-2)*sz+1+k) = Height(i,j-1) + k/sz*(Height(i,j)-Height(i,j-1));
end
end
end

% Calculate heights of inner components of each column by interpolation
for j = 1:sz*(sz-1)+1
for i = 2:sz
for k = 1:sz-1
Z2((i-2)*sz+1+k,j) = Z2((i-2)*sz+1,j) + k/sz*(Z2((i-1)*sz+1,j)-Z2((i-2)*sz+1,j));
end
end
end

% For each vertex, randomly choose to use either the splined height or the
% interpolated height
for i = 1:sz*(sz-1)+1
for j = 1:sz*(sz-1)+1
if randn(1) > 0
Z(i,j) = Z2(i,j);
end
end
end
end

if (method == 'spplus') | (method == 'rndpls') | (method == 'randsp')
% sHeight is a miniaturized version of the original height field
if method == 'rndpls'
sHeight = Height / divs * jag;
else
sHeight = Height / divs;
end

% Duplicate sHeight (divs-1) x (divs-1) times, and store the result in aug
augtemp = sHeight;
for i = 1:divs-1
augtemp = [augtemp;sHeight(2:sz,:)];
end
aug = augtemp;
for j = 1:divs-1
aug = [aug augtemp(:,2:sz)];
end

% Artificially fractalize the height field by superimposing the miniaturized
% versions of the height field on top of the original version
if method == 'spplus'
Z = Z + aug;

% For each vertex, randomly choose whether or not to artificially fractalize
% the height field by superimposing the miniaturized versions of the height
% field on top of the original version
else
for i = 1:sz*(sz-1)+1
for j = 1:sz*(sz-1)+1
if (randn(1) > 0) & (Z(i,j) > 0)
Z(i,j) = Z(i,j) + aug(i,j);
end
end
end
end
end
Z = FillWithWater(Z);

elseif method == 'fractl'
oldsz = sz;
c = 0;
dim(1) = 0;

% Set all edge vertices to 0
% Randomly set all inner vertices to either 0 or 1
Z = zeros(sz);
Z(2:sz-1,2:sz-1) = (rand(sz-2)-0.2)./sz;

% Create a fractal height field by randomly perturbing the field by smaller
% and smaller amounts
for k = 2:iters
sz = sz * 2 - 1;
Znew = zeros(sz);
oldSize = 1 / (2*sz);
for i = 1:sz
for j = 1:sz
iold = (i+1)/2;
jold = (j+1)/2;
if Odd(i) & Odd(j)
Znew(i,j) = Z(iold,jold);
elseif Odd(i) & ~Odd(j)
Znew(i,j) = (Z(iold,floor(jold)) + Z(iold,ceil(jold))) / 2;
elseif ~Odd(i) & Odd(j)
Znew(i,j) = (Z(floor(iold),jold) + Z(ceil(iold),jold)) / 2;
else
% Add in the random perturbation and calculate the resulting fractal dimension
c = c + 1;
dy = max([abs(Z(ceil(iold),ceil(jold))-Z(floor(iold),floor(jold))) abs(Z(ceil(iold),floor(jold))-Z(floor(iold),ceil(jold)))]);
ang = (pi/2) - atan(dy/(sqrt(2)*(1/((sz-1)/2))));
perturb = (rand(1)-0.5) / sz;
newSize = sqrt(oldSize^2 + (sin(ang)*perturb)^2);
dim(1) = dim(1) + log(4) / log((1/sz)/newSize);
Znew(i,j) = (Z(floor(iold),floor(jold)) + Z(floor(iold),ceil(jold)) + Z(ceil(iold),floor(jold)) + Z(ceil(iold),ceil(jold))) / 4 + perturb;
end
end
end
Z = Znew;
end
dim(1) = dim(1) / c;
set(dim(3),'String',num2str(dim(1)));

X = zeros(sz);
Y = zeros(sz);
for k = 1:sz
pos = (k-1) / (sz-1);
X(:,k) = pos;
Y(k,:) = pos;
end
Z = FillWithWater(Z);
sz = oldsz;

% Use the original, unmodified height field to create the landscape
else
X = zeros(sz);
Y = zeros(sz);
for k = 1:sz
pos = (k-1) / (sz-1);
X(:,k) = pos;
Y(k,:) = pos;
end
Z = Height;
end

% C is the colormap used for the landscape
C = zeros(100,3);
for i = 1:2
C(i,:) = [.3 .1 1];
end
for i = 3:7
C(i,:) = [.9-.02*(i-2) .76-.046*(i-2) .2];
end
for i = 8:14
C(i,:) = [.8-.1*(i-7) .53 .2];
end
for i = 15:40
C(i,:) = [0 .53+.01*(i-14) .2];
end
for i = 41:79
C(i,:) = [.02*(i-40) .8+.0025*(i-40) .2+.005*(i-40)];
end
for i = 80:96
C(i,:) = [.8+.006*(i-79) .9+.004*(i-79) .4+.033*(i-79)];
end
for i = 97:100
C(i,:) = [.908 .972 1];
end

% If water was increased from 30 percent, submerge the height field by
% the appropriate amount
maxHeight = max(max(Z));
Z = FillWithWater(Z - maxHeight * water);

% Increase the jaggedness of the higher peaks by a factor of peak
if peak > 0
exponent = peak / 6 + 1;
[M,N] = size(Z);
for i = 1:M
for j = 1:N
curHeight = 2 * Z(i,j) / maxHeight;
Z(i,j) = curHeight^exponent;
end
end
end

% Plot the landscape
surf(X,Y,Z);
colormap(C);
shading interp;
maxHeight = max(max(Z));
axis([0 1 0 1 0 maxHeight*2]);