UnrealPrototyping/Tools/Terrain/internal/uplift/painted.go

package uplift

import (
	"math"

	"salty/terrain/internal/dt"
	"salty/terrain/internal/field"
	"salty/terrain/internal/manifest"
	"salty/terrain/internal/noise"
	"salty/terrain/internal/world"
)

// The painted source, and the one decision behind all of it: paint the uplift, never the height.
//
// Docs/Terrain-Next.md 6 lists importing a painted heightmap under "do not redo these", and the reason is
// not taste. A stream-power solve handed a painted surface erodes it into something else within a few
// hundred steps, and what it produces instead has no relationship to what was drawn - while the drainage
// network, which is the entire reason this generator replaced the droplet pipeline, is thrown away and
// rebuilt from whatever the painting happened to leave behind.
//
// Painting the uplift rate instead means an author draws intent - a range here, lowlands there, a coast like
// this - and the simulation produces terrain that honours it and has real rivers, real divides and a real
// valley hierarchy, because those came out of the physics rather than out of the brush.
//
// What is left for noise to do is therefore narrow, and it is not decoration:
//
//   - The regional swell. D-49 is arithmetic: with critical_area_m2 at 0 the steady-state slope is
//     U/(K*A^m) down to a single cell, so a uniform uplift rate over a wide area gives a surface with no
//     divides at all. A painted lowland holds one rate over tens of kilometres. Without a long-wavelength
//     modulation the plains come out table-flat, the only gradient across them is the priority-flood's
//     epsilon, and the router draws the flood's traversal order as rivers. That was measured once already.
//   - The initial relief, which only breaks the symmetry. Small on purpose: the solve is what produces
//     relief, and starting it from big ridges means it spends its run tearing them down.
//
// Every noise field here is built on world coordinates through noise.WorldUV, so two regions covering the
// same physical place agree to the bit. That is rule 1 of the tiling plan.

// Pass indices for the painted path's seeded sources. They sit above the procedural path's 1..9 so that the
// two never share a stream and adding one here cannot reshuffle the other.
// Feature counts, in lattice cells per noise period. They are named because the warp amounts are derived
// from them - a warp is only meaningful as a fraction of the wavelength it is bending.
const (
	swellCells = 4  // 25 km at a 100 km period: the regional swell
	ridgeCells = 24 // 4.2 km: the initial relief
	crestCells = 64 // 1.6 km: the crest lines
	plainCells = 48 // 2.1 km: the lowland break-up
)

const (
	srcPaintSwell  = 20
	srcPaintRidges = 21
	srcPaintCrests = 22
	srcPaintPlains = 23
	srcPaintWarp   = 24
	srcPaintMassif = 25
)

// Paint is a region's painted world: which class every cell is, which cells are land, and what the legend
// says those classes mean.
type Paint struct {
	Frame world.Frame

	Class []uint8 // one legend index per frame cell
	Land  []bool  // land the region owns; everything else is water, including other regions' islands

	Rates []float32 // per class, metres a year
	Ks    []float32 // per class, the multiplier on stream-power K

	// PlainM and PlainFloor put a class's range inland: within PlainM metres of the waterline the rate ramps
	// from PlainFloor up to the class rate. Per class; zero PlainM means the class reaches the sea at its
	// full rate, which is what every class did before and still does unless an author asks otherwise.
	PlainM     []float64
	PlainFloor []float32

	// MassifFloor and MassifFraction break a class into plain and upland instead of holding it at one rate.
	// Where the fraction is zero the class is uniform, which is what every class did before this existed and
	// what a legend that asks for nothing still gets. See massif.go: the class rate is then the rate a massif
	// reaches, the floor is the plain between them, and the fraction is how much of the class stands above
	// the midpoint of the two.
	MassifFloor    []float32
	MassifFraction []float64

	// MassifCells is the fabric's wavelength in lattice cells of the noise period, from
	// manifest.Planet.MassifCells. Read only when some class asks for a massif.
	MassifCells int

	// RockCells and RockMult are the planet's lithology: the wavelength of the rock field in lattice cells,
	// and the erodibility multiplier of each rock type. LithMix is per class, how much of it shows through.
	// Zero cells, fewer than two multipliers, or every mix at zero means no rock field is built at all.
	RockCells int
	RockMult  []float64
	LithMix   []float64

	// Faults is the planet's whole fault set, in world metres. A region filters it to the traces that reach
	// into its own frame, which is why it is the planet's and not the region's: a fault crossing a region
	// boundary has to be one fault, and two decompositions of the same planet have to produce the same
	// escarpment. RunYears is how long the solve runs, which turns a total throw into a rate.
	Faults   []FaultTrace
	RunYears float64

	// ClampCeilM is the uplift rate at which a divide reaches the angle of repose, at K x1, in metres a year.
	// It bounds what a *fault* may add and nothing else: an author who paints a class past the ceiling gets
	// what they asked for and a warning from `terrain plan`, but a fault stacking on top of one is an
	// accident nobody chose. Zero switches the bound off.
	ClampCeilM float64

	// Variation is how far the swell modulates the painted rate, as a fraction. See the note above: this
	// is what gives a painted plain its divides, and it is the first thing that will be cut for time.
	Variation float64
}

// FromTemplate builds the geology inputs for one region of a painted planet.
//
// It is a sibling of Build rather than a branch inside it. Build's continent mask, percentile range band,
// normalised swell and percentile lithology split are all global operations over the grid they are given,
// and a region is not a world - two regions taking percentiles of their own extents would disagree about
// the same rock. None of them survives here; the paint replaces all four.
func FromTemplate(p Paint, m *manifest.Manifest) *Result {
	f := p.Frame
	cfg := m.Pipeline
	seed := m.Source.Seed
	w, h := f.W, f.H
	cellM := f.P.CellM

	u, v := noise.WorldUV(w, h, cellM, f.OriginXM(), f.OriginYM(), f.P.NoisePeriodM)

	// A low-frequency warp, which bends everything built on it so that ridges curve and cells are not
	// polygons. Build has one and this did not, which was a porting mistake with a very visible signature:
	// Terrain.md records that cellular crest lines without a strong enough warp "turn ranges into a honeycomb
	// of polygon walls", and that is exactly what the first painted mountains looked like - flat plates with
	// hard edges, at every uplift rate, which is how it was eventually told apart from the repose clamp.
	//
	// The warp amounts need converting rather than copying. Build works in map coordinates where 0..1 spans
	// the map once, so its 0.16 and 0.224 are fractions of a whole map; here 0..1 spans one noise period, and
	// what has to be preserved is the warp measured in the *feature's own wavelength*. Build warps the ridges
	// by 0.8 of their wavelength (0.16 against BaseCells 5) and the crests by 3.1 of theirs (0.224 against
	// BaseCells 14), so those ratios are what carry across.
	wx, wy := paintWarp(u, v, seed)

	// The regional swell: long-wavelength, so a painted lowland has hills and basins of its own rather than
	// one uniform rate across a whole continent. One turn of the planet at BaseCells 4 is a 25 km feature,
	// and four octaves take it down to about 3 km.
	ss := noise.NewSource(seed, srcPaintSwell)
	swu, swv := noise.Warp(u, v, wx, wy, 0.2/swellCells)
	swell := noise.FBMAt(swu, swv, ss, noise.Params{BaseCells: swellCells, Octaves: 4, Gain: 0.5})

	rate := field.New(w, h, cellM)
	k := field.New(w, h, cellM)
	land := field.New(w, h, cellM)
	base := make([]bool, w*h)

	// Distance from every land cell to the nearest water, for the coastal plain. One exact transform over the
	// region, computed only when some class asks for it. The region's frame is flat - it is a rectangle cut
	// out of the cylinder with water all round it - so this does not wrap, and the water it measures to is
	// this region's own coastline: anything else inside the frame is a different landmass, and a different
	// landmass is more than a margin away by construction.
	var shoreM []float32
	if wantsPlain(p.PlainM) {
		d2 := dt.Distance2(invert(p.Land), w, h, false)
		shoreM = make([]float32, len(d2))
		for i, d := range d2 {
			shoreM[i] = float32(math.Sqrt(float64(d)) * cellM)
		}
	}

	// The upland fabric, built only when a class asks for one. It is the one field here that is a cut of a
	// planet-wide measurement rather than a value read straight off a noise, which is why it lives in
	// massif.go with the note on why that measurement cannot be a percentile of the region.
	var rank *field.Field
	if anyMassif(p.MassifFraction) {
		rank = MassifRank(f.P, seed, p.MassifCells, u, v)
	}

	// The rock field, the same way and for the same reason: a quantile of the planet, never of the region.
	var rock *field.Field
	if anyMix(p.LithMix) {
		rock = RockK(f.P, seed, p.RockCells, p.RockMult, u, v)
	}

	// And the faults, which are a rate *difference* across a line rather than a field of their own. Built
	// once for the planet and filtered to this frame; nil when none of them reaches it.
	fault := FaultDelta(f, p.Faults, p.RunYears)

	maxRate := 0.0
	for _, r := range p.Rates {
		if float64(r) > maxRate {
			maxRate = float64(r)
		}
	}
	if maxRate <= 0 {
		maxRate = 1
	}

	// The painted class boundary is deliberately not smoothed. The blend rule in Docs/Terrain-Next.md 3.2
	// exists because a painted map coarser than the grid reads as blocks; here a paint pixel is 12.9 m
	// against an 8 m cell, so there is barely an upsample to soften. And a step in the uplift *rate* is a
	// step in steady-state slope, not in height: the solve grades the transition over a hillslope of its own
	// accord, which is a better answer than a blur, and blurring would have pulled the sea's zero into the
	// coastal cells - the mistake D-52 undid, where the land ending decided how fast it was rising.
	// preFault is the rate before any fault touches it: the class rate after the massif cut and the
	// coastal-plain ramp. The initial relief is scaled by it rather than by the finished rate, which is
	// D-63 and is not a detail - see the amplitude below.
	preFault := make([]float32, len(rate.Data))

	clamped := 0
	for i := range rate.Data {
		c0 := p.Class[i]
		kk := float64(p.Ks[c0])
		if rock != nil && p.LithMix[c0] > 0 {
			// The rock field multiplies the class's own erodibility rather than replacing it: `k_mult` is
			// what the author said this ground is made of, and the province is the variation within it.
			kk *= 1 + p.LithMix[c0]*(float64(rock.Data[i])-1)
		}
		k.Data[i] = float32(kk)
		if !p.Land[i] {
			base[i] = true
			continue
		}
		land.Data[i] = 1
		c := p.Class[i]
		r := float64(p.Rates[c])
		if rank != nil && p.MassifFraction[c] > 0 {
			// The class rate is the rate a massif reaches; the floor is the plain between them.
			r = MassifRate(float64(p.MassifFloor[c]), r, float64(rank.Data[i]), p.MassifFraction[c])
		}
		if shoreM != nil && p.PlainM[c] > 0 {
			// Smoothstep rather than linear, so the plain meets the range without a crease in the slope
			// field - a crease there would be a line of channel heads all starting at the same distance
			// from the sea, which is the sort of thing that reads as a contour rather than as terrain.
			t := float64(shoreM[i]) / p.PlainM[c]
			if t > 1 {
				t = 1
			}
			t = t * t * (3 - 2*t)
			floor := float64(p.PlainFloor[c])
			// Only ever downwards. Before massifs the class rate was uniform and the legend guarantees the
			// coastal floor is below it, so this could not fire; a cell of plain between two massifs now
			// sits below the coastal floor perfectly legitimately, and ramping it *up* towards the shore
			// would put a rim of hills round the edge of every continent.
			if r > floor {
				r = floor + (r-floor)*t
			}
		}
		asked := r
		preFault[i] = float32(asked)
		if fault != nil {
			// A fault is a difference in rate across a line. It adds on one side and subtracts on the other,
			// and the subtraction is what tilts the block rather than merely raising a ridge - so it is
			// allowed to take the rate down, but not below zero: subsidence is not modelled.
			if r += float64(fault[i]); r < 0 {
				r = 0
			}
			// The only ceiling in the whole painted path, and it binds on faults alone. Past
			// U = tan(talus)*K*cell the repose clamp shapes the ground instead of erosion and the surface
			// comes out as polygonal facets; an author may choose that for a class, but a fault stacking on
			// top of ground that was already near it is nobody's choice. So the bound is the ceiling *or*
			// whatever the author's own numbers asked for here, whichever is higher. The count is reported.
			if p.ClampCeilM > 0 {
				lim := p.ClampCeilM * kk
				if lim < asked {
					lim = asked
				}
				if r > lim {
					r = lim
					clamped++
				}
			}
		}
		rate.Data[i] = float32(r * (1 + p.Variation*(2*float64(swell.Data[i])-1)))
	}

	// Initial relief. The spec says 50-150 m times normalised uplift and means it; this only breaks the
	// symmetry so the solve has something to bite on.
	rs := noise.NewSource(seed, srcPaintRidges)
	ru, rv := noise.Warp(u, v, wx, wy, 0.8/ridgeCells)
	ridges := noise.FBMAt(ru, rv, rs, noise.Params{BaseCells: ridgeCells, Octaves: 6, Gain: 0.42, Ridged: true})
	cs := noise.NewSource(seed, srcPaintCrests)
	cu, cv := noise.Warp(u, v, wx, wy, 3.1/crestCells) // the stronger warp the crest lines need
	crests := noise.CellularEdges(cu, cv, cs, int(crestCells), 0.95)
	ps := noise.NewSource(seed, srcPaintPlains)
	pu, pv := noise.Warp(u, v, wx, wy, 0.5/plainCells)
	plains := noise.FBMAt(pu, pv, ps, noise.Params{BaseCells: int(plainCells), Octaves: 4, Gain: 0.45})

	ampLo := cfg.Relief.AmplitudeM.Lo()
	ampHi := cfg.Relief.AmplitudeM.Hi()
	crestW := cfg.Relief.CrestWeight

	height := field.New(w, h, cellM)
	for i := range height.Data {
		if base[i] {
			// Ocean sits at sea level for the whole solve and the coastal pass lays the floor afterwards.
			// Left at a real depth, a coastal cell drains into it and the solver cuts the river down to meet
			// it; the first run with a coast eroded the land to 174 m below sea level for exactly that.
			height.Data[i] = float32(m.SeaLevelM)
			continue
		}
		// The amplitude comes from the rate *before* the faults, and is bounded at one (D-63).
		//
		// It used to come from rate.Data, which is the finished rate with the fault delta in it and no
		// upper bound, and that coupling is half of why Bake_018's flanks came out ribbed. The initial
		// relief exists only to break the symmetry of the background so the solve has something to bite
		// on; how much noise is stamped on a hillside is not a fault's decision. With D-62's six
		// kilometre footwalls the ratio went from about 0.18 on unfaulted foreland - 39 m of relief - to
		// 1.12 on a footwall, which is 166 m, on a landmass whose whole relief is 221 m. A thousand steps
		// cannot erase initial relief the size of the landscape, so the ridged noise stopped being a
		// symmetry-breaker and became the terrain: the ribs measure 250-300 m, which is octave five of a
		// 4.2 km ridged fBm. The bound at one is a guard rather than the fix - with the fault gone the
		// rate cannot exceed the largest class rate - but it is the property worth stating.
		norm := float64(preFault[i]) / maxRate
		if norm > 1 {
			norm = 1
		} else if norm < 0 {
			norm = 0
		}
		amp := ampLo + (ampHi-ampLo)*norm
		shape := (1-crestW)*float64(ridges.Data[i]) + crestW*float64(crests.Data[i])
		height.Data[i] = float32(m.SeaLevelM + 20 + amp*shape + float64(plains.Data[i])*8)
	}

	return &Result{Rate: rate, Height: height, Land: land, K: k, Base: base, FaultClamped: clamped}
}

// anyMix reports whether any class lets the rock field through.
func anyMix(mix []float64) bool {
	for _, v := range mix {
		if v > 0 {
			return true
		}
	}
	return false
}

func wantsPlain(plain []float64) bool {
	for _, v := range plain {
		if v > 0 {
			return true
		}
	}
	return false
}

func invert(b []bool) []bool {
	out := make([]bool, len(b))
	for i, v := range b {
		out[i] = !v
	}
	return out
}