summaryrefslogtreecommitdiff
path: root/static/src/_posts/2014-01-11-diamond-square.md
blob: 528c95301f2a8938f2a33a95d585a3b13932fd7e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
---
title: Diamond Square
description: >-
    Tackling the problem of semi-realistic looking terrain generation in
    clojure.
updated: 2018-09-06
tags: tech art
---

![terrain][terrain]

I recently started looking into the diamond-square algorithm (you can find a
great article on it [here][diamondsquare]). The following is a short-ish
walkthrough of how I tackled the problem in clojure and the results. You can
find the [leiningen][lein] repo [here][repo] and follow along within that, or
simply read the code below to get an idea.

Also, Marco ported my code into clojurescript, so you can get random terrain
in your browser. [Check it out!][marco]

```clojure
(ns diamond-square.core)

; == The Goal ==
; Create a fractal terrain generator using clojure

; == The Algorithm ==
; Diamond-Square. We start with a grid of points, each with a height of 0.
;
; 1. Take each corner point of the square, average the heights, and assign that
;    to be the height of the midpoint of the square. Apply some random error to
;    the midpoint.
;
; 2. Creating a line from the midpoint to each corner we get four half-diamonds.
;    Average the heights of the points (with some random error) and assign the
;    heights to the midpoints of the diamonds.
;
; 3. We now have four square sections, start at 1 for each of them (with
;    decreasing amount of error for each iteration).
;
; This picture explains it better than I can:
; https://blog.mediocregopher.com/img/diamond-square/dsalg.png
; (http://nbickford.wordpress.com/2012/12/21/creating-fake-landscapes/dsalg/)
;
; == The Strategy ==
; We begin with a vector of vectors of numbers, and iterate over it, filling in
; spots as they become available. Our grid will have the top-left being (0,0),
; y being pointing down and x going to the right. The outermost vector
; indicating row number (y) and the inner vectors indicate the column number (x)
;
; = Utility =
; First we create some utility functions for dealing with vectors of vectors.

(defn print-m
  "Prints a grid in a nice way"
  [m]
  (doseq [n m]
    (println n)))

(defn get-m
  "Gets a value at the given x,y coordinate of the grid, with [0,0] being in the
  top left"
  [m x y]
  ((m y) x))

(defn set-m
  "Sets a value at the given x,y coordinat of the grid, with [0,0] being in the
  top left"
  [m x y v]
  (assoc m y
    (assoc (m y) x v)))

(defn add-m
  "Like set-m, but adds the given value to the current on instead of overwriting
  it"
  [m x y v]
  (set-m m x y
   (+ (get-m m x y) v)))

(defn avg
  "Returns the truncated average of all the given arguments"
  [& l]
  (int (/ (reduce + l) (count l))))

; = Grid size =
; Since we're starting with a blank grid we need to find out what sizes the
; grids can be. For convenience the size (height and width) should be odd, so we
; easily get a midpoint. And on each iteration we'll be halfing the grid, so
; whenever we do that the two resultrant grids should be odd and halfable as
; well, and so on.
;
; The algorithm that fits this is size = 2^n + 1, where 1 <= n. For the rest of
; this guide I'll be referring to n as the "degree" of the grid.


(def exp2-pre-compute
  (vec (map #(int (Math/pow 2 %)) (range 31))))

(defn exp2
  "Returns 2^n as an integer. Uses pre-computed values since we end up doing
  this so much"
  [n]
  (exp2-pre-compute n))

(def grid-sizes
  (vec (map #(inc (exp2 %)) (range 1 31))))

(defn grid-size [degree]
  (inc (exp2 degree)))

; Available grid heights/widths are as follows:
;[3 5 9 17 33 65 129 257 513 1025 2049 4097 8193 16385 32769 65537 131073
;262145 524289 1048577 2097153 4194305 8388609 16777217 33554433 67108865
;134217729 268435457 536870913 1073741825])

(defn blank-grid
  "Generates a grid of the given degree, filled in with zeros"
  [degree]
  (let [gsize (grid-size degree)]
    (vec (repeat gsize
      (vec (repeat gsize 0))))))

(comment
  (print-m (blank-grid 3))
)

; = Coordinate Pattern (The Tricky Part) =
; We now have to figure out which coordinates need to be filled in on each pass.
; A pass is defined as a square step followed by a diamond step. The next pass
; will be the square/dimaond steps on all the smaller squares generated in the
; pass. It works out that the number of passes required to fill in the grid is
; the same as the degree of the grid, where the first pass is 1.
;
; So we can easily find patterns in the coordinates for a given degree/pass,
; I've laid out below all the coordinates for each pass for a 3rd degree grid
; (which is 9x9).

; Degree 3 Pass 1 Square
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . 1 . . . .] (4,4)
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]

; Degree 3 Pass 1 Diamond
; [. . . . 2 . . . .] (4,0)
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]
; [2 . . . . . . . 2] (0,4) (8,4)
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . 2 . . . .] (4,8)

; Degree 3 Pass 2 Square
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . 3 . . . 3 . .] (2,2) (6,2)
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . . . . . . . .]
; [. . 3 . . . 3 . .] (2,6) (6,6)
; [. . . . . . . . .]
; [. . . . . . . . .]

; Degree 3 Pass 2 Diamond
; [. . 4 . . . 4 . .] (2,0) (6,0)
; [. . . . . . . . .]
; [4 . . . 4 . . . 4] (0,2) (4,2) (8,2)
; [. . . . . . . . .]
; [. . 4 . . . 4 . .] (2,4) (6,4)
; [. . . . . . . . .]
; [4 . . . 4 . . . 4] (0,6) (4,6) (8,6)
; [. . . . . . . . .]
; [. . 4 . . . 4 . .] (2,8) (6,8)

; Degree 3 Pass 3 Square
; [. . . . . . . . .]
; [. 5 . 5 . 5 . 5 .] (1,1) (3,1) (5,1) (7,1)
; [. . . . . . . . .]
; [. 5 . 5 . 5 . 5 .] (1,3) (3,3) (5,3) (7,3)
; [. . . . . . . . .]
; [. 5 . 5 . 5 . 5 .] (1,5) (3,5) (5,5) (7,5)
; [. . . . . . . . .]
; [. 5 . 5 . 5 . 5 .] (1,7) (3,7) (5,7) (7,7)
; [. . . . . . . . .]

; Degree 3 Pass 3 Square
; [. 6 . 6 . 6 . 6 .] (1,0) (3,0) (5,0) (7,0)
; [6 . 6 . 6 . 6 . 6] (0,1) (2,1) (4,1) (6,1) (8,1)
; [. 6 . 6 . 6 . 6 .] (1,2) (3,2) (5,2) (7,2)
; [6 . 6 . 6 . 6 . 6] (0,3) (2,3) (4,3) (6,3) (8,3)
; [. 6 . 6 . 6 . 6 .] (1,4) (3,4) (5,4) (7,4)
; [6 . 6 . 6 . 6 . 6] (0,5) (2,5) (4,5) (6,5) (8,5)
; [. 6 . 6 . 6 . 6 .] (1,6) (3,6) (5,6) (7,6)
; [6 . 6 . 6 . 6 . 6] (0,7) (2,7) (4,7) (6,7) (8,7)
; [. 6 . 6 . 6 . 6 .] (1,8) (3,8) (5,8) (7,8)
;
; I make two different functions, one to give the coordinates for the square
; portion of each pass and one for the diamond portion of each pass. To find the
; actual patterns it was useful to first look only at the pattern in the
; y-coordinates, and figure out how that translated into the pattern for the
; x-coordinates.

(defn grid-square-coords
  "Given a grid degree and pass number, returns all the coordinates which need
  to be computed for the square step of that pass"
  [degree pass]
  (let [gsize (grid-size degree)
        start (exp2 (- degree pass))
        interval (* 2 start)
        coords (map #(+ start (* interval %))
                (range (exp2 (dec pass))))]
    (mapcat (fn [y]
      (map #(vector % y) coords))
      coords)))
;
; (grid-square-coords 3 2)
; => ([2 2] [6 2] [2 6] [6 6])

(defn grid-diamond-coords
  "Given a grid degree and a pass number, returns all the coordinates which need
  to be computed for the diamond step of that pass"
  [degree pass]
  (let [gsize (grid-size degree)
        interval (exp2 (- degree pass))
        num-coords (grid-size pass)
        coords (map #(* interval %) (range 0 num-coords))]
    (mapcat (fn [y]
      (if (even? (/ y interval))
        (map #(vector % y) (take-nth 2 (drop 1 coords)))
        (map #(vector % y) (take-nth 2 coords))))
      coords)))

; (grid-diamond-coords 3 2)
; => ([2 0] [6 0] [0 2] [4 2] [8 2] [2 4] [6 4] [0 6] [4 6] [8 6] [2 8] [6 8])

; = Height Generation =
; We now work on functions which, given a coordinate, will return what value
; coordinate will have.

(defn avg-points
  "Given a grid and an arbitrary number of points (of the form [x y]) returns
  the average of all the given points that are on the map. Any points which are
  off the map are ignored"
  [m & coords]
  (let [grid-size (count m)]
    (apply avg
      (map #(apply get-m m %)
        (filter
          (fn [[x y]]
            (and (< -1 x) (> grid-size x)
                 (< -1 y) (> grid-size y)))
          coords)))))

(defn error
  "Returns a number between -e and e, inclusive"
  [e]
  (- (rand-int (inc (* 2 e))) e))

; The next function is a little weird. It primarily takes in a point, then
; figures out the distance from that point to the points we'll take the average
; of. The locf (locator function) is used to return back the actual points to
; use. For the square portion it'll be the points diagonal from the given one,
; for the diamond portion it'll be the points to the top/bottom/left/right from
; the given one.
;
; Once it has those points, it finds the average and applies the error. The
; error function is nothing more than a number between -interval and +interval,
; where interval is the distance between the given point and one of the averaged
; points. It is important that the error decreases the more passes you do, which
; is why the interval is used.
;
; The error function is what should be messed with primarily if you want to
; change what kind of terrain you generate (a giant mountain instead of
; hills/valleys, for example). The one we use is uniform for all intervals, so
; it generates a uniform terrain.

(defn- grid-fill-point
  [locf m degree pass x y]
  (let [interval (exp2 (- degree pass))
        leftx  (- x interval)
        rightx (+ x interval)
        upy    (- y interval)
        downy  (+ y interval)
        v      (apply avg-points m
                (locf x y leftx rightx upy downy))]
    (add-m m x y (+ v (error interval)))))

(def grid-fill-point-square
  "Given a grid, the grid's degree, the current pass number, and a point on the
  grid, fills in that point with the average (plus some error) of the
  appropriate corner points, and returns the resultant grid"
  (partial grid-fill-point
    (fn [_ _ leftx rightx upy downy]
      [[leftx upy]
       [rightx upy]
       [leftx downy]
       [rightx downy]])))

(def grid-fill-point-diamond
  "Given a grid, the grid's degree, the current pass number, and a point on the
  grid, fills in that point with the average (plus some error) of the
  appropriate edge points, and returns the resultant grid"
  (partial grid-fill-point
    (fn [x y leftx rightx upy downy]
      [[leftx y]
       [rightx y]
       [x upy]
       [x downy]])))

; = Filling in the Grid =
; We finally compose the functions we've been creating to fill in the entire
; grid

(defn- grid-fill-point-passes
  "Given a grid, a function to fill in coordinates, and a function to generate
  those coordinates, fills in all coordinates for a given pass, returning the
  resultant grid"
  [m fill-f coord-f degree pass]
  (reduce
    (fn [macc [x y]] (fill-f macc degree pass x y))
    m
    (coord-f degree pass)))

(defn grid-pass
  "Given a grid and a pass number, does the square then the diamond portion of
  the pass"
  [m degree pass]
  (-> m
    (grid-fill-point-passes
      grid-fill-point-square grid-square-coords degree pass)
    (grid-fill-point-passes
      grid-fill-point-diamond grid-diamond-coords degree pass)))

; The most important function in this guide, does all the work
(defn terrain
  "Given a grid degree, generates a uniformly random terrain on a grid of that
  degree"
  ([degree]
    (terrain (blank-grid degree) degree))
  ([m degree]
    (reduce
      #(grid-pass %1 degree %2)
      m
      (range 1 (inc degree)))))

(comment
  (print-m
    (terrain 5))
)

; == The Results ==
; We now have a generated terrain, probably. We should check it. First we'll
; create an ASCII representation. But to do that we'll need some utility
; functions.

(defn max-terrain-height
  "Returns the maximum height found in the given terrain grid"
  [m]
  (reduce max
    (map #(reduce max %) m)))

(defn min-terrain-height
  "Returns the minimum height found in the given terrain grid"
  [m]
  (reduce min
    (map #(reduce min %) m)))

(defn norm
  "Given x in the range (A,B), normalizes it into the range (0,new-height)"
  [A B new-height x]
  (int (/ (* (- x A) new-height) (- B A))))

(defn normalize-terrain
  "Given a terrain map and a number of \"steps\", normalizes the terrain so all
  heights in it are in the range (0,steps)"
  [m steps]
  (let [max-height (max-terrain-height m)
        min-height (min-terrain-height m)
        norm-f (partial norm min-height max-height steps)]
    (vec (map #(vec (map norm-f %)) m))))

; We now define which ASCII characters we want to use for which heights. The
; vector starts with the character for the lowest height and ends with the
; character for the heighest height.

(def tiles
  [\~ \~ \" \" \x \x \X \$ \% \# \@])

(defn tile-terrain
  "Given a terrain map, converts it into an ASCII tile map"
  [m]
  (vec (map #(vec (map tiles %))
    (normalize-terrain m (dec (count tiles))))))

(comment
  (print-m
    (tile-terrain
      (terrain 5)))

; [~ ~ " " x x x X % $ $ $ X X X X X X $ x x x X X X x x x x " " " ~]
; [" ~ " " x x X X $ $ $ X X X X X X X X X X X X X X x x x x " " " "]
; [" " " x x x X X % $ % $ % $ $ X X X X $ $ $ X X X X x x x x " " "]
; [" " " x x X $ % % % % % $ % $ $ X X $ $ $ $ X X x x x x x x " " x]
; [" x x x x X $ $ # % % % % % % $ X $ X X % $ % X X x x x x x x x x]
; [x x x X $ $ $ % % % % % $ % $ $ $ % % $ $ $ $ X X x x x x x x x x]
; [X X X $ % $ % % # % % $ $ % % % % $ % $ $ X $ X $ X X x x x X x x]
; [$ $ X $ $ % $ % % % % $ $ $ % # % % % X X X $ $ $ X X X x x x x x]
; [% X X % % $ % % % $ % $ % % % # @ % $ $ X $ X X $ X x X X x x x x]
; [$ $ % % $ $ % % $ $ X $ $ % % % % $ $ X $ $ X X X X X X x x x x x]
; [% % % X $ $ % $ $ X X $ $ $ $ % % $ $ X X X $ X X X x x X x x X X]
; [$ $ $ X $ $ X $ X X X $ $ $ $ % $ $ $ $ $ X $ X x X X X X X x X X]
; [$ $ $ $ X X $ X X X X X $ % % % % % $ X $ $ $ X x X X X $ X X $ $]
; [X $ $ $ $ $ X X X X X X X % $ % $ $ $ X X X X X x x X X x X X $ $]
; [$ $ X X $ X X x X $ $ X X $ % X X X X X X X X X x X X x x X X X X]
; [$ $ X X X X X X X $ $ $ $ $ X $ X X X X X X X x x x x x x x X X X]
; [% % % $ $ X $ X % X X X % $ $ X X X X X X x x x x x x x x x X X $]
; [$ % % $ $ $ X X $ $ $ $ $ $ X X X X x X x x x x " x x x " x x x x]
; [$ X % $ $ $ $ $ X X X X X $ $ X X X X X X x x " " " " " " " " x x]
; [$ X $ $ % % $ X X X $ X X X x x X X x x x x x " " " " " ~ " " " "]
; [$ $ X X % $ % X X X X X X X X x x X X X x x x " " " " " " ~ " " "]
; [$ $ X $ % $ $ X X X X X X x x x x x x x x x " " " " " " " " " ~ ~]
; [$ $ $ $ $ X X $ X X X X X x x x x x x x x " " " " " " " ~ " " " ~]
; [$ % X X $ $ $ $ X X X X x x x x x x x x x x " " " " ~ " " ~ " " ~]
; [% $ $ X $ X $ X $ X $ X x x x x x x x x x x " " " " ~ ~ ~ " ~ " ~]
; [$ X X X X $ $ $ $ $ X x x x x x x x x x x " " " " ~ ~ ~ ~ ~ ~ ~ ~]
; [X x X X x X X X X X X X X x x x x x x x x x " " " ~ ~ " " ~ ~ ~ ~]
; [x x x x x x X x X X x X X X x x x x x x x " x " " " " " ~ ~ ~ ~ ~]
; [x x x x x x x x X X X X $ X X x X x x x x x x x x " ~ ~ ~ ~ ~ ~ ~]
; [" x x x x x X x X X X X X X X X X x x x x x x " " " " ~ ~ ~ ~ ~ ~]
; [" " " x x x X X X X $ $ $ X X X X X X x x x x x x x x " " ~ ~ ~ ~]
; [" " " " x x x X X X X X $ $ X X x X X x x x x x x x " " " " " ~ ~]
; [~ " " x x x x X $ X $ X $ $ X x X x x x x x x x x x x x x " " " ~]
)

; = Pictures! =
; ASCII is cool, but pictures are better. First we import some java libraries
; that we'll need, then define the colors for each level just like we did tiles
; for the ascii representation.

(import
  'java.awt.image.BufferedImage
  'javax.imageio.ImageIO
  'java.io.File)

(def colors
  [0x1437AD 0x04859D 0x007D1C 0x007D1C 0x24913C
   0x00C12B 0x38E05D 0xA3A3A4 0x757575 0xFFFFFF])

; Finally we reduce over a BufferedImage instance to output every tile as a
; single pixel on it.

(defn img-terrain
  "Given a terrain map and a file name, outputs a png representation of the
  terrain map to that file"
  [m file]
  (let [img (BufferedImage. (count m) (count m) BufferedImage/TYPE_INT_RGB)]
    (reduce
      (fn [rown row]
        (reduce
          (fn [coln tile]
            (.setRGB img coln rown (colors tile))
            (inc coln))
          0 row)
        (inc rown))
      0 (normalize-terrain m (dec (count colors))))
    (ImageIO/write img "png" (File. file))))

(comment
  (img-terrain
    (terrain 10)
    "resources/terrain.png")

  ; https://blog.mediocregopher.com/img/diamond-square/terrain.png
)

; == Conclusion ==
; There's still a lot of work to be done. The algorithm starts taking a
; non-trivial amount of time around the 10th degree, which is only a 1025x1025px
; image. I need to profile the code and find out where the bottlenecks are. It's
; possible re-organizing the code to use pmaps instead of reduces in some places
; could help.
```

[marco]: http://marcopolo.io/diamond-square/
[terrain]: /img/diamond-square/terrain.png
[diamondsquare]: http://www.gameprogrammer.com/fractal.html
[lein]: https://github.com/technomancy/leiningen
[repo]: https://github.com/mediocregopher/diamond-square