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

package planet

import (
	"encoding/json"
	"fmt"
	"image/png"
	"math"
	"os"
	"path/filepath"
	"sync"
	"time"

	"salty/terrain/internal/detail"
	"salty/terrain/internal/field"
	"salty/terrain/internal/manifest"
	"salty/terrain/internal/noise"
	"salty/terrain/internal/overlay"
	"salty/terrain/internal/thermal"
	"salty/terrain/internal/tile"
)

// The detail bake: the geology grid becomes ground somebody can stand on, one tile at a time.
//
// It reads the heightmap a geology bake left behind rather than solving anything itself, which is what makes
// it batchable. The geology is hours; a tile is seconds, and the islands somebody cares about can be baked
// first and the rest later or never.
//
// Every pass here is local, and every hash and noise lattice in them is keyed on absolute world position, so
// a tile's interior comes out the same as it would have in one impossible whole-world run. That is measured
// rather than asserted: internal/detail's TestHowFarTheCutEdgeReachesIn is where the margin comes from.

// TileOptions steer a detail bake.
type TileOptions struct {
	In      *Inputs
	HeightM *field.Field // the geology heightmap, painted rows only, in metres
	Sea     []bool       // painted rows

	// Exposure is the coastal pass's fetch field, painted rows, 0 sheltered to 1 open water, or nil when the
	// bake predates it. A tile cannot compute this - see detail.CoastalParams - so without it the coastal
	// detail pass treats every shore as fully exposed and the run says so once.
	Exposure *field.Field

	Out     string
	Prefix  string
	Only    [4]int // x0,y0,x1,y1 in tile indices; zero means all
	OnlySet bool

	// NoDetail writes the tile as the geology upsampled and nothing else. It is a diagnostic and it earns
	// its place: when ground looks wrong at two metres, the first question is always whether the detail
	// passes did it or whether they are faithfully magnifying something the solve produced, and there is no
	// other way to ask.
	NoDetail bool

	// NoShore skips pass 11b and nothing else, for the same reason NoDetail exists one level up: when a
	// coastline looks wrong the first question is whether the shore pass did it or whether it is faithfully
	// magnifying what the geology handed it, and a diff between two runs is the only way to ask.
	NoShore bool

	Jobs int
	Log  func(string, ...any)
}

// TileRecord is one tile in the index.
type TileRecord struct {
	IX           int     `json:"ix"`
	IY           int     `json:"iy"`
	File         string  `json:"file"`
	W            int     `json:"w"`
	H            int     `json:"h"`
	OriginXM     float64 `json:"origin_x_m"`
	OriginYM     float64 `json:"origin_y_m"`
	MinM         float64 `json:"min_m"`
	MaxM         float64 `json:"max_m"`
	ClipFrac     float64 `json:"clip_fraction"`
	Droplets     int     `json:"droplets"`
	Rounds       int     `json:"rounds"`
	LargestCutM  float64 `json:"largest_cut_m"`
	LargestFillM float64 `json:"largest_fill_m"`
	Seconds      float64 `json:"seconds"`

	// Coastal is what pass 11b moved on this tile, or nil on a tile with no shore in it.
	Coastal *detail.CoastalStats `json:"coastal,omitempty"`

	// OverlayFile is the annotation mask beside this tile - one mark index a detail cell, zero for nothing -
	// or empty when the planet has no overlay. The key is in overlay.json in the same directory.
	OverlayFile string `json:"overlay_file,omitempty"`
}

// TileIndex is tiles.json: everything a consumer needs to place the tiles back into a world.
type TileIndex struct {
	When       time.Time `json:"when"`
	Prefix     string    `json:"prefix"`
	CellM      float64   `json:"cell_m"`
	TilePx     int       `json:"tile_px"`
	MarginPx   int       `json:"margin_px"`
	NX         int       `json:"nx"`
	NY         int       `json:"ny"`
	WrapX      bool      `json:"wrap_x"`
	WorldWM    float64   `json:"world_w_m"`
	WorldHM    float64   `json:"world_h_m"`
	ElevationM struct {
		Min float64 `json:"min"`
		Max float64 `json:"max"`
	} `json:"elevation_m"`
	Tiles []TileRecord `json:"tiles"`
}

// CheckBake refuses a detail bake whose geology was produced by a different manifest.
//
// The failure it exists for is silent and total. A heightmap is 16-bit samples over an elevation range, so a
// bake made under one range and decoded under another comes out shifted - and if the shift takes the land
// below sea level, every tile decides it is ocean, holds itself at sea level, and writes a flat zero. That
// happened on the first run of this command and there was nothing in the output to say why.
func CheckBake(dir string, m *manifest.Manifest, warn func(string, ...any)) error {
	raw, err := os.ReadFile(filepath.Join(dir, "meta.json"))
	if err != nil {
		return fmt.Errorf("%s: %w (run `terrain bake` first)", filepath.Join(dir, "meta.json"), err)
	}
	// Pointers, so that a field a bake did not record is distinguishable from one it recorded as zero. An
	// absent field is a bake older than this check, which is a reason to say so and carry on; a different
	// field is a reason to stop. The first version conflated the two and refused a perfectly good bake.
	var meta struct {
		Seed *int64 `json:"seed"`
		Plan struct {
			CircumferenceKm *float64 `json:"circumference_km"`
			CellM           *float64 `json:"cell_m"`
			ElevationMinM   *float64 `json:"elevation_min_m"`
			ElevationMaxM   *float64 `json:"elevation_max_m"`
		} `json:"plan"`
	}
	if err := json.Unmarshal(raw, &meta); err != nil {
		return fmt.Errorf("%s: %w", filepath.Join(dir, "meta.json"), err)
	}
	if warn == nil {
		warn = func(string, ...any) {}
	}

	unknown := 0
	for _, c := range []struct {
		what string
		was  *float64
		now  float64
	}{
		{"elevation_m.min", meta.Plan.ElevationMinM, m.ElevationM.Min},
		{"elevation_m.max", meta.Plan.ElevationMaxM, m.ElevationM.Max},
		{"circumference_km", meta.Plan.CircumferenceKm, m.Planet.CircumferenceKm},
		{"the geology cell", meta.Plan.CellM, m.GeologyCellM()},
	} {
		if c.was == nil {
			unknown++
			continue
		}
		if *c.was != c.now {
			return fmt.Errorf("%s was baked with %s %v and the manifest now says %v. The heightmap on disk "+
				"means something different from what this run would read it as; rebake, or put the manifest "+
				"back", dir, c.what, *c.was, c.now)
		}
	}
	if meta.Seed != nil && *meta.Seed != m.Source.Seed {
		return fmt.Errorf("%s was baked with seed %d and the manifest now says %d; the detail passes would "+
			"be hashing a different world from the one in the heightmap", dir, *meta.Seed, m.Source.Seed)
	}
	if unknown > 0 {
		warn("warning %s predates this check and does not record %d of the numbers it would be checked "+
			"against; if the manifest has moved since it was baked, the heights will be read as something "+
			"they are not", dir, unknown)
	}
	return nil
}

// BakeTiles runs the detail passes over a rectangle of tiles and writes them.
func BakeTiles(opt TileOptions) (*TileIndex, error) {
	log := opt.Log
	if log == nil {
		log = func(string, ...any) {}
	}
	in := opt.In
	m := in.M
	cfg := m.Pipeline

	marginPx := detail.MarginCells(cfg.Particle)
	g, err := tile.NewGrid(in.P, cfg.GeologyFactor, cfg.Detail.TilePx, marginPx)
	if err != nil {
		return nil, err
	}
	detailCellM := in.P.CellM / float64(cfg.GeologyFactor)
	log("tiles   %d x %d of %d px at %.1f m (%.2f km), margin %d px (%.0f m)",
		g.NX, g.NY, cfg.Detail.TilePx, detailCellM,
		float64(cfg.Detail.TilePx)*detailCellM/1000, g.MarginGeo*g.Factor, float64(g.MarginGeo)*in.P.CellM)
	if want := int(cfg.Detail.ClassBlendM/in.P.CellM + 0.5); want > g.MarginGeo/2 {
		log("warning class_blend_m is %.0f m, which a tile cannot reach past its own margin; it will blend "+
			"over %.0f m instead. Two passes of the blur reach twice its radius, and the margin is %.0f m",
			cfg.Detail.ClassBlendM, float64(g.MarginGeo/2)*in.P.CellM, float64(g.MarginGeo)*in.P.CellM)
	}
	if cd := cfg.CoastDetail; cd.Enabled && !opt.NoShore {
		log("shore   the coastal detail pass is on: %.0f m of surf reach is %.0f cells here, a beach below "+
			"%.0f m of backshore and a cliff above %.0f", cfg.Coast.SurfReachM,
			cfg.Coast.SurfReachM/detailCellM, cd.CliffFromM, cd.CliffToM)
	} else {
		log("shore   the coastal detail pass is off; the shore is the geology upsampled")
	}

	all := g.Tiles()
	wanted := all[:0:0]
	for _, t := range all {
		if opt.OnlySet {
			if t.IX < opt.Only[0] || t.IX > opt.Only[2] || t.IY < opt.Only[1] || t.IY > opt.Only[3] {
				continue
			}
		}
		wanted = append(wanted, t)
	}
	if len(wanted) == 0 {
		return nil, fmt.Errorf("no tiles selected; the grid is %d x %d", g.NX, g.NY)
	}
	if err := os.MkdirAll(opt.Out, 0o755); err != nil {
		return nil, err
	}

	prefix := opt.Prefix
	if prefix == "" {
		prefix = "Planet"
	}
	jobs := opt.Jobs
	if jobs <= 0 {
		jobs = 4
	}
	if jobs > len(wanted) {
		jobs = len(wanted)
	}

	recs := make([]TileRecord, len(wanted))
	errs := make([]error, len(wanted))
	var wg sync.WaitGroup
	var mu sync.Mutex
	next := make(chan int)
	go func() {
		for i := range wanted {
			next <- i
		}
		close(next)
	}()
	for w := 0; w < jobs; w++ {
		wg.Add(1)
		go func() {
			defer wg.Done()
			for i := range next {
				rec, err := bakeOneTile(g, wanted[i], opt, prefix)
				recs[i], errs[i] = rec, err
				mu.Lock()
				if err != nil {
					log("tile %s FAILED: %v", wanted[i].Name(prefix), err)
				} else {
					log("tile %s  %d x %d  %.0f..%.0f m  %.3f%% clipped  [%.1f s]",
						rec.File, rec.W, rec.H, rec.MinM, rec.MaxM, rec.ClipFrac*100, rec.Seconds)
				}
				mu.Unlock()
			}
		}()
	}
	wg.Wait()
	for _, err := range errs {
		if err != nil {
			return nil, err
		}
	}

	idx := &TileIndex{
		When: time.Now().UTC().Truncate(time.Second), Prefix: prefix,
		CellM: detailCellM, TilePx: cfg.Detail.TilePx, MarginPx: g.MarginGeo * g.Factor,
		NX: g.NX, NY: g.NY, WrapX: true,
		WorldWM: in.P.CircumferenceM(), WorldHM: in.P.HeightM(),
		Tiles: recs,
	}
	idx.ElevationM.Min, idx.ElevationM.Max = m.ElevationM.Min, m.ElevationM.Max
	data, err := json.MarshalIndent(idx, "", "  ")
	if err != nil {
		return nil, err
	}
	if err := os.WriteFile(filepath.Join(opt.Out, "tiles.json"), append(data, '\n'), 0o644); err != nil {
		return nil, err
	}
	// The key to every *_overlay.png, plus the features in world metres, written beside them so an importer
	// reads one directory rather than two. It is the same document the plan and the bake write; it is small,
	// it describes the whole planet, and a tile batch that did not carry it would be a folder of masks with
	// no legend.
	if in.OverlayDoc != nil {
		if err := in.OverlayDoc.WriteJSON(opt.Out); err != nil {
			return nil, err
		}
	}
	return idx, nil
}

// bakeOneTile is passes 8 to 12 and 14 over one tile.
func bakeOneTile(g *tile.Grid, t tile.Tile, opt TileOptions, prefix string) (TileRecord, error) {
	in := opt.In
	m := in.M
	cfg := m.Pipeline

	// Pass 8: cut with the margin and upsample. UpsampleInt is exact-factor Catmull-Rom, so every geology
	// sample lands exactly on a detail sample and there is no phase error to accumulate along a tile row.
	//
	// The sea is flattened to sea level *before* the upsample, not after, and both halves of that matter. The
	// geology raster drops from the shore to the painted ocean depth in a single cell, so a Catmull-Rom
	// upsample of it rings at every coastline - hundreds of metres of overshoot in the water and a wave of it
	// back into the land. And with the sea flat, the interpolated height crosses sea level on a smooth
	// contour, so the detail land mask can be read off the height itself; taken up from the geology mask by
	// nearest neighbour instead, the coastline comes out as a staircase of 8 m blocks and it is plainly
	// visible in a hillshade.
	//
	// It is the same invariant the fluvial solve keeps, for the same reason: with the floor in place a cell at
	// the waterline stands five hundred metres above its neighbour, and thermal weathering would find the
	// whole coastline past the angle of repose and pour it into the sea.
	geo, _, _ := g.Cut(t, opt.HeightM, 0)
	geoSea := g.CutMask(t, opt.Sea, opt.HeightM.W, 0)
	floor := geo.Clone()
	for i, isSea := range geoSea {
		if isSea {
			geo.Data[i] = float32(m.SeaLevelM)
		}
	}

	h := geo.UpsampleInt(cfg.GeologyFactor)
	land := make([]bool, len(h.Data))
	for i, v := range h.Data {
		land[i] = float64(v) > m.SeaLevelM
	}

	f := g.Frame(t)
	periodM := m.Planet.DetailNoisePeriodKm * 1000

	classes := blendedClasses(g, t, opt, geo, h)

	if opt.NoDetail {
		restoreSeaFloor(h, land, geo, floor, m, cfg.GeologyFactor)
		return finishTile(g, t, opt, prefix, h, land, nil, nil)
	}

	// Pass 9.
	detail.RunDetailNoise(h, land, detail.DetailNoiseParams{
		Cfg: cfg.Detail, Seed: m.Source.Seed, Frame: f, PeriodM: periodM, SeaLevelM: m.SeaLevelM,
		Classes: classes,
	})

	// Pass 10 feeds pass 11 rather than standing alone: strata is hardness, and hardness is what the
	// droplets scale their cutting by, which is how a hard band ends up holding a shelf on a cut face.
	hard := detail.NewHardness(f, m.Source.Seed, m.Planet.NoisePeriodKm*1000,
		cfg.Strata.PeriodM, cfg.Strata.Contrast, classes)

	// Pass 11.
	maps, _ := detail.RunParticle(h, land, detail.ParticleParams{
		Cfg: cfg.Particle, Seed: m.Source.Seed, Frame: f, SeaLevelM: m.SeaLevelM, Hardness: hard,
		Classes: classes,
	})

	// Pass 12: the same mass-conserving weathering the coarse grid gets, at the cell size where scree and a
	// cliff face are actually resolved.
	fixed := make([]bool, len(land))
	for i := range land {
		fixed[i] = !land[i]
	}
	thermal.Apply(h.Data, h.W, h.H, h.CellM, thermal.TalusFromDegrees(cfg.Thermal.TalusDeg),
		cfg.Thermal.FinePasses, fixed, nil)

	// The sea floor goes back *here*, before the shore is drawn, rather than on the way out. Pass 11b works
	// on both sides of the waterline - a foreshore is below it and a berm is above it - so a shore laid onto
	// water that is about to be overwritten would be half a shore.
	restoreSeaFloor(h, land, geo, floor, m, cfg.GeologyFactor)

	// Pass 9b: the same texture as pass 9, under water, now that there is a sea bed to put it on. It runs
	// before the shore rather than after, so the beach the shore pass draws is smooth sand over it rather
	// than sand with noise on top.
	detail.RunSeabedNoise(h, detail.DetailNoiseParams{
		Cfg: cfg.Detail, Seed: m.Source.Seed, Frame: f, PeriodM: periodM, SeaLevelM: m.SeaLevelM,
		Classes: classes,
	})

	// Pass 11b: the shore. It runs last of the detail passes because marine processes are the last thing to
	// act on a coast and they act faster than anything inland: a berm is rebuilt by every tide, while the
	// hillslope creep that pass 12 stands for takes the age of the cliff behind it. Running it before the
	// fine thermal would have that creep immediately relax the one face on the map that is meant to be
	// steeper than the angle of repose.
	var coastal *detail.CoastalStats
	if cfg.CoastDetail.Enabled && !opt.NoShore {
		st := detail.RunCoastal(h, land, detail.CoastalParams{
			Cfg: cfg.CoastDetail, Surf: cfg.Coast, Seed: m.Source.Seed, Frame: f, PeriodM: periodM,
			SeaLevelM: m.SeaLevelM, Exposure: cutExposure(g, t, opt, geo, h), Hardness: hard,
		})
		coastal = &st
	}

	return finishTile(g, t, opt, prefix, h, land, maps, coastal)
}

// cutExposure lifts the coastal pass's fetch field onto this tile's detail grid, or nil when the bake did not
// carry one. Interpolated rather than nearest: it is a smooth field and a staircase in it would put a
// staircase into the berm height along every beach.
func cutExposure(g *tile.Grid, t tile.Tile, opt TileOptions, geo, h *field.Field) []float32 {
	if opt.Exposure == nil {
		return nil
	}
	cut, _, _ := g.Cut(t, opt.Exposure, 0)
	up := cut.UpsampleInt(opt.In.M.Pipeline.GeologyFactor)
	if len(up.Data) != len(h.Data) {
		return nil
	}
	return up.Data
}

// restoreSeaFloor puts the water back after the land passes, which ran with the sea flattened to sea level.
//
// It used to be nearest neighbour, unconditionally, and the comment said why: the geology raster dropped from
// the shore to the painted ocean depth in a single cell, and interpolating a five-hundred-metre step is
// exactly what the flattening exists to avoid. The cost was a four-fold staircase over the whole sea floor,
// which nobody could see while the shore was a cliff into five hundred metres of water.
//
// D-60 changed the input. There is a continental shelf now, and a surf-cut platform, and a beach, and between
// them they carry the sea floor down from the waterline to the shelf break over kilometres rather than over
// one cell. So the shallow water is interpolated - from the *unflattened* cut, which still holds the land
// heights, so the surface runs across the waterline with no seam in it - and only the drop past the break is
// still nearest. The two are blended over a depth band rather than switched between, because a hard switch
// would put back a smaller version of the step it exists to avoid.
func restoreSeaFloor(h *field.Field, land []bool, geo, floor *field.Field, m *manifest.Manifest, factor int) {
	breakM := m.ShelfBreakM()
	if breakM <= 0 {
		breakM = 30
	}
	// The clamp sits at the *far* end of the blend band rather than at the break, so that everywhere the blend
	// is still reading the interpolation, the interpolation is of the real sea floor. Clamped at the break
	// instead, the smooth half of the blend was a flat surface at break depth while the nearest half followed
	// the slope down, and the mixture lifted the floor by up to half the band - ten metres of invented shelf
	// in exactly the strip the blend exists to make invisible.
	const bandM = 40.0
	deepest := float32(m.SeaLevelM - (breakM + bandM))

	shallow := floor.Clone()
	for i, v := range shallow.Data {
		if v < deepest {
			shallow.Data[i] = deepest
		}
	}
	smooth := shallow.UpsampleInt(factor)

	// Not smoothed, and it is worth saying why not, because the first version was.
	//
	// The interpolated sea floor comes out of a hillshade covered in dotted contour lines, which look exactly
	// like an interpolation artefact and are not: measured, an eighty by hundred patch of open water takes
	// three distinct 8-bit shade values, 94 % of them the same one. It is the hillshade's own quantisation on
	// a surface that slopes at one in three hundred, it was there before and it is in the picture rather than
	// in the ground. A box blur over the floor was tried against it and changed the tile by a fifth of a
	// height quantum on average - its only real effect was to soften genuine one-cell steps in the geology,
	// which is not what it was for.
	for y := 0; y < h.H; y++ {
		sy := y / factor
		if sy >= geo.H {
			sy = geo.H - 1
		}
		for x := 0; x < h.W; x++ {
			i := y*h.W + x
			if land[i] {
				continue
			}
			sx := x / factor
			if sx >= geo.W {
				sx = geo.W - 1
			}
			near := float64(floor.Data[sy*geo.W+sx])
			depth := m.SeaLevelM - near
			t := (depth - breakM) / bandM
			if t <= 0 {
				h.Data[i] = smooth.Data[i]
				continue
			}
			if t >= 1 {
				h.Data[i] = float32(near)
				continue
			}
			w := noise.Smoothstep(t)
			h.Data[i] = float32((1-w)*float64(smooth.Data[i]) + w*near)
		}
	}
}

// finishTile restores the sea floor, crops the margin away and writes everything out.
func finishTile(g *tile.Grid, t tile.Tile, opt TileOptions, prefix string, h *field.Field, land []bool,
	maps *detail.Maps, coastal *detail.CoastalStats) (TileRecord, error) {

	start := time.Now()
	_ = land
	m := opt.In.M
	cfg := m.Pipeline

	// Pass 14: crop the margin away and write. Everything outside the interior was only ever there so the
	// passes above had somewhere to read from.
	ix, iy := g.MarginGeo*g.Factor, g.MarginGeo*g.Factor
	iw, ih := g.DetailW(t), g.DetailH(t)
	out := crop(h, ix, iy, iw, ih)

	rec := TileRecord{
		IX: t.IX, IY: t.IY, File: t.Name(prefix) + ".png", W: iw, H: ih,
		OriginXM: g.OriginXM(t), OriginYM: g.OriginYM(t),
	}
	if coastal != nil && coastal.ShoreCells > 0 {
		rec.Coastal = coastal
	}

	// The annotation layer, sampled onto this tile's interior. It is written before the height, because it is
	// the cheap one and a failure here should not leave a heightmap with no mask beside it.
	//
	// Nothing in the detail passes read it and nothing here consults it: it is the author's layer travelling
	// through to whatever builds the level. The values are mark indices, zero for nothing, and the key is in
	// overlay.json beside tiles.json.
	if ov := opt.In.OverlayRaster; ov != nil {
		marks := ov.SampleWorld(rec.OriginXM, rec.OriginYM, h.CellM, iw, ih, opt.In.OverlayScale())
		rec.OverlayFile = t.Name(prefix) + "_overlay.png"
		if err := overlay.WriteMask(filepath.Join(opt.Out, rec.OverlayFile), iw, ih, marks); err != nil {
			return rec, err
		}
	}
	lo, hi := out.MinMax()
	rec.MinM, rec.MaxM = float64(lo), float64(hi)
	rec.ClipFrac = m.ClipFraction(out.Data)

	if err := field.WriteGray16(filepath.Join(opt.Out, rec.File), iw, ih,
		m.Encode(out.Data), png.DefaultCompression); err != nil {
		return rec, err
	}
	// A hillshade beside the heightmap, at full resolution. A 16-bit grey PNG of a hundred metres of relief
	// is a flat grey rectangle to look at, and the whole reason these passes exist is what they do to the
	// surface - which cannot be judged from a number.
	// Shaded with the water clamped at the shelf break rather than at sea level. Clamping at sea level was
	// right while the shore was a step into five hundred metres of water and there was nothing below the
	// waterline worth looking at; now there is a shore platform, a foreshore and a beach down there, and
	// they are most of what pass 11b does. The break is still clamped, because a continental slope in the
	// corner of a tile would otherwise set the whole hillshade's contrast.
	shade := out.Clone()
	shadeFloor := float32(m.SeaLevelM - m.ShelfBreakM())
	for i := range shade.Data {
		if shade.Data[i] < shadeFloor {
			shade.Data[i] = shadeFloor
		}
	}
	if err := field.WriteHillshade(filepath.Join(opt.Out, t.Name(prefix)+"_shade.png"), shade, iw, 1); err != nil {
		return rec, err
	}
	// A slice, not a map: map iteration order is randomised in Go and nothing in this generator is allowed
	// to depend on it (cross-cutting rule 12). Here it would only reorder two file writes, which is exactly
	// the kind of "it does not matter this time" that makes the rule worth keeping without exception.
	// The full-scale values are fixed constants, not percentiles of the tile.
	//
	// Field.ToUnit takes the 99th percentile of whatever it is given, which is exactly right for one map of
	// one world and exactly wrong here: it is a statistic of the tile's own extent, so two tiles would stretch
	// by different anchors and their shared valley would come out two different greys. That is the same
	// mistake the coastal pass's exposure made and had withdrawn, and the same rule - no pass computes a
	// statistic of the piece of the world it happens to be looking at.
	//
	// Flow is water-units accumulated and runs over decades, so it is log-scaled; wear and deposit are metres
	// and a metre of either is a great deal at a 2 m cell.
	flowFull := 40 * cfg.Particle.DropletsPerCell * float64(cfg.Particle.Lifetime)
	var derived []struct {
		name string
		data []float32
		full float64
		log  bool
	}
	if maps != nil {
		derived = []struct {
			name string
			data []float32
			full float64
			log  bool
		}{
			{"flow", maps.Flow, flowFull, true},
			{"wear", maps.Wear, 1.0, false},
			{"deposit", maps.Deposit, 1.0, false},
		}
	}
	for _, d := range derived {
		c := crop(&field.Field{W: h.W, H: h.H, CellM: h.CellM, Data: d.data}, ix, iy, iw, ih)
		if err := field.WriteGray8(filepath.Join(opt.Out, t.Name(prefix)+"_"+d.name+".png"), iw, ih,
			toBytes(normalise(c.Data, d.full, d.log)), png.BestSpeed); err != nil {
			return rec, err
		}
	}
	rec.Seconds = time.Since(start).Seconds()
	return rec, nil
}

// upsampleClass takes the class raster up by an integer factor, nearest. A class index is a name and not a
// quantity: interpolating one would invent a class that is neither of its neighbours.
func crop(f *field.Field, x0, y0, w, h int) *field.Field {
	out := field.New(w, h, f.CellM)
	for y := 0; y < h; y++ {
		copy(out.Data[y*w:(y+1)*w], f.Data[(y0+y)*f.W+x0:(y0+y)*f.W+x0+w])
	}
	return out
}

// normalise maps values onto 0..1 against a fixed full-scale, never a percentile of the data. See the note
// where the constants are chosen.
func normalise(data []float32, full float64, logScale bool) []float32 {
	if full <= 0 {
		full = 1
	}
	top := full
	if logScale {
		top = math.Log1p(full)
	}
	out := make([]float32, len(data))
	for i, v := range data {
		x := float64(v)
		if x < 0 {
			x = 0
		}
		if logScale {
			x = math.Log1p(x)
		}
		out[i] = float32(x / top)
	}
	return out
}

func toBytes(data []float32) []uint8 {
	out := make([]uint8, len(data))
	for i, v := range data {
		x := v
		if x < 0 {
			x = 0
		} else if x > 1 {
			x = 1
		}
		out[i] = uint8(x*255 + 0.5)
	}
	return out
}

// blendedClasses turns the painted class raster into the four numbers the detail passes read, per cell, with
// the boundaries between classes faded rather than stepped.
//
// **Why the fade.** A class is a name and a name is never interpolated - the mask that travels to whatever
// builds the level is still nearest neighbour, and it has to be. But the numbers a class stands for are
// quantities. Kept as a lookup on the class index, a desert meeting a lowland went from seven metres of dune
// amplitude to two, and from a fifth of the running water to all of it, in the width of one cell, along a
// line somebody drew with a mouse. It read as what it was: a boundary in a picture rather than a change in
// the ground. Faded over `class_blend_m`, the same boundary is a few hundred metres of one becoming the
// other, which is what the edge of a sand sea looks like from inside it.
//
// **Why at the geology grid.** The class raster is a geology-resolution field, so blurring it there costs a
// four-hundredth of blurring at detail resolution, and the upsample afterwards is the same exact-factor
// Catmull-Rom every other field gets - so the result is smoother than a blur at detail resolution would have
// been, not coarser.
//
// **Why the radius is clamped.** A blur reads outside the cell it writes, and a tile only has its margin to
// read from. Two passes of a box blur of radius r reach 2r, so r is capped at half the margin and the run
// says so once when the manifest asks for more. Past that cap a tile would be blending against its own cut
// edge and two tiles would disagree about the same ground, which is the one thing the tiling may not do.
func blendedClasses(g *tile.Grid, t tile.Tile, opt TileOptions, geo, h *field.Field) *detail.Classes {
	in := opt.In
	cfg := in.M.Pipeline
	if !in.Legend.Overrides() {
		return nil
	}
	tbl := in.Legend.DetailTables(cfg.Particle.DropletsPerCell,
		cfg.Detail.AmplitudeM.Lo(), cfg.Detail.AmplitudeM.Hi(), cfg.Strata.Contrast)
	cls := g.CutClass(t, in.Map.Class, in.P.W, in.P.PadY)

	radius := int(cfg.Detail.ClassBlendM/in.P.CellM + 0.5)
	if max := g.MarginGeo / 2; radius > max {
		radius = max
	}
	lift := func(table []float64) []float32 {
		f := field.New(geo.W, geo.H, geo.CellM)
		for i, k := range cls {
			f.Data[i] = float32(table[k])
		}
		if radius > 0 {
			field.BoxSmooth(f.Data, f.W, f.H, radius, 2)
		}
		return f.UpsampleInt(cfg.GeologyFactor).Data
	}
	return &detail.Classes{
		Droplets: lift(tbl.Droplets), AmpLo: lift(tbl.AmpLo),
		AmpHi: lift(tbl.AmpHi), Contrast: lift(tbl.Contrast),
	}
}