{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE NamedFieldPuns #-} {-# LANGUAGE TypeApplications #-} module Moonlight.Triangulation.Internal.CircleSweep ( circleSweepInsert ) where import Control.Monad (forM_, when) import Control.Monad.ST (ST) import Data.Bits (xor) import Data.STRef ( STRef , newSTRef , readSTRef , writeSTRef ) import qualified Data.Vector.Algorithms.Intro as Intro import qualified Data.Vector.Unboxed.Mutable as MUV import Data.Word (Word32) import Moonlight.Triangulation.Handles.HandleDefs (DirectedEdgeId (..)) import Moonlight.Triangulation.Insertion (insertExistingVertex) import Moonlight.Triangulation.Internal.DcelOperations ( closeOuterTurn , drainLegalization , fixHullConvexity , insertOutsideHullBetween , LegalizationLaw (..) , noStarVertex , seedGenericEdges ) import Moonlight.Triangulation.Internal.Mutable import Moonlight.Triangulation.Internal.OperationState ( Counter (..) , OperationState , addCounter , maxCounter , readScratch , writeScratch ) import Moonlight.Triangulation.Internal.Probe (Probe (..)) import Moonlight.Triangulation.Scalar (orient2dCoordinates) import Moonlight.Triangulation.Types (BuildError (..)) -- | The hull is an angular index over the DCEL outer-face cycle. The cycle -- itself already owns hull adjacency, so left and right walking is -- 'readPrevious'/'readNext' on the live mesh and this record only caches what -- topology cannot state: the pseudo-angle of each outer edge's origin and the -- bucket anchors accelerating the predecessor search. An outer edge's key is -- its origin's @(angle, x, y)@ with the edge id as the final tie-break; the -- angle lives in 'hullAngleByEdge' and the coordinates are re-read from the -- immutable origin only when two cached angles compare exactly equal. Slots -- of edges that have left the outer cycle are never read again, so the cache -- needs no invalidation, and there is no second ring beside the authoritative -- one. data Hull s = Hull { hullCenterX :: {-# UNPACK #-} !Double , hullCenterY :: {-# UNPACK #-} !Double , hullBucketCapacity :: {-# UNPACK #-} !Int , hullAngleByEdge :: !(MUV.MVector s Double) , hullActiveCount :: !(MUV.MVector s Int) , hullBuckets :: !(STRef s (MUV.MVector s Word32)) } readActiveCount :: Hull s -> ST s Int readActiveCount hull = MUV.unsafeRead (hullActiveCount hull) 0 {-# INLINE readActiveCount #-} writeActiveCount :: Hull s -> Int -> ST s () writeActiveCount hull = MUV.unsafeWrite (hullActiveCount hull) 0 {-# INLINE writeActiveCount #-} noOuterEdge :: Word32 noOuterEdge = maxBound -- | Circle sweep over one mutable DCEL, consuming one packed radial arena of -- @(squaredDistance, x, y, vertex)@ records. The arena is sorted in place and -- then read directly — no decorated freeze, no undecoration pass. Insertions -- initially close only acute hull turns. One terminal Graham pass restores -- full convexity, so construction does not repeatedly pay for global -- convexity that no intermediate observer can see. circleSweepInsert :: MutableDcel s vertex directed undirected face -> OperationState s -> MUV.MVector s (Double, Double, Double, Word32) -> ST s (Either BuildError Int) circleSweepInsert mutable operation arena | MUV.length arena == 0 = pure (Right 0) | otherwise = do Intro.sort arena let !ordered = MUV.length arena seed <- insertSeed 0 case seed of Left failure -> pure (Left failure) Right seedCount -> do faces <- faceCount mutable if faces <= 1 || seedCount >= ordered then pure (Right seedCount) else do (centerX, centerY) <- seedCentre mutable builtHull <- buildHull mutable operation centerX centerY case builtHull of Left failure -> pure (Left failure) Right hull -> do skipped <- MUV.new (ordered - seedCount) inserted <- insertRemaining hull ordered seedCount skipped 0 0 0 0 case inserted of Left failure -> pure (Left failure) Right (!skippedCount, !fastCount, !flips, !maxDepth) -> do -- The sweep counts its own hull insertions; its drains -- hand their flip and depth tallies up once. addCounter operation CounterHullInsertions fastCount -- The angular candidate leaves a point to the fallback -- whenever it lands right of, or on, the hull edge its -- own angle selected — spade's -- `is_on_right_side_or_on_line` branch, which spade's own -- source calls "very slow". Both counts ride out in -- BuildStats so the split is comparable across the two -- implementations and satisfies -- seed + fast + skipped = unique. -- -- CounterSweepFastPoints is NOT CounterHullInsertions -- renamed: the latter is also charged by -- insertOutsideHull, so it counts hull-adjacent -- insertions from either path and dominates this one -- whenever a skipped point lands outside the hull. The -- two coincide only while skipped is zero, and charging -- them independently is what makes the identity a check -- rather than a restatement. addCounter operation CounterSweepFastPoints fastCount addCounter operation CounterSweepSkippedPoints skippedCount addCounter operation CounterEdgeFlips flips maxCounter operation CounterLegalizationMaxStack maxDepth repaired <- fixHullConvexity @'ProbeOff mutable operation case repaired of Left failure -> pure (Left failure) Right (!_closures, !terminalFlips, !terminalMaxDepth) -> do addCounter operation CounterEdgeFlips terminalFlips maxCounter operation CounterLegalizationMaxStack terminalMaxDepth insertedSkipped <- insertSkipped skipped skippedCount 0 pure (seedCount <$ insertedSkipped) where insertSeed !index | index >= MUV.length arena = pure (Right index) | otherwise = do (_, _, _, raw) <- MUV.unsafeRead arena index result <- insertExistingVertex @'ProbeOff mutable operation (fromIntegral raw) case result of Left failure -> pure (Left failure) Right () -> do faces <- faceCount mutable if faces > 1 then pure (Right (index + 1)) else insertSeed (index + 1) insertRemaining !hull !ordered !index !skipped !skippedCount !fastCount !flips !maxDepth | index >= ordered = pure (Right (skippedCount, fastCount, flips, maxDepth)) | otherwise = do (_, queryXWide, queryYWide, raw) <- MUV.unsafeRead arena index let !vertex = fromIntegral raw queryX <- readPointX mutable vertex queryY <- readPointY mutable vertex let !queryAngle = pseudoAngle (hullCenterX hull) (hullCenterY hull) queryXWide queryYWide edge <- hullCandidate mutable hull queryAngle queryXWide queryYWide fromVertex <- readOrigin mutable edge toVertex <- readOrigin mutable (edge `xor` 1) fromX <- readPointX mutable fromVertex fromY <- readPointY mutable fromVertex toX <- readPointX mutable toVertex toY <- readPointY mutable toVertex if orient2dCoordinates fromX fromY toX toY queryX queryY == GT then do deferred <- insertDeferred mutable operation hull edge vertex queryAngle case deferred of Left failure -> pure (Left failure) Right (!_newClosures, !newFlips, !newMaxDepth) -> insertRemaining hull ordered (index + 1) skipped skippedCount (fastCount + 1) (flips + newFlips) (max maxDepth newMaxDepth) else do MUV.unsafeWrite skipped skippedCount raw insertRemaining hull ordered (index + 1) skipped (skippedCount + 1) fastCount flips maxDepth insertSkipped skipped !count = go where go !index | index >= count = pure (Right ()) | otherwise = do vertex <- fromIntegral <$> MUV.unsafeRead skipped index inserted <- insertExistingVertex @'ProbeOff mutable operation vertex case inserted of Left failure -> pure (Left failure) Right () -> go (index + 1) -- | The hull centre: the centroid of the first inner face, in the widened -- comparison format, computed exactly as @centroid@ states it. seedCentre :: MutableDcel s vertex directed undirected face -> ST s (Double, Double) seedCentre mutable = do (e0, e1, e2) <- faceEdges mutable 1 o0 <- readOrigin mutable e0 o1 <- readOrigin mutable e1 o2 <- readOrigin mutable e2 x0 <- readPointX mutable o0 y0 <- readPointY mutable o0 x1 <- readPointX mutable o1 y1 <- readPointY mutable o1 x2 <- readPointX mutable o2 y2 <- readPointY mutable o2 pure (x0 + (x1 - x0) / 3 + (x2 - x0) / 3, y0 + (y1 - y0) / 3 + (y2 - y0) / 3) -- | Index the authoritative outer cycle. The seed fan is star-shaped around -- the first face's centroid — a collinear chain closed by its apex — so the -- cycle is already the angular order the predecessor search assumes; what -- remains is caching each edge's angle and anchoring the buckets. buildHull :: MutableDcel s vertex directed undirected face -> OperationState s -> Double -> Double -> ST s (Either BuildError (Hull s)) buildHull mutable operation centerX centerY = do countResult <- collectOuterEdges mutable operation case countResult of Left failure -> pure (Left failure) Right count | count <= 0 -> pure (Left CircleSweepHullEmpty) | otherwise -> do let !capacity = max (count + 4) (pointCapacity mutable + 8) hullAngleByEdge <- MUV.new (halfEdgeCapacity mutable) hullActiveCount <- MUV.replicate 1 count initialBuckets <- MUV.replicate (initialBucketCount count capacity) noOuterEdge hullBuckets <- newSTRef initialBuckets let hull = Hull { hullCenterX = centerX , hullCenterY = centerY , hullBucketCapacity = capacity , hullAngleByEdge , hullActiveCount , hullBuckets } forM_ [0 .. count - 1] $ \index -> do edge <- readScratch operation index origin <- readOrigin mutable edge x <- readPointX mutable origin y <- readPointY mutable origin MUV.unsafeWrite hullAngleByEdge edge (pseudoAngle centerX centerY x y) installBucketMaximum mutable hull edge pure (Right hull) collectOuterEdges :: MutableDcel s vertex directed undirected face -> OperationState s -> ST s (Either BuildError Int) collectOuterEdges mutable operation = do start <- readFaceEdge mutable 0 if start < 0 then pure (Right 0) else do bound <- directedEdgeCount mutable go (bound + 1) start start False 0 where go !remaining !start !edge !seen !count | remaining <= 0 = pure ( Left ( OuterCycleDidNotTerminate (DirectedEdgeId (fromIntegral start)) (DirectedEdgeId (fromIntegral edge)) count ) ) | seen && edge == start = pure (Right count) | otherwise = do writeScratch operation count edge following <- readNext mutable edge go (remaining - 1) start following True (count + 1) initialBucketCount :: Int -> Int -> Int initialBucketCount active capacity = min capacity (nextPowerOfTwo (max 8 ((active + 7) `quot` 8))) nextPowerOfTwo :: Int -> Int nextPowerOfTwo requested = go 1 where target = max 1 requested go !value | value >= target = value | value > maxBound `quot` 2 = maxBound | otherwise = go (value * 2) readAngle :: Hull s -> Int -> ST s Double readAngle hull edge = MUV.unsafeRead (hullAngleByEdge hull) edge {-# INLINE readAngle #-} rebuildBuckets :: MutableDcel s vertex directed undirected face -> Hull s -> Int -> ST s () rebuildBuckets mutable hull requested = do let !count = max 1 (min (hullBucketCapacity hull) requested) buckets <- MUV.replicate count noOuterEdge writeSTRef (hullBuckets hull) buckets active <- readActiveCount hull start <- readFaceEdge mutable 0 let go !remaining !edge | remaining <= 0 = pure () | otherwise = do installBucketMaximum mutable hull edge following <- readNext mutable edge go (remaining - 1) following when (active > 0 && start >= 0) (go active start) maybeGrowBuckets :: MutableDcel s vertex directed undirected face -> Hull s -> ST s () maybeGrowBuckets mutable hull = do active <- readActiveCount hull buckets <- readSTRef (hullBuckets hull) let !current = MUV.length buckets when (active > 8 * current && current < hullBucketCapacity hull) $ rebuildBuckets mutable hull (min (hullBucketCapacity hull) (2 * current)) bucketFor :: Int -> Double -> Int bucketFor count angle = min (count - 1) (max 0 (floor (angle * fromIntegral count * 0.25))) {-# INLINE bucketFor #-} -- | Compare two outer edges by their origins' keys: cached angle first, then -- the immutable origin coordinates, then the edge id. The coordinate reads -- only happen on an exact angle tie, which is the same-ray case and no other. compareEdgeKeys :: MutableDcel s vertex directed undirected face -> Hull s -> Int -> Int -> ST s Ordering compareEdgeKeys mutable hull left right = do leftAngle <- readAngle hull left rightAngle <- readAngle hull right case compare leftAngle rightAngle of LT -> pure LT GT -> pure GT EQ -> do leftOrigin <- readOrigin mutable left rightOrigin <- readOrigin mutable right leftX <- readPointX mutable leftOrigin rightX <- readPointX mutable rightOrigin case compare leftX rightX of LT -> pure LT GT -> pure GT EQ -> do leftY <- readPointY mutable leftOrigin rightY <- readPointY mutable rightOrigin case compare leftY rightY of LT -> pure LT GT -> pure GT EQ -> pure (compare left right) -- | Whether an outer edge's key orders at or before the stated query key, -- settled field by field without materializing either key. edgeAtMost :: MutableDcel s vertex directed undirected face -> Hull s -> Int -> Double -> Double -> Double -> Int -> ST s Bool edgeAtMost mutable hull edge queryAngle queryX queryY tie = do angle <- readAngle hull edge case compare angle queryAngle of LT -> pure True GT -> pure False EQ -> do origin <- readOrigin mutable edge x <- readPointX mutable origin case compare x queryX of LT -> pure True GT -> pure False EQ -> do y <- readPointY mutable origin case compare y queryY of LT -> pure True GT -> pure False EQ -> pure (edge <= tie) installBucketMaximum :: MutableDcel s vertex directed undirected face -> Hull s -> Int -> ST s () installBucketMaximum mutable hull edge = do buckets <- readSTRef (hullBuckets hull) angle <- readAngle hull edge let !bucket = bucketFor (MUV.length buckets) angle currentRaw <- MUV.unsafeRead buckets bucket if currentRaw == noOuterEdge then MUV.unsafeWrite buckets bucket (fromIntegral edge) else do verdict <- compareEdgeKeys mutable hull edge (fromIntegral currentRaw) when (verdict == GT) (MUV.unsafeWrite buckets bucket (fromIntegral edge)) -- | Drop @edge@ from the bucket index. When it was its bucket's maximum, the -- cyclic predecessor @left@ stands in provided it still belongs to the same -- bucket — exact, because a bucket's edges form one contiguous arc of the -- angular order. The caller states @left@ explicitly: the mesh has already -- been mutated by the time the index is updated, so the retired edge no -- longer knows its own neighbour. removeBucketMaximum :: Hull s -> Int -> Int -> ST s () removeBucketMaximum hull edge left = do buckets <- readSTRef (hullBuckets hull) angle <- readAngle hull edge let !bucket = bucketFor (MUV.length buckets) angle current <- MUV.unsafeRead buckets bucket when (current == fromIntegral edge) $ do leftAngle <- readAngle hull left if left /= edge && bucketFor (MUV.length buckets) leftAngle == bucket then MUV.unsafeWrite buckets bucket (fromIntegral left) else MUV.unsafeWrite buckets bucket noOuterEdge retireEdge :: Hull s -> Int -> Int -> ST s () retireEdge hull edge left = do removeBucketMaximum hull edge left active <- readActiveCount hull writeActiveCount hull (active - 1) {-# INLINE retireEdge #-} activateEdge :: MutableDcel s vertex directed undirected face -> Hull s -> Int -> ST s () activateEdge mutable hull edge = do active <- readActiveCount hull writeActiveCount hull (active + 1) installBucketMaximum mutable hull edge maybeGrowBuckets mutable hull -- | The outer edge whose key is the greatest key at or below the query: the -- visible candidate the sweep inserts against. Bucket anchors land the walk -- near the answer and the live outer cycle carries it the rest of the way. hullCandidate :: MutableDcel s vertex directed undirected face -> Hull s -> Double -> Double -> Double -> ST s Int hullCandidate mutable hull queryAngle queryX queryY = do buckets <- readSTRef (hullBuckets hull) let !count = MUV.length buckets !bucket = bucketFor count queryAngle raw <- MUV.unsafeRead buckets bucket if raw == noOuterEdge then previousNonEmpty buckets (if bucket == 0 then count - 1 else bucket - 1) count else adjustWithinBucket buckets bucket (fromIntegral raw) where previousNonEmpty buckets !bucket !remaining | remaining <= 0 = readFaceEdge mutable 0 | otherwise = do raw <- MUV.unsafeRead buckets bucket if raw /= noOuterEdge then pure (fromIntegral raw) else previousNonEmpty buckets (if bucket == 0 then MUV.length buckets - 1 else bucket - 1) (remaining - 1) adjustWithinBucket buckets !bucket !initial = do active <- readActiveCount hull go active initial where go !remaining !edge | remaining <= 0 = pure initial | otherwise = do atMost <- edgeAtMost mutable hull edge queryAngle queryX queryY maxBound if atMost then pure edge else do left <- readPrevious mutable edge leftAngle <- readAngle hull left if bucketFor (MUV.length buckets) leftAngle /= bucket then pure left else go (remaining - 1) left insertDeferred :: MutableDcel s vertex directed undirected face -> OperationState s -> Hull s -> Int -> Int -> Double -> ST s (Either BuildError (Int, Int, Int)) insertDeferred mutable operation hull replacedEdge vertex insertedAngle = do inserted <- insertOutsideHullBetween @'ProbeOff mutable operation replacedEdge replacedEdge vertex case inserted of Left failure -> pure (Left failure) Right (firstEdge, lastEdge) -> insertDeferredBetween firstEdge lastEdge where insertDeferredBetween firstEdge lastEdge = do -- The patch proves the new keys algebraically, so no key is rediscovered -- from the mesh: replacing outer edge a->b with a->v, v->b keeps the -- replaced edge's origin for the first spoke (its cached angle carries -- over) and starts the last spoke at the inserted vertex, whose angle is the -- one the candidate search was already given. The tie-break is each fresh -- edge id itself. replacedAngle <- readAngle hull replacedEdge leftOfReplaced <- readPrevious mutable firstEdge retireEdge hull replacedEdge leftOfReplaced -- Both spokes join the outer cycle before either is indexed, and a bucket -- rebuild during activation walks the live cycle: both angle slots must be -- initialized before the first activation can read them. MUV.unsafeWrite (hullAngleByEdge hull) firstEdge replacedAngle MUV.unsafeWrite (hullAngleByEdge hull) lastEdge insertedAngle activateEdge mutable hull firstEdge activateEdge mutable hull lastEdge -- Every turn this insertion closes legalizes together in one epoch. Closure -- never deletes an edge and never touches the outer cycle, so which turns -- close does not depend on when the interior is repaired; the fan above -- keeps its own drain because its candidates are oriented against the -- inserted vertex, and that star must still be intact when they are tested. leftOutcome <- closeLeft firstEdge 0 0 case leftOutcome of Left obstruction -> pure (Left obstruction) Right (_left, closuresLeft, topLeft) -> do rightOutcome <- closeRight lastEdge 0 topLeft case rightOutcome of Left obstruction -> pure (Left obstruction) Right (_right, closuresRight, topAll) -> do (flips, maxDepth) <- drainLegalization @'ProbeOff mutable operation topAll noStarVertex ValidMesh pure (Right (closuresLeft + closuresRight, flips, maxDepth)) closeLeft !current !closures !top = do left <- readPrevious mutable current close <- shouldCloseTurn mutable hull insertedAngle left current if not close then pure (Right (current, closures, top)) else do leftAngle <- readAngle hull left closed <- closeOuterTurn mutable left case closed of Left obstruction -> pure (Left obstruction) Right replacement -> do -- Closing consecutive a->b, b->c into a->c keeps the first edge's -- origin, so the replacement inherits its angle; the tie is the fresh -- edge id. Both retired edges answer to the same left neighbour, the -- edge now preceding the replacement on the outer cycle. Seeding -- happens here, after the replacement's links exist. leftOfReplacement <- readPrevious mutable replacement retireEdge hull left leftOfReplacement retireEdge hull current leftOfReplacement MUV.unsafeWrite (hullAngleByEdge hull) replacement leftAngle activateEdge mutable hull replacement nextTop <- seedGenericEdges operation top [left, current] closeLeft replacement (closures + 1) nextTop closeRight !current !closures !top = do right <- readNext mutable current close <- shouldCloseTurn mutable hull insertedAngle current right if not close then pure (Right (current, closures, top)) else do currentAngle <- readAngle hull current closed <- closeOuterTurn mutable current case closed of Left obstruction -> pure (Left obstruction) Right replacement -> do leftOfReplacement <- readPrevious mutable replacement retireEdge hull current leftOfReplacement retireEdge hull right leftOfReplacement MUV.unsafeWrite (hullAngleByEdge hull) replacement currentAngle activateEdge mutable hull replacement nextTop <- seedGenericEdges operation top [current, right] closeRight replacement (closures + 1) nextTop {-# INLINE insertDeferred #-} shouldCloseTurn :: MutableDcel s vertex directed undirected face -> Hull s -> Double -> Int -> Int -> ST s Bool shouldCloseTurn mutable hull insertedAngle first second = do following <- readNext mutable first if following /= second then pure False else do fromVertex <- readOrigin mutable first middleVertex <- readOrigin mutable (first `xor` 1) targetVertex <- readOrigin mutable (second `xor` 1) fromX <- readPointX mutable fromVertex fromY <- readPointY mutable fromVertex middleX <- readPointX mutable middleVertex middleY <- readPointY mutable middleVertex targetX <- readPointX mutable targetVertex targetY <- readPointY mutable targetVertex if orient2dWide fromX fromY middleX middleY targetX targetY /= GT then pure False else do -- The second edge begins at the middle vertex and stands on the outer -- cycle, so its cached key already is that vertex's pseudo-angle; -- same-ray is a slot read rather than a second angle. middleAngle <- readAngle hull second pure ( middleAngle == insertedAngle || acuteAtMiddle fromX fromY middleX middleY targetX targetY ) {-# INLINE shouldCloseTurn #-} -- The deferred turn test deliberately runs the widened Binary64 predicate -- rather than the exact binary64 one; only the terminal Graham pass owns exact -- convexity. This is the class's Double instance called by name, which is -- what the boxed 'orient2d' on widened points resolved to. orient2dWide :: Double -> Double -> Double -> Double -> Double -> Double -> Ordering orient2dWide = orient2dCoordinates {-# INLINE orient2dWide #-} -- Spade's deferred-convexity rule is local: close the turn when the angle at -- the shared hull vertex is strictly below 90 degrees. Requiring the entire -- triangle to be acute leaves avoidable star-hull work for the terminal pass. acuteAtMiddle :: Double -> Double -> Double -> Double -> Double -> Double -> Bool acuteAtMiddle ax ay bx by cx cy = let !ux = ax - bx !uy = ay - by !vx = cx - bx !vy = cy - by !dot = ux * vx + uy * vy !scale = max 1 (ux * ux + uy * uy + vx * vx + vy * vy) in dot > 64 * encodeFloat 1 (-52) * scale {-# INLINE acuteAtMiddle #-} -- Clockwise pseudo-angle in [0,4), matching the orientation of the outer-face -- cycle. pseudoAngle :: Double -> Double -> Double -> Double -> Double pseudoAngle centerX centerY x y | norm == 0 = 0 | raw >= 4 = 0 | otherwise = raw where !dx = x - centerX !dy = y - centerY !norm = abs dx + abs dy !projection = dx / norm !raw = if dy > 0 then 1 + projection else 3 - projection {-# INLINE pseudoAngle #-}