package field

import "math"

// SmoothEdgePreserving relaxes a height field towards its neighbours with a weight that falls away as the
// step between them grows, so a channel wall or a ridge crest survives a pass that takes a grid-cut facet
// off. It is a port of the bilateral smooth in the World Orogen browser generator, which has one for exactly
// this reason - to blend the artefacts its own routing leaves without rounding the landforms off with them.
//
// It is a filter and not a process. It conserves nothing, it has no time in it, and running it inside the
// solve loop would act as an uncontrolled extra diffusivity: that changes the steady-state slope, which is
// U/K, which is the one knob the whole generator's relief hangs on. It runs once, after the solve, and it is
// off by default. The point of having it is that the alternative - raising diffusion_m2_yr until the
// artefacts go - is measured to smooth away the landforms too, at about 0.05.
//
// Two deviations from the reference, both about units.
//
// The weight is 1/(1 + |dh|/(d*slopeRef)) rather than 1/(1 + |dh|*sensitivity). A sensitivity in 1/m is a
// height threshold, and a height threshold means one thing on a 32 m geology cell and something four times
// as aggressive on an 8 m one, so the same painted world would come out differently at two resolutions -
// which is the property Docs/Terrain-Next.md section 4.D says the generator lives or dies by. slopeRef is a
// rise over run and carries across. Ground steeper than it is preserved; ground gentler is relaxed.
//
// And a diagonal neighbour is sqrt(2) further away, so it carries both its own distance in the slope and an
// inverse-distance geometric weight - which is what a Gaussian would give those two offsets.
//
// The waterline is a wall, not a value. A neighbour that is not land is skipped entirely rather than clamped:
// clamping to sea level would pull the shore down, and clamping the other way would drown the beach the
// coastal pass built. Sea cells are never written.
//
// scratch must be at least len(h); it is used as the destination of each pass.
func SmoothEdgePreserving(h []float32, w, hgt int, cellM float64, land []bool, passes int, slopeRef float64, scratch []float32) {
	if passes <= 0 || slopeRef <= 0 || cellM <= 0 {
		return
	}
	if passes > smoothMaxPasses {
		passes = smoothMaxPasses
	}
	tmp := scratch[:len(h)]

	// dh/(d*slopeRef) per face, folded into one reciprocal each.
	invCard := float32(1 / (cellM * slopeRef))
	invDiag := float32(1 / (cellM * math.Sqrt2 * slopeRef))
	const geomDiag = float32(1 / math.Sqrt2)

	for p := 0; p < passes; p++ {
		src := h
		Rows(hgt, func(y0, y1 int) {
			for y := y0; y < y1; y++ {
				for x := 0; x < w; x++ {
					i := y*w + x
					if !land[i] {
						tmp[i] = src[i]
						continue
					}
					c := src[i]
					var sumW, sumH float32
					face := func(nx, ny int, inv, geom float32) {
						if nx < 0 || ny < 0 || nx >= w || ny >= hgt {
							return
						}
						ni := ny*w + nx
						if !land[ni] {
							return
						}
						n := src[ni]
						d := n - c
						if d < 0 {
							d = -d
						}
						wk := geom / (1 + d*inv)
						sumW += wk
						sumH += wk * n
					}
					face(x-1, y, invCard, 1)
					face(x+1, y, invCard, 1)
					face(x, y-1, invCard, 1)
					face(x, y+1, invCard, 1)
					face(x-1, y-1, invDiag, geomDiag)
					face(x+1, y-1, invDiag, geomDiag)
					face(x-1, y+1, invDiag, geomDiag)
					face(x+1, y+1, invDiag, geomDiag)
					tmp[i] = (c + sumH) / (1 + sumW)
				}
			}
		})
		copy(h, tmp)
	}
}

// smoothMaxPasses is a hard ceiling, not a default. Past about three passes the edge weight has stopped
// protecting anything - every face inside a landform is gentler than slopeRef by then - and what is left is a
// box blur with extra steps.
const smoothMaxPasses = 4