moonlight-triangulation-0.1.0.0: src-dcel/Moonlight/Triangulation/Internal/DcelOperations/Normalize.hs
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
-- | The normalization procedure of the flip rewrite system.
module Moonlight.Triangulation.Internal.DcelOperations.Normalize
( drainLegalization
) where
import Control.Monad (when)
import Control.Monad.ST (ST)
import Data.Bits ((.&.), shiftR)
import Data.STRef (readSTRef)
import qualified Data.Vector.Unboxed.Mutable as MUV
import Moonlight.Triangulation.Internal.DcelOperations.CandidateArena
( genericCandidate
, growLegalizationArena
, starCandidate
)
import Moonlight.Triangulation.Internal.DcelOperations.FlipRewrite (applyFlip)
import Moonlight.Triangulation.Internal.DcelOperations.FlipRule
( LegalizationLaw (..)
, diagonalFires
)
import Moonlight.Triangulation.Internal.DcelOperations.Twin (reverseIndex)
import Moonlight.Triangulation.Internal.Mutable
( MutableDcel (..)
, readConstraint
, readFace
, readNext
, readOrigin
, readPointX
, readPointY
, readPrevious
)
import Moonlight.Triangulation.Internal.OperationState
( Counter (..)
, OperationState
, legalizationArena
, storeLegalizationArena
)
import Moonlight.Triangulation.Internal.PackedIndex (packIndex)
import Moonlight.Triangulation.Internal.Probe (KnownProbe (..))
-- | The normalization procedure of a confluent terminating rewrite system, and
-- one canonical legalization engine because a normalization procedure is what
-- it is.
--
-- The objects are the triangulations of a fixed point set. The single rule is
-- the Lawson flip: a locally illegal diagonal is replaced by the other
-- diagonal of its quadrilateral. A normal form is a mesh with no illegal
-- diagonal left to fire on.
--
-- /Termination/ is by the lifted-paraboloid potential. Send each point to
-- @(x, y, x² + y²)@ and read a triangulation as a piecewise-linear surface
-- over the point set; @illegalDiagonal@ is exactly the test that the flip
-- lowers that surface, so every rewrite strictly decreases it, and a finite
-- point set has finitely many triangulations. Exact cocircularity is the one
-- case where the surface does not move — the four lifted points are coplanar
-- and both diagonals give the same surface — so there the potential is the
-- diagonal's own key order, and the rule fires only downward in it. Without
-- that tie-break a cocircular quadrilateral flips forever.
--
-- /Confluence/ is Delaunay's theorem, in its strong form: a triangulation with
-- no locally illegal diagonal is globally Delaunay. Local normality is thus
-- global normality, the normal form is unique, and every rewrite order reaches
-- it. That is what licenses the arena below to be a LIFO stack rather than a
-- priority queue, and it is why callers may seed it in whatever order is
-- cheapest to produce — the fan first, or the hull turns first, or both
-- interleaved — without any of them changing the mesh that comes out.
--
-- All callers differ only in how they seed the arena; topology mutation and
-- propagation have exactly one owner.
-- The stack top, the maximum top, and the flip count are strict loop
-- variables, returned once when the drain finishes — the mesh reports nothing
-- per candidate, and 'applyFlip' reports nothing at all. The phantom
-- 'KnownProbe' parameter counts popped candidates for the instrumented lane
-- and is erased everywhere else.
--
-- A popped candidate is read once. Turning it against the star vertex, judging
-- it, and rewriting it are three questions about the same two half-edge
-- records, and the apex the turn looks for is the apex the judgement needs, so
-- one pass over the quadrilateral answers all three.
drainLegalization
:: forall p s vertex directed undirected face
. KnownProbe p
=> MutableDcel s vertex directed undirected face
-> OperationState s
-> Int
-> Int
-> LegalizationLaw
-> ST s (Int, Int)
drainLegalization mutable operation seededTop starVertex law = do
-- No constraint can appear during a drain, so a mesh holding none at entry
-- never needs the per-candidate protection read.
constrained <- readSTRef (mdConstraintCount mutable)
initialArena <- legalizationArena operation
-- A star candidate is turned so that the inserted vertex is the apex
-- opposite the diagonal, so that apex is the same vertex on every star
-- candidate the drain pops and its two coordinates are read once here
-- instead of once per candidate.
starX <- if starVertex < 0 then pure 0 else readPointX mutable starVertex
starY <- if starVertex < 0 then pure 0 else readPointY mutable starVertex
let !guarded = constrained /= 0
-- Everything the cavity fan did not create is pinned. The border loop's
-- legality is already decided by the outside triangle it keeps, and
-- testing it against a neighbourhood that is still inverted could only
-- produce a spurious verdict. 'ValidMesh' pins nothing, which is the
-- floor no undirected index falls below.
!floorPair = case law of
ValidMesh -> 0
CavityRepair floorEdge -> floorEdge
loop !arena !top !maxTop !flips !candidates
| top <= 0 = pure (flips, maxTop, candidates, arena)
| otherwise = do
let !nextTop = top - 1
packedWord <- MUV.unsafeRead arena nextTop
let !packed = fromIntegral packedWord :: Int
!rawEdge = packed `shiftR` 1
!isStar = packed .&. 1 == 0
eligible <-
if rawEdge `shiftR` 1 < floorPair
then pure False
else if guarded then not <$> readConstraint mutable rawEdge else pure True
if not eligible
then loop arena nextTop maxTop flips (probeBump @p candidates)
else do
let !rawTwin = reverseIndex rawEdge
rawFace <- readFace mutable rawEdge
rawTwinFace <- readFace mutable rawTwin
if rawFace == 0 || rawTwinFace == 0
then loop arena nextTop maxTop flips (probeBump @p candidates)
else do
rawBefore <- readPrevious mutable rawEdge
rawTwinBefore <- readPrevious mutable rawTwin
rawApex <- readOrigin mutable rawBefore
rawTwinApex <- readOrigin mutable rawTwinBefore
-- A star candidate names an undirected edge and has to be
-- turned so the inserted vertex is the apex opposite it; a
-- generic one already names the direction to test. Both
-- want the same two apexes, so the one read that settles
-- the turn is also the one that supplies the quadrilateral.
let !turned = isStar && rawApex /= starVertex
if turned && rawTwinApex /= starVertex
then loop arena nextTop maxTop flips (probeBump @p candidates)
else do
let !edge = if turned then rawTwin else rawEdge
!twin = if turned then rawEdge else rawTwin
!leftFace = if turned then rawTwinFace else rawFace
!rightFace = if turned then rawFace else rawTwinFace
!edgePrevious = if turned then rawTwinBefore else rawBefore
!twinPrevious = if turned then rawBefore else rawTwinBefore
!c = if turned then rawTwinApex else rawApex
!d = if turned then rawApex else rawTwinApex
a <- readOrigin mutable edge
b <- readOrigin mutable twin
ax <- readPointX mutable a
ay <- readPointY mutable a
bx <- readPointX mutable b
by <- readPointY mutable b
cx <- if isStar then pure starX else readPointX mutable c
cy <- if isStar then pure starY else readPointY mutable c
dx <- readPointX mutable d
dy <- readPointY mutable d
if not (diagonalFires law ax ay bx by cx cy dx dy)
then loop arena nextTop maxTop flips (probeBump @p candidates)
else do
-- The neighbourhood is complete before the rewrite
-- consumes it, and the pushes carry the same edges
-- in the same order they always did.
edgeNext <- readNext mutable edge
twinNext <- readNext mutable twin
applyFlip
mutable
edge
twin
edgeNext
edgePrevious
twinNext
twinPrevious
leftFace
rightFace
a
b
c
d
if isStar
then do
let !addedTop = nextTop + 2
grown <- growLegalizationArena arena addedTop
MUV.unsafeWrite grown nextTop (packIndex (starCandidate twinPrevious))
MUV.unsafeWrite grown (nextTop + 1) (packIndex (starCandidate twinNext))
loop
grown
addedTop
(max maxTop addedTop)
(flips + 1)
(probeBump @p candidates)
else do
let !addedTop = nextTop + 4
grown <- growLegalizationArena arena addedTop
MUV.unsafeWrite grown nextTop (packIndex (genericCandidate edgeNext))
MUV.unsafeWrite grown (nextTop + 1) (packIndex (genericCandidate edgePrevious))
MUV.unsafeWrite grown (nextTop + 2) (packIndex (genericCandidate twinNext))
MUV.unsafeWrite grown (nextTop + 3) (packIndex (genericCandidate twinPrevious))
loop
grown
addedTop
(max maxTop addedTop)
(flips + 1)
(probeBump @p candidates)
(flips, maxTop, candidates, finalArena) <-
loop initialArena seededTop seededTop 0 (probeZero @p)
when (MUV.length finalArena /= MUV.length initialArena) $
storeLegalizationArena operation finalArena
probeCharge @p operation CounterDiagLegalizationCandidates candidates
pure (flips, maxTop)
{-# INLINE drainLegalization #-}