// Package dt is the exact Euclidean distance transform, with a feature index and an optional cylinder.
//
// It lives on its own because three different things need it and two of them are nowhere near the coast:
// the coastal pass writes every one of its processes as "how far is this cell from the waterline and which
// stretch of shore does it belong to"; the region partitioner dilates the land mask to decide which
// landmasses are close enough to be solved together; and the template classifier dissolves the decorative
// stroke an artist drew by handing each of its pixels to the nearest pixel that means something.
//
// Exact, not a chamfer approximation: Felzenszwalb and Huttenlocher's transform is two 1-D passes and O(n)
// whatever the radius, so there is nothing to buy by approximating, and a chamfer's 2 % anisotropy would
// show up directly as a shelf wider along the grid axes than across them.
package dt

import (
	"math"

	"salty/terrain/internal/field"
)

// Transform returns, for every cell, the squared distance in cells to the nearest seed cell and the flat
// index of that seed. A column pass finds the nearest seed in each column; a row pass takes the lower
// envelope of the parabolas those distances define.
//
// With wrapX the row pass is periodic, so the left and right edges of the grid are neighbours. That is what
// a planet needs: a landmass straddling the seam is one landmass, and the shelf in front of it is one shelf.
//
// Cells in a column with no seed at all are given a cost above any real distance rather than an infinity, so
// the envelope arithmetic never sees a NaN; they are then never chosen unless the grid has no seeds
// anywhere, in which case every near index comes back -1.
func Transform(seed []bool, w, h int, wrapX bool) (d2 []float32, near []int32) {
	return transform(seed, w, h, wrapX, true)
}

// Distance2 is Transform without the feature index, for a caller that only wants "how far".
//
// It is a separate entry point rather than a nil argument because the saving is the point: at planet scale
// the index and the column scratch it needs are two more arrays of four bytes a cell, which is most of a
// gigabyte for an answer nobody reads. The region partitioner only asks whether a cell is within a margin
// of land.
func Distance2(seed []bool, w, h int, wrapX bool) []float32 {
	d2, _ := transform(seed, w, h, wrapX, false)
	return d2
}

func transform(seed []bool, w, h int, wrapX, wantNear bool) (d2 []float32, near []int32) {
	d2 = make([]float32, w*h)
	if wantNear {
		near = make([]int32, w*h)
	}

	bigF := float64(w*w+h*h) * 4 // above any achievable dx² + dy²
	bigD := float32(math.Sqrt(bigF))

	colD := make([]float32, w*h) // distance in cells to the nearest seed in this column
	var colN []int32             // that seed's row, or -1; only needed for the feature index
	if wantNear {
		colN = make([]int32, w*h)
	}

	field.Rows(w, func(x0, x1 int) {
		for x := x0; x < x1; x++ {
			best := -1
			for y := 0; y < h; y++ {
				i := y*w + x
				if seed[i] {
					best = y
				}
				if best < 0 {
					colD[i] = bigD
					if wantNear {
						colN[i] = -1
					}
				} else {
					colD[i] = float32(y - best)
					if wantNear {
						colN[i] = int32(best)
					}
				}
			}
			best = -1
			for y := h - 1; y >= 0; y-- {
				i := y*w + x
				if seed[i] {
					best = y
				}
				if best >= 0 {
					if d := float32(best - y); d < colD[i] {
						colD[i] = d
						if wantNear {
							colN[i] = int32(best)
						}
					}
				}
			}
		}
	})

	// The row pass. On a cylinder the row is laid out three times - one turn to the left, the row itself,
	// one turn to the right - and the answer is read out of the middle copy. From a cell in the middle copy
	// the three images of any column sit at offsets d, d-w and d+w, whose smallest absolute value is the
	// cyclic distance, so the envelope returns exactly the wrapped answer with no special cases in it.
	span := w
	off := 0
	if wrapX {
		span = 3 * w
		off = w
	}
	field.Rows(h, func(y0, y1 int) {
		f := make([]float64, span)
		v := make([]int, span)
		z := make([]float64, span+1)
		for y := y0; y < y1; y++ {
			row := y * w
			for j := 0; j < span; j++ {
				d := float64(colD[row+srcX(j, off, w)])
				f[j] = d * d
			}
			k := 0
			v[0] = 0
			z[0] = math.Inf(-1)
			z[1] = math.Inf(1)
			for q := 1; q < span; q++ {
				s := intersect(f, v[k], q)
				for s <= z[k] {
					k--
					s = intersect(f, v[k], q)
				}
				k++
				v[k] = q
				z[k] = s
				z[k+1] = math.Inf(1)
			}
			k = 0
			for q := 0; q < span; q++ {
				for z[k+1] < float64(q) {
					k++
				}
				if q < off || q >= off+w {
					continue // a replica column; only the middle copy is the answer
				}
				dx := float64(q - v[k])
				o := row + q - off
				d2[o] = float32(dx*dx + f[v[k]])
				if !wantNear {
					continue
				}
				sx := srcX(v[k], off, w)
				if n := colN[row+sx]; n < 0 {
					near[o] = -1
				} else {
					near[o] = n*int32(w) + int32(sx)
				}
			}
		}
	})
	return d2, near
}

// srcX maps a column of the (possibly replicated) row back to a real column.
func srcX(j, off, w int) int {
	x := j - off
	for x < 0 {
		x += w
	}
	for x >= w {
		x -= w
	}
	return x
}

// intersect is where the parabolas rooted at p and q cross.
func intersect(f []float64, p, q int) float64 {
	return ((f[q] + float64(q*q)) - (f[p] + float64(p*p))) / float64(2*q-2*p)
}