UnrealPrototyping/Tools/Terrain/internal/planet/crater.go

package planet

import (
	"math"

	"salty/terrain/internal/dt"
	"salty/terrain/internal/template"
)

// Impacts, stamped onto the finished terrain.
//
// A crater is not an uplift rate and it cannot be one, which is worth writing down because it is the obvious
// thing to try. A closed basin does not survive the fluvial solve: the priority-flood runs every step and
// *raises* every depression to its spill level, so a crater built out of negative uplift would be filled in
// before the run was a hundred steps old. It is also the wrong model. A crater is an event, not a rate - it
// postdates the landscape it sits in - and a pass running after the solve is what that means.
//
// The shape is derived from the painted blob rather than drawn. Distance inward from the blob's own boundary,
// normalised by its widest point, is a coordinate that runs 0 at the shore to 1 at the centre whatever size
// and shape the author painted, so one set of numbers describes every crater on the map.

// craterProfile is the height at a normalised distance t in from the shore.
//
//	t = 0        the waterline: sea level, so the island keeps the outline that was painted
//	t = RimAt    the crest
//	t = WallAt   the foot of the inner wall
//	t > WallAt   floor
//
// Both segments are smoothstepped, so the crest is a ridge rather than a corner and the floor meets the wall
// without a crease. A corner at either would be ground the detail passes then spend their time sanding off.
func craterProfile(c template.Crater, seaLevelM, t float64) float64 {
	switch {
	case t <= 0:
		return seaLevelM
	case t < c.RimAt:
		return seaLevelM + (c.RimM-seaLevelM)*smoothstep(t/c.RimAt)
	case t < c.WallAt:
		return c.RimM + (c.FloorM-c.RimM)*smoothstep((t-c.RimAt)/(c.WallAt-c.RimAt))
	default:
		return c.FloorM
	}
}

func smoothstep(t float64) float64 {
	if t <= 0 {
		return 0
	}
	if t >= 1 {
		return 1
	}
	return t * t * (3 - 2*t)
}

// CraterStats is what the pass stamped.
//
// The floor is measured over the cells that actually reached it rather than over the whole blob, and that is
// not fussiness: the profile starts at sea level on the shoreline, so the minimum over a blob is always zero
// and reporting it as the floor says nothing at all. What is worth knowing is whether the blob was wide
// enough for the profile to get there - a crater painted smaller than its own rim is a hill.
type CraterStats struct {
	Class      string
	Blobs      int
	Cells      int
	FloorCells int     // cells past wall_at, which are the ones at the floor
	RadiusM    float64 // the widest blob's inradius: what the profile is normalised by
	RimM       float64 // the highest point actually stamped
	FloorM     float64 // the lowest point among the floor cells
}

// stampCraters reshapes every blob of every crater class in the planet raster.
//
// It runs after the regions are composited and before the ocean is laid, so it sees finished land and writes
// only onto land the paint marked as crater.
func stampCraters(in *Inputs, res *Result, log func(string, ...any)) []CraterStats {
	if !in.Legend.HasCraters() {
		return nil
	}
	p := in.P
	n := p.W * p.H
	var out []CraterStats

	for ci := range in.Legend.Classes {
		c := in.Legend.Classes[ci]
		if c.Crater == nil {
			continue
		}
		mask := make([]bool, n)
		outside := make([]bool, n)
		count := 0
		for i := 0; i < n; i++ {
			if in.Map.Class[i] == uint8(ci) && !in.Map.Sea[i] {
				mask[i] = true
				count++
			} else {
				outside[i] = true
			}
		}
		if count == 0 {
			continue
		}

		// Distance inward from the blob's boundary: seed the transform with everything that is *not* this
		// class, and every cell of it then knows how far it is from the nearest edge. Wrapped, because a
		// crater on the seam is one crater.
		d2 := dt.Distance2(outside, p.W, p.H, true)
		dist := make([]float64, n)
		for i := range d2 {
			if mask[i] {
				dist[i] = math.Sqrt(float64(d2[i]))
			}
		}

		// Each blob is normalised by its own widest point, so a big crater and a small one get the same
		// shape rather than the same depth. Components are found with the same wrap-aware flood the region
		// partitioner uses; there are a handful of them and they are tiny.
		comp, maxDist, blobs := craterComponents(mask, dist, p.W, p.H, p.WrapX)

		st := CraterStats{Class: c.Name, Blobs: blobs, Cells: count, FloorM: math.Inf(1),
			RimM: math.Inf(-1)}
		for _, d := range maxDist {
			if d*p.CellM > st.RadiusM {
				st.RadiusM = d * p.CellM
			}
		}
		for i := 0; i < n; i++ {
			if !mask[i] || comp[i] < 0 {
				continue
			}
			d := maxDist[comp[i]]
			if d <= 0 {
				continue
			}
			t := dist[i] / d
			h := craterProfile(*c.Crater, in.M.SeaLevelM, t)
			res.Height.Data[i] = float32(h)
			if h > st.RimM {
				st.RimM = h
			}
			if t >= c.Crater.WallAt {
				st.FloorCells++
				if h < st.FloorM {
					st.FloorM = h
				}
			}
		}
		if st.FloorCells == 0 {
			st.FloorM = 0
		}
		out = append(out, st)
		log("crater  %s: %d blob(s), %d cells, widest %.0f m across; rim reached %.0f m, "+
			"%d cells at the floor (%.0f m)",
			st.Class, st.Blobs, st.Cells, 2*st.RadiusM, st.RimM, st.FloorCells, st.FloorM)
		if st.FloorCells == 0 {
			log("        WARNING no cell reached the floor: every blob is narrower than wall_at asks for, " +
				"so this is a hill rather than a crater. Paint it wider or lower wall_at.")
		}
	}
	return out
}

// craterComponents labels each blob and records its widest point, which is the radius the profile is
// normalised by.
func craterComponents(mask []bool, dist []float64, w, h int, wrap func(int) int) (comp []int32, maxDist []float64, n int) {
	comp = make([]int32, len(mask))
	for i := range comp {
		comp[i] = -1
	}
	var stack []int32
	for start := 0; start < len(mask); start++ {
		if !mask[start] || comp[start] >= 0 {
			continue
		}
		id := int32(len(maxDist))
		maxDist = append(maxDist, 0)
		comp[start] = id
		stack = append(stack[:0], int32(start))
		for len(stack) > 0 {
			c := stack[len(stack)-1]
			stack = stack[:len(stack)-1]
			if dist[c] > maxDist[id] {
				maxDist[id] = dist[c]
			}
			cx, cy := int(c)%w, int(c)/w
			for dy := -1; dy <= 1; dy++ {
				ny := cy + dy
				if ny < 0 || ny >= h {
					continue
				}
				base := ny * w
				for dx := -1; dx <= 1; dx++ {
					if dx == 0 && dy == 0 {
						continue
					}
					ni := int32(base + wrap(cx+dx))
					if mask[ni] && comp[ni] < 0 {
						comp[ni] = id
						stack = append(stack, ni)
					}
				}
			}
		}
	}
	return comp, maxDist, len(maxDist)
}