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

package field

// SlidingMax is the maximum over a square window, separable and O(1) a cell whatever the radius.
//
// The naive form is a loop over the window, which is what internal/stats' localRelief used to do and is fine
// on the 500 m window it uses at 8 m cells - until the map is a planet. 28 million land cells times a 63-cell
// radius is 1.1e11 comparisons, which is not a slow diagnostic, it is one nobody will ever see the end of.
// The monotonic deque is the standard answer: each index enters and leaves once, so the row pass is linear in
// the row however wide the window.
//
// wrapX makes the row pass periodic, which is what a cylinder needs; the column pass always clamps, because
// the top and bottom of the map are the poles and not each other.
func SlidingMax(f *Field, radius int, wrapX bool) *Field {
	return sliding(f, radius, wrapX, func(inDeque, arriving float32) bool { return inDeque <= arriving })
}

// SlidingMin is the same window, the other way up. The pair is what local relief is made of.
func SlidingMin(f *Field, radius int, wrapX bool) *Field {
	return sliding(f, radius, wrapX, func(inDeque, arriving float32) bool { return inDeque >= arriving })
}

// LocalRelief is max minus min over a square window: the standard field measure of how rugged a place is, and
// the one thing slope cannot tell you. A 5 m hummock and a 500 m mountainside both stand at 30 degrees.
//
// Two sliding passes and a subtract, so it costs the same as one of them twice and nothing per radius. It
// holds two fields at once at the peak, which at planet scale is 600 MB - worth saying, because the naive
// version held none and could not finish.
func LocalRelief(f *Field, radius int, wrapX bool) *Field {
	hi := SlidingMax(f, radius, wrapX)
	lo := SlidingMin(f, radius, wrapX)
	for i := range hi.Data {
		hi.Data[i] -= lo.Data[i]
	}
	return hi
}

// sliding is the shared separable pass. keep reports whether the value already at the back of the deque can
// be dropped when a new one arrives, which is the only thing that differs between the maximum and the
// minimum: the deque holds indices whose values are monotone, so its front is always the answer for the live
// window and anything the arriving value dominates can never be the answer again.
func sliding(f *Field, radius int, wrapX bool, keep func(inDeque, arriving float32) bool) *Field {
	if radius < 1 {
		return f.Clone()
	}
	w, h := f.W, f.H
	row := New(w, h, f.CellM)
	buf := make([]float32, 0, w+2*radius)
	idx := make([]int, 0, w+2*radius)

	for y := 0; y < h; y++ {
		// The row, extended by the radius at each end so the deque never has to special-case an edge.
		buf = buf[:0]
		for x := -radius; x < w+radius; x++ {
			sx := x
			if wrapX {
				sx = ((sx % w) + w) % w
			} else if sx < 0 {
				sx = 0
			} else if sx >= w {
				sx = w - 1
			}
			buf = append(buf, f.Data[y*w+sx])
		}
		slide(buf, idx[:0], 2*radius+1, keep, func(i int, v float32) {
			if i < w {
				row.Data[y*w+i] = v
			}
		})
	}

	out := New(w, h, f.CellM)
	col := make([]float32, 0, h+2*radius)
	for x := 0; x < w; x++ {
		col = col[:0]
		for y := -radius; y < h+radius; y++ {
			sy := y
			if sy < 0 {
				sy = 0
			} else if sy >= h {
				sy = h - 1
			}
			col = append(col, row.Data[sy*w+x])
		}
		slide(col, idx[:0], 2*radius+1, keep, func(i int, v float32) {
			if i < h {
				out.Data[i*w+x] = v
			}
		})
	}
	return out
}

// slide walks a padded line with a monotonic deque and reports the window's answer ending at each output
// position.
func slide(line []float32, dq []int, window int, keep func(inDeque, arriving float32) bool,
	emit func(i int, v float32)) {

	dq = dq[:0]
	for i, v := range line {
		for len(dq) > 0 && keep(line[dq[len(dq)-1]], v) {
			dq = dq[:len(dq)-1]
		}
		dq = append(dq, i)
		if dq[0] <= i-window {
			dq = dq[1:]
		}
		if out := i - window + 1; out >= 0 {
			emit(out, line[dq[0]])
		}
	}
}

// BoxSmooth blurs a field in place with `passes` of a separable box blur of the given radius, clamping at the
// edges. Two passes are near enough to a Gaussian for anything here and cost four linear sweeps.
//
// Deterministic by construction: fixed traversal order, running sums, no goroutines. It lives here rather than
// in the pass that first wanted it because two now do - the coastal detail pass smooths the signed distance to
// the shoreline, and the tile bake smooths the interpolated sea floor.
func BoxSmooth(data []float32, w, h, radius, passes int) {
	if radius < 1 || passes < 1 || len(data) < w*h {
		return
	}
	tmp := make([]float32, len(data))
	for p := 0; p < passes; p++ {
		boxRows(data, tmp, w, h, radius)
		boxCols(tmp, data, w, h, radius)
	}
}

func boxRows(src, dst []float32, w, h, radius int) {
	n := float32(2*radius + 1)
	for y := 0; y < h; y++ {
		row := y * w
		var sum float32
		for k := -radius; k <= radius; k++ {
			sum += src[row+clampIdx(k, w)]
		}
		for x := 0; x < w; x++ {
			dst[row+x] = sum / n
			sum += src[row+clampIdx(x+radius+1, w)] - src[row+clampIdx(x-radius, w)]
		}
	}
}

func boxCols(src, dst []float32, w, h, radius int) {
	n := float32(2*radius + 1)
	for x := 0; x < w; x++ {
		var sum float32
		for k := -radius; k <= radius; k++ {
			sum += src[clampIdx(k, h)*w+x]
		}
		for y := 0; y < h; y++ {
			dst[y*w+x] = sum / n
			sum += src[clampIdx(y+radius+1, h)*w+x] - src[clampIdx(y-radius, h)*w+x]
		}
	}
}

func clampIdx(i, n int) int {
	if i < 0 {
		return 0
	}
	if i >= n {
		return n - 1
	}
	return i
}

// BoxSmoothMasked is BoxSmooth restricted to the cells the mask selects: a cell outside it is neither read
// nor written, so the blur never averages across the boundary.
//
// That distinction is the whole reason it exists. The coastal detail pass damps the metre-scale texture near
// the shore, and an unmasked blur there does not damp texture, it bridges the waterline: measured on a
// fixture with forty metres of water against the land, the plain blur lifted the sea floor by twenty metres.
// The step at a shoreline is a landform, not roughness, and a filter that cannot tell them apart is the wrong
// filter.
//
// Separable and weighted: the row pass carries a running sum of values and of weights, the column pass sums
// those, and the quotient is the mean over the masked cells in the window. Deterministic, like BoxSmooth.
func BoxSmoothMasked(data []float32, mask []bool, w, h, radius, passes int) {
	if radius < 1 || passes < 1 || len(data) < w*h || len(mask) < w*h {
		return
	}
	n := w * h
	val := make([]float32, n)
	wgt := make([]float32, n)
	tv := make([]float32, n)
	tw := make([]float32, n)
	for p := 0; p < passes; p++ {
		for i := 0; i < n; i++ {
			if mask[i] {
				val[i], wgt[i] = data[i], 1
			} else {
				val[i], wgt[i] = 0, 0
			}
		}
		boxRowsSum(val, tv, w, h, radius)
		boxRowsSum(wgt, tw, w, h, radius)
		boxColsSum(tv, val, w, h, radius)
		boxColsSum(tw, wgt, w, h, radius)
		for i := 0; i < n; i++ {
			if mask[i] && wgt[i] > 0 {
				data[i] = val[i] / wgt[i]
			}
		}
	}
}

// boxRowsSum and boxColsSum are the running sums BoxSmooth uses, without the division: a masked blur needs
// the weight sum as well as the value sum, and dividing in the middle would be dividing by the wrong thing.
func boxRowsSum(src, dst []float32, w, h, radius int) {
	for y := 0; y < h; y++ {
		row := y * w
		var sum float32
		for k := -radius; k <= radius; k++ {
			sum += src[row+clampIdx(k, w)]
		}
		for x := 0; x < w; x++ {
			dst[row+x] = sum
			sum += src[row+clampIdx(x+radius+1, w)] - src[row+clampIdx(x-radius, w)]
		}
	}
}

func boxColsSum(src, dst []float32, w, h, radius int) {
	for x := 0; x < w; x++ {
		var sum float32
		for k := -radius; k <= radius; k++ {
			sum += src[clampIdx(k, h)*w+x]
		}
		for y := 0; y < h; y++ {
			dst[y*w+x] = sum
			sum += src[clampIdx(y+radius+1, h)*w+x] - src[clampIdx(y-radius, h)*w+x]
		}
	}
}