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

// Package planet is the driver for a painted world: template in, terrain out.
//
// It owns the order of operations and nothing else. The image and its legend belong to internal/template,
// the cylinder to internal/world, the cutting up to internal/region, and every physical process to the
// packages that already had it. What lives here is the sequence, the reporting, and the two things that are
// only true of a whole planet: that its open ocean is painted rather than solved, and that its statistics
// pool across regions rather than being computed per region and averaged.
package planet

import (
	"fmt"
	"math"
	"os"
	"path/filepath"
	"strconv"
	"time"

	"salty/terrain/internal/field"
	"salty/terrain/internal/manifest"
	"salty/terrain/internal/overlay"
	"salty/terrain/internal/plates"
	"salty/terrain/internal/region"
	"salty/terrain/internal/template"
	"salty/terrain/internal/uplift"
	"salty/terrain/internal/world"
)

// Inputs is everything a bake needs before a single erosion step has run: the painted map classified and
// projected onto the cylinder, and the cylinder cut into regions.
//
// It is a type of its own because it is worth looking at on its own. Two decisions can wreck an hour-long
// bake - how the legend read the paint, and how the planet was cut up - and both are settled here, in about
// a minute. That is what the plan command exists to show.
type Inputs struct {
	M      *manifest.Manifest
	P      world.Planet
	Legend *template.Legend
	Raster *template.Raster // at paint resolution
	Map    *template.Map    // at planet resolution
	Part   *region.Partition

	// Palette is how the preview is drawn. Never nil: the generator's own when the manifest names none.
	Palette *field.Palette

	PaintW, PaintH int
	Match          template.Match
	EdgeRewritten  int
	Dissolved      int
	MarginCells    int
	Elapsed        time.Duration

	// The annotation layer, or nil throughout when the manifest configures none. It is prepared here rather
	// than at export time for one reason: its coast_jitter marks have to be read *before* the waterline is
	// roughened, which is the first thing that happens to the painting, so by the time a plan exists the
	// overlay has already had its say. Everything else it carries is inert - see internal/overlay.
	Overlay       *overlay.Legend
	OverlayRaster *overlay.Raster
	OverlayMatch  overlay.Match
	OverlayDoc    *overlay.Document

	// Faults is the planet's whole fault set, in world metres, drawn once here because placing a trace needs
	// the class raster of the *whole* cylinder - and because a set drawn per region would put a different
	// fault in every one of them, which is the defect that kept the procedural path's version from being
	// portable at all. Every region reads the same slice and filters it to its own frame.
	Faults []uplift.FaultTrace

	// Plates is the tectonic model, or nil when the manifest asks for none. Drawn here for the same reason
	// the fault set is: a plate is a planet-wide object, and a partition computed per region would give the
	// same physical margin a different classification in every one of them.
	Plates *plates.Model
}

// Painting is the two images already in memory, which is what the studio has: the pictures being edited are
// the ones in the browser, so a plan run against the files on disk would answer a question nobody asked.
//
// A nil Painting, or a nil half of one, means "read what the manifest names". The two halves are separate
// because they are edited separately: moving a road does not re-roughen the coast unless a coast_jitter mark
// moved with it, and the studio's plan cache keys on them one at a time.
type Painting struct {
	Class          []uint8
	ClassW, ClassH int

	Overlay            []uint8
	OverlayAlpha       []uint8
	OverlayW, OverlayH int

	// The tectonic layer, RGB with no alpha: every pixel of it is some plate, so there is no "nothing" to
	// carry. Like the other two it is the authority while the studio is open.
	Plates           []uint8
	PlatesW, PlatesH int
}

// Prepare reads the template, classifies it, projects it onto the planet and partitions it.
//
// Nothing here erodes anything, and nothing here is expensive: the whole thing is a few image passes and
// three distance transforms.
func Prepare(m *manifest.Manifest, log func(string, ...any)) (*Inputs, error) {
	return PrepareWith(m, nil, log)
}

// PrepareWith is Prepare over paintings already in memory. See Painting; nil reads what the manifest names,
// which is what Prepare does.
func PrepareWith(m *manifest.Manifest, art *Painting, log func(string, ...any)) (*Inputs, error) {

	if !m.IsPlanet() {
		return nil, fmt.Errorf("%s has no planet block; this is the square canvas that `generate` builds",
			m.Path)
	}
	start := time.Now()
	if log == nil {
		log = func(string, ...any) {}
	}
	pb := m.Planet
	cell := m.GeologyCellM()

	lg, err := template.Load(m.LegendPath())
	if err != nil {
		return nil, err
	}
	log("legend   %s: %d classes", m.Planet.Legend, len(lg.Classes))

	var px []uint8
	var pw, ph int
	if art != nil {
		px, pw, ph = art.Class, art.ClassW, art.ClassH
	}
	if px == nil {
		var err error
		px, pw, ph, err = template.DecodeRGB(m.TemplatePath())
		if err != nil {
			return nil, err
		}
		log("template %s: %d x %d px", pb.Template, pw, ph)
	} else {
		log("template in memory: %d x %d px", pw, ph)
	}

	ras, match := lg.Classify(px, pw, ph)
	px = nil
	edge, dissolved := ras.DissolveStrokes(lg)
	speckle := ras.Despeckle()
	log("classify %s; %d px rescued at the poles, %d dissolved, %d despeckled",
		match, edge, dissolved, speckle)

	marginCells := pb.MarginCells(cell)
	p, err := world.New(pb.CircumferenceKm*1000, cell, pw, ph, marginCells, pb.NoisePeriodKm*1000)
	if err != nil {
		return nil, err
	}
	log("planet   %d x %d cells of %.1f m: %.1f x %.1f km, %.0f km2 (+%d rows of polar pad)",
		p.W, p.PaintH(), p.CellM, p.CircumferenceM()/1000, p.HeightM()/1000,
		p.CircumferenceM()*p.HeightM()/1e6, p.PadY)

	padClass := lg.FirstSea()
	if pb.PadClass != "" {
		padClass = lg.Index(pb.PadClass)
		if padClass < 0 {
			return nil, fmt.Errorf("%s: planet.pad_class %q is not a class in the legend", m.Path, pb.PadClass)
		}
		if !lg.Classes[padClass].Sea {
			return nil, fmt.Errorf("%s: planet.pad_class %q is land; the pad is the ocean a polar cap "+
				"drains into", m.Path, pb.PadClass)
		}
	}
	if padClass < 0 {
		return nil, fmt.Errorf("%s: the legend has no sea class, so there is nothing to fill the polar pad "+
			"with", m.LegendPath())
	}

	pal := field.DefaultPalette()
	if path := m.PalettePath(); path != "" {
		pal, err = field.LoadPalette(path)
		if err != nil {
			return nil, err
		}
		log("palette  %s", m.Planet.Palette)
	}

	// The annotation layer, read before the coast is roughened rather than after: its coast_jitter marks are
	// the one thing on it the generator reads, and what they decide is how far the waterline may move.
	ov, ovRas, ovMatch, err := loadOverlay(m, art, pw, ph, log)
	if err != nil {
		return nil, err
	}

	// The painted waterline is roughened before it is projected: a drawn shore is a smooth curve and a coast
	// is not. See template/coast.go. Off by default only in the sense that an amplitude of zero is the
	// painting exactly as drawn.
	coast := template.Coast{
		AmplitudePx:  pb.CoastJitterPx,
		WavelengthPx: pb.CoastJitterWavelengthPx,
		Octaves:      pb.CoastJitterOctaves,
		Gain:         pb.CoastJitterGain,
		Seed:         m.Source.Seed,
	}
	if ov != nil && ovRas != nil {
		coast.Scale = ov.CoastScale(ovRas)
	}
	if coast.Amount() {
		ras = ras.RoughenCoast(lg, p, coast)
		masked := ""
		if coast.Scale != nil {
			masked = ", masked by the overlay"
		}
		log("coast    the painted waterline roughened by up to %.0f px over %.0f px bays, %d octaves%s",
			coast.AmplitudePx, coast.WavelengthPx, coast.Octaves, masked)
	}

	pm := ras.Project(p, lg, padClass)

	part, err := region.Build(pm, marginCells, pb.MinLandCells)
	if err != nil {
		return nil, err
	}

	// The plates first, because the fault set is now partly a consequence of them.
	tect, err := buildPlates(m, p, pm, art, log)
	if err != nil {
		return nil, err
	}

	faults := buildFaultSet(pm, lg, pb.FaultGrainKm, m.Source.Seed)
	if len(faults) > 0 {
		log("faults   %d class traces over %d classes, grain %.0f km",
			len(faults), faultClasses(lg), pb.FaultGrainKm)
	}
	if belt := buildBeltFaults(p, pm, tect, pb.Plates.Faults, m.Source.Seed, log); len(belt) > 0 {
		faults = append(faults, belt...)
	}

	in := &Inputs{
		M: m, P: p, Legend: lg, Raster: ras, Map: pm, Part: part, Palette: pal,
		PaintW: pw, PaintH: ph, Match: match,
		EdgeRewritten: edge, Dissolved: dissolved, MarginCells: marginCells,
		Overlay: ov, OverlayRaster: ovRas, OverlayMatch: ovMatch,
		Faults:  faults,
		Plates:  tect,
		Elapsed: time.Since(start),
	}
	if ov != nil && ovRas != nil {
		in.OverlayDoc = ov.Describe(ovRas, ovMatch, in.OverlayScale(),
			m.Planet.Overlay, m.Planet.OverlayLegend)
	}
	return in, nil
}

// faultCandidateTarget is roughly how many strided samples of the planet the fault placement draws from. The
// sample is only ever used for a uniform draw - the *areas* are the projection's own exact counts - so what
// it has to be is dense enough that a small class still has somewhere to put a trace, not dense enough to
// measure anything. Sixty-odd thousand over a 76-million-cell planet is a stride of about 34 cells, 270 m,
// and leaves a class covering a fifth of a per cent with over a hundred candidates.
const faultCandidateTarget = 65536

// buildFaultSet draws the planet's faults, or returns nil when no class asks for any.
func buildFaultSet(pm *template.Map, lg *template.Legend, grainKm float64, seed int64) []uplift.FaultTrace {
	if !lg.HasFaults() || grainKm <= 0 {
		return nil
	}
	specs := make([]uplift.FaultSpec, len(lg.Classes))
	for i := range lg.Classes {
		f := lg.Classes[i].Faults
		if f == nil || !lg.Classes[i].Land() {
			continue
		}
		specs[i] = uplift.FaultSpec{Per1000Km2: f.Per1000Km2, ThrowM: f.ThrowM, LengthKm: f.LengthKm}
	}

	p := pm.P
	stride := int(math.Sqrt(float64(p.W)*float64(p.PaintH())/faultCandidateTarget) + 0.5)
	if stride < 1 {
		stride = 1
	}
	candidates := make([][]int32, len(lg.Classes))
	// Painted rows only. The polar pad is synthetic ocean that no class was ever painted on, and a trace
	// placed there would be a fault in scaffolding.
	for y := p.PadY; y < p.H-p.PadY; y += stride {
		for x := 0; x < p.W; x += stride {
			i := y*p.W + x
			c := pm.Class[i]
			if specs[c].Wanted() {
				candidates[c] = append(candidates[c], int32(i))
			}
		}
	}
	areaCells, _, _ := pm.Counts()
	return uplift.BuildFaults(p, seed, grainKm, specs, candidates, areaCells)
}

// buildPlates draws the planet's tectonics, or returns nil when the manifest asks for none.
//
// The land mask is read through a callback at world coordinates rather than handed over as a raster, and
// that is what keeps internal/plates ignorant of templates. What it wants from the painting is one bit per
// position - continent or ocean - and that bit is what decides whether a plate is continental, and therefore
// whether a margin between two of them is a collision or a subduction zone.
func buildPlates(mf *manifest.Manifest, p world.Planet, pm *template.Map, art *Painting,
	log func(string, ...any)) (*plates.Model, error) {

	cfg := mf.Planet.Plates
	inMemory := art != nil && art.Plates != nil
	painted := inMemory || mf.HasPaintedPlates()
	if !painted && cfg.Count <= 0 {
		return nil, nil
	}

	var m *plates.Model
	var err error
	if painted {
		m, err = paintedPlates(mf, p, pm, art, cfg, log)
	} else {
		m, err = plates.Build(p, mf.Source.Seed, cfg, landAt(p, pm))
	}
	if err != nil {
		return nil, err
	}

	continental := 0
	for i := range m.Plates {
		if m.Plates[i].Continental {
			continental++
		}
	}
	byKind := plates.LengthByKind(m.Boundaries)
	source := "from seed " + strconv.FormatInt(mf.Source.Seed, 10)
	if painted {
		source = "painted"
	}
	log("plates   %d %s (%d continental), %d boundaries, tectonic grid %.0f m",
		len(m.Plates), source, continental, len(m.Boundaries), m.GCellM)
	log("         collision %.0f km, subduction %.0f km, rift %.0f km, ridge %.0f km, transform %.0f km",
		byKind[plates.Collision]/1000, byKind[plates.Subduction]/1000, byKind[plates.Rift]/1000,
		byKind[plates.Ridge]/1000, byKind[plates.Transform]/1000)
	return m, nil
}

// landAt answers "is this world position painted land" against the projected map.
//
// A callback rather than the raster itself, because internal/plates and internal/uplift's belt placement both
// want exactly this one bit and neither should know what a template is. It is also the only thing the
// painting tells the tectonic model, which is worth being able to point at: everything else about a plate
// comes from the seed and the manifest.
func landAt(p world.Planet, pm *template.Map) func(xM, yM float64) bool {
	return func(xM, yM float64) bool {
		x := p.WrapX(int(math.Floor(xM / p.CellM)))
		y := p.ClampY(int(math.Floor(yM/p.CellM)) + p.PadY)
		return !pm.Sea[y*p.W+x]
	}
}

// buildBeltFaults places the traces that belong to the plate margins rather than to a painted class.
//
// Separate from buildFaultSet and added to the same slice, because from the solve's point of view a fault is
// a fault: both end up in Inputs.Faults, both are rasterised by the same FaultDelta, and the difference is
// only in how they were placed. Where they came from survives in the log line and in the trace's own Class,
// which is -1 for a belt fault because no painted colour asked for it.
func buildBeltFaults(p world.Planet, pm *template.Map, tect *plates.Model, cfg plates.Belt,
	seed int64, log func(string, ...any)) []uplift.FaultTrace {

	if tect == nil || !cfg.Wanted() {
		return nil
	}
	out := uplift.BuildBeltFaults(p, seed, cfg, tect.Boundaries, landAt(p, pm))
	cfg = cfg.WithDefaults()
	log("         %d belt traces, %.0f km deformation half-width at %.0f cm/yr, %.0f%% conjugate",
		len(out), cfg.ZoneKm, cfg.ReferenceCmYr, cfg.Conjugate()*100)
	return out
}

// paintedPlates reads the tectonic layer and turns it into a model.
//
// The decoding happens here rather than in internal/plates for the same reason the land mask does: that
// package deals in geometry and motion, and giving it a file path would give it an opinion about image
// formats, paths and the manifest. It is handed pixels.
func paintedPlates(mf *manifest.Manifest, p world.Planet, pm *template.Map, art *Painting,
	cfg plates.Config, log func(string, ...any)) (*plates.Model, error) {

	legendPath := mf.PlatesLegendPath()
	lg, err := plates.LoadPaintLegend(legendPath)
	if err != nil {
		return nil, err
	}

	// The sheet in memory wins when there is one, for the same reason Painting exists at all: the studio is
	// editing a picture in a browser, and a plan run against the file on disk would answer a question nobody
	// asked.
	layerPath := mf.Planet.Plates.Layer
	var px []uint8
	var pw, ph int
	if art != nil && art.Plates != nil {
		px, pw, ph = art.Plates, art.PlatesW, art.PlatesH
		layerPath = "in memory"
	} else {
		path := mf.PlatesLayerPath()
		if px, pw, ph, err = template.DecodeRGB(path); err != nil {
			return nil, fmt.Errorf("tectonic layer %s: %w", path, err)
		}
	}

	m, match, err := plates.FromPainting(p, cfg, lg, px, pw, ph, landAt(p, pm))
	if err != nil {
		return nil, fmt.Errorf("tectonic layer %s: %w", layerPath, err)
	}
	log("tectonic %s, %dx%d px, %d plates", layerPath, pw, ph, len(lg.Plates))
	if match.Far > 0 {
		// Reported rather than fatal, and loudly. Every pixel becomes the nearest plate whatever happens, so
		// a layer whose colours have drifted still produces a model - it just produces the wrong one, with
		// boundaries somewhere nobody put them.
		log("         %d of %d sampled cells are over %.0f from any plate colour (worst %.0f); the layer and "+
			"%s disagree", match.Far, match.Cells, lg.WarnDistance, match.MaxDistance,
			filepath.Base(legendPath))
	}
	return m, nil
}

// RebuildFaults draws the fault set again from the legend as it stands now.
//
// It exists for the studio's plan cache, which reuses a whole prepare when only the legend's *numbers*
// changed - and a class's `faults` block is a number that changes the set without changing a pixel of the
// raster the cache is keyed on. Cheap: a strided scan and a few dozen walks, against the six seconds the
// cache is there to avoid.
func (in *Inputs) RebuildFaults() {
	in.Faults = buildFaultSet(in.Map, in.Legend, in.M.Planet.FaultGrainKm, in.M.Source.Seed)
	// The belt set comes back too. It is not the legend's, but it is in the same slice, and a rebuild that
	// dropped it would silently unfault every margin on the planet the first time a class number changed.
	quiet := func(string, ...any) {}
	in.Faults = append(in.Faults,
		buildBeltFaults(in.P, in.Map, in.Plates, in.M.Planet.Plates.Faults, in.M.Source.Seed, quiet)...)
}

// faultClasses is how many classes asked for traces, for the log line.
func faultClasses(lg *template.Legend) int {
	n := 0
	for i := range lg.Classes {
		if lg.Classes[i].Faults != nil {
			n++
		}
	}
	return n
}

// OverlayScale is how an overlay pixel maps to world metres. The overlay is registered to the template and
// shares its frame, so this is the template's scale and not the geology grid's - a pixel is 12.9 m where a
// cell is 8.
func (in *Inputs) OverlayScale() overlay.Scale {
	w, h := in.PaintW, in.PaintH
	if in.OverlayRaster != nil {
		w, h = in.OverlayRaster.W, in.OverlayRaster.H
	}
	circ := in.P.CircumferenceM()
	return overlay.Scale{
		MetresPerPxX:   circ / float64(w),
		MetresPerPxY:   in.P.HeightM() / float64(h),
		CircumferenceM: circ,
	}
}

// loadOverlay reads and classifies the annotation layer, from memory when the studio has one and from disk
// otherwise. Returns nils all round when the manifest configures none, which is not an error anywhere.
//
// The overlay must be the same size as the template. It is registered to it - a mark means "here", and here
// is a place on the painting - so two different sizes is not something to resample past, it is an author who
// exported one of the two at the wrong scale and would otherwise find their villages drifting.
func loadOverlay(m *manifest.Manifest, art *Painting, paintW, paintH int, log func(string, ...any)) (
	*overlay.Legend, *overlay.Raster, overlay.Match, error) {

	var none overlay.Match
	if !m.HasOverlay() {
		return nil, nil, none, nil
	}
	ov, err := overlay.Load(m.OverlayLegendPath())
	if err != nil {
		return nil, nil, none, err
	}

	var px, alpha []uint8
	var w, h int
	if art != nil && art.Overlay != nil {
		px, alpha, w, h = art.Overlay, art.OverlayAlpha, art.OverlayW, art.OverlayH
	} else {
		path := m.OverlayPath()
		if path == "" {
			if ov.Image == "" {
				return nil, nil, none, fmt.Errorf("%s: planet.overlay_legend is set but neither it nor "+
					"planet.overlay names an image", m.Path)
			}
			path = filepath.Join(filepath.Dir(m.OverlayLegendPath()), ov.Image)
		}
		if _, statErr := os.Stat(path); statErr != nil {
			// A configured overlay whose image is not there yet is the state the studio starts an author in:
			// the legend is written and the sheet is blank. Worth saying, not worth failing on.
			log("overlay  %s: no image yet (%s); nothing is marked",
				m.Planet.OverlayLegend, filepath.Base(path))
			return ov, nil, none, nil
		}
		px, alpha, w, h, err = template.DecodeRGBA(path)
		if err != nil {
			return nil, nil, none, err
		}
	}
	if w != paintW || h != paintH {
		return nil, nil, none, fmt.Errorf("the overlay is %dx%d and the template is %dx%d; they are "+
			"registered to each other, so they have to be the same size", w, h, paintW, paintH)
	}
	ras, match := ov.Classify(px, alpha, w, h)
	log("overlay  %d marks: %s", len(ov.Marks), match)
	return ov, ras, match, nil
}

// SolveCells is how many cells the geology solve will actually visit, summed over regions. It is the number
// the bake time is proportional to, and it is well above the land area because every region carries water.
func (in *Inputs) SolveCells() int {
	n := 0
	for _, r := range in.Part.Regions {
		n += r.Cells()
	}
	return n
}

// LandCells is how many painted land cells the regions own.
func (in *Inputs) LandCells() int {
	n := 0
	for _, r := range in.Part.Regions {
		n += r.LandCells
	}
	return n
}

// The measurement the estimates below are scaled from, and the thing it cannot know.
//
// Docs/Terrain.md's time budget records 256 s for 3.2 M cells over 1000 steps on the development machine's
// 16 cores, which is 8.0e-8 s a cell-step. Measured again on a lowland region of this planet - 14.0 M cells,
// 200 steps, 415 s - it is 1.48e-7, very nearly twice as slow, most likely because the square canvas is two
// thirds land while a region is two thirds water: a cheap ocean cell is not a free one.
//
// What no single constant can capture is that **the cost per cell depends on the uplift rate, and by a lot**.
// Measured on the same bake at 1000 steps with four regions in flight:
//
//	lowland  0.08 mm/yr   14.0 M cells   1014 s     72 s per million cells
//	highland 0.90 mm/yr    4.0 M cells   1394 s    350 s per million cells
//	crater   1.60 mm/yr    1.4 M cells   1829 s   1278 s per million cells
//
// Eighteen-fold, and it is not the stream power. It is the hillslope: DiffuseNonlinear sub-steps to stay
// stable, the count rises with the steepest slope on the grid, and it saturates at max_hillslope_substeps -
// 24 by default. Steep ground pays all 24 every step; a plain pays one.
//
// So the estimate is calibrated on the plains and **badly under-predicts a mountainous template**. It is a
// floor rather than a forecast, the printed line says so, and the practical consequence for an author is
// that raising an uplift rate does not only change the terrain, it changes how long the bake takes.
const secondsPerCellStep = 415.0 / (14.02e6 * 200)

// EstimateSeconds is how long a region's solve should take at the manifest's step count.
func (in *Inputs) EstimateSeconds(cells int) float64 {
	return float64(cells) * float64(in.M.Pipeline.Fluvial.Steps) * secondsPerCellStep
}

// bytesPerCell is what a region costs while it is being solved: the fluvial.Grid's eight int32/float32
// arrays and three masks, plus the height, uplift and erodibility fields the solve reads. It is an estimate
// and it is labelled as one wherever it is printed.
const bytesPerCell = 35 + 12 + 16

// EstimateBytes is roughly how much memory a region's solve holds at once.
func (in *Inputs) EstimateBytes(cells int) int64 { return int64(cells) * bytesPerCell }