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 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 = 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) = 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
        (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)
    ]