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

package plates

import (
	"encoding/json"
	"fmt"
	"math"
	"os"
	"regexp"

	"salty/terrain/internal/world"
)

// The painted tectonic layer: a third painting beside the template and the overlay, where a colour is a
// plate and the legend says how that plate is moving.
//
// **You paint the cause, not the conclusion.** A colour does not say "there is a collision here" - it says
// "this piece of lithosphere is moving north-east at three centimetres a year", and where two of them meet,
// what happens is worked out from the two motions and the shape of the contact. That is the whole reason to
// paint plates rather than to paint boundary lines: a drawn line has to be told what it is, while a contact
// between two painted plates *becomes* a collision, a transform or a rift by itself, and changes character
// along its own length wherever it turns relative to the motion. The Alpide belt is a collision at the
// Himalaya and a strike-slip fault through Anatolia for exactly that reason, and no author should have to
// hand-annotate it.
//
// It also means the tracer needs no new code. Build's weighted Voronoi and this both produce the same thing -
// a plate id at every cell of the tectonic grid - and everything downstream reads that.
//
// **Registration is by extent, not by pixel.** The layer is stretched over the painted map's own rectangle,
// so it does not have to be the template's size. Paint plates at a quarter of it if you like: the tectonic
// grid is a few hundred metres a cell and a plate is tens of kilometres across, so detail below that is
// detail nothing will ever read. The polar pad has no painting under it and takes the nearest painted row,
// which is right - a plate does not stop at the top of the author's canvas.

// PaintLegend is what the colours on a tectonic layer mean.
type PaintLegend struct {
	// Comment is the legend's own note to whoever opens it next. Propose writes the conventions into it,
	// because "which way does heading 90 point" is the first thing an author needs and the last thing they
	// should have to find in a source file.
	Comment string `json:"_comment,omitempty"`

	// Image is the layer's file name, resolved beside the legend. The manifest may name one instead.
	Image string `json:"image"`

	// WarnDistance is how far, in RGB, a sampled pixel may sit from the nearest plate before the run says so.
	// It exists for the same reason the class template's does: a JPEG bleeds several units of each channel
	// across a painted edge, and a silent mismatch is a plate boundary in the wrong place.
	WarnDistance float64 `json:"warn_distance"`

	Plates []PaintPlate `json:"plates"`
}

// PaintPlate is one painted plate: a colour, and how that piece of lithosphere is moving.
type PaintPlate struct {
	Name string `json:"name"`

	// RGB is the colour on the layer. Every sampled pixel becomes the *nearest* plate in RGB, because on a
	// tectonic layer every pixel has to be some plate - the same rule the class template uses, and the
	// opposite of the overlay's, where most of the image is deliberately nothing.
	RGB [3]int `json:"rgb"`

	// SpeedCmYr and HeadingDeg are the plate's drift. The heading is a compass bearing over the map: 0 points
	// at the top of the image, 90 to the right, 180 to the bottom. Earth's plates run 1 to 10 cm/yr, and what
	// matters at a margin is the *difference* between two of these, so two plates both drifting east at 4 are
	// a boundary doing nothing at all.
	SpeedCmYr  float64 `json:"speed_cm_yr"`
	HeadingDeg float64 `json:"heading_deg"`

	// SpinDegMyr turns the plate about its own centre, in degrees per million years, positive clockwise on
	// the map.
	//
	// It is worth setting on at least one plate. A planet of plates that only drift has margins that are the
	// same all the way along, because the relative velocity is then one constant vector and the only thing
	// that varies is where the contact happens to point. A little spin is what makes one end of a margin
	// collide while the other slides - which is the Anatolia case, and the most useful thing a tectonic map
	// can give a fault set.
	SpinDegMyr float64 `json:"spin_deg_myr"`

	// Continental overrides what the painting says. Left out - which is the usual case - a plate is
	// continental when enough of its painted area is land, so the template decides and the two paintings
	// cannot contradict each other. Set it when they should: an oceanic plate carrying a chain of islands, or
	// a continental fragment currently underwater.
	Continental *bool `json:"continental,omitempty"`
}

// rgbOneLine finds an indented colour triple so MarshalLegend can put it back on one line.
var rgbOneLine = regexp.MustCompile(`"rgb": \[\s*(\d+),\s*(\d+),\s*(\d+)\s*\]`)

// MarshalLegend writes a legend as JSON somebody will want to edit.
//
// json.MarshalIndent puts every colour on five lines, because Indent reformats every array whatever a custom
// marshaller does, and a seven-plate legend then runs to ninety lines of mostly punctuation. Putting the
// triples back on one line each is cosmetic and it is worth the ten lines: this file is meant to be opened
// and changed by hand, beside the painting, and a legend nobody can read at a glance is a legend nobody
// keeps in step with the picture.
func MarshalLegend(lg *PaintLegend) ([]byte, error) {
	data, err := json.MarshalIndent(lg, "", "  ")
	if err != nil {
		return nil, err
	}
	return append(rgbOneLine.ReplaceAll(data, []byte(`"rgb": [$1, $2, $3]`)), '\n'), nil
}

// LoadPaintLegend reads a tectonic layer's legend.
func LoadPaintLegend(path string) (*PaintLegend, error) {
	data, err := os.ReadFile(path)
	if err != nil {
		return nil, err
	}
	var lg PaintLegend
	if err := json.Unmarshal(data, &lg); err != nil {
		return nil, fmt.Errorf("%s: %w", path, err)
	}
	if err := lg.validate(path); err != nil {
		return nil, err
	}
	return &lg, nil
}

func (l *PaintLegend) validate(path string) error {
	if len(l.Plates) < 2 {
		return fmt.Errorf("%s: %d plate(s); a planet in one plate has no boundaries", path, len(l.Plates))
	}
	if l.WarnDistance <= 0 {
		l.WarnDistance = 60
	}
	seen := map[[3]int]string{}
	for i := range l.Plates {
		p := &l.Plates[i]
		if p.Name == "" {
			return fmt.Errorf("%s: plate %d has no name", path, i)
		}
		for c := range 3 {
			if p.RGB[c] < 0 || p.RGB[c] > 255 {
				return fmt.Errorf("%s: plate %q has rgb %v", path, p.Name, p.RGB)
			}
		}
		if prev, dup := seen[p.RGB]; dup {
			return fmt.Errorf("%s: plates %q and %q are both rgb %v; a colour is one plate",
				path, prev, p.Name, p.RGB)
		}
		seen[p.RGB] = p.Name
		if p.SpeedCmYr < 0 {
			return fmt.Errorf("%s: plate %q moves at %v cm/yr; speed is a magnitude and the heading is "+
				"where it points", path, p.Name, p.SpeedCmYr)
		}
	}
	return nil
}

// PaintMatch is how well the painting matched the legend, reported the way the class template's match is: a
// layer whose colours have drifted is a tectonic model quietly built on the wrong plates.
type PaintMatch struct {
	Cells       int     `json:"cells"`
	Far         int     `json:"far"`
	MaxDistance float64 `json:"max_distance"`
}

// FromPainting builds a tectonic model from a painted layer instead of from a seed.
//
// px is the layer decoded to RGB triples, pw by ph. Decoding happens in the caller so that this package keeps
// knowing nothing about files or image formats - the same reason land is a callback.
func FromPainting(p world.Planet, cfg Config, lg *PaintLegend, px []uint8, pw, ph int,
	land func(xM, yM float64) bool) (*Model, PaintMatch, error) {

	var match PaintMatch
	if lg == nil || len(lg.Plates) < 2 {
		return nil, match, fmt.Errorf("a tectonic layer needs at least two plates")
	}
	if pw <= 0 || ph <= 0 || len(px) < pw*ph*3 {
		return nil, match, fmt.Errorf("the tectonic layer is %dx%d with %d bytes", pw, ph, len(px))
	}

	cfg = cfg.withDefaults()
	circ := p.CircumferenceM()
	gw := int(circ/cfg.ResolutionM + 0.5)
	if gw < 8 {
		gw = 8
	}
	gcell := circ / float64(gw)
	gh := int(float64(p.H)*p.CellM/gcell + 0.5)
	if gh < 2 {
		gh = 2
	}

	m := &Model{P: p, Cfg: cfg, GW: gw, GH: gh, GCellM: gcell, Cell: make([]int16, gw*gh)}
	m.Plates = make([]Plate, len(lg.Plates))
	for i := range m.Plates {
		m.Plates[i] = Plate{ID: i, Weight: 1}
	}

	heightM := p.HeightM()
	for gy := range gh {
		yM := m.GridYM(gy)
		// The painted map covers 0..heightM; the polar pad above and below it takes the nearest painted row.
		v := clamp01(yM / heightM)
		py := int(v * float64(ph-1))
		row := gy * gw
		for gx := range gw {
			xM := m.GridXM(gx)
			pxi := int(xM / circ * float64(pw))
			if pxi >= pw {
				pxi = pw - 1
			}
			o := (py*pw + pxi) * 3
			id, dist := nearestPlate(lg.Plates, px[o], px[o+1], px[o+2])
			m.Cell[row+gx] = int16(id)

			match.Cells++
			if dist > lg.WarnDistance {
				match.Far++
			}
			if dist > match.MaxDistance {
				match.MaxDistance = dist
			}
		}
	}

	for i := range m.Plates {
		pl := &m.Plates[i]
		src := lg.Plates[i]
		// Compass bearing over the map: 0 points at the top of the image, which is -Y, and 90 to the right.
		speed := src.SpeedCmYr / 100
		bearing := src.HeadingDeg * math.Pi / 180
		pl.TransXM = speed * math.Sin(bearing)
		pl.TransYM = -speed * math.Cos(bearing)
		// Positive spin is clockwise on the map: with Y running down the image, v = T + omega x r sends the
		// point east of the centre southwards. The centre itself comes from measure, below.
		pl.OmegaRadYr = src.SpinDegMyr * math.Pi / 180 / 1e6
	}

	// measure fills in the area, the land fraction and the centre of each plate - and the centre is the pole
	// every one of them turns about, so nothing has a usable velocity field until this has run.
	m.measure(land)
	for i := range m.Plates {
		// A painted plate has no Voronoi site. Its centre of area is the only position it has, and it is what
		// the map and the reports point at.
		m.Plates[i].SiteXM = m.Plates[i].CentroidXM
		m.Plates[i].SiteYM = m.Plates[i].CentroidYM
	}
	// The painting has the last word where it asks for one, after measure has read the template's land.
	for i := range m.Plates {
		if c := lg.Plates[i].Continental; c != nil {
			m.Plates[i].Continental = *c
		}
	}

	m.Boundaries = m.buildBoundaries()
	return m, match, nil
}

// nearestPlate is the legend entry closest to a colour, and how far away it was.
//
// Nearest rather than exact, and unlike the overlay there is no "no plate" answer: every pixel of a tectonic
// layer is some piece of lithosphere, so a colour that matches nothing is a painting mistake to report rather
// than a hole to leave. WarnDistance is what reports it.
func nearestPlate(ps []PaintPlate, r, g, b uint8) (id int, dist float64) {
	best, bestID := math.Inf(1), 0
	for i := range ps {
		dr := float64(int(r) - ps[i].RGB[0])
		dg := float64(int(g) - ps[i].RGB[1])
		db := float64(int(b) - ps[i].RGB[2])
		if d := dr*dr + dg*dg + db*db; d < best {
			best, bestID = d, i
		}
	}
	return bestID, math.Sqrt(best)
}

func clamp01(v float64) float64 {
	if v < 0 {
		return 0
	}
	if v > 1 {
		return 1
	}
	return v
}

// Propose turns a generated model into a painting and a legend to start from.
//
// An author should not face a blank canvas for this. Seven plates with plausible motions is a minute's work
// for the Voronoi and an afternoon's by hand, and what an author actually wants to do is move two of them and
// change a heading - which is editing, not authoring from nothing.
//
// The returned pixels are the layer at the given width, and the legend has one entry per plate carrying the
// motion the generator drew. Writing them out is the caller's job.
func (m *Model) Propose(width int) (px []uint8, w, h int, lg *PaintLegend) {
	if width < 64 {
		width = 64
	}
	h = int(float64(width) * m.P.HeightM() / m.P.CircumferenceM())
	if h < 1 {
		h = 1
	}
	w = width

	lg = &PaintLegend{
		Comment: "A painted tectonic layer: one colour per plate, and how that plate is moving. " +
			"heading_deg is a compass bearing over the map - 0 points at the top of the image, 90 to the " +
			"right, 180 to the bottom. speed_cm_yr is drift; what happens at a margin is the difference " +
			"between the two plates either side of it, so two plates drifting the same way are a boundary " +
			"doing nothing. spin_deg_myr turns a plate about its own centre, positive clockwise, and it is " +
			"worth setting on at least one: without it every margin is the same all the way along, and " +
			"with it one end collides while the other slides. Paint the plates, not the mountains - where " +
			"two of these meet, the collision, the belt and its faults are worked out from the motions. " +
			"Repaint the blobs freely; only the colours have to keep matching this file.",
		WarnDistance: 60,
		Plates:       make([]PaintPlate, len(m.Plates)),
	}
	colours := make([][3]uint8, len(m.Plates))
	for i := range m.Plates {
		pl := &m.Plates[i]
		// Hues walked by the golden ratio, so that neighbouring ids are not neighbouring colours and an
		// author can tell two touching plates apart at a glance.
		c := hsvBytes(math.Mod(float64(i)*0.61803398875, 1)*360, 0.62, 0.86)
		colours[i] = c

		speed := math.Hypot(pl.TransXM, pl.TransYM) * 100 // m/yr to cm/yr
		// Back to a compass bearing: 0 at the top of the image, 90 to the right.
		bearing := math.Atan2(pl.TransXM, -pl.TransYM) * 180 / math.Pi
		if bearing < 0 {
			bearing += 360
		}
		lg.Plates[i] = PaintPlate{
			Name:       fmt.Sprintf("plate_%d", i),
			RGB:        [3]int{int(c[0]), int(c[1]), int(c[2])},
			SpeedCmYr:  math.Round(speed*10) / 10,
			HeadingDeg: math.Round(bearing),
			SpinDegMyr: math.Round(pl.OmegaRadYr*180/math.Pi*1e6*100) / 100,
		}
	}

	px = make([]uint8, w*h*3)
	for y := range h {
		yM := m.P.HeightM() * (float64(y) + 0.5) / float64(h)
		for x := range w {
			xM := m.P.CircumferenceM() * (float64(x) + 0.5) / float64(w)
			c := colours[m.PlateAt(xM, yM)]
			o := (y*w + x) * 3
			px[o], px[o+1], px[o+2] = c[0], c[1], c[2]
		}
	}
	return px, w, h, lg
}

// hsvBytes is a hue in degrees, saturation and value in 0..1, as an RGB triple.
func hsvBytes(hue, sat, val float64) [3]uint8 {
	hue = math.Mod(math.Mod(hue, 360)+360, 360) / 60
	i := math.Floor(hue)
	f := hue - i
	p := val * (1 - sat)
	q := val * (1 - sat*f)
	t := val * (1 - sat*(1-f))
	var r, g, b float64
	switch int(i) % 6 {
	case 0:
		r, g, b = val, t, p
	case 1:
		r, g, b = q, val, p
	case 2:
		r, g, b = p, val, t
	case 3:
		r, g, b = p, q, val
	case 4:
		r, g, b = t, p, val
	default:
		r, g, b = val, p, q
	}
	return [3]uint8{byte(r*255 + 0.5), byte(g*255 + 0.5), byte(b*255 + 0.5)}
}