moonlight-triangulation-1.2.0.1: src-dcel/Moonlight/Triangulation/Valuation.hs
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE GADTs #-}
-- | Intrinsic valuations of exact closed cell selections and admitted planar
-- regions. Euler characteristic and area remain exact; Euclidean length is an
-- exact radical expression accompanied by outward-rounded binary64 bounds.
module Moonlight.Triangulation.Valuation
( EulerCharacteristic
, eulerCharacteristicValue
, ExactArea
, exactAreaValue
, ExactLengthTerm
, lengthCoefficient
, squaredLength
, ExactLengthExpression
, exactLengthTerms
, CertifiedInterval (..)
, ExactLengthMeasurement
, exactLengthExpression
, exactLengthBounds
, PlanarValuations
, valuationEuler
, valuationArea
, valuationIntrinsic1
, ValuationError (..)
, cellValuations
, regionValuations
, cellSetPerimeter
, regionPerimeter
) where
import Control.DeepSeq (NFData)
import Data.Bifunctor (first)
import Data.Bits (shiftL)
import Data.Foldable (foldlM)
import qualified Data.List as List
import Data.List.NonEmpty (NonEmpty (..))
import qualified Data.IntMap.Strict as IntMap
import qualified Data.IntSet as IntSet
import qualified Data.Map.Strict as Map
import Data.Maybe (catMaybes)
import qualified Data.Ratio as Ratio
import qualified Data.Set as Set
import qualified Data.Vector as V
import GHC.Float (castDoubleToWord64, castWord64ToDouble)
import GHC.Generics (Generic)
import Moonlight.Triangulation.Dcel
( faceVertices
, incidentFace
, undirectedEndpoints
)
import Moonlight.Triangulation.Exact
( ExactGeometryError
, ExactPoint
, ExactSegment
, exactOnClosedSegment
, exactPointCross
, exactPointCoordinates
, exactSegment
, exactSegmentEndpoints
)
import Moonlight.Triangulation.Handles.HandleDefs
( FaceId (..)
, UndirectedEdgeId (..)
, VertexId (..)
, directedPair
)
import Moonlight.Triangulation.Internal.CellSet
( ExactCellSet (..)
, exactCellSetIsFaceClosure
)
import Moonlight.Triangulation.Internal.BoundaryCycle
( consecutivePairs
, cyclePairs
, orderedPair
)
import Moonlight.Triangulation.Internal.Dyadic (integerBitLength)
import Moonlight.Triangulation.Internal.ExactRational
( ExactRational
, exactRationalDenominator
, exactRationalFromDyadic
, exactRationalFromFiniteDouble
, exactRationalIsZero
, exactRationalNumerator
, exactSignum
)
import Moonlight.Triangulation.Internal.ExactSegmentEvents
( ExactSegmentEvent (..)
, ExactSegmentEventObstruction
, ExactSweepSegmentId (..)
, exactSegmentEventPlan
, exactSegmentEvents
, exactSegmentSplitPoints
)
import Moonlight.Triangulation.Internal.Region.Types
( ExactLoop (..)
, PlanarRegion (..)
, PolygonComponent (..)
)
import Moonlight.Triangulation.Internal.Region.Bounds
( ExactBounds
, componentBounds
, overlappingPredecessors
)
import Moonlight.Triangulation.Internal.Representation (Triangulation)
newtype EulerCharacteristic = EulerCharacteristic Int
deriving stock (Eq, Ord, Show, Generic)
deriving anyclass (NFData)
eulerCharacteristicValue :: EulerCharacteristic -> Int
eulerCharacteristicValue (EulerCharacteristic value) = value
newtype ExactArea = ExactArea ExactRational
deriving stock (Eq, Ord, Show, Generic)
deriving anyclass (NFData)
exactAreaValue :: ExactArea -> ExactRational
exactAreaValue (ExactArea value) = value
data ExactLengthTerm = ExactLengthTerm
{ lengthCoefficient :: !ExactRational
, squaredLength :: !ExactRational
}
deriving stock (Eq, Ord, Show, Generic)
deriving anyclass (NFData)
-- | A normalized sum of rational coefficients times square roots of rational
-- squared lengths. It intentionally has no 'Eq' instance: syntactic radical
-- normalization is not algebraic-number equality.
newtype ExactLengthExpression = ExactLengthExpression [ExactLengthTerm]
deriving stock (Show, Generic)
deriving anyclass (NFData)
exactLengthTerms :: ExactLengthExpression -> [ExactLengthTerm]
exactLengthTerms (ExactLengthExpression terms) = terms
data CertifiedInterval = CertifiedInterval
{ intervalLower :: !Double
, intervalUpper :: !Double
}
deriving stock (Eq, Ord, Show, Generic)
deriving anyclass (NFData)
data ExactLengthMeasurement = ExactLengthMeasurement
{ exactLengthExpression :: !ExactLengthExpression
, exactLengthBounds :: !CertifiedInterval
}
deriving stock (Show, Generic)
deriving anyclass (NFData)
data PlanarValuations = PlanarValuations
{ valuationEuler :: !EulerCharacteristic
, valuationArea :: !ExactArea
, valuationIntrinsic1 :: !ExactLengthMeasurement
}
deriving stock (Show, Generic)
deriving anyclass (NFData)
data ValuationError
= ValuationCoordinateMissing !VertexId
| ValuationFaceArity !FaceId !Int
| ValuationInvalidRegionSegment !ExactGeometryError
| ValuationSegmentEventsInvalid !ExactSegmentEventObstruction
| ValuationSegmentMissing !ExactSweepSegmentId
| ValuationBoundaryMultiplicity !ExactPoint !ExactPoint !Int
| ValuationNegativeSquaredLength !ExactRational
| ValuationCellSetNotPureRegion
deriving stock (Eq, Show, Generic)
deriving anyclass (NFData)
cellValuations :: ExactCellSet -> Either ValuationError PlanarValuations
cellValuations (ExactCellSet triangulation points selectedEdges selectedFaces) = do
faceDoubleAreas <-
traverse
(cellFaceDoubleArea triangulation points . FaceId . fromIntegral)
(IntSet.toAscList selectedFaces)
edgeContributions <-
traverse
( cellEdgeLengthContribution triangulation points selectedFaces
. UndirectedEdgeId
. fromIntegral
)
(IntSet.toAscList selectedEdges)
assembleValuations
(IntMap.size points - IntSet.size selectedEdges + IntSet.size selectedFaces)
(List.foldl' (+) 0 faceDoubleAreas)
(normalizeLengthContributions id edgeContributions)
regionValuations :: PlanarRegion -> Either ValuationError PlanarValuations
regionValuations (PlanarRegion components) = do
componentBoundaries <- traverse componentBoundaryData components
let boundaryCover =
overlappingPredecessors componentBoundaryBounds componentBoundaries
hasPotentialBoundaryContacts = any (not . null . snd) boundaryCover
euler <- regionEuler boundaryCover
boundaryAtoms <-
normalizedRegionBoundaryAtoms
hasPotentialBoundaryContacts
componentBoundaries
let doubleArea =
List.foldl'
(\area component -> area + componentDoubleArea component)
0
components
assembleValuations
euler
doubleArea
( normalizeLengthContributions
(\(from, to) -> (oneHalf, segmentSquaredLength from to))
boundaryAtoms
)
assembleValuations
:: Int
-> ExactRational
-> ExactLengthExpression
-> Either ValuationError PlanarValuations
assembleValuations euler doubleArea lengthExpression =
PlanarValuations (EulerCharacteristic euler) (ExactArea (oneHalf * doubleArea))
<$> measureLength lengthExpression
cellSetPerimeter
:: ExactCellSet
-> Either ValuationError ExactLengthMeasurement
cellSetPerimeter cellSet
| exactCellSetIsFaceClosure cellSet =
cellValuations cellSet >>= conventionalPerimeter
| otherwise = Left ValuationCellSetNotPureRegion
regionPerimeter
:: PlanarRegion
-> Either ValuationError ExactLengthMeasurement
regionPerimeter region = regionValuations region >>= conventionalPerimeter
conventionalPerimeter
:: PlanarValuations
-> Either ValuationError ExactLengthMeasurement
conventionalPerimeter valuations =
measureLength
(scaleLengthExpression 2 (exactLengthExpression (valuationIntrinsic1 valuations)))
cellFaceDoubleArea
:: Triangulation mode vertex directed undirected face
-> IntMap.IntMap ExactPoint
-> FaceId
-> Either ValuationError ExactRational
cellFaceDoubleArea triangulation points face =
case faceVertices triangulation face of
[firstVertex, secondVertex, thirdVertex] -> do
firstPoint <- cellPoint points firstVertex
secondPoint <- cellPoint points secondVertex
thirdPoint <- cellPoint points thirdVertex
pure (triangleDoubleArea firstPoint secondPoint thirdPoint)
vertices -> Left (ValuationFaceArity face (length vertices))
cellEdgeLengthContribution
:: Triangulation mode vertex directed undirected face
-> IntMap.IntMap ExactPoint
-> IntSet.IntSet
-> UndirectedEdgeId
-> Either ValuationError (ExactRational, ExactRational)
cellEdgeLengthContribution triangulation points selectedFaces edge = do
let (fromVertex, toVertex) = undirectedEndpoints triangulation edge
(forward, backward) = directedPair edge
selected face =
let FaceId raw = face
in IntSet.member (fromIntegral raw) selectedFaces
coefficient = case (selected (incidentFace triangulation forward), selected (incidentFace triangulation backward)) of
(False, False) -> 1
(True, True) -> 0
_ -> oneHalf
from <- cellPoint points fromVertex
to <- cellPoint points toVertex
pure (coefficient, segmentSquaredLength from to)
cellPoint
:: IntMap.IntMap ExactPoint
-> VertexId
-> Either ValuationError ExactPoint
cellPoint points vertex@(VertexId raw) =
maybe
(Left (ValuationCoordinateMissing vertex))
Right
(IntMap.lookup (fromIntegral raw) points)
triangleDoubleArea :: ExactPoint -> ExactPoint -> ExactPoint -> ExactRational
triangleDoubleArea firstPoint secondPoint thirdPoint =
exactPointCross firstPoint secondPoint
+ exactPointCross secondPoint thirdPoint
+ exactPointCross thirdPoint firstPoint
componentDoubleArea :: PolygonComponent -> ExactRational
componentDoubleArea component =
List.foldl'
(\area loop -> area + loopDoubleArea loop)
0
(polygonOuterLoop component : polygonHoleLoops component)
loopDoubleArea :: ExactLoop -> ExactRational
loopDoubleArea (ExactLoop points) =
List.foldl'
(\area (from, to) -> area + exactPointCross from to)
0
(cyclePairs points)
segmentSquaredLength :: ExactPoint -> ExactPoint -> ExactRational
segmentSquaredLength from to =
let (fromX, fromY) = exactPointCoordinates from
(toX, toY) = exactPointCoordinates to
deltaX = toX - fromX
deltaY = toY - fromY
in deltaX * deltaX + deltaY * deltaY
normalizeLengthContributions
:: Foldable collection
=> (value -> (ExactRational, ExactRational))
-> collection value
-> ExactLengthExpression
normalizeLengthContributions contribution contributions =
ExactLengthExpression
[ ExactLengthTerm coefficient square
| (square, coefficient) <- Map.toAscList coefficientsBySquare
, not (exactRationalIsZero coefficient)
]
where
coefficientsBySquare =
List.foldl' accumulateContribution Map.empty contributions
accumulateContribution coefficients value =
case contribution value of
(coefficient, square)
| exactRationalIsZero coefficient -> coefficients
| otherwise -> Map.insertWith (+) square coefficient coefficients
scaleLengthExpression
:: Integer
-> ExactLengthExpression
-> ExactLengthExpression
scaleLengthExpression scalar (ExactLengthExpression terms) =
let exactScalar = fromInteger scalar
in ExactLengthExpression
[ term
{ lengthCoefficient =
exactScalar * lengthCoefficient term
}
| term <- terms
]
measureLength
:: ExactLengthExpression
-> Either ValuationError ExactLengthMeasurement
measureLength expression@(ExactLengthExpression terms) = do
(lower, upper) <-
foldlM
addTermBounds
(0, 0)
terms
pure
ExactLengthMeasurement
{ exactLengthExpression = expression
, exactLengthBounds =
CertifiedInterval
{ intervalLower = directedLowerDouble lower
, intervalUpper = directedUpperDouble upper
}
}
where
addTermBounds (lowerTotal, upperTotal) term = do
(lowerRoot, upperRoot) <- exactSquareRootBounds (squaredLength term)
let coefficient = lengthCoefficient term
pure
( lowerTotal + coefficient * lowerRoot
, upperTotal + coefficient * upperRoot
)
exactSquareRootBounds
:: ExactRational
-> Either ValuationError (ExactRational, ExactRational)
exactSquareRootBounds value =
case exactSignum value of
LT -> Left (ValuationNegativeSquaredLength value)
_ ->
let numerator = exactRationalNumerator value
denominator = exactRationalDenominator value
scale = 1 `shiftL` radicalPrecisionBits
scaledNumerator = numerator * scale * scale
root = integerSquareRoot (scaledNumerator `div` denominator)
exact = root * root * denominator == scaledNumerator
dyadicPower = negate radicalPrecisionBits
in Right
( exactRationalFromDyadic root dyadicPower
, exactRationalFromDyadic (if exact then root else root + 1) dyadicPower
)
radicalPrecisionBits :: Int
radicalPrecisionBits = 128
integerSquareRoot :: Integer -> Integer
integerSquareRoot value
| value < 2 = value
| otherwise = descend initial
where
initial = 1 `shiftL` ((integerBitLength value + 1) `div` 2)
descend estimate =
let refined = (estimate + value `div` estimate) `div` 2
in if refined >= estimate then estimate else descend refined
directedLowerDouble :: ExactRational -> Double
directedLowerDouble value =
let candidate = rationalToDouble value
in if isInfinite candidate
then maximumFiniteDouble
else
if exactRationalFromFiniteDouble candidate <= value
then candidate
else previousPositiveDouble candidate
directedUpperDouble :: ExactRational -> Double
directedUpperDouble value =
let candidate = rationalToDouble value
in if isInfinite candidate
|| exactRationalFromFiniteDouble candidate >= value
then candidate
else nextPositiveDouble candidate
rationalToDouble :: ExactRational -> Double
rationalToDouble value =
fromRational
( exactRationalNumerator value
Ratio.% exactRationalDenominator value
)
previousPositiveDouble :: Double -> Double
previousPositiveDouble value
| value <= 0 = 0
| otherwise = castWord64ToDouble (castDoubleToWord64 value - 1)
nextPositiveDouble :: Double -> Double
nextPositiveDouble value
| value == 0 = castWord64ToDouble 1
| otherwise = castWord64ToDouble (castDoubleToWord64 value + 1)
maximumFiniteDouble :: Double
maximumFiniteDouble = castWord64ToDouble 0x7fefffffffffffff
data ComponentBoundaryData = ComponentBoundaryData
{ componentBoundaryEuler :: !Int
, componentBoundaryBounds :: !ExactBounds
, componentBoundarySegments :: !(V.Vector ExactSegment)
}
componentBoundaryData
:: PolygonComponent
-> Either ValuationError ComponentBoundaryData
componentBoundaryData component = do
segments <-
V.fromList
<$> traverse
admittedSegment
( concatMap
(cyclePairs . loopPoints)
(polygonOuterLoop component : polygonHoleLoops component)
)
pure
ComponentBoundaryData
{ componentBoundaryEuler = 1 - length (polygonHoleLoops component)
, componentBoundaryBounds = componentBounds component
, componentBoundarySegments = segments
}
regionEuler
:: [(ComponentBoundaryData, [ComponentBoundaryData])]
-> Either ValuationError Int
regionEuler = foldlM attachComponent 0
where
attachComponent accumulatedEuler (current, priorCandidates) = do
let priorSegments = V.concat (map componentBoundarySegments priorCandidates)
intersectionEuler <-
boundaryIntersectionEuler
(componentBoundarySegments current)
priorSegments
pure
( accumulatedEuler
+ componentBoundaryEuler current
- intersectionEuler
)
boundaryIntersectionEuler
:: V.Vector ExactSegment
-> V.Vector ExactSegment
-> Either ValuationError Int
boundaryIntersectionEuler current prior
| V.null prior = Right 0
| otherwise = do
let currentCount = V.length current
segments = current <> prior
plan <- first ValuationSegmentEventsInvalid (exactSegmentEventPlan segments)
contacts <-
traverse
(contactFromEvent segments)
[ event
| event <- exactSegmentEvents plan
, crossPartition currentCount event
]
let contactPoints =
Set.fromList
[ point
| ContactPoint point <- contacts
]
intervals =
[ interval
| ContactInterval interval <- contacts
]
allSplitPoints =
Set.fromList
( concatMap
(exactSegmentSplitPoints plan)
[ ExactSweepSegmentId index
| index <- [0 .. V.length segments - 1]
]
)
contactEdges =
Set.fromList
[ orderedPair from to
| interval <- intervals
, let (lower, upper) = interval
points =
Set.toAscList
( Set.filter
(exactOnClosedSegment lower upper)
allSplitPoints
)
, (from, to) <- consecutivePairs points
, from /= to
]
vertices =
Set.unions
[ contactPoints
, Set.fromList
[ point
| (from, to) <- Set.toAscList contactEdges
, point <- [from, to]
]
]
pure (Set.size vertices - Set.size contactEdges)
data BoundaryContact
= ContactPoint !ExactPoint
| ContactInterval !(ExactPoint, ExactPoint)
contactFromEvent
:: V.Vector ExactSegment
-> ExactSegmentEvent
-> Either ValuationError BoundaryContact
contactFromEvent _ (ExactProperCrossing _ _ point) = Right (ContactPoint point)
contactFromEvent _ (ExactEndpointTouch _ _ point) = Right (ContactPoint point)
contactFromEvent _ (ExactSharedEndpoint _ _ point) = Right (ContactPoint point)
contactFromEvent segments (ExactDuplicateSegments leftId _) =
ContactInterval . canonicalSegmentEndpoints
<$> requireSegment segments leftId
contactFromEvent _ (ExactCollinearOverlap _ _ lower upper) =
Right (ContactInterval (orderedPair lower upper))
crossPartition :: Int -> ExactSegmentEvent -> Bool
crossPartition boundary event =
let (ExactSweepSegmentId left, ExactSweepSegmentId right) = eventIds event
in (left < boundary) /= (right < boundary)
eventIds
:: ExactSegmentEvent
-> (ExactSweepSegmentId, ExactSweepSegmentId)
eventIds (ExactProperCrossing left right _) = (left, right)
eventIds (ExactEndpointTouch left right _) = (left, right)
eventIds (ExactSharedEndpoint left right _) = (left, right)
eventIds (ExactDuplicateSegments left right) = (left, right)
eventIds (ExactCollinearOverlap left right _ _) = (left, right)
normalizedRegionBoundaryAtoms
:: Bool
-> [ComponentBoundaryData]
-> Either ValuationError (Set.Set (ExactPoint, ExactPoint))
normalizedRegionBoundaryAtoms hasPotentialBoundaryContacts boundaries
| V.null segments = Right Set.empty
| not hasPotentialBoundaryContacts =
Right
( Set.fromList
(map canonicalSegmentEndpoints (V.toList segments))
)
| otherwise = do
plan <- first ValuationSegmentEventsInvalid (exactSegmentEventPlan segments)
let orientedAtoms =
concatMap
segmentAtoms
[ exactSegmentSplitPoints plan (ExactSweepSegmentId index)
| index <- [0 .. V.length segments - 1]
]
traverseMultiplicity
(Map.toAscList (Map.fromListWith (+) orientedAtoms))
where
segments = V.concat (map componentBoundarySegments boundaries)
segmentAtoms :: [ExactPoint] -> [((ExactPoint, ExactPoint), Int)]
segmentAtoms points =
[ ( orderedPair firstPoint secondPoint
, if firstPoint <= secondPoint then 1 else -1
)
| (firstPoint, secondPoint) <- consecutivePairs points
, firstPoint /= secondPoint
]
traverseMultiplicity
:: [((ExactPoint, ExactPoint), Int)]
-> Either ValuationError (Set.Set (ExactPoint, ExactPoint))
traverseMultiplicity entries = do
retained <-
traverse
(\(edge@(from, to), multiplicity) ->
case abs multiplicity of
0 -> Right Nothing
1 -> Right (Just edge)
_ -> Left (ValuationBoundaryMultiplicity from to multiplicity))
entries
pure (Set.fromList (catMaybes retained))
admittedSegment
:: (ExactPoint, ExactPoint)
-> Either ValuationError ExactSegment
admittedSegment (from, to) = first ValuationInvalidRegionSegment (exactSegment from to)
requireSegment
:: V.Vector ExactSegment
-> ExactSweepSegmentId
-> Either ValuationError ExactSegment
requireSegment segments segmentId@(ExactSweepSegmentId index) =
maybe
(Left (ValuationSegmentMissing segmentId))
Right
(segments V.!? index)
canonicalSegmentEndpoints :: ExactSegment -> (ExactPoint, ExactPoint)
canonicalSegmentEndpoints = uncurry orderedPair . exactSegmentEndpoints
loopPoints :: ExactLoop -> NonEmpty ExactPoint
loopPoints (ExactLoop points) = points
oneHalf :: ExactRational
oneHalf = exactRationalFromDyadic 1 (-1)