{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE ViewPatterns         #-}
{-# LANGUAGE UndecidableInstances #-}

{-# OPTIONS_GHC -fno-warn-orphans #-}
-- Orphan Traced instances for Segment Closed V2 and FixedSegment V2.
-- They can't go in Traced; but they shouldn't really go in
-- Diagrams.Segment either because we only have Traced instances for
-- the special case of R2.
-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.TwoD.Segment
-- Copyright   :  (c) 2012 diagrams-lib team (see LICENSE)
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  diagrams-discuss@googlegroups.com
--
-- Segments in two dimensions are special since we may meaningfully
-- compute their point of intersection with a ray.
--
-----------------------------------------------------------------------------

module Diagrams.TwoD.Segment
  ( -- * Segment intersections

    intersectPointsS
  , intersectPointsS'

    -- * Closest point on a segment

  , closestPoint
  , closestPoint'
  , closestDistance
  , closestDistance'
  , closestParam
  , closestParam'

    -- ** Low level functions
  , segmentSegment
  , lineSegment
  )
  where

import           Control.Lens                    hiding (at, contains, transform, ( # ))
import           Data.Maybe

import           Diagrams.Core

import           Diagrams.Direction
import           Diagrams.Located
import           Diagrams.Parametric
import           Diagrams.Segment
import           Diagrams.TwoD.Points
import           Diagrams.TwoD.Segment.Bernstein
import           Diagrams.TwoD.Transform
import           Diagrams.TwoD.Types             hiding (p2)
import           Diagrams.TwoD.Vector

import           Linear.Affine
import           Linear.Metric

{- All instances of Traced should maintain the invariant that the list of
   traces is sorted in increasing order.
-}

instance OrderedField n => Traced (Segment Closed V2 n) where
  getTrace :: Segment Closed V2 n
-> Trace (V (Segment Closed V2 n)) (N (Segment Closed V2 n))
getTrace = forall a. Traced a => a -> Trace (V a) (N a)
getTrace forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall n (v :: * -> *).
(Num n, Additive v) =>
Located (Segment Closed v n) -> FixedSegment v n
mkFixedSeg forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall a. a -> Point (V a) (N a) -> Located a
`at` forall (f :: * -> *) a. (Additive f, Num a) => Point f a
origin)

instance OrderedField n => Traced (FixedSegment V2 n) where
  getTrace :: FixedSegment V2 n
-> Trace (V (FixedSegment V2 n)) (N (FixedSegment V2 n))
getTrace FixedSegment V2 n
seg = forall (v :: * -> *) n.
(Point v n -> v n -> SortedList n) -> Trace v n
mkTrace forall a b. (a -> b) -> a -> b
$ \Point (V (FixedSegment V2 n)) (N (FixedSegment V2 n))
p V (FixedSegment V2 n) (N (FixedSegment V2 n))
v ->
    forall a. Ord a => [a] -> SortedList a
mkSortedList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view forall s t a b. Field1 s t a b => Lens s t a b
_1) forall a b. (a -> b) -> a -> b
$ forall n.
OrderedField n =>
n -> Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
lineSegment forall n. Fractional n => n
defEps (V (FixedSegment V2 n) (N (FixedSegment V2 n))
v forall a. a -> Point (V a) (N a) -> Located a
`at` Point (V (FixedSegment V2 n)) (N (FixedSegment V2 n))
p) FixedSegment V2 n
seg

defEps :: Fractional n => n
defEps :: forall n. Fractional n => n
defEps = n
1e-8

-- | Compute the intersections between two fixed segments.
intersectPointsS :: OrderedField n => FixedSegment V2 n -> FixedSegment V2 n -> [P2 n]
intersectPointsS :: forall n.
OrderedField n =>
FixedSegment V2 n -> FixedSegment V2 n -> [P2 n]
intersectPointsS = forall n.
OrderedField n =>
n -> FixedSegment V2 n -> FixedSegment V2 n -> [P2 n]
intersectPointsS' forall n. Fractional n => n
defEps

-- | Compute the intersections between two segments using the given tolerance.
intersectPointsS' :: OrderedField n => n -> FixedSegment V2 n -> FixedSegment V2 n -> [P2 n]
intersectPointsS' :: forall n.
OrderedField n =>
n -> FixedSegment V2 n -> FixedSegment V2 n -> [P2 n]
intersectPointsS' n
eps FixedSegment V2 n
s1 FixedSegment V2 n
s2 = forall a b. (a -> b) -> [a] -> [b]
map (forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view forall s t a b. Field3 s t a b => Lens s t a b
_3) forall a b. (a -> b) -> a -> b
$ forall n.
OrderedField n =>
n -> FixedSegment V2 n -> FixedSegment V2 n -> [(n, n, P2 n)]
segmentSegment n
eps FixedSegment V2 n
s1 FixedSegment V2 n
s2

-- | Get the closest distance(s) from a point to a 'FixedSegment'.
closestDistance :: OrderedField n => FixedSegment V2 n -> P2 n -> [n]
closestDistance :: forall n. OrderedField n => FixedSegment V2 n -> P2 n -> [n]
closestDistance = forall n. OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestDistance' forall n. Fractional n => n
defEps

-- | Get the closest distance(s) from a point to a 'FixedSegment' within given
--   tolerance.
closestDistance' :: OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestDistance' :: forall n. OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestDistance' n
eps FixedSegment V2 n
seg P2 n
p = forall a b. (a -> b) -> [a] -> [b]
map (forall a (p :: * -> *).
(Floating a, Foldable (Diff p), Affine p) =>
p a -> p a -> a
distanceA P2 n
p) forall a b. (a -> b) -> a -> b
$ forall n.
OrderedField n =>
n -> FixedSegment V2 n -> P2 n -> [P2 n]
closestPoint' n
eps FixedSegment V2 n
seg P2 n
p

-- | Get the closest point(s) on a 'FixedSegment' from a point.
closestPoint :: OrderedField n => FixedSegment V2 n -> P2 n -> [P2 n]
closestPoint :: forall n. OrderedField n => FixedSegment V2 n -> P2 n -> [P2 n]
closestPoint = forall n.
OrderedField n =>
n -> FixedSegment V2 n -> P2 n -> [P2 n]
closestPoint' forall n. Fractional n => n
defEps

-- | Get the closest point(s) on a 'FixedSegment' from a point within given
--   tolerance.
closestPoint' :: OrderedField n => n -> FixedSegment V2 n -> P2 n -> [P2 n]
closestPoint' :: forall n.
OrderedField n =>
n -> FixedSegment V2 n -> P2 n -> [P2 n]
closestPoint' n
eps FixedSegment V2 n
seg = forall a b. (a -> b) -> [a] -> [b]
map (FixedSegment V2 n
seg forall p. Parametric p => p -> N p -> Codomain p (N p)
`atParam`) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall n. OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestParam' n
eps FixedSegment V2 n
seg

-- | Find the closest value(s) on the Bêzier to the given point.
closestParam :: OrderedField n => FixedSegment V2 n -> P2 n -> [n]
closestParam :: forall n. OrderedField n => FixedSegment V2 n -> P2 n -> [n]
closestParam = forall n. OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestParam' forall n. Fractional n => n
defEps

-- | Find the closest value(s) on the Bêzier to the given point within given
--   tolerance.
closestParam' :: OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestParam' :: forall n. OrderedField n => n -> FixedSegment V2 n -> P2 n -> [n]
closestParam' n
_ (FLinear Point V2 n
p0 Point V2 n
p1) Point V2 n
p
  | n
t forall a. Ord a => a -> a -> Bool
< n
0     = [n
0]
  | n
t forall a. Ord a => a -> a -> Bool
> n
1     = [n
1]
  | Bool
otherwise = [n
t]
  where
    vp :: Diff (Point V2) n
vp = Point V2 n
p  forall (p :: * -> *) a. (Affine p, Num a) => p a -> p a -> Diff p a
.-. Point V2 n
p0
    v :: Diff (Point V2) n
v  = Point V2 n
p1 forall (p :: * -> *) a. (Affine p, Num a) => p a -> p a -> Diff p a
.-. Point V2 n
p0
    dp :: n
dp = Diff (Point V2) n
vp forall (f :: * -> *) a. (Metric f, Num a) => f a -> f a -> a
`dot` Diff (Point V2) n
v
    t :: n
t  = n
dp forall a. Fractional a => a -> a -> a
/ forall (f :: * -> *) a. (Metric f, Num a) => f a -> a
quadrance Diff (Point V2) n
v
closestParam' n
eps FixedSegment V2 n
cb (P (V2 n
px n
py)) = forall n. OrderedField n => n -> BernsteinPoly n -> n -> n -> [n]
bezierFindRoot n
eps BernsteinPoly n
poly n
0 n
1
  where
    (BernsteinPoly n
bx, BernsteinPoly n
by) = forall n.
Fractional n =>
FixedSegment V2 n -> (BernsteinPoly n, BernsteinPoly n)
bezierToBernstein FixedSegment V2 n
cb
    bx' :: BernsteinPoly n
bx'  = forall n. Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv BernsteinPoly n
bx
    by' :: BernsteinPoly n
by'  = forall n. Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv BernsteinPoly n
by
    poly :: BernsteinPoly n
poly = (BernsteinPoly n
bx forall a. Num a => a -> a -> a
- forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein [n
px, n
px, n
px, n
px]) forall a. Num a => a -> a -> a
* BernsteinPoly n
bx'
         forall a. Num a => a -> a -> a
+ (BernsteinPoly n
by forall a. Num a => a -> a -> a
- forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein [n
py, n
py, n
py, n
py]) forall a. Num a => a -> a -> a
* BernsteinPoly n
by'

------------------------------------------------------------------------
-- Low level
------------------------------------------------------------------------

-- | Return the intersection points with the parameters at which each segment
--   intersects.
segmentSegment :: OrderedField n => n -> FixedSegment V2 n -> FixedSegment V2 n -> [(n, n, P2 n)]
segmentSegment :: forall n.
OrderedField n =>
n -> FixedSegment V2 n -> FixedSegment V2 n -> [(n, n, P2 n)]
segmentSegment n
eps FixedSegment V2 n
s1 FixedSegment V2 n
s2 =
  case (FixedSegment V2 n
s1,FixedSegment V2 n
s2) of
    (FCubic{}, FCubic{})  -> forall a b. (a -> b) -> [a] -> [b]
map (\(n
t1,n
t2) -> (n
t1,n
t2, FixedSegment V2 n
s1 forall p. Parametric p => p -> N p -> Codomain p (N p)
`atParam` n
t1))
                           forall a b. (a -> b) -> a -> b
$ forall n.
OrderedField n =>
n -> FixedSegment V2 n -> FixedSegment V2 n -> [(n, n)]
bezierClip n
eps FixedSegment V2 n
s1 FixedSegment V2 n
s2
    (FCubic{}, FLinear{}) -> forall a b. (a -> b) -> [a] -> [b]
map forall {b} {a} {c}. (b, a, c) -> (a, b, c)
flip12 forall a b. (a -> b) -> a -> b
$ Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
linearSeg (forall (v :: * -> *) n.
InSpace v n (v n) =>
FixedSegment v n -> Located (v n)
segLine FixedSegment V2 n
s2) FixedSegment V2 n
s1
    (FixedSegment V2 n, FixedSegment V2 n)
_                     -> Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
linearSeg (forall (v :: * -> *) n.
InSpace v n (v n) =>
FixedSegment v n -> Located (v n)
segLine FixedSegment V2 n
s1) FixedSegment V2 n
s2 -- s1 is linear
  where
    linearSeg :: Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
linearSeg Located (V2 n)
l FixedSegment V2 n
s  = forall a. (a -> Bool) -> [a] -> [a]
filter (forall n. (Fractional n, Ord n) => n -> Bool
inRange forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view forall s t a b. Field1 s t a b => Lens s t a b
_1) forall a b. (a -> b) -> a -> b
$ forall n.
OrderedField n =>
n -> Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
lineSegment n
eps Located (V2 n)
l FixedSegment V2 n
s
    flip12 :: (b, a, c) -> (a, b, c)
flip12 (b
a,a
b,c
c) = (a
b,b
a,c
c)

-- | Return the intersection points with the parameters at which the line and segment
--   intersect.
lineSegment :: OrderedField n => n -> Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
lineSegment :: forall n.
OrderedField n =>
n -> Located (V2 n) -> FixedSegment V2 n -> [(n, n, P2 n)]
lineSegment n
_ Located (V2 n)
l1 p :: FixedSegment V2 n
p@(FLinear Point V2 n
p0 Point V2 n
p1)
  = forall a b. (a -> b) -> [a] -> [b]
map (\(n
tl,n
tp) -> (n
tl, n
tp, FixedSegment V2 n
p forall p. Parametric p => p -> N p -> Codomain p (N p)
`atParam` n
tp))
  forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (a -> Bool) -> [a] -> [a]
filter (forall n. (Fractional n, Ord n) => n -> Bool
inRange forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Maybe a -> [a]
maybeToList forall a b. (a -> b) -> a -> b
$ forall n.
(Fractional n, Eq n) =>
Located (V2 n) -> Located (V2 n) -> Maybe (n, n)
lineLine Located (V2 n)
l1 (forall (v :: * -> *) n.
InSpace v n (v n) =>
Point v n -> Point v n -> Located (v n)
mkLine Point V2 n
p0 Point V2 n
p1)
lineSegment n
eps (forall a. Located a -> (Point (V a) (N a), a)
viewLoc -> (Point (V (V2 n)) (N (V2 n))
p,V2 n
r)) FixedSegment V2 n
cb = forall a b. (a -> b) -> [a] -> [b]
map n -> (n, n, Point V2 n)
addPoint [n]
params
  where
    params :: [n]
params = forall n. OrderedField n => n -> BernsteinPoly n -> n -> n -> [n]
bezierFindRoot n
eps (forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein forall a b. (a -> b) -> a -> b
$ FixedSegment V2 n
cb' forall s a. s -> Getting (Endo [a]) s a -> [a]
^.. forall s t a b. Each s t a b => Traversal s t a b
each forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. R2 t => Lens' (t a) a
_y) n
0 n
1
    cb' :: FixedSegment V2 n
cb'    = forall t. Transformable t => Transformation (V t) (N t) -> t -> t
transform (forall (v :: * -> *) n.
(Functor v, Num n) =>
Transformation v n -> Transformation v n
inv (forall n. OrderedField n => Direction V2 n -> T2 n
rotationTo forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) n. v n -> Direction v n
dir V2 n
r)) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. HasOrigin t => Point (V t) (N t) -> t -> t
moveOriginTo Point (V (V2 n)) (N (V2 n))
p forall a b. (a -> b) -> a -> b
$ FixedSegment V2 n
cb
    --
    addPoint :: n -> (n, n, Point V2 n)
addPoint n
bt = (n
lt, n
bt, Codomain (FixedSegment V2 n) (N (FixedSegment V2 n))
intersect)
      where
        intersect :: Codomain (FixedSegment V2 n) (N (FixedSegment V2 n))
intersect = FixedSegment V2 n
cb forall p. Parametric p => p -> N p -> Codomain p (N p)
`atParam` n
bt
        lt :: n
lt        = (FixedSegment V2 n
cb' forall p. Parametric p => p -> N p -> Codomain p (N p)
`atParam` n
bt) forall s a. s -> Getting a s a -> a
^. forall (t :: * -> *) a. R1 t => Lens' (t a) a
_x forall a. Fractional a => a -> a -> a
/ forall (f :: * -> *) a. (Metric f, Floating a) => f a -> a
norm V2 n
r

-- Adapted from from kuribas's cubicbezier package https://github.com/kuribas/cubicbezier

-- | Use the Bêzier clipping algorithm to return the parameters at which the
--   Bêzier curves intersect.
bezierClip :: OrderedField n => n -> FixedSegment V2 n -> FixedSegment V2 n -> [(n, n)]
bezierClip :: forall n.
OrderedField n =>
n -> FixedSegment V2 n -> FixedSegment V2 n -> [(n, n)]
bezierClip n
eps FixedSegment V2 n
p_ FixedSegment V2 n
q_ = forall a. (a -> Bool) -> [a] -> [a]
filter (forall s a. Getting All s a -> (a -> Bool) -> s -> Bool
allOf forall (r :: * -> * -> *) a b.
Bitraversable r =>
Traversal (r a a) (r b b) a b
both forall n. (Fractional n, Ord n) => n -> Bool
inRange) -- sometimes this returns NaN
                     forall a b. (a -> b) -> a -> b
$ FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
p_ FixedSegment V2 n
q_ n
0 n
1 n
0 n
1 n
0 Bool
False
  where
  go :: FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
p FixedSegment V2 n
q n
tmin n
tmax n
umin n
umax n
clip Bool
revCurves
    | forall a. Maybe a -> Bool
isNothing Maybe (n, n)
chopInterval = []

    -- This check happens before the subdivision
    -- test to avoid non-termination as values
    -- transition to within epsilon.
    | forall a. Ord a => a -> a -> a
max (n
umax forall a. Num a => a -> a -> a
- n
umin) (n
tmax' forall a. Num a => a -> a -> a
- n
tmin') forall a. Ord a => a -> a -> Bool
< n
eps =
      if Bool
revCurves -- return parameters in correct order
      then [ (forall a. Fractional a => a -> a -> a
avg n
umin  n
umax,  forall a. Fractional a => a -> a -> a
avg n
tmin' n
tmax') ]
      else [ (forall a. Fractional a => a -> a -> a
avg n
tmin' n
tmax', forall a. Fractional a => a -> a -> a
avg n
umin  n
umax ) ]

    -- split the curve if there isn't enough reduction
    | n
clip forall a. Ord a => a -> a -> Bool
> n
0.8 Bool -> Bool -> Bool
&& n
clip' forall a. Ord a => a -> a -> Bool
> n
0.8 =
      if n
tmax' forall a. Num a => a -> a -> a
- n
tmin' forall a. Ord a => a -> a -> Bool
> n
umax forall a. Num a => a -> a -> a
- n
umin -- split the longest segment
      then let (FixedSegment V2 n
pl, FixedSegment V2 n
pr) = FixedSegment V2 n
p' forall p. Sectionable p => p -> N p -> (p, p)
`splitAtParam` n
0.5
               tmid :: n
tmid = forall a. Fractional a => a -> a -> a
avg n
tmin' n
tmax'
           in  FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
q FixedSegment V2 n
pl n
umin n
umax n
tmin' n
tmid  n
clip' (Bool -> Bool
not Bool
revCurves) forall a. [a] -> [a] -> [a]
++
               FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
q FixedSegment V2 n
pr n
umin n
umax n
tmid  n
tmax' n
clip' (Bool -> Bool
not Bool
revCurves)
      else let (FixedSegment V2 n
ql, FixedSegment V2 n
qr) = FixedSegment V2 n
q forall p. Sectionable p => p -> N p -> (p, p)
`splitAtParam` n
0.5
               umid :: n
umid = forall a. Fractional a => a -> a -> a
avg n
umin n
umax
           in  FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
ql FixedSegment V2 n
p' n
umin n
umid n
tmin' n
tmax' n
clip' (Bool -> Bool
not Bool
revCurves) forall a. [a] -> [a] -> [a]
++
               FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
qr FixedSegment V2 n
p' n
umid n
umax n
tmin' n
tmax' n
clip' (Bool -> Bool
not Bool
revCurves)

    -- iterate with the curves reversed.
    | Bool
otherwise = FixedSegment V2 n
-> FixedSegment V2 n -> n -> n -> n -> n -> n -> Bool -> [(n, n)]
go FixedSegment V2 n
q FixedSegment V2 n
p' n
umin n
umax n
tmin' n
tmax' n
clip' (Bool -> Bool
not Bool
revCurves)
    where
      chopInterval :: Maybe (n, n)
chopInterval              = forall n.
OrderedField n =>
FixedSegment V2 n -> FixedSegment V2 n -> Maybe (n, n)
chopCubics FixedSegment V2 n
p FixedSegment V2 n
q
      Just (n
tminChop, n
tmaxChop) = Maybe (n, n)
chopInterval
      p' :: FixedSegment V2 n
p'    = forall p. Sectionable p => p -> N p -> N p -> p
section FixedSegment V2 n
p n
tminChop n
tmaxChop
      clip' :: n
clip' = n
tmaxChop forall a. Num a => a -> a -> a
- n
tminChop
      tmin' :: n
tmin' = n
tmax forall a. Num a => a -> a -> a
* n
tminChop forall a. Num a => a -> a -> a
+ n
tmin forall a. Num a => a -> a -> a
* (n
1 forall a. Num a => a -> a -> a
- n
tminChop)
      tmax' :: n
tmax' = n
tmax forall a. Num a => a -> a -> a
* n
tmaxChop forall a. Num a => a -> a -> a
+ n
tmin forall a. Num a => a -> a -> a
* (n
1 forall a. Num a => a -> a -> a
- n
tmaxChop)

-- | Find the zero of a 1D Bêzier curve of any degree.  Note that this
--   can be used as a Bernstein polynomial root solver by converting from
--   the power basis to the Bernstein basis.
bezierFindRoot :: OrderedField n
               => n   -- ^ The accuracy
               -> BernsteinPoly n -- ^ the Bernstein coefficients of the polynomial
               -> n   -- ^ The lower bound of the interval
               -> n   -- ^ The upper bound of the interval
               -> [n] -- ^ The roots found
bezierFindRoot :: forall n. OrderedField n => n -> BernsteinPoly n -> n -> n -> [n]
bezierFindRoot n
eps BernsteinPoly n
p n
tmin n
tmax
  | forall a. Maybe a -> Bool
isNothing Maybe (n, n)
chopInterval    = []
  | n
clip forall a. Ord a => a -> a -> Bool
> n
0.8 = let (BernsteinPoly n
p1, BernsteinPoly n
p2) = forall p. Sectionable p => p -> N p -> (p, p)
splitAtParam BernsteinPoly n
newP n
0.5
                     tmid :: n
tmid     = n
tmin' forall a. Num a => a -> a -> a
+ (n
tmax' forall a. Num a => a -> a -> a
- n
tmin') forall a. Fractional a => a -> a -> a
/ n
2
                 in  forall n. OrderedField n => n -> BernsteinPoly n -> n -> n -> [n]
bezierFindRoot n
eps BernsteinPoly n
p1 n
tmin' n
tmid  forall a. [a] -> [a] -> [a]
++
                     forall n. OrderedField n => n -> BernsteinPoly n -> n -> n -> [n]
bezierFindRoot n
eps BernsteinPoly n
p2 n
tmid  n
tmax'
  | n
tmax' forall a. Num a => a -> a -> a
- n
tmin' forall a. Ord a => a -> a -> Bool
< n
eps = [forall a. Fractional a => a -> a -> a
avg n
tmin' n
tmax']
  | Bool
otherwise           = forall n. OrderedField n => n -> BernsteinPoly n -> n -> n -> [n]
bezierFindRoot n
eps BernsteinPoly n
newP n
tmin' n
tmax'
  where
    chopInterval :: Maybe (n, n)
chopInterval              = forall n. OrderedField n => [n] -> Maybe (n, n)
chopYs (forall n. BernsteinPoly n -> [n]
bernsteinCoeffs BernsteinPoly n
p)
    Just (n
tminChop, n
tmaxChop) = Maybe (n, n)
chopInterval
    newP :: BernsteinPoly n
newP  = forall p. Sectionable p => p -> N p -> N p -> p
section BernsteinPoly n
p n
tminChop n
tmaxChop
    clip :: n
clip  = n
tmaxChop forall a. Num a => a -> a -> a
- n
tminChop
    tmin' :: n
tmin' = n
tmax forall a. Num a => a -> a -> a
* n
tminChop forall a. Num a => a -> a -> a
+ n
tmin forall a. Num a => a -> a -> a
* (n
1 forall a. Num a => a -> a -> a
- n
tminChop)
    tmax' :: n
tmax' = n
tmax forall a. Num a => a -> a -> a
* n
tmaxChop forall a. Num a => a -> a -> a
+ n
tmin forall a. Num a => a -> a -> a
* (n
1 forall a. Num a => a -> a -> a
- n
tmaxChop)



------------------------------------------------------------------------
-- Internal
------------------------------------------------------------------------

-- | An approximation of the fat line for a cubic Bêzier segment. Returns
--   @(0,0)@ for a linear segment.
fatLine :: OrderedField n => FixedSegment V2 n -> (n,n)
fatLine :: forall n. OrderedField n => FixedSegment V2 n -> (n, n)
fatLine (FCubic Point V2 n
p0 Point V2 n
p1 Point V2 n
p2 Point V2 n
p3)
  = case (n
d1 forall a. Ord a => a -> a -> Bool
> n
0, n
d2 forall a. Ord a => a -> a -> Bool
> n
0) of
      (Bool
True,  Bool
True)  -> (n
0,                n
0.75 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
max n
d1 n
d2)
      (Bool
False, Bool
False) -> (n
0.75 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
min n
d1 n
d2, n
0               )
      (Bool
True,  Bool
False) -> (n
4forall a. Fractional a => a -> a -> a
/n
9 forall a. Num a => a -> a -> a
* n
d2,         n
4forall a. Fractional a => a -> a -> a
/n
9 forall a. Num a => a -> a -> a
* n
d1        )
      (Bool
False, Bool
True)  -> (n
4forall a. Fractional a => a -> a -> a
/n
9 forall a. Num a => a -> a -> a
* n
d1,         n
4forall a. Fractional a => a -> a -> a
/n
9 forall a. Num a => a -> a -> a
* n
d2        )
  where
    d :: Point V2 n -> n
d = forall n. (Ord n, Floating n) => P2 n -> P2 n -> P2 n -> n
lineDistance Point V2 n
p0 Point V2 n
p3
    d1 :: n
d1 = Point V2 n -> n
d Point V2 n
p1; d2 :: n
d2 = Point V2 n -> n
d Point V2 n
p2
fatLine FixedSegment V2 n
_ = (n
0,n
0)

chopYs :: OrderedField n => [n] -> Maybe (n, n)
chopYs :: forall n. OrderedField n => [n] -> Maybe (n, n)
chopYs [n]
ds = forall n. OrderedField n => n -> n -> [P2 n] -> Maybe (n, n)
chopHull n
0 n
0 [P2 n]
points
  where
    points :: [P2 n]
points = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall n. n -> n -> P2 n
mkP2 [forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n | Int
i <- [Int
0..Int
n]] [n]
ds
    n :: Int
n      = forall (t :: * -> *) a. Foldable t => t a -> Int
length [n]
ds forall a. Num a => a -> a -> a
- Int
1

chopCubics :: OrderedField n => FixedSegment V2 n -> FixedSegment V2 n -> Maybe (n,n)
chopCubics :: forall n.
OrderedField n =>
FixedSegment V2 n -> FixedSegment V2 n -> Maybe (n, n)
chopCubics FixedSegment V2 n
p q :: FixedSegment V2 n
q@(FCubic P2 n
q0 P2 n
_ P2 n
_ P2 n
q3)
  = forall n. OrderedField n => n -> n -> [P2 n] -> Maybe (n, n)
chopHull n
dmin n
dmax [P2 n]
dps
  where
    dps :: [P2 n]
dps = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall n. n -> n -> P2 n
mkP2 [n
0, n
1forall a. Fractional a => a -> a -> a
/n
3, n
2forall a. Fractional a => a -> a -> a
/n
3, n
1] [n]
ds
    ds :: [n]
ds  = FixedSegment V2 n
p forall s a. s -> Getting (Endo [a]) s a -> [a]
^.. forall s t a b. Each s t a b => Traversal s t a b
each forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (p :: * -> * -> *) (f :: * -> *) s a.
(Profunctor p, Contravariant f) =>
(s -> a) -> Optic' p f s a
to P2 n -> n
d
    d :: P2 n -> n
d   = forall n. (Ord n, Floating n) => P2 n -> P2 n -> P2 n -> n
lineDistance P2 n
q0 P2 n
q3
    --
    (n
dmin,n
dmax) = forall n. OrderedField n => FixedSegment V2 n -> (n, n)
fatLine FixedSegment V2 n
q
chopCubics FixedSegment V2 n
_ FixedSegment V2 n
_ = forall a. Maybe a
Nothing

-- Reduce the interval which the intersection is known to lie in using the fat
-- line of one curve and convex hull of the points formed from the distance to
-- the thin line of the other
chopHull :: OrderedField n => n -> n -> [P2 n] -> Maybe (n, n)
chopHull :: forall n. OrderedField n => n -> n -> [P2 n] -> Maybe (n, n)
chopHull n
dmin n
dmax [P2 n]
dps = do
  n
tL <- [P2 n] -> Maybe n -> Maybe n
testBelow [P2 n]
upper           forall a b. (a -> b) -> a -> b
$ P2 n -> Maybe n -> Maybe n
testBetween (forall a. [a] -> a
head [P2 n]
upper) forall a b. (a -> b) -> a -> b
$ [P2 n] -> Maybe n
testAbove [P2 n]
lower
  n
tR <- [P2 n] -> Maybe n -> Maybe n
testBelow (forall a. [a] -> [a]
reverse [P2 n]
upper) forall a b. (a -> b) -> a -> b
$ P2 n -> Maybe n -> Maybe n
testBetween (forall a. [a] -> a
last [P2 n]
upper) forall a b. (a -> b) -> a -> b
$ [P2 n] -> Maybe n
testAbove (forall a. [a] -> [a]
reverse [P2 n]
lower)
  forall a. a -> Maybe a
Just (n
tL, n
tR)
    where
      ([P2 n]
upper, [P2 n]
lower) = forall n. OrderedField n => [P2 n] -> ([P2 n], [P2 n])
sortedConvexHull [P2 n]
dps

      testBelow :: [P2 n] -> Maybe n -> Maybe n
testBelow (p1 :: P2 n
p1@(P (V2 n
_ n
y1)) : p2 :: P2 n
p2@(P (V2 n
_ n
y2)) : [P2 n]
ps) Maybe n
continue
        | n
y1 forall a. Ord a => a -> a -> Bool
>= n
dmin = Maybe n
continue
        | n
y1 forall a. Ord a => a -> a -> Bool
>  n
y2   = forall a. Maybe a
Nothing
        | n
y2 forall a. Ord a => a -> a -> Bool
<  n
dmin = [P2 n] -> Maybe n -> Maybe n
testBelow (P2 n
p2forall a. a -> [a] -> [a]
:[P2 n]
ps) Maybe n
continue
        | Bool
otherwise  = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall {a}. Fractional a => a -> Point V2 a -> Point V2 a -> a
intersectPt n
dmin P2 n
p1 P2 n
p2
      testBelow [P2 n]
_ Maybe n
_  = forall a. Maybe a
Nothing

      testBetween :: P2 n -> Maybe n -> Maybe n
testBetween (P (V2 n
x n
y)) Maybe n
continue
        | n
y forall a. Ord a => a -> a -> Bool
<= n
dmax = forall a. a -> Maybe a
Just n
x
        | Bool
otherwise = Maybe n
continue

      testAbove :: [P2 n] -> Maybe n
testAbove (p1 :: P2 n
p1@(P (V2 n
_ n
y1)) : p2 :: P2 n
p2@(P (V2 n
_ n
y2)) : [P2 n]
ps)
        | n
y1 forall a. Ord a => a -> a -> Bool
< n
y2      = forall a. Maybe a
Nothing
        | n
y2 forall a. Ord a => a -> a -> Bool
> n
dmax    = [P2 n] -> Maybe n
testAbove (P2 n
p2forall a. a -> [a] -> [a]
:[P2 n]
ps)
        | n
y2 forall a. Num a => a -> a -> a
- n
y1 forall a. Eq a => a -> a -> Bool
== n
0 = forall a. Maybe a
Nothing  -- Check this condition to prevent
                                  -- division by zero in `intersectPt`.
        | Bool
otherwise    = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall {a}. Fractional a => a -> Point V2 a -> Point V2 a -> a
intersectPt n
dmax P2 n
p1 P2 n
p2
      testAbove [P2 n]
_      = forall a. Maybe a
Nothing

      -- find the x value where the line through the two points
      -- intersect the line y=d.  Note that `y2 - y1 != 0` due
      -- to checks above.
      intersectPt :: a -> Point V2 a -> Point V2 a -> a
intersectPt a
d (P (V2 a
x1 a
y1)) (P (V2 a
x2 a
y2)) =
          a
x1 forall a. Num a => a -> a -> a
+ (a
d forall a. Num a => a -> a -> a
- a
y1) forall a. Num a => a -> a -> a
* (a
x2 forall a. Num a => a -> a -> a
- a
x1) forall a. Fractional a => a -> a -> a
/ (a
y2 forall a. Num a => a -> a -> a
- a
y1)



bezierToBernstein :: Fractional n => FixedSegment V2 n -> (BernsteinPoly n, BernsteinPoly n)
bezierToBernstein :: forall n.
Fractional n =>
FixedSegment V2 n -> (BernsteinPoly n, BernsteinPoly n)
bezierToBernstein FixedSegment V2 n
seg =
    (forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view forall (t :: * -> *) a. R1 t => Lens' (t a) a
_x) [Point V2 n]
coeffs, forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view forall (t :: * -> *) a. R2 t => Lens' (t a) a
_y) [Point V2 n]
coeffs)
  where coeffs :: [Point V2 n]
coeffs = forall a s. Getting (Endo [a]) s a -> s -> [a]
toListOf forall s t a b. Each s t a b => Traversal s t a b
each FixedSegment V2 n
seg

------------------------------------------------------------------------
-- Lines
------------------------------------------------------------------------

-- Could split this into a separate module.

-- | Returns @(a, b, c, d)@ such that @ax + by + c = 0@ is the line going through
--   @p1@ and @p2@ with @(a^2)/d + (b^2)/d = 1@.  We delay the division by
--   @d@ as it may not be needed in all cases and @d@ may be zero.
lineEquation :: Floating n => P2 n -> P2 n -> (n, n, n, n)
lineEquation :: forall n. Floating n => P2 n -> P2 n -> (n, n, n, n)
lineEquation (P (V2 n
x1 n
y1)) (P (V2 n
x2 n
y2)) = (n
a, n
b, n
c, n
d)
  where
    c :: n
c  = -(n
x1forall a. Num a => a -> a -> a
*n
a forall a. Num a => a -> a -> a
+ n
y1forall a. Num a => a -> a -> a
*n
b)
    a :: n
a = n
y1 forall a. Num a => a -> a -> a
- n
y2
    b :: n
b = n
x2 forall a. Num a => a -> a -> a
- n
x1
    d :: n
d  = n
aforall a. Num a => a -> a -> a
*n
a forall a. Num a => a -> a -> a
+ n
bforall a. Num a => a -> a -> a
*n
b

-- | Return the distance from a point to the line.
lineDistance :: (Ord n, Floating n) => P2 n -> P2 n -> P2 n -> n
lineDistance :: forall n. (Ord n, Floating n) => P2 n -> P2 n -> P2 n -> n
lineDistance P2 n
p1 P2 n
p2 p3 :: P2 n
p3@(P (V2 n
x n
y))
    -- I have included the check that d' <= 0 in case
    -- there exists some d > 0 where sqrt d == 0.  I don't
    -- think this can happen as sqrt is at least recommended
    -- to be within one value of correct for sqrt and near
    -- zero values get bigger.
    | n
d forall a. Ord a => a -> a -> Bool
<= n
0 Bool -> Bool -> Bool
|| n
d' forall a. Ord a => a -> a -> Bool
<= n
0 = forall (f :: * -> *) a. (Metric f, Floating a) => f a -> a
norm (P2 n
p1 forall (p :: * -> *) a. (Affine p, Num a) => p a -> p a -> Diff p a
.-. P2 n
p3)
    | Bool
otherwise = (n
aforall a. Num a => a -> a -> a
*n
x forall a. Num a => a -> a -> a
+ n
bforall a. Num a => a -> a -> a
*n
y forall a. Num a => a -> a -> a
+ n
c) forall a. Fractional a => a -> a -> a
/ n
d'
  where
    (n
a, n
b, n
c, n
d) = forall n. Floating n => P2 n -> P2 n -> (n, n, n, n)
lineEquation P2 n
p1 P2 n
p2
    d' :: n
d' = forall a. Floating a => a -> a
sqrt n
d

-- clockwise :: (Num n, Ord n) => V2 n -> V2 n -> Bool
-- clockwise a b = a `cross2` b <= 0

avg :: Fractional n => n -> n -> n
avg :: forall a. Fractional a => a -> a -> a
avg n
a n
b = (n
a forall a. Num a => a -> a -> a
+ n
b)forall a. Fractional a => a -> a -> a
/n
2

lineLine :: (Fractional n, Eq n) => Located (V2 n) -> Located (V2 n) -> Maybe (n,n)
lineLine :: forall n.
(Fractional n, Eq n) =>
Located (V2 n) -> Located (V2 n) -> Maybe (n, n)
lineLine (forall a. Located a -> (Point (V a) (N a), a)
viewLoc -> (Point (V (V2 n)) (N (V2 n))
p,V2 n
r)) (forall a. Located a -> (Point (V a) (N a), a)
viewLoc -> (Point (V (V2 n)) (N (V2 n))
q,V2 n
s))
  | n
x1 forall a. Eq a => a -> a -> Bool
== n
0 Bool -> Bool -> Bool
&& n
x2 forall a. Eq a => a -> a -> Bool
/= n
0 = forall a. Maybe a
Nothing                 -- parallel
  | Bool
otherwise          = forall a. a -> Maybe a
Just (n
x3 forall a. Fractional a => a -> a -> a
/ n
x1, n
x2 forall a. Fractional a => a -> a -> a
/ n
x1) -- intersecting or collinear
  where
    x1 :: n
x1 = V2 n
r forall n. Num n => V2 n -> V2 n -> n
× V2 n
s
    x2 :: n
x2 = Diff (Point V2) n
v forall n. Num n => V2 n -> V2 n -> n
× V2 n
r
    x3 :: n
x3 = Diff (Point V2) n
v forall n. Num n => V2 n -> V2 n -> n
× V2 n
s
    v :: Diff (Point V2) n
v  = Point (V (V2 n)) (N (V2 n))
q forall (p :: * -> *) a. (Affine p, Num a) => p a -> p a -> Diff p a
.-. Point (V (V2 n)) (N (V2 n))
p

(×) :: Num n => V2 n -> V2 n -> n
× :: forall n. Num n => V2 n -> V2 n -> n
(×) = forall n. Num n => V2 n -> V2 n -> n
cross2

mkLine :: InSpace v n (v n) => Point v n -> Point v n -> Located (v n)
mkLine :: forall (v :: * -> *) n.
InSpace v n (v n) =>
Point v n -> Point v n -> Located (v n)
mkLine Point v n
p0 Point v n
p1 = (Point v n
p1 forall (p :: * -> *) a. (Affine p, Num a) => p a -> p a -> Diff p a
.-. Point v n
p0) forall a. a -> Point (V a) (N a) -> Located a
`at` Point v n
p0

segLine :: InSpace v n (v n) => FixedSegment v n -> Located (v n)
segLine :: forall (v :: * -> *) n.
InSpace v n (v n) =>
FixedSegment v n -> Located (v n)
segLine (FLinear Point v n
p0 Point v n
p1)    = forall (v :: * -> *) n.
InSpace v n (v n) =>
Point v n -> Point v n -> Located (v n)
mkLine Point v n
p0 Point v n
p1
segLine (FCubic Point v n
p0 Point v n
_ Point v n
_ Point v n
p3) = forall (v :: * -> *) n.
InSpace v n (v n) =>
Point v n -> Point v n -> Located (v n)
mkLine Point v n
p0 Point v n
p3

-- This function uses `defEps`, but is used in functions
-- above that take an epsilon parameter.  It would be nice
-- to clearify the meaning of each of these epsilons.
inRange :: (Fractional n, Ord n) => n -> Bool
inRange :: forall n. (Fractional n, Ord n) => n -> Bool
inRange n
x = n
x forall a. Ord a => a -> a -> Bool
< (n
1forall a. Num a => a -> a -> a
+forall n. Fractional n => n
defEps) Bool -> Bool -> Bool
&& n
x forall a. Ord a => a -> a -> Bool
> (-forall n. Fractional n => n
defEps)