Refinement.hs

{-# LANGUAGE NamedFieldPuns #-}

-- | Ruppert refinement: Steiner insertion composed after a built mesh, never a
-- second kind of mesh. Parameters are reached through checked verbs rather than
-- a raw record, so an unrealizable quality bar is a refusal.
module Moonlight.Triangulation.Refinement
  ( refine
  , refineWithinDomain
  , validateRefinementParameters
  , withMinimumAngle
  , radiusEdgeRatioForAngle
  ) where

import Control.Monad.ST (ST, runST)
import Data.Foldable (traverse_)
import qualified Data.IntSet as IntSet
import Data.List (sort)
import qualified Data.Map.Strict as Map
import Data.Maybe (fromMaybe)
import qualified Data.Set as Set
import qualified Data.Vector as V
import Moonlight.Triangulation.Dcel
  ( faceVertices
  , incidentFace
  , isConstraintEdge
  , numFaces
  , numUndirectedEdges
  , numVertices
  , undirectedEndpoints
  , vertexPoint
  )
import Moonlight.Triangulation.FloodFillIterator (facesAtEvenBarrierDepth)
import Moonlight.Triangulation.Handles
  ( FaceId (..)
  , UndirectedEdgeId (..)
  , directedPair
  )
import Moonlight.Triangulation.Internal.Mutable
import Moonlight.Triangulation.Internal.Capacity (ensureCapacity)
import Moonlight.Triangulation.Internal.OperationState (freezeBuildStats, newOperationState)
import Moonlight.Triangulation.Internal.Refinement
import Moonlight.Triangulation.Types
import Moonlight.Triangulation.Validation (validateTopology)

data RefinementExecution mode vertex directed undirected face = RefinementExecution
  { refinementExecutionResult :: !(RefinementResult mode vertex directed undirected face)
  , refinementExecutionVisitedFaces :: ![Int]
  , refinementExecutionInterfaceBoundaryReads :: {-# UNPACK #-} !Int
  , refinementExecutionBoundaryCrossingAttempts :: {-# UNPACK #-} !Int
  }

-- | The radius-edge ratio admitting a minimum angle, in degrees.
radiusEdgeRatioForAngle :: Double -> Either BuildError (Maybe Double)
radiusEdgeRatioForAngle degrees =
  case classifyNonFinite degrees of
    Just nonFinite ->
      Left (RefinementMinimumAngleNotFinite nonFinite)
    Nothing
      | degrees < 0 || degrees > 60 ->
          Left (RefinementMinimumAngleOutOfRange degrees)
      | degrees == 0 -> Right Nothing
      | otherwise ->
          let ratio = 0.5 / sin (degrees * pi / 180)
           in case classifyNonFinite ratio of
                Just nonFinite ->
                  Left (RefinementMinimumAngleDerivedRatioNotFinite nonFinite)
                Nothing -> Right (Just ratio)

-- | Set the minimum angle, in degrees.
withMinimumAngle :: Double -> RefinementParameters -> Either BuildError (RefinementParameters)
withMinimumAngle degrees parameters = do
  ratio <- radiusEdgeRatioForAngle degrees
  pure parameters{refineMaxRadiusEdgeRatio = ratio}

-- | Refine an unconstrained or constrained triangulation in the same finite
-- DCEL. The constructor supplies application payloads for Steiner vertices;
-- payloads annotate the requested geometric point and cannot reauthor it.
refine
  :: (Point -> vertex)
  -> RefinementParameters
  -> Triangulation mode vertex directed undirected face
  -> Either BuildError (RefinementResult mode vertex directed undirected face)
refine makeVertex parameters =
  fmap refinementExecutionResult
    . refineWithInitialSeed RefineEveryFace Nothing makeVertex parameters

-- | Prove that a set of active faces is separated from every protected inner
-- face by exactly the supplied immutable edge section. The witness is tied to
-- the input arena cardinalities and is revalidated before interpretation.
mkRefinementDomain
  :: Set.Set FaceId
  -> Set.Set UndirectedEdgeId
  -> Triangulation mode vertex directed undirected face
  -> Either BuildError RefinementDomain
mkRefinementDomain permittedFaces interfaceEdges triangulation = do
  permitted <- validateSeedFaces triangulation permittedFaces
  interface <- IntSet.fromList <$> traverse validateInterfaceEdge (Set.toAscList interfaceEdges)
  let expected = expectedInterface permitted
  case IntSet.lookupMin (expected IntSet.\\ interface) of
    Just pair -> Left (RefinementDomainInterfaceMissing (UndirectedEdgeId (fromIntegral pair)))
    Nothing ->
      case IntSet.lookupMin (interface IntSet.\\ expected) of
        Just pair -> Left (RefinementDomainInterfaceExtraneous (UndirectedEdgeId (fromIntegral pair)))
        Nothing ->
          Right
            RefinementDomain
              { refinementDomainPermittedFaces = permitted
              , refinementDomainInterfacePairs = interface
              , refinementDomainInputFaces =
                  Map.fromAscList
                    [ (face, faceSignature triangulation face)
                    | raw <- [1 .. numFaces triangulation - 1]
                    , let face = FaceId (fromIntegral raw)
                    ]
              , refinementDomainInputFaceCount = numFaces triangulation
              , refinementDomainInputEdgeCount = numUndirectedEdges triangulation
              }
 where
  totalEdges = numUndirectedEdges triangulation
  validateInterfaceEdge edge@(UndirectedEdgeId raw)
    | toInteger raw >= toInteger totalEdges =
        Left (RefinementDomainInterfaceEdgeNotActive edge totalEdges)
    | otherwise = Right (fromIntegral raw)

  expectedInterface permitted =
    IntSet.fromList
      [ pair
      | pair <- [0 .. totalEdges - 1]
      , let edge = UndirectedEdgeId (fromIntegral pair)
            (forward, backward) = directedPair edge
            FaceId forwardFace = incidentFace triangulation forward
            FaceId backwardFace = incidentFace triangulation backward
            forwardInner = forwardFace /= 0
            backwardInner = backwardFace /= 0
            forwardPermitted = IntSet.member (fromIntegral forwardFace) permitted
            backwardPermitted = IntSet.member (fromIntegral backwardFace) permitted
      , forwardInner && backwardInner && forwardPermitted /= backwardPermitted
      ]

-- | Refine exactly one checked local section. Interface edges are installed as
-- transaction-local legalization barriers and removed before publication;
-- every protected face signature is then compared with the input witness.
refineWithinDomain
  :: (Point -> vertex)
  -> RefinementParameters
  -> Set.Set FaceId
  -> Set.Set UndirectedEdgeId
  -> Triangulation mode vertex directed undirected face
  -> Either BuildError (RefinementDomainResult mode vertex directed undirected face)
refineWithinDomain makeVertex parameters permittedFaces interfaceEdges triangulation = do
  validateDomainParameters parameters
  domain <- mkRefinementDomain permittedFaces interfaceEdges triangulation
  execution <-
    refineWithInitialSeed
      (RefineSeededFaces (refinementDomainPermittedFaces domain))
      (Just domain)
      makeVertex
      parameters
      triangulation
  let result = refinementExecutionResult execution
      receipt =
        buildRefinementReceipt
          domain
          triangulation
          (refinedTriangulation result)
          (refinementExecutionVisitedFaces execution)
          (refinementExecutionInterfaceBoundaryReads execution)
          (refinementExecutionBoundaryCrossingAttempts execution)
  validateProtectedFaces domain triangulation (refinedTriangulation result)
  case V.toList (refinementVisitedProtectedFaces receipt) of
    protected : _ -> Left (RefinementDomainWouldRewriteProtectedFace protected)
    [] ->
      Right
        RefinementDomainResult
          { refinementDomainResult = result
          , refinementDomainReceipt = receipt
          }

validateDomainParameters :: RefinementParameters -> Either BuildError ()
validateDomainParameters parameters
  | not (refinePreserveConvexHull parameters) =
      Left RefinementDomainRequiresConvexHullPreservation
  | not (refineKeepConstraintEdges parameters) =
      Left RefinementDomainRequiresConstraintPreservation
  | refineExcludeOuterFaces parameters =
      Left RefinementDomainForbidsOuterFaceExclusion
  | otherwise = Right ()

refineWithInitialSeed
  :: RefinementInitialSeed
  -> Maybe RefinementDomain
  -> (Point -> vertex)
  -> RefinementParameters
  -> Triangulation mode vertex directed undirected face
  -> Either BuildError (RefinementExecution mode vertex directed undirected face)
refineWithInitialSeed initialSeed domain makeVertex parameters triangulation = do
  validateRefinementParameters parameters
  case validateTopology triangulation of
    violation : _ -> Left (RefinementInputTopologyInvalid violation)
    [] -> pure ()
  let originalCount = numVertices triangulation
      budget = max 0 (fromMaybe (10 * max 1 originalCount) (refineMaxAdditionalVertices parameters))
      maximumVerticesInteger = toInteger originalCount + toInteger budget
      initialExcludedFaces =
        if refineExcludeOuterFaces parameters
          then
            IntSet.fromList
              [ fromIntegral raw
              | FaceId raw <-
                  facesAtEvenBarrierDepth
                    triangulation
                    (isConstraintEdge triangulation)
              ]
          else IntSet.empty
  let maximumVertices =
        if maximumVerticesInteger > toInteger (maxBound :: Int)
          then maxBound
          else fromInteger maximumVerticesInteger
  ensureCapacity maximumVertices
  runST $ do
      mutable <-
        thawTriangulationDense
          maximumVertices
          triangulation
      operation <- newOperationState (halfEdgeCapacity mutable)
      installInterfaceBarriers mutable domain
      outcome <-
        refineMutable
          makeVertex
          mutable
          operation
          parameters
          originalCount
          initialExcludedFaces
          initialSeed
          domain
      removeInterfaceBarriers mutable triangulation domain
      case outcome of
        Left failure -> pure (Left failure)
        Right (complete, added, excluded, visited, interfaceBoundaryReads, boundaryCrossingAttempts) -> do
          frozenOutcome <- freezeTriangulation mutable
          case frozenOutcome of
            Left obstruction -> pure (Left obstruction)
            Right frozen -> do
              stats <- freezeBuildStats operation
              pure
                ( Right
                    RefinementExecution
                      { refinementExecutionResult =
                          RefinementResult
                            { refinedTriangulation = frozen
                            , refinementStats = stats
                            , refinementAddedVertices = added
                            , refinementComplete = complete
                            , refinementExcludedFaces = V.fromList (map (FaceId . fromIntegral) excluded)
                            }
                      , refinementExecutionVisitedFaces = visited
                      , refinementExecutionInterfaceBoundaryReads = interfaceBoundaryReads
                      , refinementExecutionBoundaryCrossingAttempts = boundaryCrossingAttempts
                      }
                )

installInterfaceBarriers
  :: MutableDcel s vertex directed undirected face
  -> Maybe RefinementDomain
  -> ST s ()
installInterfaceBarriers mutable =
  traverse_
    (\pair -> setConstraint mutable (2 * pair) >> pure ())
    . maybe [] (IntSet.toAscList . refinementDomainInterfacePairs)

removeInterfaceBarriers
  :: MutableDcel s vertex directed undirected face
  -> Triangulation mode vertex directed undirected face
  -> Maybe RefinementDomain
  -> ST s ()
removeInterfaceBarriers mutable input =
  traverse_
    (\pair ->
       let edge = UndirectedEdgeId (fromIntegral pair)
        in if isConstraintEdge input edge
             then pure ()
             else clearConstraint mutable (2 * pair) >> pure ()
    )
    . maybe [] (IntSet.toAscList . refinementDomainInterfacePairs)

validateProtectedFaces
  :: RefinementDomain
  -> Triangulation mode vertex directed undirected face
  -> Triangulation mode vertex directed undirected face
  -> Either BuildError ()
validateProtectedFaces domain before after =
  traverse_ validateProtected protectedFaces
 where
  protectedFaces =
    [ FaceId (fromIntegral face)
    | face <- [1 .. refinementDomainInputFaceCount domain - 1]
    , IntSet.notMember face (refinementDomainPermittedFaces domain)
    ]
  validateProtected face
    | Map.lookup face (refinementDomainInputFaces domain)
        == Just (faceSignature before face)
        && Map.lookup face (refinementDomainInputFaces domain)
          == Just (faceSignature after face) = Right ()
    | otherwise = Left (RefinementDomainProtectedFaceChanged face)

faceSignature
  :: Triangulation mode vertex directed undirected face
  -> FaceId
  -> [Point]
faceSignature triangulation =
  sort . fmap (vertexPoint triangulation) . faceVertices triangulation

type EdgeSignature = (Point, Point, [FaceId])

edgeSignature
  :: Triangulation mode vertex directed undirected face
  -> UndirectedEdgeId
  -> EdgeSignature
edgeSignature triangulation edge =
  (min fromPoint toPoint, max fromPoint toPoint, sort [forwardFace, backwardFace])
 where
  (fromVertex, toVertex) = undirectedEndpoints triangulation edge
  fromPoint = vertexPoint triangulation fromVertex
  toPoint = vertexPoint triangulation toVertex
  (forward, backward) = directedPair edge
  forwardFace = incidentFace triangulation forward
  backwardFace = incidentFace triangulation backward

buildRefinementReceipt
  :: RefinementDomain
  -> Triangulation mode vertex directed undirected face
  -> Triangulation mode vertex directed undirected face
  -> [Int]
  -> Int
  -> Int
  -> RefinementReceipt
buildRefinementReceipt domain before after visited interfaceBoundaryReads boundaryCrossingAttempts =
  RefinementReceipt
    { refinementVisitedJoinFaces = V.fromList (fmap toFace visitedJoin)
    , refinementVisitedProtectedFaces = V.fromList (fmap toFace visitedProtected)
    , refinementCreatedFaces =
        V.fromList
          [ targetFace
          | face <- [1 .. numFaces after - 1]
          , let targetFace = toFace face
          , face >= numFaces before
              || faceSignature before targetFace /= faceSignature after targetFace
          ]
    , refinementTouchedEdges = V.fromList touched
    , refinementRemovedEdges = V.fromList removed
    , refinementInterfaceBoundaryReads = interfaceBoundaryReads
    , refinementAttemptedBoundaryCrossings = boundaryCrossingAttempts
    }
 where
  permitted = refinementDomainPermittedFaces domain
  (visitedProtected, visitedJoin) =
    foldr
      (\face (protected, join) ->
         if IntSet.member face permitted || face >= refinementDomainInputFaceCount domain
           then (protected, face : join)
           else (face : protected, join)
      )
      ([], [])
      visited
  toFace = FaceId . fromIntegral
  beforeEdges = numUndirectedEdges before
  afterEdges = numUndirectedEdges after
  touched =
    [ edge
    | raw <- [0 .. afterEdges - 1]
    , let edge = UndirectedEdgeId (fromIntegral raw)
    , raw >= beforeEdges || edgeSignature before edge /= edgeSignature after edge
    ]
  removed =
    [ edge
    | raw <- [0 .. beforeEdges - 1]
    , let edge = UndirectedEdgeId (fromIntegral raw)
    , raw >= afterEdges || edgeSignature before edge /= edgeSignature after edge
    ]

-- | Refuse outer or absent faces rather than silently treating an invalid
-- topology witness as an empty repair. Face zero is the outer face and has no
-- active refinement equation.
validateSeedFaces
  :: Triangulation mode vertex directed undirected face
  -> Set.Set FaceId
  -> Either BuildError IntSet.IntSet
validateSeedFaces triangulation requested =
  IntSet.fromList <$> traverse validateFace (Set.toAscList requested)
 where
  totalFaces = numFaces triangulation

  validateFace face@(FaceId raw)
    | raw == 0 || toInteger raw >= toInteger totalFaces =
        Left (RefinementSeedFaceNotActive face totalFaces)
    | otherwise = Right (fromIntegral raw)

-- | Refuse parameters no mesh can satisfy.
validateRefinementParameters :: RefinementParameters -> Either BuildError ()
validateRefinementParameters RefinementParameters{refineMaxAdditionalVertices, refineMinArea, refineMaxArea, refineMaxRadiusEdgeRatio} = do
  case refineMaxAdditionalVertices of
    Just value
      | value < 0 ->
          Left (RefinementMaximumAdditionalVerticesNegative value)
    _ -> Right ()
  validateOptionalRefinementParameter
    RefinementMinimumAreaNotFinite
    RefinementMinimumAreaNegative
    (>= 0)
    refineMinArea
  validateOptionalRefinementParameter
    RefinementMaximumAreaNotFinite
    RefinementMaximumAreaNotPositive
    (> 0)
    refineMaxArea
  validateOptionalRefinementParameter
    RefinementMaximumRadiusEdgeRatioNotFinite
    RefinementMaximumRadiusEdgeRatioNotPositive
    (> 0)
    refineMaxRadiusEdgeRatio
  case (refineMinArea, refineMaxArea) of
    (Just minimumArea, Just maximumArea)
      | minimumArea > maximumArea ->
          Left
            ( RefinementMinimumAreaExceedsMaximum
                minimumArea
                maximumArea
            )
    _ -> Right ()

validateOptionalRefinementParameter
  :: (NonFiniteValue -> BuildError)
  -> (Double -> BuildError)
  -> (Double -> Bool)
  -> Maybe Double
  -> Either BuildError ()
validateOptionalRefinementParameter nonFinite outsideRange predicate value =
  case value of
    Nothing -> Right ()
    Just number ->
      case classifyNonFinite number of
        Just nonFiniteValue ->
          Left (nonFinite nonFiniteValue)
        Nothing
          | predicate number -> Right ()
          | otherwise ->
              Left (outsideRange number)