UnrealPrototyping/Tools/Terrain/internal/region/region.go

// Package region cuts a planet into the pieces the geology solve runs on.
//
// Docs/Terrain-Next.md 3.3 says the fluvial solve cannot be tiled, and that is right: drainage area is an
// integral over the whole upstream catchment and the priority-flood needs global connectivity, so a river
// crossing a tile boundary would need the next tile's catchment to know how big it is.
//
// It can be decomposed per landmass, though, and that is a different statement. Ocean cells are held fixed
// at sea level for the entire solve - fluvial.ComputeReceivers makes every outlet its own receiver, so a
// receiver chain starting on land terminates the moment it steps into water, and StreamPower, both
// diffusions, the repose clamp and thermal all skip a fixed cell. No flow path crosses open water, and every
// basin is contained in one eight-connected land component. So solving a landmass in a box of its own is not
// an approximation of solving the planet whole: on land it is the same answer.
//
// What that buys is memory. The whole planet at once is a fluvial.Grid of about 35 bytes a cell plus the
// dozen full-size fields uplift builds, which at 78 million cells is several gigabytes before anything has
// been eroded. Landmasses plus a thin margin are a fraction of that area and are solved one at a time.
//
// What it costs is that the decomposition becomes part of the world's identity: the priority-flood's epsilon
// ladder across a flat depends on the flood's traversal order, which depends on the box it is flooding. The
// seed alone no longer names a world - the seed and the margin do - so the margin lives in the manifest and
// is recorded in meta.json.
//
// Note what is NOT decomposed. The coastal pass runs once on the whole cylinder, because it is cheap (tens
// of nanoseconds a cell, against tens of nanoseconds a cell *per step* for the solve) and because cutting it
// up would truncate the fetch across every strait, split the sediment budget whose conservation is the one
// thing in that pass not derived from something already measured, and leave the shoreline length and the
// exposure percentiles as statistics that do not pool. Decompose the solve, not the map.
package region

import (
	"fmt"

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

// Region is one piece of the planet: a landmass, or a cluster of landmasses close enough that they shelter
// each other, plus a margin of ocean on every side.
type Region struct {
	ID    int
	Frame world.Frame

	LandCells int  // painted land cells this region owns
	SetCells  int  // cells in the dilated set, which the frame is the bounding box of
	Seam      bool // the frame straddles x = 0
}

// Cells is the size of the grid the solve will run on, margin included.
func (r Region) Cells() int { return r.Frame.Cells() }

// Partition is a planet cut into regions, and the map from planet cell to owning region.
type Partition struct {
	P           world.Planet
	MarginCells int

	// Owner is the region id for every planet cell, or -1 for water that belongs to no region. A cell
	// inside one region's frame may be owned by another region or by nobody, which is what keeps two
	// regions from both solving the same island.
	Owner   []int32
	Regions []Region

	// Dropped counts the specks: components with less painted land than the minimum, returned to the sea.
	DroppedRegions, DroppedCells int
}

// Build partitions a classified planet.
//
// The land mask is dilated by the margin with one exact distance transform, and the connected components of
// the dilated mask are the regions. That is the whole rule, and it is deliberately not a bounding-box
// overlap test: dilated boxes are transitively closed and one long thin landmass has an enormous box, so on
// a real template box clustering collapses most of the map into a single region. Dilating the mask itself
// groups exactly those landmasses that come within a margin of each other.
//
// The bounding box of a dilated component is the region's frame, and its edges are ocean by construction: a
// land cell dilates to reach margin cells further out, so the outermost column and row of the dilated set
// are at least margin cells from any land in that component. That is the invariant TestBorderIsAlwaysOcean
// asserts about the square canvas, and the solve depends on it - a border cell is an outlet, and land
// sitting on one would freeze at its initial relief while the interior eroded out from under it.
func Build(m *template.Map, marginCells, minLandCells int) (*Partition, error) {
	p := m.P
	n := p.W * p.H
	if marginCells < 1 {
		return nil, fmt.Errorf("margin is %d cells; a region needs at least one ring of ocean", marginCells)
	}
	if p.PadY < marginCells {
		return nil, fmt.Errorf("the polar pad is %d rows against a %d cell margin; a cap touching the top "+
			"of the painted map would not get a full margin of ocean", p.PadY, marginCells)
	}

	part := &Partition{P: p, MarginCells: marginCells, Owner: minusOne(n)}

	land := make([]bool, n)
	anyLand := false
	for i := range m.Sea {
		land[i] = !m.Sea[i]
		anyLand = anyLand || land[i]
	}
	if !anyLand {
		return part, nil
	}

	near := dilate(land, p, marginCells)

	comp := make([]int32, n)
	for i := range comp {
		comp[i] = -1
	}
	var regionOfComp []int32 // one entry per component: the region index, or -1 when it was dropped
	var stack []int32
	cols := make([]bool, p.W)

	for start := 0; start < n; start++ {
		if !near[start] || comp[start] >= 0 {
			continue
		}
		id := int32(len(regionOfComp))
		regionOfComp = append(regionOfComp, -1)
		comp[start] = id
		stack = append(stack[:0], int32(start))

		for i := range cols {
			cols[i] = false
		}
		minY, maxY := p.H, -1
		setCells, landCells := 0, 0

		for len(stack) > 0 {
			c := stack[len(stack)-1]
			stack = stack[:len(stack)-1]
			cx, cy := int(c)%p.W, int(c)/p.W
			setCells++
			cols[cx] = true
			if cy < minY {
				minY = cy
			}
			if cy > maxY {
				maxY = cy
			}
			if land[c] {
				landCells++
			}
			for dy := -1; dy <= 1; dy++ {
				ny := cy + dy
				if ny < 0 || ny >= p.H {
					continue
				}
				base := ny * p.W
				for dx := -1; dx <= 1; dx++ {
					if dx == 0 && dy == 0 {
						continue
					}
					ni := int32(base + p.WrapX(cx+dx))
					if near[ni] && comp[ni] < 0 {
						comp[ni] = id
						stack = append(stack, ni)
					}
				}
			}
		}

		if landCells < minLandCells {
			// A speck: a stray paint pixel, or a lone cell the classifier left behind. Solving it would
			// spend a whole region on a rock, so it goes back to the sea and is counted.
			part.DroppedRegions++
			part.DroppedCells += landCells
			continue
		}

		x0, width := span(cols, p.W)
		if width >= p.W {
			return nil, fmt.Errorf("a landmass reaches all the way round the planet: %d of %d columns once "+
				"the %d cell margin is added. It cannot be flattened into a rectangle with ocean on both "+
				"sides, and the solve needs that, because a grid edge is an outlet. Break it with a strait, "+
				"or reduce the margin", width, p.W, marginCells)
		}

		regionOfComp[id] = int32(len(part.Regions))
		part.Regions = append(part.Regions, Region{
			ID:        len(part.Regions),
			Frame:     world.Frame{P: p, X0: x0, Y0: minY, W: width, H: maxY - minY + 1},
			LandCells: landCells,
			SetCells:  setCells,
			Seam:      x0+width > p.W,
		})
	}

	for i, c := range comp {
		if c >= 0 {
			part.Owner[i] = regionOfComp[c]
		}
	}
	return part, nil
}

// dilate marks every cell within margin cells of a seed, on the cylinder.
func dilate(seed []bool, p world.Planet, margin int) []bool {
	d2 := dt.Distance2(seed, p.W, p.H, true)
	reach := float32(margin * margin)
	out := make([]bool, len(d2))
	for i, d := range d2 {
		out[i] = d <= reach
	}
	return out
}

// Cut is the region's own view of the world: the class raster and the land mask for its frame, with every
// cell belonging to another region - or to no region - forced to sea.
//
// Forcing them is right rather than convenient. A neighbouring island inside this frame is a separate
// landmass with its own basins, and no flow path connects the two, so leaving it as land would solve it
// twice and let its relief leak into this region's statistics. As water it is exactly what it is to this
// region's rivers: base level.
func (p *Partition) Cut(m *template.Map, r Region) (class []uint8, land []bool) {
	sea := uint8(0)
	if i := m.L.FirstSea(); i >= 0 {
		sea = uint8(i)
	}
	class = make([]uint8, r.Frame.Cells())
	land = make([]bool, r.Frame.Cells())
	for y := 0; y < r.Frame.H; y++ {
		for x := 0; x < r.Frame.W; x++ {
			pi := r.Frame.PlanetIdx(x, y)
			o := y*r.Frame.W + x
			mine := p.Owner[pi] == int32(r.ID)
			if mine && !m.Sea[pi] {
				class[o] = m.Class[pi]
				land[o] = true
				continue
			}
			if m.Sea[pi] {
				class[o] = m.Class[pi] // keep the painted water class: its depth is read later
			} else {
				class[o] = sea // somebody else's land, which to this region is open water
			}
		}
	}
	return class, land
}

// Composite writes a region's solved land back into the planet raster.
//
// Only cells the region owns and that are painted land are written. Everything else in the frame is water,
// and the sea floor is the planetary coastal pass's to lay afterwards - a region must not write it, or two
// overlapping frames would disagree about the same stretch of shelf.
func (p *Partition) Composite(dst []float32, m *template.Map, r Region, src []float32) int {
	written := 0
	for y := 0; y < r.Frame.H; y++ {
		for x := 0; x < r.Frame.W; x++ {
			pi := r.Frame.PlanetIdx(x, y)
			if p.Owner[pi] != int32(r.ID) || m.Sea[pi] {
				continue
			}
			dst[pi] = src[y*r.Frame.W+x]
			written++
		}
	}
	return written
}

// span finds the shortest run of columns covering every occupied one, going round the cylinder. The largest
// gap decides: the run starts just after it.
func span(cols []bool, w int) (x0, width int) {
	occupied := make([]int, 0, w)
	for x, on := range cols {
		if on {
			occupied = append(occupied, x)
		}
	}
	if len(occupied) == 0 {
		return 0, 0
	}
	if len(occupied) == w {
		return 0, w
	}
	bestGap, bestAt := -1, 0
	for i := range occupied {
		var gap int
		if i == len(occupied)-1 {
			gap = occupied[0] + w - occupied[i]
		} else {
			gap = occupied[i+1] - occupied[i]
		}
		if gap > bestGap {
			bestGap, bestAt = gap, (i+1)%len(occupied)
		}
	}
	return occupied[bestAt], w - bestGap + 1
}

func minusOne(n int) []int32 {
	out := make([]int32, n)
	for i := range out {
		out[i] = -1
	}
	return out
}