420 lines
14 KiB
Go
420 lines
14 KiB
Go
package detail
|
|
|
|
import (
|
|
"math"
|
|
|
|
"salty/terrain/internal/field"
|
|
"salty/terrain/internal/manifest"
|
|
"salty/terrain/internal/world"
|
|
)
|
|
|
|
// brush is the 3x3 kernel a droplet's cut goes through, weights summing to one.
|
|
//
|
|
// A one-cell footprint leaves every droplet path as a rill one cell wide, which reads across the lowlands as
|
|
// brush strokes. Deposits are *not* spread through it and land on the droplet's own bilinear cell instead:
|
|
// spread through the brush, a pit's rim rises faster than its floor, the pit never fills, and every droplet
|
|
// that drains into it adds to the rim until there is a mound.
|
|
var brush = [9]struct {
|
|
dx, dy int
|
|
w float64
|
|
}{
|
|
{0, 0, 0.36},
|
|
{0, 1, 0.12}, {0, -1, 0.12}, {1, 0, 0.12}, {-1, 0, 0.12},
|
|
{1, 1, 0.04}, {1, -1, 0.04}, {-1, 1, 0.04}, {-1, -1, 0.04},
|
|
}
|
|
|
|
// Maps are the derivative fields the droplets leave behind: how much water passed, how much bedrock was
|
|
// scraped, how much sediment was laid. The layer rules read them - scraped bedrock and convex ridges paint as
|
|
// rock, sediment fans and basins as meadow.
|
|
type Maps struct {
|
|
Flow, Wear, Deposit []float32
|
|
}
|
|
|
|
func newMaps(n int) *Maps {
|
|
return &Maps{Flow: make([]float32, n), Wear: make([]float32, n), Deposit: make([]float32, n)}
|
|
}
|
|
|
|
// ParticleParams is one particle pass over one tile.
|
|
type ParticleParams struct {
|
|
Cfg manifest.Particle
|
|
Seed int64
|
|
Frame world.Frame // the tile's cut, at detail resolution: what the hashes are keyed on
|
|
|
|
SeaLevelM float64
|
|
Hardness *Hardness
|
|
|
|
// Classes gives each cell its own droplet density, which is the difference between a rain-fed landscape
|
|
// and an arid one: drop it and the dendritic gully network thins to isolated channels.
|
|
Classes *Classes
|
|
}
|
|
|
|
// ParticleStats is what the pass moved, in metres.
|
|
type ParticleStats struct {
|
|
Droplets int
|
|
Rounds int
|
|
LargestCut, LargestFill float64
|
|
}
|
|
|
|
// RunParticle erodes a tile in place with hydraulic droplets.
|
|
//
|
|
// h is in metres and land marks the cells droplets may spawn on. Everything inside works in *cell heights* -
|
|
// metres over the cell size - so a slope of 1 is 45 degrees and every constant in the manifest means the same
|
|
// thing at any resolution, which is how the numpy was tuned and why the numbers carry across.
|
|
//
|
|
// Determinism, which is the part that is not a port. The numpy draws spawn cells from an RNG stream; that is
|
|
// index-dependent, so a cell would get different droplets depending on which tile it fell in and every seam
|
|
// would show. Here a cell's droplet count and every one of their choices is a hash of (seed, world position),
|
|
// so a droplet spawned in a tile's interior is bit-identical to the one spawned when that cell falls inside a
|
|
// neighbour's margin.
|
|
//
|
|
// The pass runs in rounds, which is the numpy's batching kept deliberately rather than inherited: droplets
|
|
// within a round read the height as it was when the round began and scatter their deltas into per-band
|
|
// buffers summed afterwards in band order, so two droplets in one cell in one round do not see each other and
|
|
// the result does not depend on which goroutine ran. Feedback - a channel deepening as more water follows it -
|
|
// comes from the rounds, not from within one.
|
|
func RunParticle(h *field.Field, land []bool, p ParticleParams) (*Maps, ParticleStats) {
|
|
var st ParticleStats
|
|
cellM := h.CellM
|
|
w, ht := h.W, h.H
|
|
maps := newMaps(w * ht)
|
|
cfg := p.Cfg
|
|
if cfg.Lifetime <= 0 || (cfg.DropletsPerCell <= 0 && p.Classes == nil) {
|
|
return maps, st
|
|
}
|
|
|
|
// Into cell heights, and back at the end.
|
|
hc := make([]float64, w*ht)
|
|
inv := 1 / cellM
|
|
for i, v := range h.Data {
|
|
hc[i] = float64(v) * inv
|
|
}
|
|
// The numpy spawns on land standing at least two metres clear of the water, which keeps droplets out of
|
|
// the surf zone where they would only churn the beach the coastal pass laid.
|
|
spawnAbove := (p.SeaLevelM + 2) / cellM
|
|
|
|
lifetime := cfg.Lifetime
|
|
inertia := cfg.Inertia
|
|
capacityF := cfg.Capacity
|
|
minSlope := cfg.MinSlope
|
|
depositRate := cfg.DepositRate
|
|
erodeRate := cfg.ErodeRate * orOne(cfg.Scale)
|
|
maxChange := cfg.MaxChange * orOne(cfg.Scale)
|
|
evaporation := cfg.Evaporation
|
|
gravity := cfg.Gravity
|
|
maxSpeed := cfg.MaxSpeed
|
|
maxLoad := cfg.MaxLoad
|
|
minErode := math.Max(cfg.MinErodeSlope, 1e-6)
|
|
limit := float64(w) - 2.001
|
|
limitY := float64(ht) - 2.001
|
|
|
|
// How many droplets each cell spawns, and therefore how many rounds. Counting first costs one pass over
|
|
// the tile and makes the round count a property of the world rather than of the loop.
|
|
total := 0
|
|
for y := 0; y < ht; y++ {
|
|
for x := 0; x < w; x++ {
|
|
i := y*w + x
|
|
if !land[i] || hc[i] <= spawnAbove {
|
|
continue
|
|
}
|
|
wx, wy := p.Frame.PlanetXY(x, y)
|
|
total += int(p.Classes.droplets(i, cfg.DropletsPerCell) + hashXY(p.Seed, wx, wy, 0))
|
|
}
|
|
}
|
|
if total == 0 {
|
|
return maps, st
|
|
}
|
|
// Rounds comes from the manifest and *not* from the droplet count, which is the one place this departs
|
|
// from the numpy on purpose. Derived from the count it would depend on how big a piece of the world was
|
|
// being worked on, so a droplet would land in a different round in a tile than in the whole map and the
|
|
// seams would not close.
|
|
rounds := cfg.Rounds
|
|
if rounds < 1 {
|
|
rounds = 1
|
|
}
|
|
st.Droplets, st.Rounds = total, rounds
|
|
|
|
reach := lifetime + 2 // a droplet steps one cell at a time; the brush adds one more
|
|
// A fixed band size, not one per core: a cell's contributions are summed band by band and floating-point
|
|
// addition is not associative, so a partition that moved with GOMAXPROCS would move the last bit with it.
|
|
const bandRows = 64
|
|
bands := field.FixedBandCount(ht, bandRows)
|
|
type buf struct {
|
|
y0, y1 int // the rows this band may touch
|
|
dh []float64
|
|
flow, wear, dep []float32
|
|
}
|
|
bufs := make([]buf, bands)
|
|
|
|
for round := 0; round < rounds; round++ {
|
|
field.FixedBands(ht, bandRows, func(b, y0, y1 int) {
|
|
lo := y0 - reach
|
|
if lo < 0 {
|
|
lo = 0
|
|
}
|
|
hi := y1 + reach
|
|
if hi > ht {
|
|
hi = ht
|
|
}
|
|
n := (hi - lo) * w
|
|
bf := &bufs[b]
|
|
if len(bf.dh) != n {
|
|
bf.dh = make([]float64, n)
|
|
bf.flow = make([]float32, n)
|
|
bf.wear = make([]float32, n)
|
|
bf.dep = make([]float32, n)
|
|
} else {
|
|
clear(bf.dh)
|
|
clear(bf.flow)
|
|
clear(bf.wear)
|
|
clear(bf.dep)
|
|
}
|
|
bf.y0, bf.y1 = lo, hi
|
|
|
|
add := func(x, y int, dh, flow, wear, dep float64) {
|
|
if y < lo || y >= hi || x < 0 || x >= w {
|
|
return
|
|
}
|
|
j := (y-lo)*w + x
|
|
bf.dh[j] += dh
|
|
bf.flow[j] += float32(flow)
|
|
bf.wear[j] += float32(wear)
|
|
bf.dep[j] += float32(dep)
|
|
}
|
|
|
|
for y := y0; y < y1; y++ {
|
|
for x := 0; x < w; x++ {
|
|
i := y*w + x
|
|
if !land[i] || hc[i] <= spawnAbove {
|
|
continue
|
|
}
|
|
wx, wy := p.Frame.PlanetXY(x, y)
|
|
count := int(p.Classes.droplets(i, cfg.DropletsPerCell) + hashXY(p.Seed, wx, wy, 0))
|
|
for j := 0; j < count; j++ {
|
|
if int(hashXY(p.Seed, wx, wy, int32(100+j))*float64(rounds)) != round {
|
|
continue
|
|
}
|
|
px := clampF(float64(x)+hashXY(p.Seed, wx, wy, int32(3*j+1)), 1, limit)
|
|
py := clampF(float64(y)+hashXY(p.Seed, wx, wy, int32(3*j+2)), 1, limitY)
|
|
runDroplet(hc, land, w, px, py, dropletConst{
|
|
lifetime: lifetime, inertia: inertia, capacityF: capacityF,
|
|
minSlope: minSlope, depositRate: depositRate, erodeRate: erodeRate,
|
|
maxChange: maxChange, evaporation: evaporation, gravity: gravity,
|
|
maxSpeed: maxSpeed, maxLoad: maxLoad, minErode: minErode,
|
|
limitX: limit, limitY: limitY,
|
|
}, p.Hardness, add)
|
|
}
|
|
}
|
|
}
|
|
})
|
|
|
|
// Summed in band order, never drained from a channel: the result must not depend on which goroutine
|
|
// finished first (cross-cutting rule 12).
|
|
for b := range bufs {
|
|
bf := &bufs[b]
|
|
if bf.dh == nil {
|
|
continue
|
|
}
|
|
for y := bf.y0; y < bf.y1; y++ {
|
|
src := (y - bf.y0) * w
|
|
dst := y * w
|
|
for x := 0; x < w; x++ {
|
|
hc[dst+x] += bf.dh[src+x]
|
|
maps.Flow[dst+x] += bf.flow[src+x]
|
|
maps.Wear[dst+x] += bf.wear[src+x]
|
|
maps.Deposit[dst+x] += bf.dep[src+x]
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for i := range h.Data {
|
|
after := float32(hc[i] * cellM)
|
|
if d := float64(after - h.Data[i]); d < st.LargestCut {
|
|
st.LargestCut = d
|
|
} else if d > st.LargestFill {
|
|
st.LargestFill = d
|
|
}
|
|
h.Data[i] = after
|
|
}
|
|
// Wear and deposit are in cell heights; report them in metres like everything else.
|
|
for i := range maps.Wear {
|
|
maps.Wear[i] = float32(float64(maps.Wear[i]) * cellM)
|
|
maps.Deposit[i] = float32(float64(maps.Deposit[i]) * cellM)
|
|
}
|
|
st.LargestCut = -st.LargestCut
|
|
return maps, st
|
|
}
|
|
|
|
type dropletConst struct {
|
|
lifetime int
|
|
inertia, capacityF, minSlope float64
|
|
depositRate, erodeRate, maxChange float64
|
|
evaporation, gravity, maxSpeed float64
|
|
maxLoad, minErode float64
|
|
limitX, limitY float64
|
|
}
|
|
|
|
// runDroplet is one droplet's whole life. It reads the height as it was at the start of the round and reports
|
|
// what it moved through add; it never writes to the shared map itself.
|
|
func runDroplet(h []float64, land []bool, w int, px, py float64, c dropletConst, hard *Hardness,
|
|
add func(x, y int, dh, flow, wear, dep float64)) {
|
|
|
|
dx, dy := 0.0, 0.0
|
|
speed, water, sediment := 1.0, 1.0, 0.0
|
|
|
|
for step := 0; step < c.lifetime; step++ {
|
|
hcv, gx, gy, x0, y0, fx, fy := sampleBilinear(h, w, px, py)
|
|
|
|
dx = dx*c.inertia - gx*(1-c.inertia)
|
|
dy = dy*c.inertia - gy*(1-c.inertia)
|
|
length := math.Hypot(dx, dy)
|
|
if length <= 1e-9 {
|
|
return // standing water: it cannot pick a direction, so it stops
|
|
}
|
|
dx /= length
|
|
dy /= length
|
|
|
|
nx, ny := px+dx, py+dy
|
|
inside := nx >= 1 && nx <= c.limitX && ny >= 1 && ny <= c.limitY
|
|
hn, _, _, _, _, _, _ := sampleBilinear(h, w, clampF(nx, 1, c.limitX), clampF(ny, 1, c.limitY))
|
|
dh := 0.0
|
|
if inside {
|
|
dh = hn - hcv
|
|
}
|
|
|
|
slope := math.Max(-dh, c.minSlope)
|
|
capacity := math.Min(slope*speed*water*c.capacityF, c.maxLoad)
|
|
hardness := 0.0
|
|
if hard != nil {
|
|
hardness = hard.At(y0*w+x0, hcv)
|
|
}
|
|
// Flat ground resists cutting. The gate has to sit well above the median lowland slope or the
|
|
// meadows come out brushed with rills, which is the lesson 0.25 encodes.
|
|
holds := math.Hypot(gx, gy) / c.minErode
|
|
if holds > 1 {
|
|
holds = 1
|
|
}
|
|
holds *= holds
|
|
|
|
deposit, erode := 0.0, 0.0
|
|
if dh > 0 {
|
|
deposit = math.Min(dh, sediment) // uphill: fill the pit it is climbing out of
|
|
} else if sediment > capacity {
|
|
deposit = (sediment - capacity) * c.depositRate
|
|
}
|
|
if dh <= 0 && sediment <= capacity {
|
|
erode = math.Min((capacity-sediment)*c.erodeRate, -dh) * (1 - hardness) * holds
|
|
}
|
|
|
|
// The sea is a sink: the droplet drops its whole load at the mouth, which is what makes a fan. It is
|
|
// the land mask that decides, not a height comparison - the sea floor is held at sea level while the
|
|
// detail passes run (the same invariant the solve keeps), so there is no depth to compare against.
|
|
intoSea := false
|
|
if inside {
|
|
nxi, nyi := int(nx+0.5), int(ny+0.5)
|
|
if nxi >= 0 && nxi < w && nyi >= 0 && nyi*w+nxi < len(land) {
|
|
intoSea = !land[nyi*w+nxi]
|
|
}
|
|
}
|
|
if intoSea {
|
|
deposit, erode = sediment, 0
|
|
} else {
|
|
deposit = math.Min(deposit, c.maxChange)
|
|
erode = math.Min(erode, c.maxChange)
|
|
}
|
|
|
|
// Neither the cut nor the deposit may touch water. Both stencils straddle the waterline whenever a
|
|
// droplet is within a cell of it, and the sea floor is held at sea level here and put back afterwards,
|
|
// so anything written there would be silently thrown away - sediment that should have built a beach,
|
|
// quietly deleted. The cut is simply skipped, because cutting a sea floor that is a placeholder means
|
|
// nothing; the deposit is given to the droplet's own cell, which is land for as long as it is alive.
|
|
onLand := func(x, y int) bool {
|
|
if x < 0 || x >= w || y < 0 {
|
|
return false
|
|
}
|
|
i := y*w + x
|
|
return i < len(land) && land[i]
|
|
}
|
|
if erode > 0 {
|
|
for _, b := range brush {
|
|
if onLand(x0+b.dx, y0+b.dy) {
|
|
add(x0+b.dx, y0+b.dy, -erode*b.w, 0, 0, 0)
|
|
}
|
|
}
|
|
}
|
|
if deposit > 0 {
|
|
put := func(x, y int, amount float64) {
|
|
if !onLand(x, y) {
|
|
x, y = x0, y0
|
|
}
|
|
add(x, y, amount, 0, 0, 0)
|
|
}
|
|
put(x0, y0, deposit*(1-fx)*(1-fy))
|
|
put(x0+1, y0, deposit*fx*(1-fy))
|
|
put(x0, y0+1, deposit*(1-fx)*fy)
|
|
put(x0+1, y0+1, deposit*fx*fy)
|
|
}
|
|
add(x0, y0, 0, water, erode, deposit)
|
|
|
|
sediment += erode - deposit
|
|
speed = math.Min(math.Sqrt(math.Max(0, speed*speed-dh*c.gravity)), c.maxSpeed)
|
|
water *= 1 - c.evaporation
|
|
|
|
if !inside || intoSea || water <= 0.001 {
|
|
return
|
|
}
|
|
px, py = nx, ny
|
|
}
|
|
}
|
|
|
|
// sampleBilinear is the height and its gradient at a float position, with the integer cell and the
|
|
// fractions the caller needs to scatter back. The caller keeps the position inside [1, size-2].
|
|
func sampleBilinear(h []float64, w int, px, py float64) (hc, gx, gy float64, x0, y0 int, fx, fy float64) {
|
|
x0 = int(px)
|
|
y0 = int(py)
|
|
fx = px - float64(x0)
|
|
fy = py - float64(y0)
|
|
i := y0*w + x0
|
|
h00 := h[i]
|
|
h10 := h[i+1]
|
|
h01 := h[i+w]
|
|
h11 := h[i+w+1]
|
|
gx = (h10-h00)*(1-fy) + (h11-h01)*fy
|
|
gy = (h01-h00)*(1-fx) + (h11-h10)*fx
|
|
hc = h00*(1-fx)*(1-fy) + h10*fx*(1-fy) + h01*(1-fx)*fy + h11*fx*fy
|
|
return hc, gx, gy, x0, y0, fx, fy
|
|
}
|
|
|
|
func clampF(v, lo, hi float64) float64 {
|
|
if v < lo {
|
|
return lo
|
|
}
|
|
if v > hi {
|
|
return hi
|
|
}
|
|
return v
|
|
}
|
|
|
|
func orOne(v float64) float64 {
|
|
if v <= 0 {
|
|
return 1
|
|
}
|
|
return v
|
|
}
|
|
|
|
// hashXY is splitmix64's finaliser over the seed and a world position, in [0, 1). The same arithmetic as the
|
|
// router's jitter and for the same reason: everything random has to be a hash of where a thing is, never of
|
|
// the order it was visited in.
|
|
func hashXY(seed int64, x, y int, k int32) float64 {
|
|
h := uint64(seed)*0x9e3779b97f4a7c15 + 0x243f6a8885a308d3
|
|
h ^= uint64(uint32(int32(x)))*0x9e3779b97f4a7c15 +
|
|
uint64(uint32(int32(y)))*0xc2b2ae3d27d4eb4f +
|
|
uint64(uint32(k))*0x165667b19e3779f9
|
|
h ^= h >> 30
|
|
h *= 0xbf58476d1ce4e5b9
|
|
h ^= h >> 27
|
|
h *= 0x94d049bb133111eb
|
|
h ^= h >> 31
|
|
return float64(h>>11) / float64(uint64(1)<<53)
|
|
}
|