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

package field

import "math"

// The vertex convention, which every resample here obeys: a field of N samples a side spans N-1 quads, so
// sample i sits at parameter i/(N-1) and the four corners are fixed points of any resize. Getting this wrong
// shifts the whole map by half a cell per resize and the error compounds over a pipeline.

// Resample returns the field at a new resolution: block means when shrinking by an exact integer factor
// (which is what the geology grid wants, and what preserves mass), bilinear otherwise. Ported from
// heightmap_io.resample, which chose the same two paths for the same reasons.
func (f *Field) Resample(w, h int) *Field {
	if w == f.W && h == f.H {
		return f.Clone()
	}
	cell := f.CellM * float64(f.W-1) / float64(w-1)
	if w < f.W && (f.W-1)%(w-1) == 0 && (f.H-1)%(h-1) == 0 && (f.W-1)/(w-1) == (f.H-1)/(h-1) {
		return f.blockMean((f.W-1)/(w-1), w, h, cell)
	}
	return f.bilinear(w, h, cell)
}

// blockMean averages each factor x factor block of quads onto one output sample. The last row and column are
// half-blocks under the vertex convention, which is why the accumulation counts what it actually summed.
func (f *Field) blockMean(factor, w, h int, cell float64) *Field {
	out := New(w, h, cell)
	Rows(h, func(y0, y1 int) {
		for oy := y0; oy < y1; oy++ {
			for ox := 0; ox < w; ox++ {
				var sum float64
				var n int
				for dy := 0; dy < factor; dy++ {
					sy := oy*factor + dy - factor/2
					if sy < 0 || sy >= f.H {
						continue
					}
					for dx := 0; dx < factor; dx++ {
						sx := ox*factor + dx - factor/2
						if sx < 0 || sx >= f.W {
							continue
						}
						sum += float64(f.At(sx, sy))
						n++
					}
				}
				if n > 0 {
					out.Data[out.Idx(ox, oy)] = float32(sum / float64(n))
				}
			}
		}
	})
	return out
}

func (f *Field) bilinear(w, h int, cell float64) *Field {
	out := New(w, h, cell)
	sx := float64(f.W-1) / float64(w-1)
	sy := float64(f.H-1) / float64(h-1)
	Rows(h, func(y0, y1 int) {
		for oy := y0; oy < y1; oy++ {
			fy := float64(oy) * sy
			iy := int(fy)
			ty := float32(fy - float64(iy))
			for ox := 0; ox < w; ox++ {
				fx := float64(ox) * sx
				ix := int(fx)
				tx := float32(fx - float64(ix))
				a := f.AtClamped(ix, iy)
				b := f.AtClamped(ix+1, iy)
				c := f.AtClamped(ix, iy+1)
				d := f.AtClamped(ix+1, iy+1)
				top := a + (b-a)*tx
				bot := c + (d-c)*tx
				out.Data[out.Idx(ox, oy)] = top + (bot-top)*ty
			}
		}
	})
	return out
}

// UpsampleInt is the geology-to-detail step: an exact integer factor on the quad count, so 1786 at factor 4
// becomes (1786-1)*4+1 = 7141 with every source sample landing exactly on an output sample and no resample
// phase error at all. Catmull-Rom between them, which is the bicubic the spec asks for and does not overshoot
// into ringing the way a plain cubic does on a ridge.
func (f *Field) UpsampleInt(factor int) *Field {
	if factor <= 1 {
		return f.Clone()
	}
	w := (f.W-1)*factor + 1
	h := (f.H-1)*factor + 1
	out := New(w, h, f.CellM/float64(factor))
	inv := 1.0 / float64(factor)
	Rows(h, func(y0, y1 int) {
		for oy := y0; oy < y1; oy++ {
			sy := oy / factor
			ty := float64(oy%factor) * inv
			for ox := 0; ox < w; ox++ {
				sx := ox / factor
				tx := float64(ox%factor) * inv
				var col [4]float64
				for k := 0; k < 4; k++ {
					col[k] = catmullRom(
						float64(f.AtClamped(sx-1, sy-1+k)),
						float64(f.AtClamped(sx, sy-1+k)),
						float64(f.AtClamped(sx+1, sy-1+k)),
						float64(f.AtClamped(sx+2, sy-1+k)), tx)
				}
				out.Data[out.Idx(ox, oy)] = float32(catmullRom(col[0], col[1], col[2], col[3], ty))
			}
		}
	})
	return out
}

func catmullRom(p0, p1, p2, p3, t float64) float64 {
	t2 := t * t
	t3 := t2 * t
	return 0.5 * ((2 * p1) +
		(-p0+p2)*t +
		(2*p0-5*p1+4*p2-p3)*t2 +
		(-p0+3*p1-3*p2+p3)*t3)
}

// ToUnit squashes a field into [0, 1] against a percentile, optionally through log1p first: what the four
// derivative maps (flow, wear, deposit) need before they become 8-bit PNGs. Ported from
// heightmap_erosion.to_unit.
func (f *Field) ToUnit(percentile float64, logScale bool) *Field {
	out := NewLike(f)
	for i, v := range f.Data {
		x := float64(v)
		if x < 0 {
			x = 0
		}
		if logScale {
			x = math.Log1p(x)
		}
		out.Data[i] = float32(x)
	}
	top := float64(out.Percentile(percentile))
	if top < 1e-6 {
		top = 1e-6
	}
	for i, v := range out.Data {
		x := float64(v) / top
		if x > 1 {
			x = 1
		}
		out.Data[i] = float32(x)
	}
	return out
}

// Sub extracts a sub-rectangle given in map coordinates (x0, y0, x1, y1 in 0..1), at the source resolution.
// Used by the preview to look at a piece of the map closely, which is the only way to judge whether hill
// country reads as hill country rather than as small mountains.
func (f *Field) Sub(crop [4]float64) *Field {
	clamp := func(v float64) float64 {
		if v < 0 {
			return 0
		}
		if v > 1 {
			return 1
		}
		return v
	}
	x0 := int(clamp(crop[0]) * float64(f.W-1))
	y0 := int(clamp(crop[1]) * float64(f.H-1))
	x1 := int(clamp(crop[2]) * float64(f.W-1))
	y1 := int(clamp(crop[3]) * float64(f.H-1))
	if x1 <= x0 {
		x1 = x0 + 1
	}
	if y1 <= y0 {
		y1 = y0 + 1
	}
	w, h := x1-x0+1, y1-y0+1
	out := New(w, h, f.CellM)
	for y := 0; y < h; y++ {
		copy(out.Data[y*w:(y+1)*w], f.Data[(y0+y)*f.W+x0:(y0+y)*f.W+x0+w])
	}
	return out
}