Fractal Landscapes

Jim Mahoney, November 2004

Google on "Fractal landscape" or "Fractal terrain" to find references.
A few are
     www.gameprogrammer.com/fractal.hml 
     www.midnightbeach.com/jon/pubs/Fractal_Landscapes.html 
     www.fractal-landscapes.co.uk
The first two have local copies in cs/docs/fractal_terrain/

There are lots of variations and refinements on this scheme - including 
different ways to add color and "render" the images.

The program below uses the "squares and diamonds" approach with linear interpolation.
My first attempt used recursion but failed

Code

Off[General::spell1];

(* Set up global variables and draw the button. *)
init[nn_]:= Module[{},
	n = nn;											(* size of map to generate; should be a power of 2 plus 1 *)
	initialHeight = 1.0 (n-1);						(* arbitray choice *)
	terrain = Table[initialHeight,{n},{n}];	        (* blank n by n map *)
	];

(* Return a random number related to size *)
rand[size_]:= 2.0 Random[Real,{-size, size}];		(* ***** not sure what function this sould be *****  *)
	
(* Return the map value at the given point; if past map edges return the value at the edge *)
getTerrain[xx_, yy_]:= Module[{x=xx, y=yy},
	If[ x < 1 || x > n || y < 1 || y > n, Return[initialHeight]];
	Return[ terrain[[x,y]] ];
	];

(* average of 4 cells in a square; input is bottom left and size of square *)
squareAverage[x_, y_, size_]:=
 ( terrain[[x,y]] + terrain[[x+size,y]] + 
   terrain[[x,y+size]] + terrain[[x+size, y+size]] )/4.0;

(* average of 4 cells in a diamond; input is center and distance from center to vertex *)
diamondAverage[x_, y_, halfSize_]:=
 ( getTerrain[x, y+halfSize] + getTerrain[x, y-halfSize] +
   getTerrain[x+halfSize, y] + getTerrain[x-halfSize, y] )/4.0;

(* given the coordinates of the center of the diamond and the halfSize, set the terrain center accordingly *)
setDiamondCenter[x_, y_, halfSize_]:= terrain[[x,y]] = diamondAverage[x, y, halfSize] + rand[1.414 halfSize];

(* This was my second attempt:
	explicit loops, 1st over squares, then over diamods, at smaller and smaller scales. *)
doLoops := Module[{x,y,size,halfSize},
	size = n-1;
	While[ size > 1,
		halfSize = size/2;
		Do[
			(* Print[" x=",x, "; y=",y, "; size=",size]; *)   (* debugging *)
			terrain[[ x+halfSize, y+halfSize ]] = squareAverage[x,y,size] + rand[size],
			{x,1,n-1,size},
			{y,1,n-1,size}
			];
	
		Do[
			setDiamondCenter[x+halfSize, y, halfSize];
			setDiamondCenter[x, y+halfSize, halfSize];
			setDiamondCenter[x+size, y+halfSize, halfSize];
			setDiamondCenter[x+halfSize, y+size, halfSize],
			{x,1,n-1,size},
			{y,1,n-1,size}
			];
		size = halfSize;
		];
	];

(* This was my original attempt - which failed...
	Generate fractal terrain in terrain[[]] starting at (x,y) for square with given size.
   This routine is called recursively to generate the terrain at smaller and smaller scales.
   BUT IT GIVES ARTIFACTS ! *)
fractalTerrain[x_, y_, size_] := Module[{halfSize}, 
	If[ size==1, Return[] ];
	halfSize = size/2;
	terrain[[ x+halfSize, y+halfSize ]] = squareAverage[x,y,size] + rand[size];
	setDiamondCenter[x+halfSize, y, halfSize];
	setDiamondCenter[x, y+halfSize, halfSize];
	setDiamondCenter[x+size, y+halfSize, halfSize];
	setDiamondCenter[x+halfSize, y+size, halfSize];
	fractalTerrain[x, y, halfSize];
	fractalTerrain[x+halfSize, y, halfSize];
	fractalTerrain[x, y+halfSize, halfSize];
	fractalTerrain[x+halfSize, y+halfSize, halfSize];
	];	
	
(* Show me what you got... *)
displayTerrain := Module[{},
	terrain = terrain / Max[terrain];
	Show[Graphics[Raster[terrain]], AspectRatio->1]; ListPlot3D[terrain];
	];

(* THIS ONE DOESN'T WORK *)
(* fractalTerrain[1, 1, n-1]; *)

(* THIS ONE WORKS! *)
(* doLoops; *)

Results

Worked

big, n=129

n = 129
Using the explicit doLoop.
Looks like that fixed the artifacts - this one looks pretty smooth.       Not bad...

init[129]; doLoops;  displayTerrain;

small, n=65

init[65]; doLoops; displayTerrain;

Didn't work

n = 129, 
Using the recursive version.
Problem : horizontal and diagonal artifacts
I think the problem may be in the order of evaluation - perhaps doing an explicit loop rather than recursive would work better.
This way, I think we're going deep and looking for pieces beyond the edge of current square before doing all the squares at a given level.

init[129];  fractalTerrain[1, 1, 128];  displayTerrain;

Triangles

Here's an idea that I haven't finished.

Start with a trianglular polygon.
Find three new points by taking the midpoints of each line, plus a random amount.
Use these six points to make four smaller triangular polygons.
Repeat on smaller and smaller scales, and VOILA!
The details get a bit tricky.
Here's some playing around in that direction.

triangle={{0. , 0. , 0. }, {0. , 1. , 1. }, {1. , 
          0. ,  .5}}

Show[Graphics3D[Polygon[triangle]]];

If we have a list of such triangles, then we can use the Map function to apply Polygon to each one and plot 'em all.

another =  a/2 + 1;     
Show[
      Graphics3D[Map[Polygon, {triangle, another}]]];

(* 
      This routine takes one triangle (e . g .  
              a list of 3 points) as input,  
and returns a list of four triangles as output. *) 
morphTriangle[p_] := 
        Module[{one, two, three, onetwo, twothree, 
            threeone}, 
one = 
            p[[1]]; 
two = 
            p[[2]]; 
three = 
            p[[3]]; 
onetwo  =  (p[[1]] + 
                    p[[2]])/2  +  {0, 0, 
                Random[Real, {0, 
                    1}]}; 
twothree  =  (p[[2]] + 
                    p[[3]])/2  +  {0, 0, 
                Random[Real, {0, 
                    1}]}; 
threeone  =  (p[[3]] + 
                    p[[1]])/2  +  {0, 0, 
                Random[Real, {0, 1}]}; 
Return[{{one, 
                onetwo, threeone}, {two, twothree, onetwo}, {three, threeone, 
                twothree}, {onetwo, twothree, 
                threeone}}];
];

Show[
        Graphics3D[Map[Polygon,  morphTriangle[triangle]]]];

triangleList = 
      morphTriangle[triangle]

Dimensions[triangleList]

Dimensions[
      Flatten[Map[morphTriangle, triangleList], 1]]

triangleList = 
        Flatten[Map[morphTriangle, triangleList], 1];

Show[
        Graphics3D[Map[Polygon, triangleList]]];

Yeah, it's a mess - the randomness is too big, and always positive.  Not good.
So all we'd need do is change the size of the random numbers on each pass, and do this repeatedly...

(* 
      Putting all that together : *) 
scale = 1.0;

    firstTriangle  =  {{0. , 0. , 0. }, {1. , 0. , 0. }, {0. , 1. , 
            0. }};

    theList  =  {firstTriangle};

    listFromTriangle[p_] := 
        Module[{one, two, three, onetwo, twothree, 
            threeone}, 
 
          scale = scale  0.6; 
one = 
            p[[1]]; 
two = 
            p[[2]]; 
three = 
            p[[3]]; 
onetwo  =  (p[[1]] + 
                    p[[2]])/2  +  {0, 0, 
                Random[Real, {-scale, 
                    scale}]}; 
twothree  =  (p[[2]] + p[[3]])/2  +  {0, 0, 
                Random[Real, {-scale, 
                    scale}]}; 
threeone  =  (p[[3]] + p[[1]])/2  +  {0, 0, 
                Random[Real, {-scale, 
                    scale}]}; 
Return[{{one, onetwo, 
                threeone}, {two, twothree, onetwo}, {three, threeone, 
                twothree}, {onetwo, twothree, 
                threeone}}];
];

    Do[ 
        theList  =  
          Flatten[Map[listFromTriangle, theList], 
            1], {4}];

    Show[Graphics3D[Map[Polygon, theList]], 
        ViewPoint -> {-1.190,  -1.485,  
            1.327}];

This has several flaws.
	* The random number generator isn't scaling all that well.  
	* We're pulling the edges away from each other - one of the articles I reference at the top of the page discusses this.

Still, this is the kind of approach that will recursively do smaller triangles.