Classical elements in NetLogo: Earth
In this first post we will see one of the most classical algorithms for the formation of realistic landscapes: the mid point displacement algorithm. This algorithm handles only generate a map of heights above a defined level, then apply a flooding process, color and shading to create more realistic results. At the end of this post you will be able to generate things like this and use them for your 3D agent models:
The Mid Point Displacement Algorithm
The mid point displacement algorithm is a subdivision algorithm that generates a terrain iteratively, in each iteration the level of detail increases. From a computational point of view this algorithm is very inexpensive when compared to other terrain generation algorithms, and it provides a very impressive results. It is also known as a Plasma Algorithm, a set of algorithms that produce liquid effects when changing the colors in a cyclic over a prefixed palette of colors (probably, the first algorithms of this type were created by demo coders in the 90's by using some specific features of the display hardware from that time).
The general process is simple, and we will start with a 1D version to exemplfy the algorithm. In this simple case we start with an initial segment an then apply the next algorithm:
- For Each Segment
- Calculate the middle point of the segment
- Add a random value to the mid point
- Reduce the random range
- Repeat from step 1 until we reach the desired number of points / segment size
In NetLogo we can easily do this by using links and iterating on them in every step:
globals [ disp ; Displacement value for the mid points ] ; Create the initial conditions for the algorithm to start ca ; White background ask patches [set pcolor white] ; Create the 2 first points (in the borders), and a link (segment) ; between them crt 2 [ setxy (2 * who - 1) * max-pxcor random-ycor / 2 create-links-with other turtles [set color black] set color blue set shape "circle" ] ; Initial displacement value set disp max-pycor / 2 end ; Procedure to create the final curve: only a single loop to go repeat number-iterations [ iter ] ; For a better visualization, we hide the turtles at the end ask turtles [ht] end
The iteration procedure will pass through all the links/segments, calculate their mid points and displace them. Then, remove the original links. Note that we have to use a foreach iteration, and not only a ask links [...], because the creation of new turtles can't be done inside a link context. By convenience, we define an auxiliary report to return a random number in a symmetrical interval.
to iter foreach sort links [ ask (midpoint ?) [ set ycor ycor + ran disp] ask ? [die] ] ; Remove the original link after division
; In every iteratio the displacement is reduced according to roghness factor
set disp disp * 2 ^ (-1 * roughness) end ; Report that returns a random between -x and x to-report ran [x] report 2 * (random-float x) - x end
The midpoint of a link will return a turtle in the middle of the segment. As secondary effect, this new turtle will connect with the extremes of the segment.
to-report midpoint [s] let x 0 ; xcor of midpoint let y 0 ; ycor of midpoint ; Fill this values from the link, and remove it ask s [ set x ([xcor] of end1 + [xcor] of end2) / 2 set y ([ycor] of end1 + [ycor] of end2) / 2 ] ; Create one turtle in the middle point of e1 - e2 and report it let rep 0 crt 1 [ setxy x y ; Locate in the mid point create-links-with [both-ends] of s [; link it to e1 set color black ] set color blue set shape "circle" set rep self ] report rep end
Here you can try the model (or download it, if you want):
The same idea works on building 2D terrains. In this case we will need to split the iteration in two stages, in the first one we need to calculate the central point of every square (diammond phase), and later, we use this central point and the corners to calculate the midle points again (square phase):
This iteration produces a deeper resolution in the terrain formation:
Now, we will use a different approach. Instead of using a iterative segmentation of a simple square to get a grid of points, we prefix the points of the grid and fill them with the same method (note that, the 2D coordinates of the points are prefixed in any case). In order to make easier this iterative process, we will consider that the size of the world will be of the form \((2^s + 1) \times (2^s + 1)\), where \(s\) will be the number of iterations. In this way, with only 8 iterations we will obtain a terrain with \(257\times 257\) heights generated on a grid of patches (hence, be careful with the number of iterations or NetLogo will complain about the memory).
Also, we will have a roughness parameter that determines the look of the final terrain, from extremelly rough (0 value) to absolutely smooth (\(\infty\) value). Of course, we will keep this value between bounded values, for example between 0 and 3.
The code to perform this algorithm is given now. It starts with a square of maximum size \((2s+1)x(2s+1)\) points (\(2s\times 2s\) length), and in every iteration it will work on squares half sized. As we have seen in the previous pictures, we need two stages: a diamond phase to calculate the centers of current squares, and a square phase to interpolate the corners for the new squares. After every iteration, the displacement amount is changed by the roughness factor. More details can be found in the comments of the code.
After the filling of heights of patches is completed, we normalize them in order to have values in \([0,1]\).
patches-own [ height ; Where the height of the landscape will be stored, in [0,1] ] ; The Generate procedure follows the "Random midpoint displacement method" to generate ca ; Half size of the world, initial step size let s 2 ^ Pow-size ; The size of the world is (2s+1)x(2s+1) resize-world (- s) s (- s) s set-patch-size 400 / world-width ; Initial displacement let disp 1 ; Ratio of change of displacement in every iteration let disp-ratio 2 ^ (- r) ; Diamond-set and square-set will contain, respectively, the cross and corners ; for computing the points in the "Random midpoint displacement method" let diamond-set 0 let square-set 0 ; Initial diamod set: Central point of the world set diamond-set (patch-set patch 0 0) ; In order to decrease the time, the plasma will be shown after the complete ; calculation no-display ; In every iteration the step size will be half sized, the process is ; iterated while s >= 1 (the grid is integer indexed with patches positions) while [s >= 1] [ let -s (- s) ; relCorners stores the relative positions of the vertex for the squares ; in this iteration (the same size for all the squares) let relCorners (list (list -s s) (list s s) (list s -s) (list -s -s) ) ; relCross stores the relative positions of the vertex for the diamonds ; in this iteration (the same size for all the diamonds) let relCross (list (list 0 -s) (list s 0) (list 0 s) (list -s 0) ) ; First, the squares are computed for the diamonds set square-set (patch-set [patches at-points relcorners] of diamond-set) ; and there heights filled (by randomly displace the average of their vertex) ask square-set [ set height ram disp + mean [ height ] of patches at-points relcorners set pcolor scale-color white height -.5 1.5 ] ; Next, the diamonds are computed for the squares set diamond-set (patch-set [ patches at-points relCross ] of square-set) ; and there heights filled (by randomly displace the average of their vertex) ask diamond-set [ set height ram disp + mean [ height ] of patches at-points relCross set pcolor scale-color white height -.5 1.5 ] ; Finally, jump step and maximum displacement are adjuts for next iteration set s s / 2 set disp disp * disp-ratio ] ; Since the heights can be out of [0,1] interval, the heights must be normalized normalize ; After the new plas has been generated, it will be displayed display end ; Normalization procedure to put all the heights in [0,1] to normalize let hMax max [height] of patches let hMin min [height] of patches ask patches [ set height (height - hMin) / (hMax - hMin) set pcolor scale-color white height 0 1] end ; Generate a random float number in (-x,x) to-report ram [ x ] report random-float (2 * x) - x end
In this moment we have a system to generate height maps that can be used as a landscape generator. To provide a more final version, in the next procedures we will get more realistic lookings for this landscapes.
The first step is to flood the landscape to add some water. A simple procedure will give blue color to those patches under the desired water level. Of course, we could give more complex procedures, to create lakes or rivers at several heights, and they are easily programmed from the provided code.
to flood ask patches with [height <= water-level] [ set pcolor blue] ask patches with [height > water-level] [ set pcolor scale-color white height -.5 1.5] end
The next step is to give ome realistic color depending on the heights, to differentiate low lands, high lands and water. To get this effect we will use an extension, gradient, a little bit stronger than the scale-color function that is built-in NetLogo and that allows to define colors in a similar way (the color depending a numeric variable) but using different colors as intermediate points. Again, it is easy to change this gradient to get another landscape looking. In these processes, we use RGB colors rather than the original linear schema from NetLogo. (Note: Remember to add the line extensions [gradient] at the beginning of the code)
to real-color no-display ask patches [ ; Allow some mix in colors (diffuse borders) let v height + mix-colors * (random-float 1 - .5) / 10 ; Use one palette for water, and other for land ifelse height <= water-level [ set pcolor extract-rgb scale-color blue v 0 1 ] [ set pcolor gradient:scale (list extract-rgb (yellow + 2) extract-rgb (yellow - 2) extract-rgb (lime - 1.5) extract-rgb (green - 2) extract-rgb (brown - 1.5) extract-rgb (gray - 2) extract-rgb (gray + 4) ) v water-level 1 ] ] display end
Finally, with a simple trick we can add much more realistic looking by using shades on the previous colorization. To get this shadow effect, we decide the position of the sun and dark the color of patches that get obscured by the adjacent patch in the direction of the sun.
to shade ; First, compute the higher slope angle in the landscape: ; Difference between heights in adjacent patches let M 0 ask patches [
; We reduce the number of comparissons ask neighbors with [pxcor >= [pxcor] of myself] [ let dif abs (height - [height] of myself) if dif > M [set M dif] ] ] ; Scale will contain the amount of shade to be considered let scale intensity * 100 / M ; Add shadows to patches having higher neighbors in the direction of the sun ask patches with [abs pycor < max-pycor and abs pxcor < max-pxcor] [ let h [height] of dir-sun if h > height [set pcolor map [cut (? - scale * (h - height) )] pcolor] ; In a similar way you can light the patches that are higher than the adjacents ; in the sun direction ;if h < height [set pcolor map [cut (? - scale / 4 * (h - height) )] pcolor] ;set pcolor map [cut (? - scale * (h - height) )] pcolor ] end ; Cut [0, 255] to-report cut [x] ifelse x > 255 [report 255][ifelse x < 0 [report 0][report x]] end ; Set sun direction to-report dir-sun ifelse sun = "up" [report patch-at 0 1] [ ifelse sun = "down" [report patch-at 0 -1] [ ifelse sun = "right" [report patch-at 1 0] [report patch-at -1 0] ]] end
You can download this last model from here. With some changes and improvements the model can be adjusted to work with periodical border conditions working on a torus.
All this process can be made (or imported the results) in NetLogo3D and get some real 3D landscape to be used in other simulations with agents. The only fact to be in mind is that the performance of the render can vary depending on the way we represent the landscape: patches, 3D turtles (as cubes or spheres), or 2D turtles (squares). Also, usually the resolution we can achieve in NetLogo3D is lower than that for 2D.
In future posts we will see how to model the other three elements: water, fire and air. Hope to see you there.