moonlight-triangulation-1.2.0.1: src-dcel/Moonlight/Triangulation/FloodFillIterator.hs
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE FlexibleInstances #-}
-- | Shape queries and face flood fills over immutable triangulations.
module Moonlight.Triangulation.FloodFillIterator
( DistanceMetric (..)
, CircleMetric
, CircleMetricError (..)
, RadiusSquaredError (..)
, RectangleMetric
, RectangleMetricError (..)
, circleMetric
, rectangleMetric
, edgesInShape
, verticesInShape
, edgesInCircle
, verticesInCircle
, edgesInRectangle
, verticesInRectangle
, floodFillFaces
, outerFaceFloodFill
, facesAtEvenBarrierDepth
, FaceComponent
, faceComponentFaces
, BoundaryLoop
, boundaryLoopVertices
, RegionBoundary
, regionBoundaryOuterLoop
, regionBoundaryHoleLoops
, BoundaryObstruction (..)
, faceComponents
, labelledRegionBoundaries
, componentBoundary
, RadiusSquared
, mkRadiusSquared
, alphaShapeContainsFace
) where
import Control.DeepSeq (NFData)
import Data.Bifunctor (first)
import qualified Data.IntMap.Strict as IntMap
import qualified Data.IntSet as IntSet
import Data.List (partition, unfoldr)
import qualified Data.List as List
import Data.List.NonEmpty (NonEmpty (..))
import qualified Data.Sequence as Seq
import qualified Data.Vector as V
import GHC.Generics (Generic)
import Moonlight.Triangulation.Dcel
import Moonlight.Triangulation.Handles.HandleDefs
import Moonlight.Triangulation.Handles.Iterators.FixedIterators (undirectedEdges)
import Moonlight.Triangulation.Internal.BoundaryCycle (simplifyBoundaryCycle)
import Moonlight.Triangulation.Math
import Moonlight.Triangulation.PointLocation
import Moonlight.Triangulation.Scalar (circumradiusSquaredWithinCoordinates)
import Moonlight.Triangulation.Types
-- | One non-empty connected set of equally labelled bounded face indices.
newtype FaceComponent = FaceComponent IntSet.IntSet
deriving stock (Eq, Show, Generic)
deriving anyclass (NFData)
-- | Component faces in ascending DCEL order.
faceComponentFaces :: FaceComponent -> [FaceId]
faceComponentFaces (FaceComponent faces) =
fmap (FaceId . fromIntegral) (IntSet.toAscList faces)
-- | One non-empty simple boundary loop. Outer loops are counter-clockwise and
-- hole loops clockwise.
newtype BoundaryLoop = BoundaryLoop
{ boundaryLoopVertices :: NonEmpty VertexId
}
deriving stock (Eq, Show, Generic)
deriving anyclass (NFData)
-- | The authoritative boundary of one face component.
data RegionBoundary = RegionBoundary
{ regionBoundaryOuterLoop :: !BoundaryLoop
, regionBoundaryHoleLoops :: ![BoundaryLoop]
}
deriving stock (Eq, Show, Generic)
deriving anyclass (NFData)
-- | Typed failure to descend boundary half-edges into simple polygon loops.
data BoundaryObstruction
= BoundaryComponentFaceOutOfRange !FaceId {-# UNPACK #-} !Int
| BoundaryPinch !VertexId !DirectedEdgeId !DirectedEdgeId
| BoundaryCycleDidNotClose !DirectedEdgeId !DirectedEdgeId
| BoundaryLoopDegenerate ![VertexId]
| BoundaryOuterLoopCardinality !Int
deriving stock (Eq, Show, Generic)
deriving anyclass (NFData)
-- | An admitted finite, non-negative squared radius.
newtype RadiusSquared = RadiusSquared Double
deriving stock (Eq, Ord, Show, Generic)
deriving anyclass (NFData)
-- | Typed refusal shared by every squared-radius constructor.
data RadiusSquaredError
= NonFiniteRadiusSquared !NonFiniteValue
| NegativeRadiusSquared !Double
deriving stock (Eq, Ord, Show)
-- | Admit a finite, non-negative squared radius.
mkRadiusSquared :: Double -> Either RadiusSquaredError RadiusSquared
mkRadiusSquared value =
case classifyNonFinite value of
Just nonFinite -> Left (NonFiniteRadiusSquared nonFinite)
Nothing
| value < 0 -> Left (NegativeRadiusSquared value)
| otherwise -> Right (RadiusSquared value)
-- | A query shape that can admit points, test edges, and supply a location
-- seed.
class DistanceMetric metric where
metricContainsPoint :: metric -> Point -> Bool
metricIntersectsEdge :: metric -> Point -> Point -> Bool
metricStartPoint :: metric -> QueryPoint
-- | An admitted center and squared radius.
data CircleMetric = CircleMetric !QueryPoint !RadiusSquared
deriving stock (Eq, Ord, Show)
-- | Typed refusal for an invalid circle query.
data CircleMetricError
= InvalidCircleCenter !PointValidationError
| InvalidCircleRadius !RadiusSquaredError
deriving stock (Eq, Ord, Show)
-- | Admitted lower corner, upper corner, and center of an axis-aligned box.
data RectangleMetric = RectangleMetric !(QueryPoint) !(QueryPoint) !(QueryPoint)
deriving stock (Eq, Ord, Show)
-- | Typed refusal for an invalid rectangle query.
data RectangleMetricError
= InvalidRectangleLower !PointValidationError
| InvalidRectangleUpper !PointValidationError
| InvalidRectangleCenter !PointValidationError
deriving stock (Eq, Ord, Show)
-- | A circle metric, or why the radius is unusable.
circleMetric :: Point -> Double -> Either CircleMetricError CircleMetric
circleMetric center radiusSquared =
CircleMetric
<$> first InvalidCircleCenter (mkQueryPoint center)
<*> first InvalidCircleRadius (mkRadiusSquared radiusSquared)
-- | An axis-aligned rectangle metric, or why the corners are unusable.
rectangleMetric :: Point -> Point -> Either RectangleMetricError RectangleMetric
rectangleMetric lower@(Point lowerX lowerY) upper@(Point upperX upperY) = do
queryLower <- either (Left . InvalidRectangleLower) Right (mkQueryPoint lower)
queryUpper <- either (Left . InvalidRectangleUpper) Right (mkQueryPoint upper)
queryCenter <-
either
(Left . InvalidRectangleCenter)
Right
(mkQueryPoint (Point ((lowerX + upperX) * 0.5) ((lowerY + upperY) * 0.5)))
Right (RectangleMetric queryLower queryUpper queryCenter)
instance DistanceMetric CircleMetric where
metricContainsPoint (CircleMetric center (RadiusSquared radiusSquared)) point =
squaredDistanceWide (queryPointValue center) point <= radiusSquared
metricIntersectsEdge (CircleMetric center (RadiusSquared radiusSquared)) from to =
segmentDistanceSquaredWide from to (queryPointValue center) <= radiusSquared
metricStartPoint (CircleMetric center _) = center
instance DistanceMetric RectangleMetric where
metricContainsPoint (RectangleMetric lower upper _) (Point x y) =
lowerX <= upperX && lowerY <= upperY && x >= lowerX && x <= upperX && y >= lowerY && y <= upperY
where
Point lowerX lowerY = queryPointValue lower
Point upperX upperY = queryPointValue upper
metricIntersectsEdge rectangle from to =
metricContainsPoint rectangle from
|| metricContainsPoint rectangle to
|| segmentRectangleIntersection rectangle from to
metricStartPoint (RectangleMetric _ _ center) = center
-- | Edges meeting a circle.
edgesInCircle :: Triangulation mode vertex directed undirected face -> Point -> Double -> Either CircleMetricError [UndirectedEdgeId]
edgesInCircle triangulation center radiusSquared =
edgesInShape triangulation <$> circleMetric center radiusSquared
-- | Vertices inside a circle.
verticesInCircle :: Triangulation mode vertex directed undirected face -> Point -> Double -> Either CircleMetricError [VertexId]
verticesInCircle triangulation center radiusSquared =
verticesInShape triangulation <$> circleMetric center radiusSquared
-- | Edges meeting an axis-aligned rectangle.
edgesInRectangle :: Triangulation mode vertex directed undirected face -> Point -> Point -> Either RectangleMetricError [UndirectedEdgeId]
edgesInRectangle triangulation lower upper = edgesInShape triangulation <$> rectangleMetric lower upper
-- | Vertices inside an axis-aligned rectangle.
verticesInRectangle :: Triangulation mode vertex directed undirected face -> Point -> Point -> Either RectangleMetricError [VertexId]
verticesInRectangle triangulation lower upper = verticesInShape triangulation <$> rectangleMetric lower upper
-- | Edges meeting any metric shape.
edgesInShape :: DistanceMetric metric => Triangulation mode vertex directed undirected face -> metric -> [UndirectedEdgeId]
edgesInShape triangulation metric
| numVertices triangulation <= 1 = []
| numInnerFaces triangulation == 0 =
[edge | edge <- undirectedEdges triangulation, edgeInside edge]
| otherwise =
let starts = shapeStartFaces triangulation metric
(_, accepted) = floodFillFacesWithEdges triangulation starts edgeInside
in map (UndirectedEdgeId . fromIntegral) (IntSet.toAscList accepted)
where
edgeInside edge =
let (fromVertex, toVertex) = undirectedEndpoints triangulation edge
in metricIntersectsEdge metric (vertexPoint triangulation fromVertex) (vertexPoint triangulation toVertex)
-- | Vertices inside any metric shape.
verticesInShape :: DistanceMetric metric => Triangulation mode vertex directed undirected face -> metric -> [VertexId]
verticesInShape triangulation metric =
[ vertex
| vertex <- candidateVertices
, metricContainsPoint metric (vertexPoint triangulation vertex)
]
where
edges = edgesInShape triangulation metric
set = List.foldl' addEndpoints IntSet.empty edges
addEndpoints acc edge =
let (VertexId from, VertexId to) = undirectedEndpoints triangulation edge
in IntSet.insert (fromIntegral from) (IntSet.insert (fromIntegral to) acc)
candidateVertices
| numVertices triangulation == 1 = [VertexId 0]
| otherwise = map (VertexId . fromIntegral) (IntSet.toAscList set)
-- | Reach inner faces from the supplied seeds by crossing only admitted edges.
floodFillFaces
:: Triangulation mode vertex directed undirected face -> [FaceId]
-> (UndirectedEdgeId -> Bool)
-> [FaceId]
floodFillFaces triangulation starts canCross =
fst (floodFillFacesWithEdges triangulation starts canCross)
floodFillFacesWithEdges
:: Triangulation mode vertex directed undirected face -> [FaceId]
-> (UndirectedEdgeId -> Bool)
-> ([FaceId], IntSet.IntSet)
floodFillFacesWithEdges triangulation starts canCross =
let (faces, accepted, _) = go initialStack initialVisited IntSet.empty IntSet.empty []
in (reverse faces, accepted)
where
valid face@(FaceId value) = face /= outerFace && fromIntegral value < numFaces triangulation
(initialStack, initialVisited) = List.foldl' enqueueStart ([], IntSet.empty) starts
enqueueStart state face
| valid face = enqueue face state
| otherwise = state
go [] _ accepted rejected result = (result, accepted, rejected)
go (face : stack) visited accepted rejected result =
let (stack', visited', accepted', rejected') =
foldFaceDirectedEdges'
triangulation
face
expand
(stack, visited, accepted, rejected)
in go stack' visited' accepted' rejected' (face : result)
expand (stack, visited, accepted, rejected) edge =
let undirected@(UndirectedEdgeId raw) = asUndirected edge
edgeIndex = fromIntegral raw
adjacent = incidentFace triangulation (reverseEdge edge)
edgeAdmission
| IntSet.member edgeIndex accepted = (True, accepted, rejected)
| IntSet.member edgeIndex rejected = (False, accepted, rejected)
| canCross undirected = (True, IntSet.insert edgeIndex accepted, rejected)
| otherwise = (False, accepted, IntSet.insert edgeIndex rejected)
(crosses, accepted', rejected') = edgeAdmission
(stack', visited') =
if crosses && valid adjacent
then enqueue adjacent (stack, visited)
else (stack, visited)
in (stack', visited', accepted', rejected')
enqueue face@(FaceId value) (stack, visited)
| IntSet.member index visited = (stack, visited)
| otherwise = (face : stack, IntSet.insert index visited)
where
index = fromIntegral value
-- | Connected components of equally labelled bounded faces. Labels are
-- evaluated once; the component carrier is the same 'IntSet' used by descent.
faceComponents
:: Eq label
=> Triangulation mode vertex directed undirected face
-> (FaceId -> label)
-> [(label, FaceComponent)]
faceComponents triangulation labelFace = unfoldr descend initialUnvisited
where
labels =
V.generate
(numInnerFaces triangulation)
(\index -> labelFace (FaceId (fromIntegral (index + 1))))
initialUnvisited =
IntSet.fromRange (1, numFaces triangulation - 1)
labelAt faceIndex = labels V.!? (faceIndex - 1)
descend remaining =
case IntSet.minView remaining of
Nothing -> Nothing
Just (seedIndex, unseeded) ->
case labelAt seedIndex of
Nothing -> descend unseeded
Just componentLabel ->
let unvisited =
collectComponent componentLabel (Seq.singleton seedIndex) unseeded
componentFaces = IntSet.difference remaining unvisited
in Just ((componentLabel, FaceComponent componentFaces), unvisited)
collectComponent componentLabel queued unvisited =
case Seq.viewl queued of
Seq.EmptyL -> unvisited
faceIndex Seq.:< remainingQueue ->
let face = FaceId (fromIntegral faceIndex)
(expandedQueue, remainingUnvisited) =
foldFaceDirectedEdges'
triangulation
face
(admitAdjacent componentLabel)
(remainingQueue, unvisited)
in collectComponent
componentLabel
expandedQueue
remainingUnvisited
admitAdjacent componentLabel (queued, unvisited) edge =
if
adjacent /= outerFace
&& IntSet.member adjacentIndex unvisited
&& labelAt adjacentIndex == Just componentLabel
then
( queued Seq.|> adjacentIndex
, IntSet.delete adjacentIndex unvisited
)
else (queued, unvisited)
where
adjacent@(FaceId adjacentRaw) =
incidentFace triangulation (reverseEdge edge)
adjacentIndex = fromIntegral adjacentRaw
-- | Descend every equally labelled bounded-face component through the one
-- authoritative boundary tracer. Components are converted independently;
-- callers may group equal labels only after this descent has succeeded.
labelledRegionBoundaries
:: Eq label
=> Triangulation mode vertex directed undirected face
-> (FaceId -> label)
-> Either BoundaryObstruction [(label, RegionBoundary)]
labelledRegionBoundaries triangulation labelFace =
traverse
(\(label, component) -> (label,) <$> componentBoundary triangulation component)
(faceComponents triangulation labelFace)
-- | Extract all simple boundary loops of a component from the DCEL. Boundary
-- half-edges retain their incident component face on the left; this gives the
-- outer loop counter-clockwise and holes clockwise without a later guess.
componentBoundary
:: Triangulation mode vertex directed undirected face
-> FaceComponent
-> Either BoundaryObstruction RegionBoundary
componentBoundary triangulation (FaceComponent componentFaces) = do
case IntSet.lookupGE (numFaces triangulation) componentFaces of
Just invalid ->
Left
( BoundaryComponentFaceOutOfRange
(FaceId (fromIntegral invalid))
(numFaces triangulation)
)
Nothing -> Right ()
(boundaryEdges, outgoingSuccessor) <-
IntSet.foldl'
collectFaceBoundaryEdges
(Right (IntSet.empty, IntMap.empty))
componentFaces
orientedLoops <- traceBoundaryLoops triangulation outgoingSuccessor boundaryEdges
let (outer, holes) = partition ((== GT) . fst) orientedLoops
outerLoops = fmap snd outer
holeLoops = fmap snd holes
case outerLoops of
[outerLoop] ->
Right
RegionBoundary
{ regionBoundaryOuterLoop = outerLoop
, regionBoundaryHoleLoops = holeLoops
}
_ -> Left (BoundaryOuterLoopCardinality (length outerLoops))
where
collectFaceBoundaryEdges boundaryGraph face =
foldFaceDirectedEdges'
triangulation
(FaceId (fromIntegral face))
insertBoundaryEdge
boundaryGraph
insertBoundaryEdge outcome edge = do
graph@(edges, successors) <- outcome
let FaceId adjacent = incidentFace triangulation (reverseEdge edge)
if IntSet.member (fromIntegral adjacent) componentFaces
then Right graph
else
let vertex@(VertexId rawVertex) = origin triangulation edge
vertexIndex = fromIntegral rawVertex
DirectedEdgeId rawEdge = edge
in case IntMap.lookup vertexIndex successors of
Just previousEdge -> Left (BoundaryPinch vertex previousEdge edge)
Nothing ->
Right
( IntSet.insert (fromIntegral rawEdge) edges
, IntMap.insert vertexIndex edge successors
)
traceBoundaryLoops
:: Triangulation mode vertex directed undirected face
-> IntMap.IntMap DirectedEdgeId
-> IntSet.IntSet
-> Either BoundaryObstruction [(Ordering, BoundaryLoop)]
traceBoundaryLoops triangulation outgoingByVertex = descend []
where
descend loops unvisited =
case IntSet.minView unvisited of
Nothing -> Right (reverse loops)
Just (rawStart, _) -> do
let start = DirectedEdgeId (fromIntegral rawStart)
(vertices, remaining) <- traceCycle start start unvisited []
oriented <- simplifyBoundaryLoop triangulation vertices
descend (oriented : loops) remaining
traceCycle start current unvisited reversedVertices =
let DirectedEdgeId rawCurrent = current
remaining = IntSet.delete (fromIntegral rawCurrent) unvisited
accumulated = origin triangulation current : reversedVertices
VertexId rawTarget = destination triangulation current
in case IntMap.lookup (fromIntegral rawTarget) outgoingByVertex of
Just successor
| successor == start -> Right (reverse accumulated, remaining)
| let DirectedEdgeId rawSuccessor = successor
, IntSet.member (fromIntegral rawSuccessor) remaining ->
traceCycle start successor remaining accumulated
| otherwise -> Left (BoundaryCycleDidNotClose start successor)
Nothing -> Left (BoundaryCycleDidNotClose start current)
simplifyBoundaryLoop
:: Triangulation mode vertex directed undirected face
-> [VertexId]
-> Either BoundaryObstruction (Ordering, BoundaryLoop)
simplifyBoundaryLoop triangulation =
fmap (\(orientation, vertices) -> (orientation, BoundaryLoop vertices))
. simplifyBoundaryCycle BoundaryLoopDegenerate redundant winding key
where
point vertex = vertexPoint triangulation vertex
redundant previousVertex current nextVertex =
orient2d (point previousVertex) (point current) (point nextVertex) == EQ
&& onClosedSegment (point previousVertex) (point nextVertex) (point current)
winding previousVertex current nextVertex =
orient2d (point previousVertex) (point current) (point nextVertex)
key vertex = (point vertex, vertex)
-- | Membership of a bounded face in the closed alpha shape. Exact dyadic
-- comparison makes equality independent of circumcenter rounding.
alphaShapeContainsFace
:: RadiusSquared
-> Triangulation 'Unconstrained vertex directed undirected face
-> FaceId
-> Bool
alphaShapeContainsFace (RadiusSquared threshold) triangulation =
maybe False withinRadius . innerFaceVertices triangulation
where
withinRadius (firstVertex, secondVertex, thirdVertex) =
let Point ax ay = vertexPoint triangulation firstVertex
Point bx by = vertexPoint triangulation secondVertex
Point cx cy = vertexPoint triangulation thirdVertex
in circumradiusSquaredWithinCoordinates threshold ax ay bx by cx cy
-- | Inner faces separated from the outer face by an even minimum number of
-- barriers. A 0–1 BFS floods freely within one depth before crossing a barrier,
-- so a free-ended barrier can be walked around at depth zero while nested
-- closed barriers alternate outside and inside.
facesAtEvenBarrierDepth
:: Triangulation mode vertex directed undirected face
-> (UndirectedEdgeId -> Bool)
-> [FaceId]
facesAtEvenBarrierDepth triangulation isBarrier =
concat (evenLayers (barrierDepthLayers triangulation isBarrier))
where
evenLayers :: [[FaceId]] -> [[FaceId]]
evenLayers (outsideLayer : _insideLayer : deeper) =
outsideLayer : evenLayers deeper
evenLayers shallow = shallow
barrierDepthLayers
:: Triangulation mode vertex directed undirected face
-> (UndirectedEdgeId -> Bool)
-> [[FaceId]]
barrierDepthLayers triangulation isBarrier =
map (filter (/= outerFace)) (layers IntSet.empty [outerFace])
where
known (FaceId value) = fromIntegral value < numFaces triangulation
key :: FaceId -> Int
key (FaceId value) = fromIntegral value
layers visited frontier = case flood visited [] frontier of
([], _) -> []
(layer, visited') -> layer : layers visited' (concatMap (neighbours isBarrier) layer)
flood visited acc [] = (reverse acc, visited)
flood visited acc (face : rest)
| not (known face) || IntSet.member (key face) visited = flood visited acc rest
| otherwise =
flood
(IntSet.insert (key face) visited)
(face : acc)
(neighbours (not . isBarrier) face <> rest)
neighbours admit face =
[ incidentFace triangulation (reverseEdge edge)
| edge <- faceDirectedEdges triangulation face
, admit (asUndirected edge)
]
-- | Faces reachable from the outer face without crossing a barrier edge.
outerFaceFloodFill :: Triangulation mode vertex directed undirected face -> (UndirectedEdgeId -> Bool) -> [FaceId]
outerFaceFloodFill triangulation canCross = floodFillFaces triangulation starts canCross
where
starts =
[ face
| outerEdge <- faceDirectedEdges triangulation outerFace
, let edge = asUndirected outerEdge
, canCross edge
, let face = incidentFace triangulation (reverseEdge outerEdge)
, face /= outerFace
]
-- | The faces a shape's start point lands in.
shapeStartFaces :: DistanceMetric metric => Triangulation mode vertex directed undirected face -> metric -> [FaceId]
shapeStartFaces triangulation metric =
case locatePoint triangulation (metricStartPoint metric) of
InFace face -> [face]
OnEdge edge -> filter (/= outerFace) [incidentFace triangulation edge, incidentFace triangulation (reverseEdge edge)]
OnVertex vertex ->
intSetToFaces
(List.foldl' (\set edge -> let FaceId value = incidentFace triangulation edge in if value == 0 then set else IntSet.insert (fromIntegral value) set) IntSet.empty (vertexOutgoingEdges triangulation vertex))
OutsideConvexHull _ ->
[ incidentFace triangulation (reverseEdge edge)
| edge <- faceDirectedEdges triangulation outerFace
, let from = vertexPoint triangulation (origin triangulation edge)
, let to = vertexPoint triangulation (destination triangulation edge)
, metricIntersectsEdge metric from to
, incidentFace triangulation (reverseEdge edge) /= outerFace
]
EmptyTriangulation -> []
where
intSetToFaces = map (FaceId . fromIntegral) . IntSet.toAscList
segmentRectangleIntersection :: RectangleMetric -> Point -> Point -> Bool
segmentRectangleIntersection (RectangleMetric lowerQuery upperQuery _) from to
| lx > ux || ly > uy = False
| lower == upper = onClosedSegment from to lower
| otherwise = any (uncurry (segmentsIntersect from to)) boundaries
where
lower@(Point lx ly) = queryPointValue lowerQuery
upper@(Point ux uy) = queryPointValue upperQuery
boundaries =
[ (Point lx ly, Point lx uy)
, (Point lx uy, Point ux uy)
, (Point ux uy, Point ux ly)
, (Point ux ly, Point lx ly)
]