UnrealPrototyping/Tools/Terrain/internal/overlay/generate.go

package overlay

import (
	"fmt"
	"math"
	"sort"

	"salty/terrain/internal/field"
	"salty/terrain/internal/noise"
)

// Filling an overlay in from a baked world, so an author starts from something rather than from nothing.
//
// The annotation layer is hand-painted and it starts blank, which is the right default and a bad starting
// point. Where a forest can grow, where a town would actually stand and what a road between two towns would
// follow are all *consequences of the terrain* - of slope, of where the rivers run, of how far the sea is -
// and the terrain is the one thing an author cannot see while painting, because the solve has not happened
// yet when they are painting classes and the heightmap is 29 million pixels when it has. So the generator
// reads a finished bake and proposes the marks the terrain implies. The author then moves them.
//
// Four rules, and they are the whole design:
//
//   - **A painted pixel is never touched.** Generation fills blank pixels only. An author who has drawn the
//     capital exactly where they want it can regenerate everything else around it as often as they like, and
//     the two halves compose rather than competing. This is what makes the feature safe to re-run.
//   - **It is opt-in per mark.** A mark generates only if it carries a `generate` block. A legend written
//     before this existed produces exactly the blank sheet it always did, and a mark the author wants to own
//     completely simply says nothing.
//   - **It runs at the template's resolution**, which is the overlay's own. Generating on the 8 m geology
//     grid and downsampling would smear a road across two colours, and the classifier reads exact colours -
//     a blended pixel is dropped or becomes a different mark. Nothing here antialiases anything, for the
//     same reason the studio's brush does not.
//   - **It proposes, it does not decide.** These are starting points. The numbers below are chosen to put
//     something plausible on the sheet, not to be a settlement model.
//
// What it is emphatically not: a simulation. There is no economy, no history and no climate here - the Go
// generator has no climate model at all - so "where would a city be" is answered with drainage, slope and
// distance to the sea, which is the part of the question the terrain can actually answer.

// Generate kinds. A mark's `generate.kind` picks one.
const (
	// GenForest fills ground that could carry trees: shallow enough, below the treeline, and broken up by a
	// noise field so it reads as woodland rather than as a contour band.
	GenForest = "forest"

	// GenSettlement places discs at scored sites - rivers, flat ground, the coast - with a minimum spacing,
	// largest tier first. Several marks may use it; they share one spacing rule, so a village never lands
	// inside a city.
	GenSettlement = "settlement"

	// GenRoad joins the settlements that were placed, along least-cost paths over the terrain. Water is
	// impassable, so roads never swim: an island group comes out as one road network per island.
	GenRoad = "road"

	// GenCoast bands the waterline. It is the one kind whose mark usually carries `coast_jitter`, which is
	// the only overlay property any pass reads.
	GenCoast = "coast"
)

// GenSpec is a mark's `generate` block: what to put where, and the few numbers worth varying. Every zero
// field takes a default that is derived from the world being generated rather than from a constant, because
// a treeline in metres means nothing until you know how high the land got.
type GenSpec struct {
	Kind string `json:"kind"`

	// MaxSlopeDeg is the steepest ground this mark will be put on. Forests stop at cliffs, towns stand on
	// flat ground, and roads climb but grudgingly.
	MaxSlopeDeg float64 `json:"max_slope_deg"`

	// MinHeightM and MaxHeightM bound the elevation band. MaxHeightM zero means "derive a treeline from the
	// land's own height distribution", which is the only honest default on a world whose relief is unknown
	// until it is baked.
	MinHeightM float64 `json:"min_height_m"`
	MaxHeightM float64 `json:"max_height_m"`

	// Cover is roughly the fraction of the eligible ground this mark should take, for area kinds. It is a
	// quantile of the noise field rather than a count, so it means the same thing on any size of world.
	Cover float64 `json:"cover"`

	// WavelengthKm is how big the patches are, for area kinds.
	WavelengthKm float64 `json:"wavelength_km"`

	// Count is how many of this mark to place, for settlements.
	Count int `json:"count"`

	// MinSpacingKm is how far apart settlements must stand. Shared across every settlement mark, taken from
	// the largest that sets one.
	MinSpacingKm float64 `json:"min_spacing_km"`

	// RadiusM is how big the painted blob is. Zero derives one from the mark's own min_area_px, so the blob
	// this writes is never one the feature reducer would then discard as a speck.
	RadiusM float64 `json:"radius_m"`

	// WidthM is how wide a band or a road is painted. For a road the legend's own width_m is used when this
	// is zero, because that is the same number said once.
	WidthM float64 `json:"width_m"`

	// CoastKm is how far inland a coast band reaches, and how close to the sea a settlement wants to be for
	// its coastal bonus.
	CoastKm float64 `json:"coast_km"`

	// OnlyClasses and NotClasses restrict a mark to, or bar it from, ground painted with named classes from
	// the *class* legend.
	//
	// They exist because height and slope cannot tell an ice cap from a meadow. The first run of this
	// generator grew woodland across both polar caps: the caps are flat, they are below the treeline, and
	// nothing the terrain knows says otherwise - the only thing that does is the colour the author painted
	// there. A class name that is not in the legend is an error rather than an empty filter, because a
	// misspelt exclusion is a forest on an ice cap that nobody notices.
	OnlyClasses []string `json:"only_classes"`
	NotClasses  []string `json:"not_classes"`

	// Resolved forms of the two lists above, as class indices. Filled in by Generate.
	onlyIdx map[int]bool
	notIdx  map[int]bool
}

// resolveClasses turns the class names into indices against the class legend that was actually loaded.
func (g *GenSpec) resolveClasses(markName string, names []string) error {
	find := func(list []string) (map[int]bool, error) {
		if len(list) == 0 {
			return nil, nil
		}
		if len(names) == 0 {
			return nil, fmt.Errorf("mark %q names classes, but no class legend was handed to the generator",
				markName)
		}
		out := map[int]bool{}
		for _, want := range list {
			found := -1
			for i, n := range names {
				if n == want {
					found = i
					break
				}
			}
			if found < 0 {
				return nil, fmt.Errorf("mark %q names the class %q, which is not in the class legend",
					markName, want)
			}
			out[found] = true
		}
		return out, nil
	}
	var err error
	if g.onlyIdx, err = find(g.OnlyClasses); err != nil {
		return err
	}
	g.notIdx, err = find(g.NotClasses)
	return err
}

func (g *GenSpec) validate(markName string) error {
	switch g.Kind {
	case GenForest, GenSettlement, GenRoad, GenCoast:
	default:
		return fmt.Errorf("mark %q: generate.kind %q is not one of %q, %q, %q, %q",
			markName, g.Kind, GenForest, GenSettlement, GenRoad, GenCoast)
	}
	if g.Cover < 0 || g.Cover > 1 {
		return fmt.Errorf("mark %q: generate.cover is %v, outside 0..1", markName, g.Cover)
	}
	if g.Count < 0 {
		return fmt.Errorf("mark %q: generate.count is %d", markName, g.Count)
	}
	return nil
}

// GenInputs is the baked world the marks are read off, at the overlay's own resolution.
type GenInputs struct {
	W, H  int
	CellM float64 // metres per overlay pixel

	// HeightM is the surface in metres and Sea is which cells are under water, both over the painted rows
	// only - the polar pad is scaffolding and has no marks on it.
	HeightM []float32
	Sea     []bool

	// FlowM2 is drainage area in square metres. Nil is allowed: rivers then contribute nothing to a
	// settlement's score, which is worth saying out loud rather than silently scoring zero everywhere.
	FlowM2 []float32

	// ClassAt is the class legend's index per cell, and ClassNames the names those indices mean. Both are
	// optional together: without them only_classes and not_classes cannot be honoured, and asking for one is
	// then an error rather than a filter that quietly does nothing.
	ClassAt    []uint8
	ClassNames []string

	Seed int64

	// Existing is the overlay as it stands. Its painted pixels are preserved exactly and generation fills
	// around them. Nil is a blank sheet.
	Existing *Raster
}

// GenReport is what was placed, for the run summary.
type GenReport struct {
	Marks      []GenMarkReport
	Kept       int // pixels that were already painted and were left alone
	Painted    int // pixels this generation filled
	TreelineM  float64
	Settlement []Placed
}

// GenMarkReport is one mark's share of a generation.
type GenMarkReport struct {
	Name   string
	Kind   string
	Cells  int
	Pieces int // settlements placed, or roads traced

	// Wanted is how many were asked for, when that is a number the legend gave. Reported separately from
	// Pieces so a run that could not fit them all says so: the spacing and the amount of flat ground are
	// what ration settlements, and an author who asked for forty and got eighteen needs to be told, not left
	// to count the dots.
	Wanted int
}

// Placed is one settlement, kept so the roads can be run between them and so the summary can say where they
// went.
type Placed struct {
	Mark   int // raster index
	X, Y   int
	Score  float64
	RadPx  int
	Region int // which connected landmass, so roads never try to cross open water
}

// wantedFor is how many of a mark the legend asked for, or zero when it is not a counted kind.
func wantedFor(m *Mark) int {
	if m.Generate == nil {
		return 0
	}
	return m.Generate.Count
}

// Generate fills the blank parts of an overlay from a baked world.
func (l *Legend) Generate(in GenInputs) (*Raster, GenReport, error) {
	var rep GenReport
	if in.W <= 0 || in.H <= 0 {
		return nil, rep, fmt.Errorf("overlay generation needs a size, got %dx%d", in.W, in.H)
	}
	if len(in.HeightM) != in.W*in.H || len(in.Sea) != in.W*in.H {
		return nil, rep, fmt.Errorf("overlay generation: height and sea must be %d cells", in.W*in.H)
	}
	for i := range l.Marks {
		if g := l.Marks[i].Generate; g != nil {
			if err := g.validate(l.Marks[i].Name); err != nil {
				return nil, rep, err
			}
			if err := g.resolveClasses(l.Marks[i].Name, in.ClassNames); err != nil {
				return nil, rep, err
			}
		}
	}
	if in.ClassAt != nil && len(in.ClassAt) != in.W*in.H {
		return nil, rep, fmt.Errorf("overlay generation: the class raster is %d cells and the grid is %d",
			len(in.ClassAt), in.W*in.H)
	}

	out := &Raster{W: in.W, H: in.H, Mark: make([]uint8, in.W*in.H)}
	// What was on the sheet before this run, kept separately from what is on it now. The distinction is the
	// whole layering rule: a hand-painted pixel is never touched, while a mark this run has just put down
	// may be built over by a later one - a road through generated woodland is a road, and a town on it is a
	// town. Without the two being different, whichever kind painted first would block every kind after it,
	// which is exactly what happened on the first run: a coastal band claimed a fifth of the world and the
	// settlements and roads placed inside it painted nothing at all.
	protectedPx := make([]bool, in.W*in.H)
	if in.Existing != nil {
		if in.Existing.W != in.W || in.Existing.H != in.H {
			return nil, rep, fmt.Errorf("the overlay on disk is %dx%d and the generator is working at %dx%d",
				in.Existing.W, in.Existing.H, in.W, in.H)
		}
		copy(out.Mark, in.Existing.Mark)
		for i, m := range out.Mark {
			if m != Blank {
				protectedPx[i] = true
				rep.Kept++
			}
		}
	}

	d := newGenData(in)
	d.protected = protectedPx
	rep.TreelineM = d.treelineM

	// Painting order is coarse to fine: the coastal band, then woodland, then the roads across it, then the
	// settlements the roads run between.
	//
	// Two area marks never overwrite each other - the first in legend order claims the overlap, because
	// deciding that a forest beats a coastline or the reverse is an authoring judgement and not one a
	// generator should make silently. Roads and settlements do overwrite generated areas, because they are
	// the thing being placed and the area is the ground it stands on.
	order := []string{GenCoast, GenForest, GenRoad, GenSettlement}
	byKind := map[string][]int{}
	for i := range l.Marks {
		if g := l.Marks[i].Generate; g != nil {
			byKind[g.Kind] = append(byKind[g.Kind], i)
		}
	}

	// Settlements are placed before the roads are drawn even though they are painted after, because the
	// roads are the paths between them and cannot be traced until they exist.
	if len(byKind[GenSettlement]) > 0 {
		rep.Settlement = l.placeSettlements(byKind[GenSettlement], d)
	}

	for _, kind := range order {
		for _, mi := range byKind[kind] {
			m := &l.Marks[mi]
			idx := uint8(mi + 1)
			var cells, pieces int
			switch kind {
			case GenCoast:
				cells = l.paintCoastBand(m, d, out, idx)
			case GenForest:
				cells = l.paintForest(m, d, out, idx)
			case GenRoad:
				cells, pieces = l.paintRoads(m, d, out, idx, rep.Settlement)
			case GenSettlement:
				cells, pieces = paintSettlements(d, out, idx, rep.Settlement, l.MinArea(m))
			}
			rep.Marks = append(rep.Marks, GenMarkReport{
				Name: m.Name, Kind: kind, Cells: cells, Pieces: pieces, Wanted: wantedFor(m),
			})
			rep.Painted += cells
		}
	}
	return out, rep, nil
}

// genData is everything derived once and shared by the kinds: slope, distance to the sea, the treeline and
// the landmass labels.
type genData struct {
	in GenInputs

	// protected marks the pixels that were already painted when this run started. Nothing here may write to
	// one, whatever kind it is.
	protected []bool

	slopeDeg  []float32
	coastKm   []float32 // distance to the nearest sea cell, kilometres; land only
	landID    []int32   // connected landmass, -1 at sea
	treelineM float64
	landMaxM  float64
	flowLog   []float32 // log10 of drainage area, normalised 0..1 over the land
}

func newGenData(in GenInputs) *genData {
	d := &genData{in: in}
	d.slopeDeg = slopeField(in.HeightM, in.W, in.H, in.CellM)
	d.coastKm = coastDistanceKm(in.Sea, in.W, in.H, in.CellM)
	d.landID = labelLandmasses(in.Sea, in.W, in.H)

	// The treeline is a quantile of the land's own heights rather than a number in metres, because a metre
	// means nothing until the world is baked: the same legend over a 47 m plain and a 2800 m range has to
	// put trees on both. Two thirds of the way up leaves the summits bare on a world that has summits and
	// takes almost nothing off a world that does not - which is correct, a lowland has no treeline.
	var hs []float32
	for i, s := range in.Sea {
		if !s {
			hs = append(hs, in.HeightM[i])
		}
	}
	if len(hs) > 0 {
		sort.Slice(hs, func(a, b int) bool { return hs[a] < hs[b] })
		d.landMaxM = float64(hs[len(hs)-1])
		d.treelineM = float64(hs[int(float64(len(hs)-1)*0.94)])
	}

	if in.FlowM2 != nil && len(in.FlowM2) == in.W*in.H {
		d.flowLog = make([]float32, in.W*in.H)
		cell := in.CellM * in.CellM
		// Normalised against a trunk river's catchment rather than the map's largest, so one enormous basin
		// cannot flatten every other river to nothing.
		hi := math.Log10(math.Max(cell*4, 5e7))
		lo := math.Log10(math.Max(cell, 1))
		for i, f := range in.FlowM2 {
			if in.Sea[i] || f <= 0 {
				continue
			}
			t := (math.Log10(float64(f)) - lo) / (hi - lo)
			d.flowLog[i] = float32(math.Max(0, math.Min(1, t)))
		}
	}
	return d
}

// slopeField is the surface gradient in degrees, central differences, X wrapped because the world is a
// cylinder and Y clamped because it is not a sphere.
func slopeField(h []float32, w, hgt int, cellM float64) []float32 {
	out := make([]float32, w*hgt)
	field.Rows(hgt, func(y0, y1 int) {
		for y := y0; y < y1; y++ {
			ym := y - 1
			if ym < 0 {
				ym = 0
			}
			yp := y + 1
			if yp >= hgt {
				yp = hgt - 1
			}
			for x := 0; x < w; x++ {
				xm := (x - 1 + w) % w
				xp := (x + 1) % w
				dzdx := float64(h[y*w+xp]-h[y*w+xm]) / (2 * cellM)
				dzdy := float64(h[yp*w+x]-h[ym*w+x]) / (2 * cellM)
				out[y*w+x] = float32(math.Atan(math.Hypot(dzdx, dzdy)) * 180 / math.Pi)
			}
		}
	})
	return out
}

// coastDistanceKm is how far each land cell is from the sea, by a multi-source breadth-first walk over the
// eight neighbours with X wrapped. Hop distance rather than Euclidean: it is a score input, and a BFS over
// 29 million cells costs one pass where a distance transform costs several.
func coastDistanceKm(sea []bool, w, h int, cellM float64) []float32 {
	out := make([]float32, w*h)
	for i := range out {
		out[i] = -1
	}
	queue := make([]int32, 0, w*8)
	for i, s := range sea {
		if s {
			continue
		}
		x, y := i%w, i/w
		if touchesSea(sea, w, h, x, y) {
			out[i] = 0
			queue = append(queue, int32(i))
		}
	}
	hop := float32(cellM / 1000)
	for head := 0; head < len(queue); head++ {
		c := int(queue[head])
		cx, cy := c%w, c/w
		d := out[c] + hop
		for _, o := range neighbours8 {
			nx := (cx + o[0] + w) % w
			ny := cy + o[1]
			if ny < 0 || ny >= h {
				continue
			}
			n := ny*w + nx
			if sea[n] || out[n] >= 0 {
				continue
			}
			out[n] = d
			queue = append(queue, int32(n))
		}
	}
	return out
}

var neighbours8 = [8][2]int{{-1, -1}, {0, -1}, {1, -1}, {-1, 0}, {1, 0}, {-1, 1}, {0, 1}, {1, 1}}

func touchesSea(sea []bool, w, h, x, y int) bool {
	for _, o := range neighbours8 {
		nx := (x + o[0] + w) % w
		ny := y + o[1]
		if ny < 0 || ny >= h {
			continue
		}
		if sea[ny*w+nx] {
			return true
		}
	}
	return false
}

// labelLandmasses numbers the connected land components, X wrapped, so a landmass across the seam is one
// landmass. Roads are built per component, which is what stops them crossing open water.
func labelLandmasses(sea []bool, w, h int) []int32 {
	out := make([]int32, w*h)
	for i := range out {
		out[i] = -1
	}
	var stack []int32
	next := int32(0)
	for start := range sea {
		if sea[start] || out[start] >= 0 {
			continue
		}
		id := next
		next++
		out[start] = id
		stack = append(stack[:0], int32(start))
		for len(stack) > 0 {
			c := int(stack[len(stack)-1])
			stack = stack[:len(stack)-1]
			cx, cy := c%w, c/w
			for _, o := range neighbours8 {
				nx := (cx + o[0] + w) % w
				ny := cy + o[1]
				if ny < 0 || ny >= h {
					continue
				}
				n := ny*w + nx
				if sea[n] || out[n] >= 0 {
					continue
				}
				out[n] = id
				stack = append(stack, int32(n))
			}
		}
	}
	return out
}

// eligible is the shared test every kind starts from: on land, not too steep, inside the height band.
func (d *genData) eligible(i int, g *GenSpec, maxDefault float64) bool {
	if d.in.Sea[i] {
		return false
	}
	maxSlope := g.MaxSlopeDeg
	if maxSlope <= 0 {
		maxSlope = maxDefault
	}
	if float64(d.slopeDeg[i]) > maxSlope {
		return false
	}
	if !d.classAllows(i, g) {
		return false
	}
	hm := float64(d.in.HeightM[i])
	if hm < g.MinHeightM {
		return false
	}
	top := g.MaxHeightM
	if top <= 0 {
		top = d.treelineM
	}
	return top <= 0 || hm <= top
}

// classAllows applies a mark's only_classes and not_classes to one cell.
func (d *genData) classAllows(i int, g *GenSpec) bool {
	if d.in.ClassAt == nil || (g.onlyIdx == nil && g.notIdx == nil) {
		return true
	}
	c := int(d.in.ClassAt[i])
	if g.notIdx != nil && g.notIdx[c] {
		return false
	}
	if g.onlyIdx != nil && !g.onlyIdx[c] {
		return false
	}
	return true
}

// paintForest fills eligible ground where a noise field stands above a quantile, so woodland has an outline
// rather than a contour edge. The field is indexed by world position (cross-cutting rule 1), so the same
// ground gets the same trees whatever else changes.
func (l *Legend) paintForest(m *Mark, d *genData, out *Raster, idx uint8) int {
	g := m.Generate
	cover := g.Cover
	if cover <= 0 {
		cover = 0.45
	}
	wavelengthKm := g.WavelengthKm
	if wavelengthKm <= 0 {
		wavelengthKm = 6
	}
	in := d.in
	circM := float64(in.W) * in.CellM
	u, v := noise.WorldUV(in.W, in.H, in.CellM, 0, 0, math.Max(circM, 1))
	cells := math.Max(1, math.Round(circM/(wavelengthKm*1000)))
	f := noise.FBMAt(u, v, noise.NewSource(in.Seed, srcOverlayForest),
		noise.Params{BaseCells: int(cells), Octaves: 4, Gain: 0.5})

	// The threshold is a quantile of the noise *over the eligible ground*, so `cover` means what it says on a
	// world whose eligible ground is a thin strip as much as on one where it is everything.
	// The quantile is taken over the ground this mark can actually take - eligible and not already claimed -
	// so `cover` means the same fraction whether or not another area mark got there first.
	var vals []float32
	for i := range out.Mark {
		if out.Mark[i] == Blank && d.eligible(i, g, 25) {
			vals = append(vals, f.Data[i])
		}
	}
	if len(vals) == 0 {
		return 0
	}
	sort.Slice(vals, func(a, b int) bool { return vals[a] < vals[b] })
	cut := vals[int(float64(len(vals)-1)*(1-cover))]

	n := 0
	for i := range out.Mark {
		if out.Mark[i] != Blank || f.Data[i] < cut || !d.eligible(i, g, 25) {
			continue
		}
		out.Mark[i] = idx
		n++
	}
	return n
}

// paintCoastBand marks a strip inland of the waterline. Its usual purpose is to carry `coast_jitter`, so it
// deliberately follows the shore rather than any other feature.
func (l *Legend) paintCoastBand(m *Mark, d *genData, out *Raster, idx uint8) int {
	g := m.Generate
	reachKm := g.CoastKm
	if reachKm <= 0 {
		if g.WidthM > 0 {
			reachKm = g.WidthM / 1000
		} else {
			reachKm = 1.5
		}
	}
	n := 0
	for i := range out.Mark {
		if out.Mark[i] != Blank || d.in.Sea[i] {
			continue
		}
		if c := d.coastKm[i]; c >= 0 && float64(c) <= reachKm {
			out.Mark[i] = idx
			n++
		}
	}
	return n
}

// placeSettlements scores the land and takes the best sites, largest tier first, with one spacing rule
// shared by every settlement mark so a village never lands inside a city.
//
// The score is the part of "where would a town be" that terrain can answer: fresh water, flat ground, and
// the sea. Everything else about a settlement - trade, history, who won a war - is the author's, which is
// why these are proposals in an editable sheet rather than a placement the bake bakes in.
func (l *Legend) placeSettlements(marks []int, d *genData) []Placed {
	in := d.in
	spacingKm := 0.0
	for _, mi := range marks {
		if s := l.Marks[mi].Generate.MinSpacingKm; s > spacingKm {
			spacingKm = s
		}
	}
	if spacingKm <= 0 {
		spacingKm = 4
	}
	spacingPx := math.Max(2, spacingKm*1000/in.CellM)

	// Tiers in the order the legend lists them, which is how an author already writes them: city, town,
	// village. The first listed takes the best sites.
	type tier struct {
		mi    int
		g     *GenSpec
		radPx int
	}
	var tiers []tier
	for _, mi := range marks {
		g := l.Marks[mi].Generate
		radM := g.RadiusM
		if radM <= 0 {
			// Big enough that the feature reducer will not drop it as a speck. The margin is generous on
			// purpose: a disc loses area wherever it meets ground that is already painted, and a settlement
			// that came out just under its own min_area_px would be placed, reported, and then silently
			// dropped by the feature pass - which is what happened to a city on the first real run. 1.6
			// linear is 2.6x the area, so it survives losing more than half of itself.
			minArea := float64(l.MinArea(&l.Marks[mi]))
			radM = math.Sqrt(minArea/math.Pi) * in.CellM * 1.6
		}
		radPx := int(math.Max(1, math.Round(radM/in.CellM)))
		tiers = append(tiers, tier{mi: mi, g: g, radPx: radPx})
	}

	// Candidates are taken on a stride rather than from every cell: two sites a quarter of the spacing apart
	// are the same site, and sorting 29 million scores to throw away all but fifty is work for nothing.
	stride := int(math.Max(1, math.Floor(spacingPx/4)))
	type cand struct {
		i     int
		score float64
	}
	var cands []cand
	for y := 0; y < in.H; y += stride {
		for x := 0; x < in.W; x += stride {
			i := y*in.W + x
			s := d.settlementScore(i)
			if s > 0 {
				// The seed picks among the plausible sites; the terrain decides which sites are plausible at
				// all. Without this the score is a pure function of the ground, so every press of the
				// studio's generate button proposes exactly the same towns and a re-roll re-rolls nothing.
				// A third either way reshuffles the ranking among comparable ground while still leaving a
				// river mouth on a plain beating a hillside.
				s *= 1 + settlementJitter*(hash01(uint64(i), uint64(in.Seed))-0.5)
				cands = append(cands, cand{i: i, score: s})
			}
		}
	}
	// Sorted by score, ties broken by index so the result does not depend on the sort's stability.
	sort.Slice(cands, func(a, b int) bool {
		if cands[a].score != cands[b].score {
			return cands[a].score > cands[b].score
		}
		return cands[a].i < cands[b].i
	})

	var placed []Placed
	taken := make([][2]int, 0, 64)
	sp2 := spacingPx * spacingPx
	farEnough := func(x, y int) bool {
		for _, t := range taken {
			dx := float64(wrapDelta(x-t[0], in.W))
			dy := float64(y - t[1])
			if dx*dx+dy*dy < sp2 {
				return false
			}
		}
		return true
	}
	for _, t := range tiers {
		want := t.g.Count
		if want <= 0 {
			continue
		}
		got := 0
		maxSlope := t.g.MaxSlopeDeg
		if maxSlope <= 0 {
			maxSlope = 8
		}
		for _, c := range cands {
			if got >= want {
				break
			}
			if float64(d.slopeDeg[c.i]) > maxSlope {
				continue
			}
			x, y := c.i%in.W, c.i/in.W
			if !farEnough(x, y) {
				continue
			}
			taken = append(taken, [2]int{x, y})
			placed = append(placed, Placed{
				Mark: t.mi + 1, X: x, Y: y, Score: c.score, RadPx: t.radPx,
				Region: int(d.landID[c.i]),
			})
			got++
		}
	}
	return placed
}

// settlementJitter is how far the seed may move a site's score, as a fraction. Large enough that the
// ranking among comparable ground genuinely reshuffles between presses, small enough that a site three times
// better than its neighbour still wins every time.
const settlementJitter = 0.65

// hash01 is a deterministic value in [0,1) from two integers: splitmix64 finalised. Not a stream, so it does
// not matter which order the cells are visited in, which is cross-cutting rule 12.
func hash01(a, b uint64) float64 {
	x := a*0x9e3779b97f4a7c15 + b*0xbf58476d1ce4e5b9
	x ^= x >> 30
	x *= 0xbf58476d1ce4e5b9
	x ^= x >> 27
	x *= 0x94d049bb133111eb
	x ^= x >> 31
	return float64(x>>11) / float64(1<<53)
}

// settlementScore is 0 where nobody would build and rises with the three things the terrain knows.
//
// Ground that is already painted scores zero, which is not a judgement about the ground: a site there cannot
// be stamped, because nothing may overwrite a hand-painted pixel. Scoring it anyway is how a settlement gets
// placed, counted and reported and then paints nothing at all - measured on the shipped template, one city of
// three and six villages of eighteen came out as empty blobs that the feature pass then dropped, so the run
// summary and `terrain plan` disagreed with each other and neither was wrong.
func (d *genData) settlementScore(i int) float64 {
	if d.in.Sea[i] || d.protected[i] {
		return 0
	}
	slope := float64(d.slopeDeg[i])
	if slope > 12 {
		return 0
	}
	flat := 1 - slope/12

	river := 0.0
	if d.flowLog != nil {
		river = float64(d.flowLog[i])
	}

	// A coast bonus that falls off over a few kilometres: a harbour is worth a great deal, being forty
	// kilometres inland is worth nothing either way.
	coast := 0.0
	if c := d.coastKm[i]; c >= 0 {
		coast = math.Max(0, 1-float64(c)/5)
	}

	// Flat ground is a precondition rather than an attraction, so it multiplies; water and the sea are the
	// reasons to be here, so they add.
	return flat * (0.15 + 1.5*river + 1.0*coast)
}

// paintSettlements stamps each placed site, growing the disc until the blob is big enough to survive the
// feature pass.
//
// The growth loop is not a flourish. A disc loses whatever part of itself falls on a coastline somebody has
// already painted, or on the sea, and settlements are scored *towards* the coast, so the loss is routine
// rather than rare. Without it the generator places a town, reports it, writes it, and the feature reducer
// then drops it as a speck - so `terrain plan` lists fewer settlements than the run said it made, with
// nothing anywhere to explain the difference. Measured on the shipped template: three cities placed and two
// reported, eighteen villages placed and twelve reported.
//
// It gives up after a few tries rather than growing without limit: a site hemmed in on every side is telling
// you it is a bad site, and a village the size of a county is worse than a missing one.
func paintSettlements(d *genData, out *Raster, idx uint8, placed []Placed, minArea int) (int, int) {
	n, pieces := 0, 0
	for _, p := range placed {
		if uint8(p.Mark) != idx {
			continue
		}
		pieces++
		got, r := 0, p.RadPx
		for try := 0; try < 4; try++ {
			// Re-stamping a larger disc only adds the new ring, because the cells already taken carry this
			// mark, so the area accumulates rather than being recounted.
			got += stampDisc(out, d, p.X, p.Y, r, idx, true)
			if got >= minArea {
				break
			}
			r = int(math.Ceil(float64(r) * 1.5))
		}
		n += got
	}
	return n, pieces
}

// stampDisc paints a filled circle, wrapping in X.
//
// overArea says whether this mark may cover ground another generated mark has already taken. A hand-painted
// pixel is never covered either way, which is what keeps a drawn stroke intact underneath a generated town.
func stampDisc(out *Raster, d *genData, cx, cy, r int, idx uint8, overArea bool) int {
	n := 0
	r2 := r * r
	for dy := -r; dy <= r; dy++ {
		y := cy + dy
		if y < 0 || y >= out.H {
			continue
		}
		for dx := -r; dx <= r; dx++ {
			if dx*dx+dy*dy > r2 {
				continue
			}
			x := ((cx+dx)%out.W + out.W) % out.W
			i := y*out.W + x
			if d.in.Sea[i] || d.protected[i] || (!overArea && out.Mark[i] != Blank) {
				continue
			}
			out.Mark[i] = idx
			n++
		}
	}
	return n
}

func wrapDelta(d, w int) int {
	if d > w/2 {
		d -= w
	} else if d < -w/2 {
		d += w
	}
	return d
}

// srcOverlayForest is this pass's seeded noise stream. It sits above the detail passes' 40s and the
// tectonic 50s so that adding one here cannot reshuffle any existing field.
const srcOverlayForest = 60