Math.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DerivingStrategies #-}

-- | Robust planar predicates and derived Euclidean constructions.
module Moonlight.Triangulation.Math
  ( orient2d
  , sideQuery
  , inCircle
  , orientDetApprox
  , inCircleDetApprox
  , onClosedSegment
  , SegmentRelation (..)
  , allSegmentRelations
  , segmentRelation
  , segmentsProperlyCross
  , segmentsIntersect
  , squaredDistance
  , squaredDistanceWide
  , segmentDistanceSquared
  , segmentDistanceSquaredWide
  , distance
  , midpoint
  , centroid
  , triangleArea
  , triangleRadiusEdgeRatio
  , triangleRadiusEdgeRatioSquaredWithArea
  , circumcenter
  , barycentricCoordinates
  , inDiametralCircle
  , projectionFactor
  , canonicalPoint
  , canonicalCoordinate
  , validateCoordinate
  , mkQueryPoint
  , validatePoint
  , mitigateUnderflow
  , isFinite
  ) where

import Control.DeepSeq (NFData)
import GHC.Generics (Generic)
import Moonlight.Triangulation.Internal.Dyadic
  ( exactBarycentricDeterminants
  , exactDiametralDot
  , integerRatioToDouble
  )
import Moonlight.Triangulation.LineSideInfo (LineSideInfo, fromOrdering)
import Moonlight.Triangulation.Scalar
  ( canonicalScalarZero
  , inCircleCoordinates
  , maximumAllowedCoordinate
  , minimumAllowedCoordinate
  , orient2dCoordinates
  , scalarCcwErrorBound
  )
import Moonlight.Triangulation.Internal.Types
  ( BuildError (..)
  , CoordinateError (..)
  , Point (..)
  , PointValidationError (..)
  , QueryPoint (..)
  )

-- | The complete exact-predicate relation between two closed segments. There
-- is one vocabulary owner; traversal and constrained-union consumers derive
-- their booleans and obstruction policy from it rather than cloning slightly
-- different orientation formulae.
data SegmentRelation
  = SegmentsDisjoint
  | SegmentsDuplicate
  | SegmentsShareEndpoint
  | SegmentsProperlyCross
  | SegmentEndpointTouchesInterior
  | SegmentsCollinearlyOverlap
  deriving stock (Bounded, Enum, Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | Every segment relation in constructor order.
allSegmentRelations :: [SegmentRelation]
allSegmentRelations = [minBound .. maxBound]

-- | Exact relation between two closed segments.
segmentRelation
  :: Point
  -> Point
  -> Point
  -> Point
  -> SegmentRelation
segmentRelation a b c d
  | sameUndirectedSegment = SegmentsDuplicate
  | sharesEndpoint = SegmentsShareEndpoint
  | opposite abC abD && opposite cdA cdB = SegmentsProperlyCross
  | abC == EQ && abD == EQ && cdA == EQ && cdB == EQ = collinearRelation
  | endpointTouches = SegmentEndpointTouchesInterior
  | otherwise = SegmentsDisjoint
 where
  !abC = orient2d a b c
  !abD = orient2d a b d
  !cdA = orient2d c d a
  !cdB = orient2d c d b
  sameUndirectedSegment = (a == c && b == d) || (a == d && b == c)
  sharesEndpoint = a == c || a == d || b == c || b == d
  endpointTouches =
    (abC == EQ && onClosedSegment a b c)
      || (abD == EQ && onClosedSegment a b d)
      || (cdA == EQ && onClosedSegment c d a)
      || (cdB == EQ && onClosedSegment c d b)
  collinearRelation =
    let !overlapLower = max (min a b) (min c d)
        !overlapUpper = min (max a b) (max c d)
     in case compare overlapLower overlapUpper of
          LT -> SegmentsCollinearlyOverlap
          EQ -> SegmentEndpointTouchesInterior
          GT -> SegmentsDisjoint
  opposite left right =
    (left == LT && right == GT) || (left == GT && right == LT)

-- | Whether two closed segments share any point.
segmentsIntersect
  :: Point
  -> Point
  -> Point
  -> Point
  -> Bool
segmentsIntersect a b c d = segmentRelation a b c d /= SegmentsDisjoint

-- | The proper-crossing section of 'segmentRelation'. Consumers which reject
-- only that constructor need not compute the collinear and endpoint-touch
-- distinctions required by the complete ADT after either side already proves
-- separation.
segmentsProperlyCross
  :: Point
  -> Point
  -> Point
  -> Point
  -> Bool
segmentsProperlyCross a b c d =
  opposite (orient2d a b c) (orient2d a b d)
    && opposite (orient2d c d a) (orient2d c d b)
 where
  opposite LT GT = True
  opposite GT LT = True
  opposite _ _ = False

-- | Classify a coordinate outside the exact-predicate input domain.
validateCoordinate :: Double -> Maybe CoordinateError
validateCoordinate value
  | isNaN value = Just CoordinateNaN
  | isInfinite value = Just CoordinateInfinite
  | value /= 0 && abs value < minimumAllowedCoordinate = Just CoordinateTooSmall
  | abs value > maximumAllowedCoordinate = Just CoordinateTooLarge
  | otherwise = Nothing

-- | Admit and normalize a finite point for read-only geometric queries.
mkQueryPoint :: Point -> Either PointValidationError (QueryPoint)
mkQueryPoint point@(Point x y) = do
  maybe (Right ()) (Left . InvalidPointX) (validateCoordinate x)
  maybe (Right ()) (Left . InvalidPointY) (validateCoordinate y)
  Right (QueryPoint (canonicalPoint point))

-- | Validate a construction point while retaining its optional input slot.
validatePoint :: Maybe Int -> Point -> Either BuildError (QueryPoint)
validatePoint slot point@(Point x y) =
  case mkQueryPoint point of
    Left (InvalidPointX reason) -> Left (InvalidCoordinate slot x reason)
    Left (InvalidPointY reason) -> Left (InvalidCoordinate slot y reason)
    Right queryPoint -> Right queryPoint

-- | Round coordinates below the robust-predicate input floor toward zero.
-- The operation never changes a coordinate already accepted by
-- 'validateCoordinate'.
mitigateUnderflow :: Point -> Point
mitigateUnderflow (Point x y) = Point (mitigate x) (mitigate y)
 where
  mitigate :: Double -> Double
  mitigate value
    | value /= 0 && abs value < minimumAllowedCoordinate = 0
    | otherwise = value

-- | Canonicalize both coordinate components for point identity.
canonicalPoint :: Point -> Point
canonicalPoint (Point x y) = Point (canonicalCoordinate x) (canonicalCoordinate y)
{-# INLINE canonicalPoint #-}

-- | Round a signed zero to the canonical zero. The law that makes two points
-- at the same position compare equal lives here; @canonicalPoint@ is its
-- component-wise form and coordinate-carrying callers use it directly so no
-- t'Point' is built only to be taken apart again.
canonicalCoordinate :: Double -> Double
canonicalCoordinate = canonicalScalarZero
{-# INLINE canonicalCoordinate #-}

-- base's 'isNaN' and 'isInfinite' are FFI calls in this GHC (static
-- ccall to isDoubleNaN/isDoubleInfinite), and a circumcentre pays four of
-- them per call. @value - value == 0@ is the same predicate in pure Double
-- arithmetic: finite values subtract to zero, while NaN and the infinities
-- subtract to NaN, which never compares equal.
-- | Whether a scalar is neither infinite nor NaN.
isFinite :: Double -> Bool
isFinite value = value - value == 0
{-# INLINE isFinite #-}

-- | Fast approximate signed orientation determinant.
orientDetApprox :: Point -> Point -> Point -> Double
orientDetApprox (Point ax ay) (Point bx by) (Point cx cy) =
  (ax - cx) * (by - cy) - (ay - cy) * (bx - cx)
{-# INLINE orientDetApprox #-}

-- | Exact orientation ordering of three points.
orient2d :: Point -> Point -> Point -> Ordering
orient2d (Point ax ay) (Point bx by) (Point cx cy) =
  orient2dCoordinates ax ay bx by cx cy
{-# INLINE orient2d #-}

-- | Exact side of an oriented line.
sideQuery :: Point -> Point -> Point -> LineSideInfo
sideQuery a b point = fromOrdering (orient2d a b point)
{-# INLINE sideQuery #-}

-- | Fast approximate oriented in-circle determinant.
inCircleDetApprox
  :: Point -> Point -> Point -> Point -> Double
inCircleDetApprox
  (Point ax ay)
  (Point bx by)
  (Point cx cy)
  (Point dx dy) =
    alift * bcdet + blift * cadet + clift * abdet
 where
  !adx = ax - dx
  !ady = ay - dy
  !bdx = bx - dx
  !bdy = by - dy
  !cdx = cx - dx
  !cdy = cy - dy
  !abdet = adx * bdy - bdx * ady
  !bcdet = bdx * cdy - cdx * bdy
  !cadet = cdx * ady - adx * cdy
  !alift = adx * adx + ady * ady
  !blift = bdx * bdx + bdy * bdy
  !clift = cdx * cdx + cdy * cdy
{-# INLINE inCircleDetApprox #-}

-- | Ordering of the oriented incircle determinant. For a counter-clockwise
-- triangle, 'GT' means the fourth point lies strictly inside its circumcircle.
inCircle
  :: Point -> Point -> Point -> Point -> Ordering
inCircle
  (Point ax ay)
  (Point bx by)
  (Point cx cy)
  (Point dx dy) =
    inCircleCoordinates ax ay bx by cx cy dx dy
{-# INLINE inCircle #-}

-- | Whether a point lies on a closed segment.
onClosedSegment :: Point -> Point -> Point -> Bool
onClosedSegment a@(Point ax ay) b@(Point bx by) query@(Point qx qy) =
  orient2d a b query == EQ
    && qx >= min ax bx
    && qx <= max ax bx
    && qy >= min ay by
    && qy <= max ay by
{-# INLINE onClosedSegment #-}

-- | Squared Euclidean distance.
squaredDistance :: Point -> Point -> Double
squaredDistance (Point ax ay) (Point bx by) =
  let !dx = ax - bx
      !dy = ay - by
   in dx * dx + dy * dy
{-# INLINE squaredDistance #-}

-- | Squared Euclidean distance in the mesh's Binary64 coordinate domain.
squaredDistanceWide :: Point -> Point -> Double
squaredDistanceWide = squaredDistance
{-# INLINE squaredDistanceWide #-}

-- | Squared distance from a point to a closed segment.
segmentDistanceSquared
  :: Point -> Point -> Point -> Double
segmentDistanceSquared from@(Point ax ay) to@(Point bx by) point@(Point px py)
  | lengthSquared == 0 = squaredDistance from point
  | factor <= 0 = squaredDistance from point
  | factor >= 1 = squaredDistance to point
  | otherwise = squaredDistance point (Point (ax + factor * dx) (ay + factor * dy))
 where
  !dx = bx - ax
  !dy = by - ay
  !lengthSquared = dx * dx + dy * dy
  !factor = ((px - ax) * dx + (py - ay) * dy) / lengthSquared
{-# INLINE segmentDistanceSquared #-}

-- | Comparison form retained beside 'segmentDistanceSquared' for callers that
-- state metric intent explicitly.
segmentDistanceSquaredWide
  :: Point -> Point -> Point -> Double
segmentDistanceSquaredWide = segmentDistanceSquared
{-# INLINE segmentDistanceSquaredWide #-}

-- | Euclidean distance.
distance :: Point -> Point -> Double
distance left right = sqrt (squaredDistance left right)
{-# INLINE distance #-}

-- | Midpoint of two points.
midpoint :: Point -> Point -> Point
midpoint (Point ax ay) (Point bx by) = Point (0.5 * ax + 0.5 * bx) (0.5 * ay + 0.5 * by)
{-# INLINE midpoint #-}

-- | Centroid of three points.
centroid :: Point -> Point -> Point -> Point
centroid (Point ax ay) (Point bx by) (Point cx cy) =
  Point (ax + (bx - ax) / 3 + (cx - ax) / 3) (ay + (by - ay) / 3 + (cy - ay) / 3)
{-# INLINE centroid #-}

-- | Unsigned area of a triangle.
triangleArea :: Point -> Point -> Point -> Double
triangleArea a b c = 0.5 * abs (orientDetApprox a b c)
{-# INLINE triangleArea #-}

-- | Circumradius divided by shortest edge length.
triangleRadiusEdgeRatio
  :: Point -> Point -> Point -> Maybe Double
triangleRadiusEdgeRatio p0 p1 p2
  | area <= 0 || shortest <= 0 = Nothing
  | not (isFinite ratio) = Nothing
  | otherwise = Just ratio
 where
  !area = triangleArea p0 p1 p2
  !side01 = distance p0 p1
  !side12 = distance p1 p2
  !side20 = distance p2 p0
  !shortest = min side01 (min side12 side20)
  !otherProduct
    | side01 <= side12 && side01 <= side20 = side12 * side20
    | side12 <= side20 = side20 * side01
    | otherwise = side01 * side12
  !ratio = otherProduct / (4 * area)

-- | The square of 'triangleRadiusEdgeRatio', for a triangle whose area the
-- caller already has.
--
-- The ratio is only ever compared against a bound, and both sides are
-- non-negative, so the comparison can be made between squares. That is the
-- whole reason to have this: it settles the same question without the three
-- square roots the lengths would need, on the path taken by every face
-- refinement considers.
--
-- The area is a parameter and a degenerate triangle answers with an infinity
-- rather than an absence, because the caller on that path has already computed
-- the area to ask the area question and does nothing with the absence but
-- compare an infinity in its place.
triangleRadiusEdgeRatioSquaredWithArea
  :: Double -> Point -> Point -> Point -> Double
triangleRadiusEdgeRatioSquaredWithArea area p0 p1 p2
  | area <= 0 || shortest <= 0 = 1 / 0
  | not (isFinite ratio) = 1 / 0
  | otherwise = ratio
 where
  !side01 = squaredDistance p0 p1
  !side12 = squaredDistance p1 p2
  !side20 = squaredDistance p2 p0
  !shortest = min side01 (min side12 side20)
  !otherProduct
    | side01 <= side12 && side01 <= side20 = side12 * side20
    | side12 <= side20 = side20 * side01
    | otherwise = side01 * side12
  !ratio = otherProduct / (16 * area * area)

-- The scale the determinants are divided by cancels out of the quotient
-- exactly, so the computation works on the unscaled differences and divides
-- once per coordinate. Scaling would only matter against overflow, and the
-- validated coordinate domain (|x| <= 3.3e60) keeps every intermediate below
-- 1e183, five orders below the Double ceiling; the four divisions it cost
-- are the circumcentre's hot-path price. Identical points answer through the
-- denominator, which is exactly zero exactly when they are collinear.
-- | Circumcenter of a nondegenerate triangle.
circumcenter
  :: Point -> Point -> Point -> Maybe (Point)
circumcenter (Point ax ay) (Point bx by) (Point cx cy)
  | denominator == 0 = Nothing
  | not (isFinite resultX && isFinite resultY) = Nothing
  | otherwise = Just (canonicalPoint (Point resultX resultY))
 where
  !bax = bx - ax
  !bay = by - ay
  !cax = cx - ax
  !cay = cy - ay
  !bLength = bax * bax + bay * bay
  !cLength = cax * cax + cay * cay
  !denominator = 2 * (bax * cay - bay * cax)
  !offsetX = (cay * bLength - bay * cLength) / denominator
  !offsetY = (bax * cLength - cax * bLength) / denominator
  !resultX = ax + offsetX
  !resultY = ay + offsetY

-- | Barycentric coordinates of a point in a nondegenerate triangle.
barycentricCoordinates
  :: Point -> Point -> Point -> Point
  -> Maybe (Double, Double, Double)
barycentricCoordinates a@(Point ax ay) b@(Point bx by) c@(Point cx cy) query@(Point qx qy)
  | all reliable [denominatorInfo, weightAInfo, weightBInfo, weightCInfo] =
      if denominator == 0
        then Nothing
        else Just (weightA / denominator, weightB / denominator, weightC / denominator)
  | exactDenominator == 0 = Nothing
  | otherwise =
      Just
        ( integerRatioToDouble exactWeightA exactDenominator
        , integerRatioToDouble exactWeightB exactDenominator
        , integerRatioToDouble exactWeightC exactDenominator
        )
 where
  !denominatorInfo@(denominator, _) = determinantInfo a b c
  !weightAInfo@(weightA, _) = determinantInfo query b c
  !weightBInfo@(weightB, _) = determinantInfo a query c
  !weightCInfo@(weightC, _) = determinantInfo a b query
  (!exactDenominator, !exactWeightA, !exactWeightB, !exactWeightC) =
    exactBarycentricDeterminants ax ay bx by cx cy qx qy

  reliable (determinant, determinantSum) =
    isFinite determinant && abs determinant > scalarCcwErrorBound * determinantSum

  determinantInfo :: Point -> Point -> Point -> (Double, Double)
  determinantInfo (Point px py) (Point rx ry) (Point sx sy) =
    let !left = (px - sx) * (ry - sy)
        !right = (py - sy) * (rx - sx)
     in (left - right, abs left + abs right)

-- | Whether a point lies in the closed diametral disk of a segment, decided by
-- the sign of @(a-p)·(b-p)@.
--
-- The approximation is two coordinate differences and one product per term
-- combined additively, which is the arithmetic shape 'orient2d' is analysed
-- under: the forward error of @fl(fl(a⊖b) ⊗ fl(c⊖d))@ combined by one rounded
-- addition does not depend on the sign of that combination. The orientation
-- coefficient therefore transfers unchanged, and the exact dot product remains
-- the oracle for the uncertain band.
inDiametralCircle :: Point -> Point -> Point -> Bool
inDiametralCircle (Point ax ay) (Point bx by) (Point px py)
  | isFinite dot && abs dot > scalarCcwErrorBound * dotSum = dot < 0
  | otherwise = exactDiametralDot ax ay bx by px py <= 0
 where
  !left = (ax - px) * (bx - px)
  !right = (ay - py) * (by - py)
  !dot = left + right
  !dotSum = abs left + abs right
{-# INLINE inDiametralCircle #-}

-- | Projection parameter of a point onto an oriented segment line.
projectionFactor :: Point -> Point -> Point -> Double
projectionFactor (Point ax ay) (Point bx by) (Point qx qy)
  | lengthSquared == 0 = 0
  | otherwise = ((qx - ax) * dx + (qy - ay) * dy) / lengthSquared
 where
  !dx = bx - ax
  !dy = by - ay
  !lengthSquared = dx * dx + dy * dy
{-# INLINE projectionFactor #-}