Orientation.hs

{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Math.Orientation
  ( Orientation(..), reverseOrientation
  , convexOrientation, splineOrientation, splineTangents
  , between
  )
  where

-- base
import Control.Monad
  ( guard )
import Data.Fixed
  ( mod' )

-- acts
import Data.Act
  ( Torsor((-->)) )

-- containers
import Data.Sequence
  ( Seq(..) )

-- generic-lens
import Data.Generics.Product.Typed
  ( HasType(typed) )
import Data.Generics.Internal.VL
  ( view )

-- brush-strokes
import Math.Epsilon
  ( nearZero )
import Math.Module
  ( (×) )
import Math.Bezier.Spline
  ( Spline(..), Curves(..), Curve(..), NextPoint(..)
  , SplineType(..), KnownSplineType(..), SSplineType(..)
  , ssplineType
  )
import Math.Linear
  ( ℝ(..), T(..) )

--------------------------------------------------------------------------------

-- | An orientation in the plane: counter-clockwise or clockwise.
data Orientation = CCW | CW
  deriving stock ( Show, Eq, Ord )

-- | Reverse an orientation, turning counter-clockwise into clockwise and vice-versa.
reverseOrientation :: Orientation -> Orientation
reverseOrientation CCW = CW
reverseOrientation CW  = CCW

-- | Compute an orientation from a sequence of tangent vectors (assumed to have monotone angle).
convexOrientation :: [ T ( ℝ 2 ) ] -> Orientation
convexOrientation ( v1 : v2 : vs )
  | nearZero crossProduct
  = convexOrientation ( v2 : vs )
  | crossProduct < 0
  = CW
  | otherwise
  = CCW
  where
    crossProduct :: Double
    crossProduct = v1 × v2
convexOrientation _ = CCW -- default

-- | Compute the orientation of a spline, assuming tangent vectors have a monotone angle.
splineOrientation
  :: forall clo crvData ptData
  .  ( KnownSplineType clo, HasType ( ℝ 2 ) ptData )
  => Spline clo crvData ptData
  -> Orientation
splineOrientation = convexOrientation . splineTangents

-- | Compute the sequence of tangent vectors given by the control points of a Bézier spline.
splineTangents
  :: forall clo crvData ptData
  .  ( KnownSplineType clo, HasType ( ℝ 2 ) ptData )
  => Spline clo crvData ptData
  -> [ T ( ℝ 2 ) ]
splineTangents spline@( Spline { splineStart = sp0, splineCurves = curves } )
  | let
      p0 :: ℝ 2
      p0 = view typed sp0
  = case ssplineType @clo of
      SOpen
        | OpenCurves { openCurves = cs } <- curves
        -> go p0 cs
      SClosed
        | OpenCurves cs@( c :<| _ ) <- splineCurves $ adjustSplineType @Open spline
        -> go p0 ( cs :|> c )
      _ -> []
  where
    go :: ℝ 2 -> Seq ( Curve Open crvData ptData ) -> [ T ( ℝ 2 ) ]
    go _ Empty = []
    go p ( crv :<| crvs ) =
      case crv of
        LineTo { curveEnd = NextPoint sq }
          | let
              q :: ℝ 2
              q = view typed sq
          -> ( p --> q ) : go q crvs
        Bezier2To { controlPoint = scp, curveEnd = NextPoint sq }
          | let
              cp, q :: ℝ 2
              cp = view typed scp
              q  = view typed sq
          -> ( p --> cp ) : ( cp --> q ) : go q crvs
        Bezier3To { controlPoint1 = scp1, controlPoint2 = scp2, curveEnd = NextPoint sq }
          | let
              cp1, cp2, q :: ℝ 2
              cp1 = view typed scp1
              cp2 = view typed scp2
              q   = view typed sq
          -> ( p --> cp1 ) : ( cp1 --> cp2 ) : ( cp2 --> q ) : go q crvs

-- | Checks whether a 2D vector lies "in between" two other vectors according to a given orientation,
-- i.e. whether the angle of the query vector lies in between the angles of the start and end vectors.
--
-- The input vectors are expected to be within π radians of each other (as we are dealing
-- with convex shapes).
--
-- Returns the proportion of the angle the vector is in between, or @Nothing@ if the query vector
-- is not in between.
--
-- >>> between CCW ( V2 1 0 ) ( V2 -1 1 ) ( V2 1 1 )
-- Just 0.3333333333333333
between
  :: Orientation
  -> T ( ℝ 2 ) -- ^ start vector
  -> T ( ℝ 2 ) -- ^ end vector
  -> T ( ℝ 2 ) -- ^ query vector: is it in between the start and end vectors w.r.t. the provided orientation?
  -> Maybe Double
between CCW ( V2 x1 y1 ) ( V2 x2 y2 ) ( V2 a b ) =
  let
    τ, η, φ, θ :: Double
    τ = 2 * pi
    η = atan2 y1 x1
    φ = ( atan2 y2 x2 - η ) `mod'` τ
    θ = ( atan2 b  a  - η ) `mod'` τ

  in do
    guard ( φ < pi ) -- deals with possible floating-point inaccuracies,
                     -- which can cause the angles of the two reference vectors
                     -- to be flipped, e.g. arg(v2) = arg(v1) - 0.000000001
    guard ( θ <= φ )
    pure  ( θ / φ )
between CW v1 v2 u = ( 1 - ) <$> between CCW v2 v1 u