% 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]);
