114 lines
4.6 KiB
Go
114 lines
4.6 KiB
Go
package fluvial
|
|
|
|
// A monotone bucket priority queue, which is what priority-flood actually needs.
|
|
//
|
|
// The flood pops cells in non-decreasing elevation and never pushes anything below the cell it just popped:
|
|
// a neighbour lower than the current front is raised to it and goes to the FIFO instead. That "monotone"
|
|
// property is exactly the condition under which a bucket queue beats a binary heap, because the read cursor
|
|
// only ever moves forward and both operations become an append and a scan. The heap was costing about
|
|
// log2(3.2M) = 22 comparisons and as many cache misses per operation, on two thirds of the solve's runtime.
|
|
//
|
|
// Elevations are quantised into fixed-width buckets. Cells inside one bucket pop in an arbitrary but
|
|
// deterministic order (last in, first out), so a spill point can be wrong by at most one bucket width. At a
|
|
// centimetre against a 2 km elevation range that is far below the millimetre-per-cell epsilon the flood adds
|
|
// anyway, and it is the same approximation an integer-elevation priority-flood makes by construction.
|
|
type bucketPQ struct {
|
|
lo float64
|
|
width float64
|
|
buckets [][]int32
|
|
cur int
|
|
count int
|
|
}
|
|
|
|
const bucketWidthM = 0.01
|
|
|
|
// reposeOrderBuckets is how far ClampToRepose scatters a cell from its own bucket, in buckets, either way.
|
|
// Sixteen buckets is sixteen centimetres: enough that a tie spreads over about thirty of them and the pop
|
|
// order stops tracking the raster, far below the metres the clamp's ordering could ever depend on. See
|
|
// pushJittered.
|
|
const reposeOrderBuckets = 16
|
|
|
|
func newBucketPQ(loM, hiM float64) *bucketPQ {
|
|
if hiM <= loM {
|
|
hiM = loM + 1
|
|
}
|
|
// Headroom above the top: the flood raises cells by epsilon as it fills, so the highest key pushed can
|
|
// sit slightly above the terrain's own maximum.
|
|
n := int((hiM-loM)/bucketWidthM) + 64
|
|
return &bucketPQ{lo: loM, width: bucketWidthM, buckets: make([][]int32, n)}
|
|
}
|
|
|
|
func (q *bucketPQ) reset() {
|
|
for i := range q.buckets {
|
|
q.buckets[i] = q.buckets[i][:0]
|
|
}
|
|
q.cur = 0
|
|
q.count = 0
|
|
}
|
|
|
|
func (q *bucketPQ) len() int { return q.count }
|
|
|
|
func (q *bucketPQ) push(elev float32, idx int32) {
|
|
b := int((float64(elev) - q.lo) / q.width)
|
|
if b < q.cur {
|
|
b = q.cur // monotone: never behind the cursor, whatever rounding says
|
|
}
|
|
if b >= len(q.buckets) {
|
|
b = len(q.buckets) - 1
|
|
}
|
|
q.buckets[b] = append(q.buckets[b], idx)
|
|
q.count++
|
|
}
|
|
|
|
// pushJittered is push with the bucket chosen from the elevation plus j buckets, so cells at the same
|
|
// elevation land in different buckets instead of in one and pop in hash order rather than in reverse raster
|
|
// order.
|
|
//
|
|
// Only ClampToRepose wants this. The flood already scatters its own order through the epsilon it adds to the
|
|
// height, so its cells rarely share a bucket; the clamp reads the raw surface, where flat ground puts every
|
|
// cell in one bucket and the LIFO below then processes them bottom-right to top-left, every time, everywhere.
|
|
//
|
|
// Half a bucket is not enough - it splits a tie across two buckets and halves the correlation instead of
|
|
// removing it - so the caller scatters over reposeOrderBuckets, and that is safe for a reason worth writing
|
|
// down. The clamp's order can only matter between two cells whose heights differ by about the talus
|
|
// allowance, which is metres: a cell popped early is marked closed and never lowered again, but every cell
|
|
// popping after it stands within the jitter width of it, so its limit is the other cell's height plus the
|
|
// full allowance and cannot bind. Reordering cells that are centimetres apart therefore cannot break a
|
|
// constraint that only bites metres apart. The bound is talus*cell/2, which is 2.8 m at 35 degrees on an 8 m
|
|
// cell; the constant below is two orders of magnitude inside it.
|
|
func (q *bucketPQ) pushJittered(elev float32, idx int32, j float32) {
|
|
b := int((float64(elev)-q.lo)/q.width + float64(j))
|
|
if b < q.cur {
|
|
b = q.cur
|
|
}
|
|
if b >= len(q.buckets) {
|
|
b = len(q.buckets) - 1
|
|
}
|
|
q.buckets[b] = append(q.buckets[b], idx)
|
|
q.count++
|
|
}
|
|
|
|
// pop returns the lowest cell. The cursor only moves forward, so the total scan cost over a whole flood is
|
|
// the number of buckets, not the number of pops.
|
|
func (q *bucketPQ) pop() int32 {
|
|
for q.cur < len(q.buckets) && len(q.buckets[q.cur]) == 0 {
|
|
q.cur++
|
|
}
|
|
if q.cur >= len(q.buckets) {
|
|
return -1
|
|
}
|
|
b := q.buckets[q.cur]
|
|
v := b[len(b)-1]
|
|
q.buckets[q.cur] = b[:len(b)-1]
|
|
q.count--
|
|
return v
|
|
}
|
|
|
|
// frontElev is the elevation the cursor is at, which the FIFO compares itself against.
|
|
func (q *bucketPQ) frontElev() float32 {
|
|
for q.cur < len(q.buckets) && len(q.buckets[q.cur]) == 0 {
|
|
q.cur++
|
|
}
|
|
return float32(q.lo + float64(q.cur)*q.width)
|
|
}
|