3-D Terrain Generation

Taran Besant & Spencer Eanes

Modern computers make the task of procedurally generating realistic 3-D terrain feasible, and many different methods have been created to do this. Procedural generation techniques for graphics are commonly used in applications like video game graphics, movie graphics, and more. The goal of this project is to implement and compare different methods of procedurally generating 3-D terrain. Specifically, we will look at the midpoint method, and the

Helper Functions

Make Array

This function takes a dimension, and creates a dim x dim dimension array of three-tuples with z value 0 and x and y values spaced evenly on x and y axis at intervals of 1.

makeArr[dim_]:=Module[{data,i,j,ctr=0},
data = ConstantArray[0,{dim*dim,3}];
Do[
data[[j]][[1]]=Mod[ctr,dim];
data[[j]][[2]]=Floor[ctr/dim];
ctr++,
{j,1,dim*dim}
];
data
]

dim =3;
data =makeArr[dim]
ListPointPlot3D[data]

So, what happens if we assign the points values randomly? 
This function takes arguments of the square input array and the dimension of the array. It runs a Do-loop and assigns the third value of the three-tuple to a random real in the range [0,1];
It returns the array once all modifications have been made.

randomZ[arr_,dim_]:=Module[{data=arr,j},
Do[
data[[j]][[3]]=RandomReal[],
{j,1,dim*dim}
];
data
]

Plot Heat Map

This function takes the input array, assuming it is an array of three-tuples and creates a 3d plot with a heat map surface to easier visualize comparative heights.

plotHeatMap[arr_]:=Module[{},
ListPlot3D[arr,Filling->Bottom, ColorFunction->(ColorData[{"TemperatureMap","Reverse"},#3]&),Boxed->False,PlotRange->Automatic]
]

Use randomZ function and plot it with plotHeatMap to see what the result looks like.

dim = 5;
data=makeArr[dim];
randomData = randomZ[data,dim]
plotHeatMap[randomData]

This does work, however it is not ideal, because points have no relation to their neighbors. For instance, a local maxima could be next to a local minima, which is not very realistic. On a larger scale, this issue becomes more obvious. It is clear that nearby points have no relation to their neighbors, and this creates a surface that is simply not useful (unless maybe you’re trying to model grass).

data2 = makeArr[25];
data2 = randomZ[data2,25];
ListPointPlot3D[data2,Filling->Bottom]

Okay, so let’s implement some methods and see how much better they are.

Assign Array

Given two arrays with the same number of elements, assign the values of one to the third value of the other (to make a plot-able structure).

assignArr[threeTuples_,values_,dim_]:=Module[{toDraw=threeTuples,ctr=1},
Do[
Do[
toDraw[[ctr]][[3]]=values[[i]][[j]];
ctr++,
{j,1,dim}
],
{i,1,dim}
];
toDraw
]

Main Methods

MidPoint Method

This method sets the corners randomly, then sets the middles and the middle edge pieces to the average of their surrounding pieces, and then recursively applies this to all sub squares.
Note that this will only work nicely if square dimensions are of the form (2^n)+1.

Next, we will set the midpoints of the edges to the averages of the two corners that they are between.

Introduction

This method sets the corners randomly, then sets the middles and the middle edge pieces to the average of their surrounding pieces, and then recursively applies this to all sub squares.
Note that this will only work nicely if square dimensions are of the form (2^n)+1.


The above image shows how the corners start with random values, which are then used to set the center. Then the same is done with a smaller subsquare of dimension (n+1)/2 where n is the dimension of the previous, larger square. This subsquare recurses (2^n)*(2^n) times over the array, where n is the exponent in the original (2^n)+1 equation used to generate the size of the array. This is why the midpoint method only works with arrays of dimension (2^n)+1. On a higher order array, a 9x9 for example, the recursion looks like this:


Where the first “pass” is the midpoints of the corners and the very middle of the array being set, and consequent passes are the subsquares doing the same in smaller and smaller increments until a 3x3 square fills in the final values.

First Attempt

@arg arr - A square dim x dim array
@arg dim - The dimension of the square array
@arg xpon - The exponent of the (2^n)+1 equation that generates the square matrix

midDisp[arr_,dim_,xpon_]:=Module[{data=arr,i=dim,topMid,leftMid,rightMid,bottomMid,veryMid,count=0,j,k=1,n=1,m=0,l=1,o=1, jitter},
(*First set the corners to random real numbers, and have a jitter function to add random noise for later*)
jitter=(RandomReal[]-.5)*.15;
data[[1]][[1]]=RandomReal[];
data[[1]][[i]]=RandomReal[];
data[[i]][[1]]=RandomReal[];
data[[i]][[i]]=RandomReal[];
(*Then set the midpoints to the averages of the corners they are between*)
topMid=(data[[1]][[1]]+data[[1]][[i]])/2;
leftMid=(data[[1]][[1]]+data[[i]][[1]])/2;
rightMid=(data[[1]][[i]]+data[[i]][[i]])/2;
bottomMid=(data[[i]][[1]]+data[[i]][[i]])/2;
data[[1]][[(i+1)/2]]=topMid;
data[[(i+1)/2]][[1]]=leftMid;
data[[(i+1)/2]][[i]]=rightMid;
data[[i]][[(i+1)/2]]=bottomMid;
(*Do the same for the very center, but averaging the midpoints*)
veryMid=(topMid+leftMid+rightMid+bottomMid)/4;
data[[(i+1)/2]][[(i+1)/2]]=veryMid;
(*Iterate across every other point in progressively smaller squares, making sure not to set the same point several times throughout the iteration*)
For[j=0,n< xpon,n++,
k=1;
i=(i+1)/2;
m=0;
o=1;
While[o≤ 2^(n),
j=0;

While[l≤ 2^(n),
topMid=(data[[k+m]][[k+j]]+data[[1]][[i]])/2+jitter;
leftMid=(data[[k+m]][[k+j]]+data[[i]][[1]])/2+jitter;
rightMid=(data[[k+m]][[i+j]]+data[[i]][[i]])/2+jitter;
bottomMid=(data[[i+m]][[k+j]]+data[[i]][[i]])/2+jitter;
veryMid=(topMid+leftMid+rightMid+bottomMid)/4+jitter;
If[data[[k+m]][[(i+1)/2+j]]==0,
data[[k+m]][[(i+1)/2+j]]=topMid;,
];(*If #1, these if statements prevent the function from setting the same midpoint value several times throughout iteration, which was causing problems with oddly large mountains or very shallow valleys*)
If[data[[(i+1)/2+m]][[k+j]]==0,
data[[(i+1)/2+m]][[k+j]]=leftMid;,
];(*If #2*)
If[data[[(i+1)/2+m]][[i+j]]==0,
data[[(i+1)/2+m]][[i+j]]=rightMid;,
];(*If #3*)
If[data[[i+m]][[(i+1)/2+j]]==0,
data[[i+m]][[(i+1)/2+j]]=bottomMid;,
];(*If #4*)
If[data[[(i+1)/2+m]][[(i+1)/2+j]]==0,
data[[(i+1)/2+m]][[(i+1)/2+j]]=veryMid;,
];(*If #5*)
l++;
j=j+i-1;
];(*While loop, this loop effectively shifted the subsquares to the right across the array until the edges were set*)
o++;
l=1;
m=m+i-1;
];(*While loop, this loop shifted the subsquared down across the array and reset the value that was shifting it to the right every time the subsquare set the rightmost edge*)
];(*For loop, this loop made the subsquares smaller and resets the variables counting how many times the while loops need to iterate*)

data
]

dim=5;
midArray=ConstantArray[0,{dim,dim}];
firstTest=midDisp[midArray,5,2];
ListPlot3D[firstTest,Filling->Bottom, ColorFunction->(ColorData[{"TemperatureMap","Reverse"},#3]&),PlotRange->{0,1}]

This tends to generate fine-looking, if plain and smooth, terrain at small dimensions like 5x5 and 9x9, but at higher dimensions something went wrong and the terrain began to look... unrealistic.

dim=17;
midArray=ConstantArray[0,{dim,dim}];
firstTest=midDisp[midArray,17,4];
ListPlot3D[firstTest,Filling->Bottom, ColorFunction->(ColorData[{"TemperatureMap","Reverse"},#3]&)]

Random “potholes” of very low values and similar instances of repreated higher values appeared. Troubleshooting became a nightmare so a new approach seemed prudent.

Mid Edges

@arg arr - the input array
@arg dim - the dimensionality of the input array
@return - the modified array

midEdge[arr_,dim_]:=Module[{toCalc,harr=arr},
toCalc=Floor[(dim+1)/2];
harr[[1]][[toCalc]]=Mean[{harr[[1]][[1]],harr[[1]][[dim]]}]; (*top*)
harr[[toCalc]][[1]]=Mean[{harr[[1]][[1]],harr[[dim]][[1]]}];(*left*)
harr[[toCalc]][[dim]]=Mean[{harr[[1]][[dim]],harr[[dim]][[dim]]}];(*right*)
harr[[dim]][[toCalc]]=Mean[{harr[[dim]][[1]],harr[[1dim]][[dim]]}]; (*bottom*)

Return@harr;
]

Middle

Same params and return as above, + 
@arg jitter - coefficient to scale the jitter amount.

middle[arr_,dim_,jitter_]:=Module[{midPoint,dist,tupVal,harr=arr,b},
b=(RandomReal[]-.5)*jitter;
(*b=RandomVariate[NormalDistribution[]];*)(*we would be interested in looking more into normal distibution jitter in the fiture*)
tupVal=Floor[(dim+1)/2];
dist = tupVal-1;
(*take the mean of the four surrounding points and set the middle point to that value plus some jitter*)
harr[[tupVal]][[tupVal]]=Mean[{harr[[tupVal-dist]][[tupVal]],harr[[tupVal+dist]][[tupVal]],harr[[tupVal]][[tupVal-dist]],harr[[tupVal]][[tupVal+dist]]}]+b;
Return@harr
]

caller

Recursive function calls midEdge and middle on the sub-square arrays of the larger array.
The base case is when dim = 3, once this is completed, the whole thing is completed. 
To accomplish this, we apply midEdge and middle to the input array, then recurse on the four subarrays which correspond to the four quadrants of the bigger array.
Once each step of recursion is complete, we recombine the four quadrants back into the new filled array. 
Since it is recursive, we have to modify the input argument, arr. To do this, we give the function the HoldAll attribute, which stops the input arguments from being evaulated immediately. This allows us to modify them later in the code. 

@param arr - the input array
@param dim - the dimension of the array (x and y dimensions are always the same, so we only need one argument for this).
@param jitter - the jitter coefficient to pass to the middle function.
@return - none, because the input array arr itself is modified.

SetAttributes[caller,HoldAll];
caller[arr_,dim_,jitter_]:=Module[{smlrDim,tempArr1,tempArr2,tempArr3,tempArr4},

smlrDim=Floor[(dim+1)/2];

arr=midEdge[arr,dim];
arr=middle[arr,dim,jitter];
(*Print["smlrDim=",smlrDim];
Print[arr//MatrixForm];*)

If[smlrDim≥3,
(*we have to use tempArr's to avoid over indexing arr, seems to be a quirk of the way this is working (just trust us on this, this was a big headache to figure out) *)
(*upper left subarray*)
tempArr1 = arr[[1;;smlrDim,1;;smlrDim]];
caller[tempArr1,smlrDim,jitter];
(*upper right subarray*)
tempArr2 = arr[[1;;smlrDim,smlrDim;;dim]];
caller[tempArr2,smlrDim,jitter];
(*lower left subarray*)
tempArr3 = arr[[smlrDim;;dim,1;;smlrDim]];
caller[tempArr3,smlrDim,jitter];
(*lower right subarray*)
tempArr4 = arr[[smlrDim;;dim,smlrDim;;dim]];
caller[tempArr4,smlrDim,jitter];

(*recombine the arrays*)
arr[[1;;smlrDim,1;;smlrDim]]=tempArr1;
arr[[1;;smlrDim,smlrDim;;dim]]=tempArr2;
arr[[smlrDim;;dim,1;;smlrDim]]=tempArr3 ;
arr[[smlrDim;;dim,smlrDim;;dim]]=tempArr4;
];

]

midPoint

The method that we actually call to perform the midpoint method. 
@param n - the power that is used to calculate the size of the matrix
@param jitter - the jitter coefficient to be passed along the middle method.
@return - the filled array.

midPoint[n_,jitter_]:=Module[{dim=(2^n)+1,arr,tupVal,sml},
(*Dim x Dim array of 0's*)
arr=ConstantArray[0,{dim,dim}];

(*set corners*)
arr[[1]][[1]]=RandomReal[]; (*upper left*)
arr[[1]][[dim]]=RandomReal[]; (*upper right*)
arr[[dim]][[1]]=RandomReal[];(*lower left*)
arr[[dim]][[dim]]=RandomReal[];(*lower right*)

(*initiate caller, which is recursive*)
caller[arr,dim,jitter];

Return@arr
]
midPoint[3]//MatrixForm

Diamond-Square Method

Introduction

The diamond-square method is similar to the midpoint method, but it calculates every point by taking an average of four (or three) surround points. This is better than the midpoint method because having only two points can lead to “seams” in the output which are obvious. This method works by first setting the corners of the array randomly, then performing diamond and square steps iteratively.

In the terminology of the above image, the diamond step sets the midpoint of a square as the average of the four corners plus some small amount of jitter. The square step sets the edge pieces as the average of their three or four surrounding pieces. 
Please note that our functions are named opposite of how they are in the above image taken from Wikipedia because that made more sense to us. That is, our midSquare method is the same as the above “diamond step”, and our diamondAround is the same as the “square step” above.

Mid Square

Calculate the mid point of the square with given dimensions, fill it into the given array and return the array.
@arg arr - a square dim x dim input array
@arg dim - the dimension of the square input array
@arg pos - a 2-tuple for the position in the square that we are computing the method on
@arg dist - how far away from pos to take the averages on.
	for instance, if we were computing midSquare in the first step for the very middle of the square, pos = {Floor[(dim+1)/2],Floor[(dim+1)/2]} and dist = Floor[dim+1];
@arg jitter - the jitter coefficient that scales the random amount of jitter. Typically should be between 0 and 1 (treated as a percent), but could be outside those ranges.
@return the modified array once the midpoint has been changed.

midSquare[arr_,dim_,pos_,dist_,jitter_]:=Module[{harr=arr,a,b},
(*mean of the surrounding four points*)
a=Mean[{arr[[pos[[2]]-dist]][[pos[[1]]-dist]],arr[[pos[[2]]+dist]][[pos[[1]]-dist]],arr[[pos[[2]]-dist]][[pos[[1]]+dist]],arr[[pos[[2]]+dist]][[pos[[1]]+dist]]}];
(*random jitter values*)
b=(RandomReal[]-.5)*jitter;(*random jitter*)
harr[[pos[[1]]]][[pos[[2]]]] =a+b;
Return@harr
]

Mid Diamond

Calculate the mid point of each “diamond” shape using three or four surrounding points. If any point is out of bounds of the matrix, it is not used.
The if statements make sure that if a point outside of the matrix is requested, it is not used.
Note that this is only one fourth of the diamondStep implemented in diamondAround. 
Arguments are return the same as midsquare.

midDiamond[arr_,dim_,pos_,dist_]:=Module[{harr=arr,b=(RandomReal[]-.5)*.2,p={}},
(*left*)
If[pos[[1]]>dist,
AppendTo[p,harr[[pos[[2]]]][[pos[[1]]-dist]]];
];
(*right*)
If[dim≥pos[[1]]+dist,
AppendTo[p,harr[[pos[[2]]]][[pos[[1]]+dist]]];
];
(*top*)
If[pos[[2]]>dist,
AppendTo[p,harr[[pos[[2]]-dist]][[pos[[1]]]]];
];
(*bottom*)
If[dim≥pos[[2]]+dist,
AppendTo[p,harr[[pos[[2]]+dist]][[pos[[1]]]]];
];
(*Print["p=",p];
Print[Mean[p]];*)
harr[[pos[[2]]]][[pos[[1]]]]=Mean[p];
Return@harr
]

Diamond Around

Calls midDiamond four times for the four surround points that are dist away from the point as pos.
Arguments and return the same as the previous two.

diamondAround[arr_,dim_,pos_,dist_]:=Module[{harr=arr},
(*bottom*)
harr=midDiamond[harr,dim,{pos[[1]]+dist,pos[[2]]},dist];
(*top*)
harr=midDiamond[harr,dim,{pos[[1]]-dist,pos[[2]]},dist];
(*right*)
harr=midDiamond[harr,dim,{pos[[1]],pos[[2]]+dist},dist];
(*left*)
harr=midDiamond[harr,dim,{pos[[1]],pos[[2]]-dist},dist];
Return@harr
]

diamondSquare

The full implementation of the diamondSquare method for an array of size (2^n)+1. We force this array size to make sure that all subarrays we calculate for are of odd size so they all have a mid point.
@arg n - the power of 2, so the overall array output is of size (2^n)+1
@return the fully filled array

diamondSqure[n_,jitter_]:=Module[{arr,dim=(2^n)+1,dist,twople,tempTup,tupVal,count},
arr=ConstantArray[0,{dim,dim}];(*create 2d array of 0's*)
(*Assign the four corners randomly*)
arr[[1]][[1]]=RandomReal[]; (*upper left*)
arr[[1]][[dim]]=RandomReal[]; (*upper right*)
arr[[dim]][[1]]=RandomReal[];(*lower left*)
arr[[dim]][[dim]]=RandomReal[];(*lower right*)

(*setup dist and twople*)
dist = Floor[(dim+1)/2]; (*input for dist arg in midSquare and diamondAround methods*)
tupVal=Floor[(dim+1) /2]; (*coordinate value for twople*)
twople={tupVal,tupVal};(*the starting point we calculate around (the midpoint of subsquare we iterate over)*)
count = tupVal; (*how far to step forward from each coord of twople*)

(*iterate...*)
(*the while loop has this bound because once we calculate on twople={2,2} with dist = 2 (aka all 3x3 subsquare), the method is complete*)
(*it calls the inner loops for progressively smaller subsquares until every space has been filled*)
While[tupVal>1,

(*first double loop to loop over x and y values to add to twople*)
(*can't combine the two loops because all midSquare values must be computed before any diamondAround values can be calculated*)
Do[
Do[
arr=midSquare[arr,dim,{twople[[1]]+y,twople[[2]]+x},tupVal-1,jitter],
{x,0,dim-twople[[2]],count}
],
{y,0,dim-twople[[1]],count}
];

(*second double loop, same bounds and ragne as the previous, but this is for diamondAround*)
Do[
Do[
arr=diamondAround[arr,dim,{twople[[1]]+y,twople[[2]]+x},tupVal-1],
{x,0,dim-twople[[2]]-1,count}
],
{y,0,dim-twople[[1]]-1,count}
];

(*increment variables properly*)
count = tupVal-1;
tupVal=Floor[(tupVal+1)/2];
twople[[1]]=tupVal;
twople[[2]]=tupVal;
];

Return[arr]
];

Test to make sure that the method fills every possible point.

diamondSqure[3]//MatrixForm

Use it

Methods

Diamond Square

val=5;(*size of output array*)
dim = (2^val)+1
arr=makeArr[dim];
toDraw=diamondSqure[val,.25];(*change the second arg to adjust jitter. vals between 0 and 1 are best*)
toDraw//MatrixForm;
toDraw=assignArr[arr,toDraw,dim];
plotHeatMap[toDraw]

Midpoint

val=5;
dim = (2^val)+1
arr=makeArr[dim];
toDraw=midPoint[val,.25];
toDraw//MatrixForm;
toDraw=assignArr[arr,toDraw,dim];
plotHeatMap[toDraw]

Conclusions

As can be easily seen, these are much nicer than the original random assignment, because points affect their neighbors to create a nice structure, which some random jitter. Try messing with the jitter variable to see different results

Timing

midpointTiming={};
diamondSquareTiming={};
Do[
AppendTo[midpointTiming,AbsoluteTiming[midPoint[n,.25]][[1]]];
AppendTo[diamondSquareTiming,AbsoluteTiming[diamondSqure[n,.25]][[1]]],
{n,1,8}
]
midpointTiming

ListPlot[{midpointTiming,diamondSquareTiming}]

diamondSquareTiming

In the end, the diamondSquare method is slower, but generates “better” terrain. However, in our currently implementation both appear to take exponential time, so super large plots are not very feasible.