tower / files / src / terrain / mountains.ts
mountains.ts
| 1 | /** |
| 2 | * Mountains: summits, the ridges between them, and the valleys water cuts into them. |
| 3 | * See world-generation.md §4.1. |
| 4 | * |
| 5 | * v1 made a mountain as a dome with ridged noise laid over it. That gives one summit |
| 6 | * somewhere near the middle of a blister and nothing a person would call a ridge, |
| 7 | * a saddle or a valley. This builds a range the way the landform is actually read: |
| 8 | * |
| 9 | * 1. **Summits.** The authored peaks, plus generated ones scattered through the range, |
| 10 | * each lower than the named summit it stands near. |
| 11 | * 2. **A ridge graph.** The summits are triangulated, and the minimum spanning tree of |
| 12 | * the triangulation — with a few of the remaining short edges as spurs — becomes the |
| 13 | * ridge lines. Each ridge dips to a saddle between its two summits. |
| 14 | * 3. **An uplift field.** Height falls away from every ridge along a concave flank, |
| 15 | * blended between neighbouring ridges rather than cut off where they meet. |
| 16 | * 4. **Erosion to steady state.** Treating that field as the rate the ground is being |
| 17 | * raised, the stream-power law has a steady state that can be computed directly |
| 18 | * over a drainage network — each cell stands above the cell it drains to by just |
| 19 | * enough slope to carry away what rises — and that state is what a range worn by |
| 20 | * water for a long time looks like: branching valleys, sharp divides. It is then |
| 21 | * scaled back to the massif's authored shape, so the Girdle stays 2,900 m high. |
| 22 | * 5. **Pinning.** Every authored summit is put back at exactly its elevation. |
| 23 | * |
| 24 | * Ranges that overlap are one group, eroded together, so the Girdle and the four named |
| 25 | * massifs standing on it make one drainage rather than five that meet at a seam. |
| 26 | * |
| 27 | * The result is a grid of relief, plain data in the compiled world. Tiles sample it and |
| 28 | * add fine detail on top (field.ts), so no tile ever needs to know about another. |
| 29 | */ |
| 30 | import Delaunator from 'delaunator' |
| 31 | import { CIRCUMFERENCE_M, RADIUS_M } from '../canon/constants.js' |
| 32 | import type { Mountains } from '../world/schema.js' |
| 33 | import { type Grid, Heap, blur, drainage, forCells, sampleGrid } from './grid.js' |
| 34 | import { hash3, mulberry32, simplex } from './noise.js' |
| 35 | import { type Box, type Shape, nearestA, shapeBox, shapeOf, signedDistance } from './shapes.js' |
| 36 | |
| 37 | export { type Grid, sampleGrid } from './grid.js' |
| 38 | |
| 39 | const C = CIRCUMFERENCE_M |
| 40 | |
| 41 | export interface CPeak { |
| 42 | z: number |
| 43 | a: number |
| 44 | elevationM: number |
| 45 | /** True for a summit the author placed; generated summits are free to move. */ |
| 46 | authored: boolean |
| 47 | } |
| 48 | |
| 49 | export interface CMountainGroup { |
| 50 | /** The features eroded together. The first is the one the status dialog names. */ |
| 51 | ids: string[] |
| 52 | name: string |
| 53 | /** True when the group rings the world, so its grid wraps around theta. */ |
| 54 | ring: boolean |
| 55 | /** Grid origin, spacing and size. Row `i` is z, column `j` is a. */ |
| 56 | z0: number |
| 57 | a0: number |
| 58 | cell: number |
| 59 | nz: number |
| 60 | na: number |
| 61 | /** Relief above the land under the range, metres, one per cell. */ |
| 62 | relief: Float32Array |
| 63 | /** Every summit, authored and generated. */ |
| 64 | peaks: CPeak[] |
| 65 | /** |
| 66 | * Just the authored ones, with `elevationM` as height above the land under them. Fine |
| 67 | * detail is kept off them so they stay at their height. |
| 68 | */ |
| 69 | summits: CPeak[] |
| 70 | ruggedness: number |
| 71 | snowlineM: number |
| 72 | seed: number |
| 73 | box: Box |
| 74 | } |
| 75 | |
| 76 | /** Grid spacing. Finer than the 264 m v1 raster; tiles add the detail below it. */ |
| 77 | export const MOUNTAIN_CELL_M = 200 |
| 78 | |
| 79 | interface Member { |
| 80 | f: Mountains |
| 81 | shape: Shape |
| 82 | box: Box |
| 83 | maxInside: number |
| 84 | } |
| 85 | |
| 86 | // -------------------------------------------------------------------- grouping |
| 87 | |
| 88 | /** Overlapping ranges, as groups of members. A ring overlaps anything it spans in z. */ |
| 89 | function groups(members: Member[]): Member[][] { |
| 90 | const parent = members.map((_, i) => i) |
| 91 | const find = (i: number): number => (parent[i] === i ? i : (parent[i] = find(parent[i]!))) |
| 92 | const overlaps = (p: Box, q: Box): boolean => { |
| 93 | if (p.z1 < q.z0 || q.z1 < p.z0) return false |
| 94 | if (!Number.isFinite(p.a0) || !Number.isFinite(q.a0)) return true |
| 95 | const shift = nearestA((q.a0 + q.a1) / 2, (p.a0 + p.a1) / 2) - (q.a0 + q.a1) / 2 |
| 96 | return !(p.a1 < q.a0 + shift || q.a1 + shift < p.a0) |
| 97 | } |
| 98 | for (let i = 0; i < members.length; i++) { |
| 99 | for (let j = i + 1; j < members.length; j++) { |
| 100 | if (overlaps(members[i]!.box, members[j]!.box)) parent[find(i)] = find(j) |
| 101 | } |
| 102 | } |
| 103 | const byRoot = new Map<number, Member[]>() |
| 104 | members.forEach((m, i) => { |
| 105 | const r = find(i) |
| 106 | byRoot.set(r, [...(byRoot.get(r) ?? []), m]) |
| 107 | }) |
| 108 | return [...byRoot.values()] |
| 109 | } |
| 110 | |
| 111 | /** The greatest distance from any interior point to the edge, sampled. */ |
| 112 | function maxInside(shape: Shape): number { |
| 113 | if (shape.type === 'band') return (shape.zTo - shape.zFrom) / 2 |
| 114 | const b = shape.box |
| 115 | const N = 48 |
| 116 | let best = 0 |
| 117 | for (let i = 0; i < N; i++) { |
| 118 | for (let j = 0; j < N; j++) { |
| 119 | best = Math.max(best, signedDistance(shape, |
| 120 | b.z0 + ((b.z1 - b.z0) * (i + 0.5)) / N, b.a0 + ((b.a1 - b.a0) * (j + 0.5)) / N)) |
| 121 | } |
| 122 | } |
| 123 | return best |
| 124 | } |
| 125 | |
| 126 | /** |
| 127 | * How much a member claims a point, 0..1: full well inside, fading out past its edge. |
| 128 | * A ring's claim is its distance from the band's centre line, since a ring range is a |
| 129 | * ridge along it rather than a plateau filling it. |
| 130 | */ |
| 131 | function claim(m: Member, z: number, a: number): number { |
| 132 | if (m.shape.type === 'band') { |
| 133 | const half = (m.shape.zTo - m.shape.zFrom) / 2 |
| 134 | const d = Math.abs(z - (m.shape.zFrom + m.shape.zTo) / 2) |
| 135 | return smoothstep(half, half * 0.35, d) |
| 136 | } |
| 137 | return smoothstep(-2000, Math.min(6000, m.maxInside * 0.4), signedDistance(m.shape, z, a)) |
| 138 | } |
| 139 | |
| 140 | const smoothstep = (e0: number, e1: number, x: number): number => { |
| 141 | const t = Math.min(1, Math.max(0, (x - e0) / (e1 - e0))) |
| 142 | return t * t * (3 - 2 * t) |
| 143 | } |
| 144 | |
| 145 | // ---------------------------------------------------------------------- summits |
| 146 | |
| 147 | /** |
| 148 | * Generated summits for one range, on a jittered grid so they are spread without |
| 149 | * clumping. A ring's run along its centre line, a little off it either side; a massif's |
| 150 | * fill its outline and stand taller toward its middle. |
| 151 | */ |
| 152 | function generatedPeaks(m: Member, rand: () => number): CPeak[] { |
| 153 | const out: CPeak[] = [] |
| 154 | const density = m.f.generatedPeaks / 100e6 // per m² |
| 155 | if (density <= 0) return out |
| 156 | const spacing = 1 / Math.sqrt(density) |
| 157 | if (m.shape.type === 'band') { |
| 158 | const zc = (m.shape.zFrom + m.shape.zTo) / 2 |
| 159 | const half = (m.shape.zTo - m.shape.zFrom) / 2 |
| 160 | // Along a ring the summits sit roughly a spacing apart around the circumference. |
| 161 | const count = Math.max(3, Math.round(C / spacing)) |
| 162 | for (let k = 0; k < count; k++) { |
| 163 | const a = ((k + 0.2 + 0.6 * rand()) / count) * C |
| 164 | const z = zc + (rand() - 0.5) * half * 0.5 |
| 165 | out.push({ z, a, elevationM: m.f.heightM * (0.55 + 0.45 * rand()), authored: false }) |
| 166 | } |
| 167 | return out |
| 168 | } |
| 169 | const b = m.shape.box |
| 170 | for (let z = b.z0; z < b.z1; z += spacing) { |
| 171 | for (let a = b.a0; a < b.a1; a += spacing) { |
| 172 | const pz = z + rand() * spacing |
| 173 | const pa = a + rand() * spacing |
| 174 | const sd = signedDistance(m.shape, pz, pa) |
| 175 | if (sd < m.maxInside * 0.2) continue |
| 176 | const depth = Math.min(1, sd / m.maxInside) |
| 177 | out.push({ z: pz, a: pa, elevationM: m.f.heightM * (0.35 + 0.45 * depth) * (0.75 + 0.25 * rand()), authored: false }) |
| 178 | } |
| 179 | } |
| 180 | return out |
| 181 | } |
| 182 | |
| 183 | // ----------------------------------------------------------------- ridge graph |
| 184 | |
| 185 | interface Ridge { p: number; q: number; saddle: number } |
| 186 | |
| 187 | /** |
| 188 | * Ridges between summits: the minimum spanning tree of their Delaunay triangulation, |
| 189 | * which connects every summit by its shortest links, plus the shortest of the leftover |
| 190 | * edges as spurs, so a range has branches and the odd enclosed cirque rather than one |
| 191 | * long line. Rings are triangulated with a copy of the seam's summits a circumference |
| 192 | * away, so ridges cross theta = 0 like anywhere else. |
| 193 | */ |
| 194 | function ridges(peaks: CPeak[], ring: boolean, rand: () => number, ruggedness: number): Ridge[] { |
| 195 | const n = peaks.length |
| 196 | if (n < 2) return [] |
| 197 | const coords: number[] = [] |
| 198 | const owner: number[] = [] |
| 199 | peaks.forEach((p, i) => { coords.push(p.z, p.a); owner.push(i) }) |
| 200 | if (ring) { |
| 201 | peaks.forEach((p, i) => { |
| 202 | if (p.a < C * 0.25) { coords.push(p.z, p.a + C); owner.push(i) } |
| 203 | if (p.a > C * 0.75) { coords.push(p.z, p.a - C); owner.push(i) } |
| 204 | }) |
| 205 | } |
| 206 | const d = new Delaunator(Float64Array.from(coords)) |
| 207 | const edges = new Map<string, { p: number; q: number; len: number }>() |
| 208 | for (let e = 0; e < d.triangles.length; e++) { |
| 209 | const i = d.triangles[e]! |
| 210 | const j = d.triangles[e % 3 === 2 ? e - 2 : e + 1]! |
| 211 | const p = owner[i]!, q = owner[j]! |
| 212 | if (p === q) continue |
| 213 | const key = p < q ? `${p}:${q}` : `${q}:${p}` |
| 214 | const len = Math.hypot(coords[2 * i]! - coords[2 * j]!, coords[2 * i + 1]! - coords[2 * j + 1]!) |
| 215 | const had = edges.get(key) |
| 216 | if (!had || len < had.len) edges.set(key, { p: Math.min(p, q), q: Math.max(p, q), len }) |
| 217 | } |
| 218 | // Kruskal. Links between two high summits are cheaper, so ridges follow the high |
| 219 | // ground rather than wandering down to a low summit and back up. |
| 220 | const sorted = [...edges.values()].sort((x, y) => cost(x) - cost(y)) |
| 221 | function cost(e: { p: number; q: number; len: number }): number { |
| 222 | return e.len / (0.5 + Math.min(peaks[e.p]!.elevationM, peaks[e.q]!.elevationM) / 4000) |
| 223 | } |
| 224 | const parent = peaks.map((_, i) => i) |
| 225 | const find = (i: number): number => (parent[i] === i ? i : (parent[i] = find(parent[i]!))) |
| 226 | const out: Ridge[] = [] |
| 227 | const tree = new Set<string>() |
| 228 | for (const e of sorted) { |
| 229 | const a = find(e.p), b = find(e.q) |
| 230 | if (a === b) continue |
| 231 | parent[a] = b |
| 232 | tree.add(`${e.p}:${e.q}`) |
| 233 | out.push(ridge(e)) |
| 234 | } |
| 235 | const treeLens = out.length ? sorted.filter((e) => tree.has(`${e.p}:${e.q}`)).map((e) => e.len) : [0] |
| 236 | const median = [...treeLens].sort((x, y) => x - y)[Math.floor(treeLens.length / 2)]! |
| 237 | for (const e of sorted) { |
| 238 | if (tree.has(`${e.p}:${e.q}`)) continue |
| 239 | if (e.len < median * 1.25 && rand() < 0.25 + 0.35 * ruggedness) out.push(ridge(e)) |
| 240 | } |
| 241 | return out |
| 242 | |
| 243 | function ridge(e: { p: number; q: number; len: number }): Ridge { |
| 244 | const low = Math.min(peaks[e.p]!.elevationM, peaks[e.q]!.elevationM) |
| 245 | // A saddle is deeper on a long ridge; never so deep that the ridge disappears. |
| 246 | const drop = Math.min(low * 0.45, e.len * (0.015 + 0.02 * rand())) |
| 247 | return { p: e.p, q: e.q, saddle: low - drop } |
| 248 | } |
| 249 | } |
| 250 | |
| 251 | // ---------------------------------------------------------------- uplift field |
| 252 | |
| 253 | /** |
| 254 | * Height a ridge or summit gives a cell: the crest height at the nearest point along it, |
| 255 | * falling away across a concave flank — steep under the crest, easing out at the foot. |
| 256 | * The flank's reach scales with the crest's height, so high ridges have broad roots. |
| 257 | */ |
| 258 | function flank(crest: number, d: number): number { |
| 259 | const reach = crest / FLANK_SLOPE |
| 260 | if (d >= reach) return 0 |
| 261 | return crest * (1 - d / reach) ** 1.2 |
| 262 | } |
| 263 | /** |
| 264 | * Crest height over flank reach: the mean slope from crest to foot, about 9°. A 4,300 m |
| 265 | * summit has roots 27 km out and still stands 2,000 m high 12 km from its top, which is |
| 266 | * the size of v1's massifs and where their tarns were. The flanks are only gently |
| 267 | * concave; the steepness a range actually has comes from erosion cutting into them. |
| 268 | */ |
| 269 | const FLANK_SLOPE = 0.16 |
| 270 | |
| 271 | function upliftField( |
| 272 | grid: Grid, peaks: CPeak[], edges: Ridge[], members: Member[], rand: () => number, seed: number, |
| 273 | ): Float32Array { |
| 274 | /** |
| 275 | * The ridges are combined by a soft maximum — the fourth root of the sum of fourth |
| 276 | * powers — rather than a hard one. A hard maximum of cone-like flanks meets its |
| 277 | * neighbour along a straight crease, and a range drawn that way reads as a folded |
| 278 | * paper model. The distance to each ridge is also nudged by noise, so flanks wander. |
| 279 | */ |
| 280 | const acc = new Float32Array(grid.nz * grid.na) |
| 281 | const wobble = simplex(hash3(seed, 17, 0)) |
| 282 | const splat = (za: number, aa: number, zb: number, ab: number, h: (t: number) => number, crestMax: number): void => { |
| 283 | const reach = crestMax / FLANK_SLOPE |
| 284 | const zLo = Math.min(za, zb) - reach, zHi = Math.max(za, zb) + reach |
| 285 | const aLo = Math.min(aa, ab) - reach, aHi = Math.max(aa, ab) + reach |
| 286 | const i0 = Math.max(0, Math.floor((zLo - grid.z0) / grid.cell)) |
| 287 | const i1 = Math.min(grid.nz - 1, Math.ceil((zHi - grid.z0) / grid.cell)) |
| 288 | const j0 = Math.floor((aLo - grid.a0) / grid.cell) |
| 289 | const j1 = Math.ceil((aHi - grid.a0) / grid.cell) |
| 290 | const ez = zb - za, ea = ab - aa |
| 291 | const len2 = ez * ez + ea * ea |
| 292 | for (let i = i0; i <= i1; i++) { |
| 293 | const z = grid.z0 + i * grid.cell |
| 294 | for (let jj = j0; jj <= j1; jj++) { |
| 295 | let j = jj |
| 296 | if (grid.ring) j = ((jj % grid.na) + grid.na) % grid.na |
| 297 | else if (jj < 0 || jj >= grid.na) continue |
| 298 | const a = grid.a0 + jj * grid.cell |
| 299 | const t = len2 > 0 ? Math.max(0, Math.min(1, ((z - za) * ez + (a - aa) * ea) / len2)) : 0 |
| 300 | const d = Math.hypot(z - (za + ez * t), a - (aa + ea * t)) |
| 301 | const jitter = 1 + 0.3 * wobble(z / 9000, a / 9000, 0) |
| 302 | const v = flank(h(t), d * jitter) |
| 303 | if (v <= 0) continue |
| 304 | const k = i * grid.na + j |
| 305 | const v2 = v * v |
| 306 | acc[k]! += v2 * v2 |
| 307 | } |
| 308 | } |
| 309 | } |
| 310 | |
| 311 | for (const r of edges) { |
| 312 | const p = peaks[r.p]!, q = peaks[r.q]! |
| 313 | // Unwrap across the seam so a ridge takes the short way round. |
| 314 | const qa = grid.ring ? nearestA(q.a, p.a) : q.a |
| 315 | // A little asymmetry in where the saddle falls, so ridges are not all symmetric. |
| 316 | const mid = 0.35 + 0.3 * rand() |
| 317 | const h = (t: number): number => { |
| 318 | const base = p.elevationM + (q.elevationM - p.elevationM) * t |
| 319 | const dip = t < mid ? Math.sin((t / mid) * (Math.PI / 2)) : Math.sin(((1 - t) / (1 - mid)) * (Math.PI / 2)) |
| 320 | const saddleBase = p.elevationM + (q.elevationM - p.elevationM) * mid |
| 321 | return base - (saddleBase - r.saddle) * dip * dip |
| 322 | } |
| 323 | splat(p.z, p.a, q.z, qa, h, Math.max(p.elevationM, q.elevationM)) |
| 324 | } |
| 325 | // Every summit stands on its own, too, so a lone peak and a ridge's end both have a top. |
| 326 | for (const p of peaks) splat(p.z, p.a, p.z, p.a, () => p.elevationM, p.elevationM) |
| 327 | |
| 328 | const out = new Float32Array(grid.nz * grid.na) |
| 329 | for (let k = 0; k < out.length; k++) out[k] = Math.sqrt(Math.sqrt(acc[k]!)) |
| 330 | |
| 331 | // Fade out past the edge of every member's claim. Each member only visits its own box. |
| 332 | const claimed = new Float32Array(grid.nz * grid.na) |
| 333 | for (const m of members) { |
| 334 | const i0 = Math.max(0, Math.floor((m.box.z0 - grid.z0) / grid.cell)) |
| 335 | const i1 = Math.min(grid.nz - 1, Math.ceil((m.box.z1 - grid.z0) / grid.cell)) |
| 336 | let j0 = 0, j1 = grid.na - 1 |
| 337 | if (Number.isFinite(m.box.a0) && !grid.ring) { |
| 338 | const shift = nearestA((m.box.a0 + m.box.a1) / 2, grid.a0 + (grid.na * grid.cell) / 2) - (m.box.a0 + m.box.a1) / 2 |
| 339 | j0 = Math.max(0, Math.floor((m.box.a0 + shift - grid.a0) / grid.cell)) |
| 340 | j1 = Math.min(grid.na - 1, Math.ceil((m.box.a1 + shift - grid.a0) / grid.cell)) |
| 341 | } |
| 342 | const ring = grid.ring && Number.isFinite(m.box.a0) |
| 343 | const jStart = ring ? Math.floor((m.box.a0 - grid.a0) / grid.cell) : j0 |
| 344 | const jEnd = ring ? Math.ceil((m.box.a1 - grid.a0) / grid.cell) : j1 |
| 345 | for (let i = i0; i <= i1; i++) { |
| 346 | const z = grid.z0 + i * grid.cell |
| 347 | for (let jj = jStart; jj <= jEnd; jj++) { |
| 348 | const j = ring ? ((jj % grid.na) + grid.na) % grid.na : jj |
| 349 | const k = i * grid.na + j |
| 350 | const c = claim(m, z, grid.a0 + jj * grid.cell) |
| 351 | if (c > claimed[k]!) claimed[k] = c |
| 352 | } |
| 353 | } |
| 354 | } |
| 355 | for (let k = 0; k < out.length; k++) out[k]! *= claimed[k]! |
| 356 | return out |
| 357 | } |
| 358 | |
| 359 | // --------------------------------------------------------------------- erosion |
| 360 | |
| 361 | /** Drainage area past which ground is a channel rather than hillslope: 1 km². */ |
| 362 | const CHANNEL_AREA_M2 = 1e6 |
| 363 | /** Hillslope gradient at full uplift. The rescaling to the authored shape sets heights. */ |
| 364 | const STEADY_SLOPE = 0.6 |
| 365 | |
| 366 | /** |
| 367 | * Erode to the stream-power steady state, and scale it back to the authored shape. |
| 368 | * |
| 369 | * Under E = K·A^m·S with n = 1, ground raised at rate U settles where erosion keeps |
| 370 | * pace: slope S = U / (K·A^m). Over a drainage network that gives every cell's height |
| 371 | * directly — its receiver's height plus the run between them times that slope — in one |
| 372 | * pass up from the outlets. Divides drain almost nothing, so their slope is steep and |
| 373 | * they stand high; valleys drain a great deal and run nearly level. That is the branching |
| 374 | * the uplift field alone never has. The network depends on the heights, so the pass is |
| 375 | * repeated a few times with the network re-found by priority flood each time |
| 376 | * (Cordonnier et al., 2016). |
| 377 | * |
| 378 | * The steady state's heights have their own scale. So it is scaled, cell by cell, by the |
| 379 | * ratio of the authored field to the steady state, both blurred over a few kilometres: |
| 380 | * the broad shape comes back — the Girdle 2,900 m high, Snowcrown's massif 4,300 m — and |
| 381 | * the ridges and valleys within it are the water's. `strength` blends the result with the |
| 382 | * uncut field. |
| 383 | */ |
| 384 | function erode(grid: Grid, field: Float32Array, outlet: Uint8Array, strength: number, seed: number): Float32Array { |
| 385 | const n = grid.nz * grid.na |
| 386 | if (strength <= 0) return field |
| 387 | const work = { |
| 388 | receiver: new Int32Array(n), stack: new Int32Array(n), reached: new Uint8Array(n), |
| 389 | key: new Float64Array(n), heap: new Heap(n), |
| 390 | } |
| 391 | const area = new Float64Array(n) |
| 392 | /** |
| 393 | * Start from the field roughened at valley scales. Water crossing a smooth plane takes |
| 394 | * one direction everywhere and cuts parallel gullies along the grid's diagonals; a |
| 395 | * few kilometres of unevenness is what lets streams find each other and join. |
| 396 | */ |
| 397 | const rough = simplex(hash3(seed, 29, 0)) |
| 398 | const h = new Float32Array(n) |
| 399 | for (let i = 0; i < grid.nz; i++) { |
| 400 | for (let j = 0; j < grid.na; j++) { |
| 401 | const k = i * grid.na + j |
| 402 | const z = i * grid.cell, a = j * grid.cell |
| 403 | const wob = rough(z / 6000, a / 6000, 1.7) * 0.6 + rough(z / 2200, a / 2200, 5.3) * 0.4 |
| 404 | h[k] = field[k]! * (1 + 0.18 * wob) |
| 405 | } |
| 406 | } |
| 407 | const jitter = new Float32Array(n) |
| 408 | const cellArea = grid.cell * grid.cell |
| 409 | // Uplift rate, scaled to 0..1 by the highest field value. The absolute rate only sets |
| 410 | // the steady state's overall scale, which the rescaling below replaces anyway. |
| 411 | let peak = 0 |
| 412 | for (let k = 0; k < n; k++) peak = Math.max(peak, field[k]!) |
| 413 | if (peak <= 0) return field |
| 414 | |
| 415 | for (let pass = 0; pass < 8; pass++) { |
| 416 | // A different small jitter each pass, so ties on flat ground break a new way and the |
| 417 | // network never settles onto the grid's own axes. |
| 418 | const jr = mulberry32(hash3(seed, 31, pass)) |
| 419 | for (let k = 0; k < n; k++) jitter[k] = h[k]! + (jr() - 0.5) * grid.cell * 0.15 |
| 420 | const top = drainage(grid, jitter, outlet, work) |
| 421 | const { receiver, stack } = work |
| 422 | for (let s = 0; s < top; s++) area[stack[s]!] = cellArea |
| 423 | for (let s = top - 1; s >= 0; s--) { |
| 424 | const c = stack[s]!, r = receiver[c]! |
| 425 | if (r >= 0) area[r]! += area[c]! |
| 426 | } |
| 427 | for (let s = 0; s < top; s++) { |
| 428 | const c = stack[s]!, r = receiver[c]! |
| 429 | if (r < 0) { h[c] = 0; continue } |
| 430 | const ci = Math.floor(c / grid.na), ri = Math.floor(r / grid.na) |
| 431 | let dj = (r % grid.na) - (c % grid.na) |
| 432 | if (grid.ring && Math.abs(dj) > 1) dj = -Math.sign(dj) |
| 433 | const run = grid.cell * Math.hypot(ri - ci, dj) |
| 434 | const u = field[c]! / peak |
| 435 | /** |
| 436 | * Channels only begin past a critical drainage area. Below it the ground is |
| 437 | * hillslope, and runs at the slope the threshold gives rather than steepening to |
| 438 | * every single-cell catchment — without it each cell is its own gully and the |
| 439 | * whole range is ribbed in parallel grooves one cell wide. |
| 440 | */ |
| 441 | const a = Math.max(area[c]!, CHANNEL_AREA_M2) |
| 442 | h[c] = h[r]! + run * u * STEADY_SLOPE * Math.sqrt(CHANNEL_AREA_M2 / a) |
| 443 | } |
| 444 | // Hillslope diffusion: soil creep rounds what water leaves sharp. |
| 445 | const soft = blur(grid, h, 1) |
| 446 | for (let k = 0; k < n; k++) if (!outlet[k]) h[k] = h[k]! * 0.5 + soft[k]! * 0.5 |
| 447 | } |
| 448 | |
| 449 | // Over 8 km: the massifs' shape is the author's; the ridges and valleys within it are the water's. |
| 450 | const blurField = blur(grid, field, 40) |
| 451 | const blurSteady = blur(grid, h, 40) |
| 452 | const out = new Float32Array(n) |
| 453 | for (let k = 0; k < n; k++) { |
| 454 | const scaled = blurSteady[k]! > 1e-3 ? h[k]! * (blurField[k]! / blurSteady[k]!) : 0 |
| 455 | out[k] = field[k]! + (scaled - field[k]!) * strength |
| 456 | } |
| 457 | return out |
| 458 | } |
| 459 | |
| 460 | /** Separable box blur, `radius` cells, three passes each way — close to a Gaussian. */ |
| 461 | // ------------------------------------------------------------------ the compile |
| 462 | |
| 463 | /** A place a range must keep a shelf of ground: a cirque for an authored tarn. */ |
| 464 | export interface Bench { z: number; a: number; radiusM: number } |
| 465 | |
| 466 | /** |
| 467 | * `landAt` is the height of the ground a range stands on, before any range: authored |
| 468 | * summit elevations are above the datum, like every other height in the world, so each |
| 469 | * summit is raised to its elevation less the land under it. |
| 470 | */ |
| 471 | export function compileMountains( |
| 472 | features: Mountains[], worldSeed: number, seedOf: (f: Mountains) => number, |
| 473 | onStep: (label: string) => void, benches: Bench[] = [], landAt: (z: number, a: number) => number = () => 0, |
| 474 | ): CMountainGroup[] { |
| 475 | const members: Member[] = [] |
| 476 | for (const f of features) { |
| 477 | onStep(`Raising ${f.name}`) |
| 478 | const shape = shapeOf(f.geometry) |
| 479 | if (!shape) continue |
| 480 | const inside = maxInside(shape) |
| 481 | const reach = shape.type === 'band' ? 0 : Math.max(f.heightM, ...f.peaks.map((p) => p.elevationM)) / FLANK_SLOPE * 0.3 |
| 482 | const box = shapeBox(shape) |
| 483 | members.push({ f, shape, maxInside: inside, box: { z0: box.z0 - reach, z1: box.z1 + reach, a0: box.a0 - reach, a1: box.a1 + reach } }) |
| 484 | } |
| 485 | |
| 486 | return groups(members).map((ms) => { |
| 487 | // The largest range names the group, so the Girdle and not a summit on it is eroded. |
| 488 | ms.sort((x, y) => area(y) - area(x)) |
| 489 | const lead = ms[0]! |
| 490 | const seed = hash3(worldSeed, seedOf(lead.f), 0x4d7) |
| 491 | const rand = mulberry32(seed) |
| 492 | const ring = ms.some((m) => m.shape.type === 'band') |
| 493 | const cell = MOUNTAIN_CELL_M |
| 494 | |
| 495 | const z0 = Math.max(0, Math.min(...ms.map((m) => m.box.z0)) - cell * 2) |
| 496 | const z1 = Math.max(...ms.map((m) => m.box.z1)) + cell * 2 |
| 497 | let a0 = 0, a1 = C |
| 498 | if (!ring) { |
| 499 | // Unwrap every box to the lead's side of the seam before taking the extent. |
| 500 | const ref = (lead.box.a0 + lead.box.a1) / 2 |
| 501 | const spans = ms.map((m) => { |
| 502 | const shift = nearestA((m.box.a0 + m.box.a1) / 2, ref) - (m.box.a0 + m.box.a1) / 2 |
| 503 | return [m.box.a0 + shift, m.box.a1 + shift] as const |
| 504 | }) |
| 505 | a0 = Math.min(...spans.map((s) => s[0])) - cell * 2 |
| 506 | a1 = Math.max(...spans.map((s) => s[1])) + cell * 2 |
| 507 | } |
| 508 | const nz = Math.ceil((z1 - z0) / cell) + 1 |
| 509 | const na = ring ? Math.round(C / cell) : Math.ceil((a1 - a0) / cell) + 1 |
| 510 | // A ring's cells are stretched a hair so a whole number of them closes the circle. |
| 511 | const grid: Grid = { z0, a0, cell: ring ? C / na : cell, nz, na, ring } |
| 512 | |
| 513 | // Summits: authored first, then generated ones kept below the named summits near them. |
| 514 | // Summits carry their height above the land under them from here on: that is what |
| 515 | // the relief grid holds. |
| 516 | const authored: CPeak[] = ms.flatMap((m) => m.f.peaks.map((p) => ({ |
| 517 | z: p.z, a: ring ? ((p.theta * RADIUS_M) % C + C) % C : nearestA(p.theta * RADIUS_M, (a0 + a1) / 2), |
| 518 | elevationM: p.elevationM - landAt(p.z, p.theta * RADIUS_M), authored: true, |
| 519 | }))) |
| 520 | const generated = ms.flatMap((m) => generatedPeaks(m, rand)).map((p) => ({ |
| 521 | ...p, a: ring ? p.a : nearestA(p.a, (a0 + a1) / 2), |
| 522 | })).filter((p) => { |
| 523 | for (const q of authored) { |
| 524 | const d = Math.hypot(p.z - q.z, nearestA(p.a, q.a) - q.a) |
| 525 | if (d < 1500) return false |
| 526 | const radius = q.elevationM / FLANK_SLOPE |
| 527 | if (d < radius) p.elevationM = Math.min(p.elevationM, q.elevationM * (0.5 + 0.4 * (d / radius))) |
| 528 | } |
| 529 | return true |
| 530 | }) |
| 531 | const peaks = [...authored, ...generated] |
| 532 | const rugged = Math.max(...ms.map((m) => m.f.ruggedness)) |
| 533 | const edges = ridges(peaks, ring, rand, rugged) |
| 534 | |
| 535 | const field = upliftField(grid, peaks, edges, ms, rand, seed) |
| 536 | |
| 537 | onStep(`Letting rain find its way down ${lead.f.name}`) |
| 538 | // Outlets: the grid's z edges, a bounded grid's a edges, and anywhere the range has |
| 539 | // faded to nothing — water leaves the range there. |
| 540 | const outlet = new Uint8Array(nz * na) |
| 541 | for (let i = 0; i < nz; i++) { |
| 542 | for (let j = 0; j < na; j++) { |
| 543 | const k = i * na + j |
| 544 | if (i === 0 || i === nz - 1 || (!ring && (j === 0 || j === na - 1)) || field[k]! < 1) outlet[k] = 1 |
| 545 | } |
| 546 | } |
| 547 | const strength = Math.max(...ms.map((m) => m.f.erosion)) |
| 548 | const relief = erode(grid, field, outlet, strength, seed) |
| 549 | keepBenches(grid, field, relief, benches) |
| 550 | crown(grid, relief, authored) |
| 551 | pin(grid, relief, authored) |
| 552 | summitCone(grid, relief, authored, landAt) |
| 553 | pin(grid, relief, authored) |
| 554 | |
| 555 | return { |
| 556 | ids: ms.map((m) => m.f.id), |
| 557 | name: lead.f.name, |
| 558 | ring, |
| 559 | z0: grid.z0, a0: grid.a0, cell: grid.cell, nz, na, |
| 560 | relief, |
| 561 | peaks, |
| 562 | summits: authored, |
| 563 | ruggedness: rugged, |
| 564 | snowlineM: Math.min(...ms.map((m) => m.f.snowlineM)), |
| 565 | seed, |
| 566 | box: { z0: grid.z0, z1: grid.z0 + (nz - 1) * grid.cell, a0: ring ? -Infinity : grid.a0, a1: ring ? Infinity : grid.a0 + (na - 1) * grid.cell }, |
| 567 | } |
| 568 | }) |
| 569 | |
| 570 | function area(m: Member): number { |
| 571 | if (m.shape.type === 'band') return (m.shape.zTo - m.shape.zFrom) * C |
| 572 | return (m.shape.box.z1 - m.shape.box.z0) * (m.shape.box.a1 - m.shape.box.a0) |
| 573 | } |
| 574 | } |
| 575 | |
| 576 | /** |
| 577 | * Keep a shelf where an authored tarn sits. A cirque is ground the ice hollowed but the |
| 578 | * rivers did not cut down, so it stands at the height of the uncut range around it — |
| 579 | * which is where v1's tarns and the ski stations beside them were, around 2,100 m |
| 580 | * (design.md §1.11). Without this the steady state cuts a valley through every one. |
| 581 | * The lake itself is carved into the shelf afterwards (field.ts). |
| 582 | */ |
| 583 | function keepBenches(grid: Grid, field: Float32Array, relief: Float32Array, benches: Bench[]): void { |
| 584 | for (const b of benches) { |
| 585 | const shelf = sampleGrid(grid, field, b.z, b.a) |
| 586 | if (shelf <= 0) continue |
| 587 | const inner = b.radiusM + 600, outer = inner + 2500 |
| 588 | forCells(grid, b.z, b.a, outer, (k, d) => { |
| 589 | const w = d <= inner ? 1 : 1 - smoothstep(inner, outer, d) |
| 590 | relief[k] = Math.max(relief[k]!, relief[k]! + (shelf - relief[k]!) * w) |
| 591 | }) |
| 592 | } |
| 593 | } |
| 594 | |
| 595 | /** |
| 596 | * Nothing close to an authored summit may stand above it. The steady state has its own |
| 597 | * idea of where the high ground is, and a named summit is named because it is the top. |
| 598 | * So where anything within a few kilometres is higher, the ground around the summit is |
| 599 | * scaled down — by the full ratio out to 2.5 km, easing to none at 7 km — which keeps its |
| 600 | * ridges and valleys, where capping it flat would leave a disc with a rim. |
| 601 | */ |
| 602 | function crown(grid: Grid, relief: Float32Array, peaks: CPeak[]): void { |
| 603 | const INNER = 2500, OUTER = 7000 |
| 604 | for (const p of peaks) { |
| 605 | let highest = 0 |
| 606 | forCells(grid, p.z, p.a, INNER, (k) => { highest = Math.max(highest, relief[k]!) }) |
| 607 | if (highest <= p.elevationM) continue |
| 608 | const ratio = p.elevationM / highest |
| 609 | forCells(grid, p.z, p.a, OUTER, (k, d) => { |
| 610 | const w = d <= INNER ? 1 : 1 - smoothstep(INNER, OUTER, d) |
| 611 | relief[k]! *= 1 - w * (1 - ratio) |
| 612 | }) |
| 613 | } |
| 614 | } |
| 615 | |
| 616 | /** |
| 617 | * The last word on a summit being the top. Pinning raises a summit and, a little less, |
| 618 | * the cells right around it, which can leave a neighbour a few metres above the summit |
| 619 | * it was lowered to meet. So within 2.5 km nothing stands above a cone falling 80 m per |
| 620 | * kilometre, 15 m lower beyond the summit's own cell. After `crown` only a handful of cells are above it, so this trims them |
| 621 | * without leaving a flat disc; the summit is pinned once more after. |
| 622 | */ |
| 623 | function summitCone(grid: Grid, relief: Float32Array, peaks: CPeak[], landAt: (z: number, a: number) => number): void { |
| 624 | for (const p of peaks) { |
| 625 | // In total height, since the land under the range is not flat: the summit's land, |
| 626 | // plus its relief, against each cell's land plus its relief. |
| 627 | const summitLand = landAt(p.z, p.a) |
| 628 | // Past the summit's own cell, 15 m more: room for the final pin's small correction, |
| 629 | // which spreads to the neighbours too. |
| 630 | forCells(grid, p.z, p.a, 2500, (k, d, z, a) => { |
| 631 | const allowed = summitLand + p.elevationM - 0.08 * d - (d > grid.cell ? 15 : 0) - landAt(z, a) |
| 632 | relief[k] = Math.min(relief[k]!, allowed) |
| 633 | }) |
| 634 | } |
| 635 | } |
| 636 | |
| 637 | /** Visit every cell within `radius` of (z, a), with its distance. */ |
| 638 | /** The pinning kernel's radius: how far around a summit the correction is spread. */ |
| 639 | export const PIN_RADIUS_M = 1500 |
| 640 | |
| 641 | /** |
| 642 | * Put every authored summit back at exactly its elevation, as a smooth bump spread over |
| 643 | * the cells around it. Repeated, because the grid is sampled between its nodes: one |
| 644 | * pass lands within a few percent of the correction, and three are exact to a millimetre. |
| 645 | */ |
| 646 | function pin(grid: Grid, relief: Float32Array, peaks: CPeak[]): void { |
| 647 | for (let pass = 0; pass < 4; pass++) { |
| 648 | for (const p of peaks) { |
| 649 | const delta = p.elevationM - sampleGrid(grid, relief, p.z, p.a) |
| 650 | if (Math.abs(delta) < 1e-3) continue |
| 651 | forCells(grid, p.z, p.a, PIN_RADIUS_M, (k, d) => { relief[k]! += delta * pinKernel(d) }) |
| 652 | } |
| 653 | } |
| 654 | } |
| 655 | |
| 656 | /** 1 at the summit, easing to 0 at PIN_RADIUS_M. Shared with field.ts's detail mask. */ |
| 657 | export function pinKernel(d: number): number { |
| 658 | if (d >= PIN_RADIUS_M) return 0 |
| 659 | const u = 1 - d / PIN_RADIUS_M |
| 660 | return u * u * (3 - 2 * u) |
| 661 | } |