package main
import (
"math"
"sync"
"github.com/paulmach/orb"
"github.com/paulmach/orb/quadtree"
)
// Using MaxUint32 as scaling factor
const mu32 = float64(math.MaxUint32)
// idFromCoord creates a unique 64-bit identifier from 2D coordinates
// The high 32 bits store x, low 32 bits store y
func idFromCoord(x, y float64) uint64 {
x0, y0 := uint64(x*mu32), uint64(y*mu32)
return (x0 << 32) | (y0 & 0xFFFFFFFF)
}
// coordFromID reverses the encoding process to extract original coordinates
func coordFromID(id uint64) (float64, float64) {
x, y := float64(id>>32), float64(id&0xFFFFFFFF)
return x / mu32, y / mu32
}
// point represents a site in the Voronoi diagram with spatial position and the
// (conservative) radius to re-compute if the point is deleted. orb.Point doesn't
// support adding metadata, so we must wrap orb.Point.
type point struct {
p orb.Point
d float64
}
// Point implements the orb.Pointer interface
func (p point) Point() orb.Point {
return p.p
}
// newPoint creates a new point with given coordinates and region radius
func newPoint(x, y, d float64) point {
return point{p: orb.Point{x, y}, d: d}
}
// VoronoiGraph maintains a dynamic set of points and supports on-the-fly Voronoi cell computation
type VoronoiGraph struct {
tree *quadtree.Quadtree
size int
mut sync.Mutex
}
// NewVoronoiGraph creates a new Voronoi graph within the specified spatial bounds
func NewVoronoiGraph(bounds orb.Bound) *VoronoiGraph {
return &VoronoiGraph{quadtree.New(bounds), 0, sync.Mutex{}}
}
// Add inserts a new point into the Voronoi graph and returns its id, the radius 'd' is used to aid in deletions
func (vg *VoronoiGraph) Add(x, y, d float64) uint64 {
id := idFromCoord(x, y)
pt := newPoint(x, y, d)
// quadtree isn't goroutine safe, operations must be atomic
vg.mut.Lock()
defer vg.mut.Unlock()
vg.tree.Add(pt)
vg.size++
return id
}
// Delete removes a point from the Voronoi graph and optionally recomputes affected cells
// The recompute flag trades performance for accuracy in dynamic scenarios
func (vg *VoronoiGraph) Delete(id uint64, recompute bool) bool {
x, y := coordFromID(id)
ptPtr := vg.tree.Find(orb.Point{x, y})
vg.mut.Lock()
defer vg.mut.Unlock()
p := ptPtr.Point()
if x == p.X() && y == p.Y() {
vg.size--
vg.tree.Remove(ptPtr, func(_ orb.Pointer) bool {
return true
})
}
// Recompute affected neighbors when deleting. This is quite expensive but
// is the only way to maintain a full cover of the vg.tree bounds.
if recompute {
pz, _ := ptPtr.(point)
d := pz.d
bb := orb.Bound{
orb.Point{pz.p[0] - d, pz.p[1] - d},
orb.Point{pz.p[0] + d, pz.p[1] + d},
}
for _, neighbor := range vg.tree.InBound(nil, bb) {
nn := neighbor.Point()
id := idFromCoord(nn[0], nn[1])
_ = vg.ApproxVoronoiCell(id, 2.0, false, nil)
}
}
return true
}
// ApproxVoronoiCell computes an approximate Voronoi cell using k-nearest neighbors
// Trades exact computation for speed by only considering nearby points
func (vg *VoronoiGraph) ApproxVoronoiCell(id uint64, conf float64, update bool, arr []orb.Pointer) orb.Ring {
// Start with unit square, progressively clipped by perpendicular bisectors, the
// number of neighbors scales logarithmically w. size
x, y := coordFromID(id)
ptPtr := vg.tree.Find(orb.Point{x, y})
pt := ptPtr.Point()
cell := orb.Ring{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}}
n := int(conf * math.Ceil(math.Log(float64(vg.size+1))))
for _, nn := range vg.tree.KNearest(arr, pt, n) {
neighbor := nn.Point()
// Calculate perpendicular bisector between pt and neighbor
mx, my := (pt[0]+neighbor[0])/2, (pt[1]+neighbor[1])/2
dx, dy := neighbor[0]-pt[0], neighbor[1]-pt[1]
px, py := -1e6*dy, 1e6*dx
// Require normal vector points toward half-space containing query point
if (pt[0]-mx)*dx+(pt[1]-my)*dy < 0 {
px, py = -px, -py
}
t := -py*(mx-px) + px*(my-py)
cell = clipPolygon(cell, -py, px, t)
}
// After computing the voronoi region, set the point's site radius.
// This (AFAIK) must be an add and delete :(
if update {
bb := cell.Bound()
l, h := bb.Min, bb.Max
dx := l[0] - h[0]
dy := l[1] - h[1]
d := math.Max(dx, dy)
vg.Delete(id, false)
vg.Add(x, y, 0.5*d)
}
return cell
}
// clipPolygon clips convex polygon against half-space using Sutherland-Hodgman algorithm
func clipPolygon(ring orb.Ring, a, b, t float64) orb.Ring {
if len(ring) < 2 {
return nil
}
// Fast path: check if all vertices are inside the half-space and return
// the current ring. Avoiding extra allocations from a new `output` is
// a big chunk of the perf gains...
n := len(ring) - 1
for i := 0; i < n; i++ {
if a*ring[i][0]+b*ring[i][1] > t {
goto needsClipping
}
}
return ring
needsClipping: // Sutherland-Hodgman clipping: four cases based on vertex positions
output := make(orb.Ring, 0, n+2)
for i := 0; i < n; i++ {
curr, next := ring[i], ring[(i+1)%n]
currInside := a*curr[0]+b*curr[1] <= t
nextInside := a*next[0]+b*next[1] <= t
if currInside {
if nextInside {
output = append(output, next)
} else {
inter := intersectLineSegment(curr, next, a, b, t)
output = append(output, inter)
}
} else if nextInside {
inter := intersectLineSegment(curr, next, a, b, t)
output = append(output, inter, next)
}
}
if len(output) > 0 && output[0] != output[len(output)-1] {
output = append(output, output[0])
}
return output
}
// intersectLineSegment finds where a line segment crosses a half-space boundary
func intersectLineSegment(p1, p2 orb.Point, a, b, c float64) orb.Point {
dx, dy := p2[0]-p1[0], p2[1]-p1[1]
denom := a*dx + b*dy
if denom == 0 {
return p1 // degenerate parallel case
}
t := (c - a*p1[0] - b*p1[1]) / denom
return orb.Point{p1[0] + t*dx, p1[1] + t*dy}
}
// ParallelVoronoiProcessor allows parallel computation of voronoi regions
type ParallelVoronoiProcessor struct {
workerCount int
vg *VoronoiGraph
}
func NewParallelVoronoiProcessor(nSites int, vg *VoronoiGraph, workerCount int) *ParallelVoronoiProcessor {
return &ParallelVoronoiProcessor{
workerCount: workerCount,
vg: vg,
}
}
func (pvp *ParallelVoronoiProcessor) ProcessParallel(jobs chan uint64, wg sync.WaitGroup) {
for i := 0; i < pvp.workerCount; i++ {
go pvp.worker(&wg, jobs)
}
}
// worker processes jobs from the shared channel with a pre-allocated array
func (pvp *ParallelVoronoiProcessor) worker(wg *sync.WaitGroup, jobs <-chan uint64) {
wg.Add(1)
defer wg.Done()
n := int(3.0 * math.Ceil(math.Log(float64(vg.size + 1))))
arr := make([]orb.Pointer, n)
for j := range jobs {
_ = pvp.vg.ApproxVoronoiCell(j, 3.0, true, arr)
}
}