UnrealPrototyping/Tools/Terrain/internal/field/preview.go

package field

import (
	"bufio"
	"image"
	"image/color"
	"image/png"
	"math"
	"os"
	"path/filepath"
	"sort"
)

// A colour preview of a height field: hypsometric tint, hillshade, and the drainage network drawn on top.
//
// The grey thumbnail is nearly useless for judging this generator, which is a problem, because the thing it
// exists to produce is a drainage network and a flat grey ramp is exactly what hides one. Rivers are drawn
// from the flow accumulation with a width that grows with drainage area, so a glance says whether the network
// branches like a river system or like noise.

type PreviewOptions struct {
	// Flow is drainage area per cell, m². Optional; without it no rivers are drawn.
	Flow *Field
	// Sea marks cells below sea level. Optional.
	Sea       []bool
	SeaLevelM float64
	// Snow marks land that is permanently under ice. Optional, and it exists because the hypsometric ramp
	// tops out at snow by *elevation*: a polar cap fifty metres above the water therefore comes out the same
	// green as a meadow, and an ice sheet that reads as a meadow is a map lying about the one thing it is
	// for. No height is touched; only the colour.
	Snow []bool
	// RiverKm2 is the drainage area at which a channel starts being drawn.
	RiverKm2 float64
	// Size is the output width in pixels. The height follows the field's own aspect, so a 2:1 planet comes
	// out 2:1 rather than squashed into a square; on the square canvas the two are the same number and
	// nothing changes.
	Size int
	// Crop is a sub-rectangle in map coordinates (x0, y0, x1, y1 in 0..1), rendered at full resolution.
	// A whole continent at 1500 px puts ten kilometres into a hundred pixels, which is enough to see that
	// there is drainage and not nearly enough to see whether it is the right *kind* of drainage. Judging
	// hill country against real hill country needs a crop.
	Crop [4]float64
	// Palette is how the picture is drawn: the ramp, the water, the rivers, the ice and the light. Nil is
	// the generator's own, which is what every caller wanted before this was a file.
	Palette *Palette
	// Hillshade exaggerates the vertical before shading. Lowland relief is a few tens of metres over
	// kilometres and disappears at true scale, which is the same reason every printed relief map lies.
	Exaggeration float64
}

// rgb is the palette's colour type under the name the drawing code uses.
type rgb = RGB

// WritePreview renders the field at opt.Size and writes an RGB PNG.
//
// It returns the height the hypsometric ramp topped out at, in metres, which a caller is expected to print.
// The ramp is relative by default and a relative picture is only honest when the reader is told so: without
// that line, a 47 m plain drawn with snow on its hills is indistinguishable from an alpine one.
func WritePreview(path string, h *Field, opt PreviewOptions) (topM float64, err error) {
	size := opt.Size
	if size <= 0 {
		size = 1024
	}
	if size > h.W {
		size = h.W
	}
	if opt.Crop[2] > opt.Crop[0] && opt.Crop[3] > opt.Crop[1] {
		fullW, fullH := h.W, h.H
		h = h.Sub(opt.Crop)
		if opt.Flow != nil {
			opt.Flow = opt.Flow.Sub(opt.Crop)
		}
		if opt.Sea != nil {
			opt.Sea = subMask(opt.Sea, fullW, fullH, opt.Crop)
		}
		if size > h.W {
			size = h.W
		}
	}
	sizeH := aspectH(h, size)
	small := h.Resample(size, sizeH)
	exag := opt.Exaggeration
	if exag <= 0 {
		exag = 1
	}

	// Land elevations only: letting the sea floor into the range squashes the whole land ramp.
	//
	// And the top of the ramp is a high percentile, not the maximum. One 2800 m summit over a continent whose
	// land is mostly under 300 m puts every other cell into the bottom tenth of the ramp, and the map reads as
	// uniform green with a white dot on it — which says far more about one pixel than about the terrain. The
	// percentile lets the tint span the distribution that is actually there; the few cells above it clamp to
	// snow, which is what they should look like anyway.
	pal := opt.Palette
	if pal == nil {
		pal = DefaultPalette()
	}
	sea := resampleMask(opt.Sea, h.W, h.H, size, sizeH)
	snow := resampleMask(opt.Snow, h.W, h.H, size, sizeH)
	landVals := make([]float64, 0, len(small.Data))
	for i, v := range small.Data {
		if sea != nil && sea[i] {
			continue
		}
		landVals = append(landVals, float64(v))
	}
	landMax := 1.0
	if pal.LandTopM > 0 {
		landMax = pal.LandTopM
	} else if len(landVals) > 0 {
		sort.Float64s(landVals)
		landMax = landVals[int(pal.LandTopPercentile/100*float64(len(landVals)-1))]
	}
	if landMax <= 0 {
		landMax = 1
	}
	var seaMin float64
	for i, v := range small.Data {
		if sea != nil && sea[i] && float64(v) < seaMin {
			seaMin = float64(v)
		}
	}

	var flow *Field
	riverA := opt.RiverKm2 * 1e6
	if opt.Flow != nil && riverA > 0 {
		flow = opt.Flow.Resample(size, sizeH)
	}

	img := image.NewRGBA(image.Rect(0, 0, size, sizeH))
	for y := 0; y < sizeH; y++ {
		for x := 0; x < size; x++ {
			i := y*size + x
			elev := float64(small.Data[i])
			var c rgb

			if sea != nil && sea[i] {
				d := 0.0
				if seaMin < 0 {
					d = math.Min(1, (opt.SeaLevelM-elev)/(opt.SeaLevelM-seaMin))
				}
				c = rgb{
					pal.SeaShallow[0] + (pal.SeaDeep[0]-pal.SeaShallow[0])*d,
					pal.SeaShallow[1] + (pal.SeaDeep[1]-pal.SeaShallow[1])*d,
					pal.SeaShallow[2] + (pal.SeaDeep[2]-pal.SeaShallow[2])*d,
				}
			} else {
				c = pal.ramp(math.Min(1, math.Max(0, elev)/landMax))
				if snow != nil && snow[i] {
					// Ice, whatever height it stands at. It still takes the hillshade below rather than
					// being stamped flat, so a dome and the valleys cut into it still read.
					c = pal.Ice
				}
				// Hillshade from the north-west at 45 degrees, the DEM convention. Applied to land only;
				// shading the sea floor would draw attention to bathymetry nobody will ever see.
				gx := float64(small.AtClamped(x+1, y)-small.AtClamped(x-1, y)) * exag
				gy := float64(small.AtClamped(x, y+1)-small.AtClamped(x, y-1)) * exag
				slope := math.Atan(math.Hypot(gx, gy) / (2 * small.CellM))
				aspect := math.Atan2(gy, -gx)
				alt := pal.SunAltitudeDeg * math.Pi / 180
				// Azimuth is clockwise from north; the shading wants the direction the light comes *from*
				// measured the way Atan2 returns it, which is this quarter turn away.
				az := (90 - pal.SunAzimuthDeg) * math.Pi / 180
				lum := math.Cos(slope)*math.Sin(alt) +
					math.Sin(slope)*math.Cos(alt)*math.Cos(az-aspect)
				lum = pal.Ambient + pal.Gain*math.Max(0, lum)
				for k := range c {
					c[k] *= lum
				}
			}

			// Rivers on top, their strength growing with the log of drainage area so a trunk reads darker
			// than a headwater without needing a width in pixels.
			if flow != nil {
				if a := float64(flow.Data[i]); a >= riverA {
					w := math.Min(1, math.Log10(a/riverA)/2.2)
					blend := 0.45 + 0.55*w
					for k := range c {
						c[k] = c[k]*(1-blend) + pal.River[k]*blend
					}
				}
			}

			img.Set(x, y, color.RGBA{clamp8(c[0]), clamp8(c[1]), clamp8(c[2]), 255})
		}
	}

	if err := os.MkdirAll(filepath.Dir(path), 0o755); err != nil {
		return landMax, err
	}
	f, err := os.Create(path)
	if err != nil {
		return landMax, err
	}
	defer f.Close()
	bw := bufio.NewWriterSize(f, 1<<20)
	enc := png.Encoder{CompressionLevel: png.DefaultCompression}
	if err := enc.Encode(bw, img); err != nil {
		return landMax, err
	}
	return landMax, bw.Flush()
}

// resampleMask takes a boolean mask down to the preview size by nearest neighbour; a mask has no meaningful
// average.
func resampleMask(mask []bool, w, h, sw, sh int) []bool {
	if mask == nil {
		return nil
	}
	out := make([]bool, sw*sh)
	for y := 0; y < sh; y++ {
		sy := 0
		if sh > 1 {
			sy = y * (h - 1) / (sh - 1)
		}
		for x := 0; x < sw; x++ {
			sx := 0
			if sw > 1 {
				sx = x * (w - 1) / (sw - 1)
			}
			out[y*sw+x] = mask[sy*w+sx]
		}
	}
	return out
}

// aspectH is the output height that keeps a field's shape. Every writer in this package uses it, so a
// rectangular world is never silently squashed into a square image.
func aspectH(f *Field, w int) int {
	h := int(float64(w)*float64(f.H)/float64(f.W) + 0.5)
	if h < 1 {
		h = 1
	}
	return h
}

func clamp8(v float64) uint8 {
	if v <= 0 {
		return 0
	}
	if v >= 255 {
		return 255
	}
	return uint8(v + 0.5)
}

// subMask is Sub for a boolean mask.
func subMask(mask []bool, w, h int, crop [4]float64) []bool {
	clamp := func(v float64) float64 { return math.Min(1, math.Max(0, v)) }
	x0 := int(clamp(crop[0]) * float64(w-1))
	y0 := int(clamp(crop[1]) * float64(h-1))
	x1 := int(clamp(crop[2]) * float64(w-1))
	y1 := int(clamp(crop[3]) * float64(h-1))
	if x1 <= x0 {
		x1 = x0 + 1
	}
	if y1 <= y0 {
		y1 = y0 + 1
	}
	cw, ch := x1-x0+1, y1-y0+1
	out := make([]bool, cw*ch)
	for y := 0; y < ch; y++ {
		copy(out[y*cw:(y+1)*cw], mask[(y0+y)*w+x0:(y0+y)*w+x0+cw])
	}
	return out
}