UnrealPrototyping/Tools/Terrain/internal/plates/boundary.go

package plates

import (
	"math"
	"sort"
)

// What a boundary does, which is the whole point of the package: "two plates hit each other" is one of these
// five and the other four are what happens when they do something else.
//
// The kind is per *vertex*, not per boundary. A margin whose plates are rotating as well as translating
// closes at one end and slides at the other - that is why the pole is in the map plane at all - so a single
// label for the whole line would throw away the thing the model was built to produce.
type Kind uint8

const (
	// Transform: the relative motion is along the line rather than across it. Little uplift, a strike-slip
	// fault, and a restraining bend that pops a range up where the line curves into the motion.
	Transform Kind = iota

	// Collision: convergent, both sides continental. Neither can subduct, so the crust thickens and the
	// result is a wide doubly-vergent belt - the thing an author means when they paint a mountain range.
	Collision

	// Subduction: convergent with at least one oceanic side. The ocean floor goes under, and the uplift is
	// an arc on the *overriding* plate, set back from the trench rather than centred on the line.
	Subduction

	// Rift: divergent, both sides continental. The axis drops and the shoulders stand up - the East African
	// pattern, and the one kind of boundary that lowers ground rather than raising it.
	Rift

	// Ridge: divergent with an oceanic side. A bathymetric ridge under water; on land it is a rift that has
	// already opened.
	Ridge
)

func (k Kind) String() string {
	switch k {
	case Collision:
		return "collision"
	case Subduction:
		return "subduction"
	case Rift:
		return "rift"
	case Ridge:
		return "ridge"
	default:
		return "transform"
	}
}

// Convergent reports whether this kind is two plates closing on each other.
func (k Kind) Convergent() bool { return k == Collision || k == Subduction }

// Divergent reports whether this kind is two plates separating.
func (k Kind) Divergent() bool { return k == Rift || k == Ridge }

// Vertex is one point on a boundary and everything a later pass reads off it.
type Vertex struct {
	XM float64 `json:"x_m"`
	YM float64 `json:"y_m"`

	// NX, NY is the unit normal, pointing out of plate A and into plate B. Every sign in this package is
	// measured against it, so "which side goes up" has one definition rather than one per consumer.
	NX float64 `json:"nx"`
	NY float64 `json:"ny"`

	// ClosingMYr is the relative velocity's component along the normal, in metres a year: positive closing,
	// negative opening. This is the number an uplift rate is a function of - "when two plates hit each other
	// they create mountains" is this field and nothing else.
	ClosingMYr float64 `json:"closing_m_yr"`

	// SlipMYr is the component along the line, signed in the polyline's own direction.
	SlipMYr float64 `json:"slip_m_yr"`

	Kind Kind `json:"kind"`

	// Over is the overriding plate at a subduction margin - the side the arc is built on - and -1 anywhere
	// else.
	Over int `json:"over"`
}

// Boundary is one continuous stretch of contact between two plates.
//
// X is **unwrapped**, exactly as uplift.FaultTrace is and for exactly the same reason: a boundary that
// crosses the seam has X running past the circumference or below zero rather than jumping, so every segment
// is a straight line between neighbouring points and no consumer has to special-case the meridian.
type Boundary struct {
	A int `json:"a"`
	B int `json:"b"`

	V []Vertex `json:"vertices"`
}

// LengthM is how long the boundary is, following the line.
func (b Boundary) LengthM() float64 {
	total := 0.0
	for i := 0; i+1 < len(b.V); i++ {
		total += math.Hypot(b.V[i+1].XM-b.V[i].XM, b.V[i+1].YM-b.V[i].YM)
	}
	return total
}

// Dominant is the kind most of this boundary's length is, which is the one word to print for it.
func (b Boundary) Dominant() Kind {
	var byKind [5]float64
	for i := 0; i+1 < len(b.V); i++ {
		d := math.Hypot(b.V[i+1].XM-b.V[i].XM, b.V[i+1].YM-b.V[i].YM)
		byKind[b.V[i].Kind] += d
	}
	best, bestK := -1.0, Transform
	for k, d := range byKind {
		if d > best {
			best, bestK = d, Kind(k)
		}
	}
	return bestK
}

// LengthByKind totals the planet's boundary length in each kind, in metres: the summary a run prints and the
// one number that says whether a seed produced a world with mountains in it.
func LengthByKind(bs []Boundary) [5]float64 {
	var out [5]float64
	for _, b := range bs {
		for i := 0; i+1 < len(b.V); i++ {
			d := math.Hypot(b.V[i+1].XM-b.V[i].XM, b.V[i+1].YM-b.V[i].YM)
			out[b.V[i].Kind] += d
		}
	}
	return out
}

// sample is one crossing of the boundary on the tectonic grid: the midpoint of two adjacent cells that
// belong to different plates.
type sample struct {
	xM, yM float64
	// dx, dy is the step from the plate-A cell towards the plate-B cell, which is what fixes the normal's
	// sign once the chain has a tangent to make it perpendicular to.
	dx, dy float64
	a, b   int
}

// minChainSamples is how short a chain is allowed to be before it is dropped. Triple junctions leave stubs
// of two or three cells that are a corner of the partition rather than a margin, and a stub cannot be given
// a meaningful tangent.
const minChainSamples = 6

// maxGapCells is how far apart two samples may be and still be the same line. Along a straight run they are
// one cell apart and on a staircase 0.71, so 1.6 chains both without reaching a parallel strand.
const maxGapCells = 1.6

// smoothPasses is how many times the chained polyline is averaged with its own neighbours.
//
// It is not cosmetic. A chain straight off the grid is a staircase, so its tangent alternates between two
// axis-aligned directions from vertex to vertex - and since the normal is the tangent's perpendicular and
// every classification is a dot product with the normal, an unsmoothed margin flickers between convergent
// and transform along its whole length. Two passes of a three-tap average cost a fraction of a grid cell in
// position and give a tangent that means something.
const smoothPasses = 2

// buildBoundaries finds every stretch of contact between two plates and says what each one is doing.
func (m *Model) buildBoundaries() []Boundary {
	groups := m.collect()

	// Sorted by pair, so the set is in the same order on every run: a planet's tectonics must not depend on
	// Go's map iteration order, or two runs of the same seed would write different meta.json files.
	keys := make([][2]int, 0, len(groups))
	for k := range groups {
		keys = append(keys, k)
	}
	sort.Slice(keys, func(i, j int) bool {
		if keys[i][0] != keys[j][0] {
			return keys[i][0] < keys[j][0]
		}
		return keys[i][1] < keys[j][1]
	})

	circ := m.P.CircumferenceM()
	maxGap := maxGapCells * m.GCellM

	var out []Boundary
	for _, k := range keys {
		for _, chain := range chainSamples(groups[k], circ, maxGap) {
			b := m.classify(k[0], k[1], chain, circ)
			if len(b.V) >= minChainSamples {
				out = append(out, b)
			}
		}
	}
	return out
}

// collect walks the tectonic grid once and records every cell edge whose two sides belong to different
// plates.
//
// East and south only. Testing all four neighbours would record each edge twice, and a chain built from
// duplicated points walks on the spot.
func (m *Model) collect() map[[2]int][]sample {
	half := m.GCellM / 2
	out := make(map[[2]int][]sample)
	add := func(a, b int, xM, yM, dx, dy float64) {
		if a == b {
			return
		}
		key := [2]int{a, b}
		if a > b {
			key = [2]int{b, a}
			dx, dy = -dx, -dy
		}
		out[key] = append(out[key], sample{xM: xM, yM: yM, dx: dx, dy: dy, a: key[0], b: key[1]})
	}

	for gy := 0; gy < m.GH; gy++ {
		row := gy * m.GW
		for gx := 0; gx < m.GW; gx++ {
			here := int(m.Cell[row+gx])
			east := int(m.Cell[m.GridIdx(gx+1, gy)])
			add(here, east, m.GridXM(gx)+half, m.GridYM(gy), 1, 0)
			if gy+1 < m.GH {
				south := int(m.Cell[m.GridIdx(gx, gy+1)])
				add(here, south, m.GridXM(gx), m.GridYM(gy)+half, 0, 1)
			}
		}
	}
	return out
}

// chainSamples orders a pair's scattered crossings into one or more polylines.
//
// A greedy nearest-unused walk rather than a proper contour tracer. The set it is given is one cell wide by
// construction, so the nearest unused neighbour is the next point along the line in every case except a
// triple junction, where the walk takes one branch and the other becomes a chain of its own - which is the
// right answer, because two plates meeting a third meet it on two different margins.
//
// O(n squared) on purpose. n is a few hundred, because the tectonic grid is a quarter of a kilometre and a
// boundary is a few tens of kilometres; a spatial index here would be more code than the thing it indexes.
func chainSamples(ss []sample, circ, maxGap float64) [][]sample {
	used := make([]bool, len(ss))
	var out [][]sample
	for {
		seed := pickEnd(ss, used, circ, maxGap)
		if seed < 0 {
			break
		}
		used[seed] = true
		fwd := walk(ss, used, seed, circ, maxGap)
		back := walk(ss, used, seed, circ, maxGap)

		chain := make([]sample, 0, len(fwd)+len(back)+1)
		for i := len(back) - 1; i >= 0; i-- {
			chain = append(chain, ss[back[i]])
		}
		chain = append(chain, ss[seed])
		for _, i := range fwd {
			chain = append(chain, ss[i])
		}
		if len(chain) >= minChainSamples {
			out = append(out, chain)
		}
	}
	return out
}

// pickEnd chooses where to start a chain: a sample with at most one unused neighbour, which is an end of the
// line. Starting in the middle would give two half-chains walked in opposite directions and joined at a
// point, which is the same line with a kink in the tangent at its centre.
func pickEnd(ss []sample, used []bool, circ, maxGap float64) int {
	best, bestDeg := -1, 1<<30
	for i := range ss {
		if used[i] {
			continue
		}
		deg := 0
		for j := range ss {
			if i == j || used[j] {
				continue
			}
			if dist(ss[i], ss[j], circ) <= maxGap {
				deg++
			}
		}
		if deg <= 1 {
			return i
		}
		if deg < bestDeg {
			best, bestDeg = i, deg
		}
	}
	return best
}

// walk steps from a sample to its nearest unused neighbour until there is none in reach.
func walk(ss []sample, used []bool, from int, circ, maxGap float64) []int {
	var out []int
	cur := from
	for {
		best, bestD := -1, maxGap
		for j := range ss {
			if used[j] {
				continue
			}
			if d := dist(ss[cur], ss[j], circ); d <= bestD {
				best, bestD = j, d
			}
		}
		if best < 0 {
			return out
		}
		used[best] = true
		out = append(out, best)
		cur = best
	}
}

func dist(a, b sample, circ float64) float64 {
	return math.Hypot(wrapDelta(a.xM-b.xM, circ), a.yM-b.yM)
}

// classify turns a chain of crossings into a boundary: unwrapped, smoothed, and with the relative motion
// resolved into a closing rate and a slip rate at every vertex.
func (m *Model) classify(a, b int, chain []sample, circ float64) Boundary {
	xs := make([]float64, len(chain))
	ys := make([]float64, len(chain))
	xs[0], ys[0] = chain[0].xM, chain[0].yM
	// Unwrap as the chain is copied: each point is put within half a circumference of the one before it, so
	// a margin crossing the seam comes out as a straight run of increasing X rather than a jump.
	for i := 1; i < len(chain); i++ {
		xs[i] = xs[i-1] + wrapDelta(chain[i].xM-xs[i-1], circ)
		ys[i] = chain[i].yM
	}
	smooth(xs, ys)

	obliqueRad := m.Cfg.ObliqueDeg * math.Pi / 180
	over := m.overriding(a, b)

	out := Boundary{A: a, B: b, V: make([]Vertex, len(chain))}
	for i := range chain {
		tx, ty := tangent(xs, ys, i)
		// The normal is the tangent's perpendicular, and the crossing itself says which of the two
		// perpendiculars points into plate B.
		nx, ny := -ty, tx
		if nx*chain[i].dx+ny*chain[i].dy < 0 {
			nx, ny = ty, -tx
		}

		vax, vay := m.Plates[a].VelocityAt(m.P, xs[i], ys[i])
		vbx, vby := m.Plates[b].VelocityAt(m.P, xs[i], ys[i])
		rx, ry := vax-vbx, vay-vby

		closing := rx*nx + ry*ny
		slip := rx*tx + ry*ty

		out.V[i] = Vertex{
			XM: xs[i], YM: ys[i], NX: nx, NY: ny,
			ClosingMYr: closing, SlipMYr: slip,
			Kind: kindOf(closing, slip, obliqueRad,
				m.Plates[a].Continental && m.Plates[b].Continental),
			Over: -1,
		}
		if out.V[i].Kind == Subduction {
			out.V[i].Over = over
		}
	}
	return out
}

// kindOf is the classification itself, and it is one comparison: is the relative motion more across the line
// or more along it, and if across, which way.
func kindOf(closing, slip, obliqueRad float64, bothContinental bool) Kind {
	if math.Atan2(math.Abs(slip), math.Abs(closing)) > obliqueRad {
		return Transform
	}
	if closing > 0 {
		if bothContinental {
			return Collision
		}
		return Subduction
	}
	if bothContinental {
		return Rift
	}
	return Ridge
}

// overriding is which of two plates ends up on top when they converge.
//
// The continental one, when exactly one is: continental crust is too buoyant to go down, which is why the
// Andes are on South America and not on the Nazca plate. When both sides are oceanic it is the larger, as a
// stand-in for the older and therefore colder and denser slab being the one that sinks.
func (m *Model) overriding(a, b int) int {
	ca, cb := m.Plates[a].Continental, m.Plates[b].Continental
	switch {
	case ca && !cb:
		return a
	case cb && !ca:
		return b
	case m.Plates[a].AreaCells >= m.Plates[b].AreaCells:
		return a
	default:
		return b
	}
}

// tangent is the local direction of the line, as a unit vector, from a central difference.
func tangent(xs, ys []float64, i int) (tx, ty float64) {
	lo, hi := i-1, i+1
	if lo < 0 {
		lo = 0
	}
	if hi >= len(xs) {
		hi = len(xs) - 1
	}
	tx, ty = xs[hi]-xs[lo], ys[hi]-ys[lo]
	if d := math.Hypot(tx, ty); d > 0 {
		return tx / d, ty / d
	}
	return 1, 0
}

// smooth averages the polyline with its own neighbours, in place, with the ends pinned. See smoothPasses for
// why an unsmoothed chain is unusable rather than merely ugly.
func smooth(xs, ys []float64) {
	if len(xs) < 3 {
		return
	}
	bx := make([]float64, len(xs))
	by := make([]float64, len(ys))
	for pass := 0; pass < smoothPasses; pass++ {
		copy(bx, xs)
		copy(by, ys)
		for i := 1; i < len(xs)-1; i++ {
			xs[i] = (bx[i-1] + 2*bx[i] + bx[i+1]) / 4
			ys[i] = (by[i-1] + 2*by[i] + by[i+1]) / 4
		}
	}
}