{-# LANGUAGE BangPatterns, MultiWayIf #-}
-- | Intersection routines using Bezier Clipping.  Provides also functions for finding the roots of onedimensional bezier curves.  This can be used as a general polynomial root solver by converting from the power basis to the bernstein basis.
module Geom2D.CubicBezier.Intersection
       (bezierIntersection, bezierLineIntersections, bezierFindRoot)
       where
import Geom2D
import Geom2D.CubicBezier.Basic
import Math.BernsteinPoly
import Data.Maybe
import Geom2D.CubicBezier.Numeric
import qualified Data.Vector.Unboxed as V

-- find the convex hull by comparing the angles of the vectors with
-- the cross product and backtracking if necessary.
findOuter' :: Bool -> DPoint -> DPoint -> [DPoint] -> Either [DPoint] [DPoint]
findOuter' :: Bool
-> Point Double
-> Point Double
-> [Point Double]
-> Either [Point Double] [Point Double]
findOuter' !Bool
upper !Point Double
dir !Point Double
p1 l :: [Point Double]
l@(Point Double
p2:[Point Double]
rest)
  -- backtrack if the direction is outward
  | if Bool
upper
    then Point Double
dir forall a. Num a => Point a -> Point a -> a
`vectorCross` (Point Double
p2forall v. AdditiveGroup v => v -> v -> v
^-^Point Double
p1) forall a. Ord a => a -> a -> Bool
> Double
0 -- left turn
    else Point Double
dir forall a. Num a => Point a -> Point a -> a
`vectorCross` (Point Double
p2forall v. AdditiveGroup v => v -> v -> v
^-^Point Double
p1) forall a. Ord a => a -> a -> Bool
< Double
0 = forall a b. a -> Either a b
Left [Point Double]
l
  -- succeed
  | Bool
otherwise = case Bool
-> Point Double
-> Point Double
-> [Point Double]
-> Either [Point Double] [Point Double]
findOuter' Bool
upper (Point Double
p2forall v. AdditiveGroup v => v -> v -> v
^-^Point Double
p1) Point Double
p2 [Point Double]
rest of
    Left [Point Double]
m -> Bool
-> Point Double
-> Point Double
-> [Point Double]
-> Either [Point Double] [Point Double]
findOuter' Bool
upper Point Double
dir Point Double
p1 [Point Double]
m
    Right [Point Double]
m -> forall a b. b -> Either a b
Right (Point Double
p1forall a. a -> [a] -> [a]
:[Point Double]
m)

findOuter' Bool
_ Point Double
_ Point Double
p1 [Point Double]
p = forall a b. b -> Either a b
Right (Point Double
p1forall a. a -> [a] -> [a]
:[Point Double]
p)

-- find the outermost point.  It doesn't look at the x values.
findOuter :: Bool -> [DPoint] -> [DPoint]
findOuter :: Bool -> [Point Double] -> [Point Double]
findOuter Bool
upper (Point Double
p1:Point Double
p2:[Point Double]
rest) =
  case Bool
-> Point Double
-> Point Double
-> [Point Double]
-> Either [Point Double] [Point Double]
findOuter' Bool
upper (Point Double
p2forall v. AdditiveGroup v => v -> v -> v
^-^Point Double
p1) Point Double
p2 [Point Double]
rest of
    Right [Point Double]
l -> Point Double
p1forall a. a -> [a] -> [a]
:[Point Double]
l
    Left [Point Double]
l -> Bool -> [Point Double] -> [Point Double]
findOuter Bool
upper (Point Double
p1forall a. a -> [a] -> [a]
:[Point Double]
l)
findOuter Bool
_ [Point Double]
l = [Point Double]
l    

-- take the y values and turn it in into a convex hull with upper en
-- lower points separated.
makeHull :: [Double] -> ([DPoint], [DPoint])
makeHull :: [Double] -> ([Point Double], [Point Double])
makeHull [Double]
ds =
  let n :: Double
n      = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ds forall a. Num a => a -> a -> a
- Int
1
      points :: [Point Double]
points = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. a -> a -> Point a
Point [Double
iforall a. Fractional a => a -> a -> a
/Double
n | Double
i <- [Double
0..Double
n]] [Double]
ds
  in (Bool -> [Point Double] -> [Point Double]
findOuter Bool
True [Point Double]
points,
      Bool -> [Point Double] -> [Point Double]
findOuter Bool
False [Point Double]
points)

-- test if the chords cross the fat line
-- return the continuation if above the line
testBelow :: Double -> [DPoint] -> Maybe Double -> Maybe Double
testBelow :: Double -> [Point Double] -> Maybe Double -> Maybe Double
testBelow Double
_    [] Maybe Double
_ = forall a. Maybe a
Nothing
testBelow Double
_    [Point Double
_] Maybe Double
_ = forall a. Maybe a
Nothing
testBelow !Double
dmin (Point Double
p:Point Double
q:[Point Double]
rest) Maybe Double
cont
  | forall a. Point a -> a
pointY Point Double
p forall a. Ord a => a -> a -> Bool
>= Double
dmin = Maybe Double
cont
  | forall a. Point a -> a
pointY Point Double
p forall a. Ord a => a -> a -> Bool
> forall a. Point a -> a
pointY Point Double
q = forall a. Maybe a
Nothing
  | forall a. Point a -> a
pointY Point Double
q forall a. Ord a => a -> a -> Bool
< Double
dmin = Double -> [Point Double] -> Maybe Double -> Maybe Double
testBelow Double
dmin (Point Double
qforall a. a -> [a] -> [a]
:[Point Double]
rest) Maybe Double
cont
  | Bool
otherwise = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Double -> Point Double -> Point Double -> Double
intersectPt Double
dmin Point Double
p Point Double
q

testBetween :: Double -> DPoint -> Maybe Double -> Maybe Double
testBetween :: Double -> Point Double -> Maybe Double -> Maybe Double
testBetween !Double
dmax (Point !Double
x !Double
y) Maybe Double
cont
  | Double
y forall a. Ord a => a -> a -> Bool
<= Double
dmax = forall a. a -> Maybe a
Just Double
x
  | Bool
otherwise = Maybe Double
cont

-- test if the chords cross the line y=dmax somewhere
testAbove :: Double -> [DPoint] -> Maybe Double
testAbove :: Double -> [Point Double] -> Maybe Double
testAbove Double
_    [] = forall a. Maybe a
Nothing
testAbove Double
_    [Point Double
_] = forall a. Maybe a
Nothing
testAbove Double
dmax (Point Double
p:Point Double
q:[Point Double]
rest)
  | forall a. Point a -> a
pointY Point Double
p forall a. Ord a => a -> a -> Bool
< forall a. Point a -> a
pointY Point Double
q = forall a. Maybe a
Nothing
  | forall a. Point a -> a
pointY Point Double
q forall a. Ord a => a -> a -> Bool
> Double
dmax = Double -> [Point Double] -> Maybe Double
testAbove Double
dmax (Point Double
qforall a. a -> [a] -> [a]
:[Point Double]
rest)
  | Bool
otherwise = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Double -> Point Double -> Point Double -> Double
intersectPt Double
dmax Point Double
p Point Double
q

-- find the x value where the line through the two points
-- intersect the line y=d
intersectPt :: Double -> DPoint -> DPoint -> Double
intersectPt :: Double -> Point Double -> Point Double -> Double
intersectPt Double
d (Point Double
x1 Double
y1) (Point Double
x2 Double
y2)
  | Double
y1 forall a. Eq a => a -> a -> Bool
== Double
y2 = Double
x1
  | Bool
otherwise =
    Double
x1 forall a. Num a => a -> a -> a
+ (Double
d forall a. Num a => a -> a -> a
- Double
y1) forall a. Num a => a -> a -> a
* (Double
x2 forall a. Num a => a -> a -> a
- Double
x1) forall a. Fractional a => a -> a -> a
/ (Double
y2 forall a. Num a => a -> a -> a
- Double
y1)

-- make a hull and test over which interval the
-- curve is garuanteed to lie inside the fat line
chopHull :: Double -> Double -> [Double] -> Maybe (Double, Double)
chopHull :: Double -> Double -> [Double] -> Maybe (Double, Double)
chopHull !Double
dmin !Double
dmax [Double]
ds = do
  let ([Point Double]
upper, [Point Double]
lower) = [Double] -> ([Point Double], [Point Double])
makeHull [Double]
ds
  Double
left_t <- Double -> [Point Double] -> Maybe Double -> Maybe Double
testBelow Double
dmin [Point Double]
upper forall a b. (a -> b) -> a -> b
$
            Double -> Point Double -> Maybe Double -> Maybe Double
testBetween Double
dmax (forall a. [a] -> a
head [Point Double]
upper) forall a b. (a -> b) -> a -> b
$
            Double -> [Point Double] -> Maybe Double
testAbove Double
dmax [Point Double]
lower
  Double
right_t <- Double -> [Point Double] -> Maybe Double -> Maybe Double
testBelow Double
dmin (forall a. [a] -> [a]
reverse [Point Double]
upper) forall a b. (a -> b) -> a -> b
$
             Double -> Point Double -> Maybe Double -> Maybe Double
testBetween Double
dmax (forall a. [a] -> a
last [Point Double]
upper) forall a b. (a -> b) -> a -> b
$
             Double -> [Point Double] -> Maybe Double
testAbove Double
dmax (forall a. [a] -> [a]
reverse [Point Double]
lower)
  forall a. a -> Maybe a
Just (Double
left_t, Double
right_t)

bezierClip :: CubicBezier Double -> CubicBezier Double -> Double -> Double
           -> Double -> Double -> Double -> Double -> Double -> Bool
           -> [(Double, Double)]
bezierClip :: CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip p :: CubicBezier Double
p@(CubicBezier !Point Double
p0 !Point Double
p1 !Point Double
p2 !Point Double
p3) q :: CubicBezier Double
q@(CubicBezier !Point Double
q0 !Point Double
q1 !Point Double
q2 !Point Double
q3)
  Double
tmin Double
tmax Double
umin Double
umax Double
prevClip Double
pEps Double
vEps Bool
revCurves = forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either forall a. a -> a
id forall a. a -> a
id forall a b. (a -> b) -> a -> b
$ do
  Point Double
q3' <- if | forall a. Floating a => Point a -> Point a -> a
vectorDistance Point Double
q0 Point Double
q3 forall a. Ord a => a -> a -> Bool
> forall a. Ord a => a -> a -> a
max (forall a. Floating a => Point a -> a
vectorMag Point Double
q0) (forall a. Floating a => Point a -> a
vectorMag Point Double
q3) forall a. Fractional a => a -> a -> a
/ (Double
2forall a. Floating a => a -> a -> a
**Double
30) -> forall a b. b -> Either a b
Right Point Double
q3
            | forall a. Floating a => Point a -> Point a -> a
vectorDistance Point Double
q0 Point Double
q1 forall a. Ord a => a -> a -> Bool
> forall a. Ord a => a -> a -> a
max (forall a. Floating a => Point a -> a
vectorMag Point Double
q0) (forall a. Floating a => Point a -> a
vectorMag Point Double
q1) forall a. Fractional a => a -> a -> a
/ (Double
2forall a. Floating a => a -> a -> a
**Double
30) -> forall a b. b -> Either a b
Right Point Double
q1
            | forall a. Floating a => Point a -> Point a -> a
vectorDistance Point Double
q0 Point Double
q2 forall a. Ord a => a -> a -> Bool
> forall a. Ord a => a -> a -> a
max (forall a. Floating a => Point a -> a
vectorMag Point Double
q0) (forall a. Floating a => Point a -> a
vectorMag Point Double
q2) forall a. Fractional a => a -> a -> a
/ (Double
2forall a. Floating a => a -> a -> a
**Double
30) -> forall a b. b -> Either a b
Right Point Double
q2
            | Bool
otherwise -> forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$
              let t :: Double
t = CubicBezier Double -> Point Double -> Double -> Double
closest CubicBezier Double
p Point Double
q0 Double
vEps
                  newT :: Double
newT = Double
tmin forall a. Num a => a -> a -> a
* (Double
1forall a. Num a => a -> a -> a
-Double
t) forall a. Num a => a -> a -> a
+ Double
tmax forall a. Num a => a -> a -> a
* Double
t
                  umid :: Double
umid | Double
umax forall a. Ord a => a -> a -> Bool
>= Double
0.5 = Double
umax
                       | Bool
otherwise = Double
umin
              in if | forall a. Floating a => Point a -> Point a -> a
vectorDistance (forall (b :: * -> *) a.
(GenericBezier b, Unbox a, Fractional a) =>
b a -> a -> Point a
evalBezier CubicBezier Double
p Double
t) (forall (b :: * -> *) a.
(GenericBezier b, Unbox a, Fractional a) =>
b a -> a -> Point a
evalBezier CubicBezier Double
q Double
0.5) forall a. Ord a => a -> a -> Bool
> Double
vEps
                      -> []
                    | Bool
revCurves -> [(Double
umid, Double
newT)]
                    | Bool
otherwise -> [(Double
newT, Double
umid)]
  let d :: Point Double -> Double
d = forall a. Floating a => Line a -> Point a -> a
lineDistance (forall a. Point a -> Point a -> Line a
Line Point Double
q0 Point Double
q3')
      d1 :: Double
d1 = Point Double -> Double
d Point Double
q1
      d2 :: Double
d2 = Point Double -> Double
d Point Double
q2
      (Double
dmin, Double
dmax) | Double
d1forall a. Num a => a -> a -> a
*Double
d2 forall a. Ord a => a -> a -> Bool
> Double
0 = (Double
3forall a. Fractional a => a -> a -> a
/Double
4 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
min Double
0 (forall a. Ord a => a -> a -> a
min Double
d1 Double
d2),
                                  Double
3forall a. Fractional a => a -> a -> a
/Double
4 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
max Double
0 (forall a. Ord a => a -> a -> a
max Double
d1 Double
d2))
                   | Bool
otherwise = (Double
4forall a. Fractional a => a -> a -> a
/Double
9 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
min Double
0 (forall a. Ord a => a -> a -> a
min Double
d1 Double
d2),
                                  Double
4forall a. Fractional a => a -> a -> a
/Double
9 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
max Double
0 (forall a. Ord a => a -> a -> a
max Double
d1 Double
d2))
  (Double
chop_tmin, Double
chop_tmax) <- forall b a. b -> (a -> b) -> Maybe a -> b
maybe (forall a b. a -> Either a b
Left []) forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$
                            Double -> Double -> [Double] -> Maybe (Double, Double)
chopHull Double
dmin Double
dmax forall a b. (a -> b) -> a -> b
$
                            forall a b. (a -> b) -> [a] -> [b]
map Point Double -> Double
d [Point Double
p0, Point Double
p1, Point Double
p2, Point Double
p3]
  let newP :: CubicBezier Double
newP = forall a (b :: * -> *).
(Ord a, Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> a -> b a
bezierSubsegment CubicBezier Double
p Double
chop_tmin Double
chop_tmax
      newClip :: Double
newClip = Double
chop_tmax forall a. Num a => a -> a -> a
- Double
chop_tmin
      new_tmin :: Double
new_tmin = Double
tmax forall a. Num a => a -> a -> a
* Double
chop_tmin forall a. Num a => a -> a -> a
+ Double
tmin forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
chop_tmin)
      new_tmax :: Double
new_tmax = Double
tmax forall a. Num a => a -> a -> a
* Double
chop_tmax forall a. Num a => a -> a -> a
+ Double
tmin forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
chop_tmax)
  if | -- within tolerance      
       (forall a. Ord a => a -> a -> a
max (Double
umax forall a. Num a => a -> a -> a
- Double
umin) (Double
new_tmax forall a. Num a => a -> a -> a
- Double
new_tmin))forall a. Num a => a -> a -> a
*Double
4 forall a. Ord a => a -> a -> Bool
< Double
pEps ->
       let newu :: Double
newu | Double
umax forall a. Ord a => a -> a -> Bool
>= Double
0.5 = Double
umax
                | Bool
otherwise = Double
umin
           newt :: Double
newt | Double
new_tmax forall a. Ord a => a -> a -> Bool
>= Double
0.5 = Double
new_tmax
                | Bool
otherwise = Double
new_tmin
       in if Bool
revCurves
       then forall a b. b -> Either a b
Right [(Double
newu, Double
newt)]
       else forall a b. b -> Either a b
Right [(Double
newt, Double
newu)]
           -- not enough reduction, so split the curve in case we have
           -- multiple intersections
     | Double
prevClip forall a. Ord a => a -> a -> Bool
> Double
0.8 Bool -> Bool -> Bool
&& Double
newClip forall a. Ord a => a -> a -> Bool
> Double
0.8 ->
             if | (Double
new_tmaxforall a. Num a => a -> a -> a
-Double
new_tmin) forall a. Num a => a -> a -> a
* (Double
new_tmaxforall a. Num a => a -> a -> a
-Double
new_tmin) forall a. Num a => a -> a -> a
* forall a. Floating a => Point a -> a
vectorMagSquare (Point Double
p3 forall v. AdditiveGroup v => v -> v -> v
^-^ Point Double
p0) forall a. Ord a => a -> a -> Bool
>
                  (Double
umaxforall a. Num a => a -> a -> a
-Double
umin) forall a. Num a => a -> a -> a
* (Double
umaxforall a. Num a => a -> a -> a
-Double
umin) forall a. Num a => a -> a -> a
* forall a. Floating a => Point a -> a
vectorMagSquare (Point Double
q3 forall v. AdditiveGroup v => v -> v -> v
^-^ Point Double
q0) ->
                    -- split the longest segment
                  let (CubicBezier Double
pl, CubicBezier Double
pr) = forall a (b :: * -> *).
(Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> (b a, b a)
splitBezier CubicBezier Double
newP Double
0.5
                      half_t :: Double
half_t = Double
new_tmin forall a. Num a => a -> a -> a
+ (Double
new_tmax forall a. Num a => a -> a -> a
- Double
new_tmin) forall a. Fractional a => a -> a -> a
/ Double
2
                  in forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip CubicBezier Double
q CubicBezier Double
pl Double
umin Double
umax Double
new_tmin Double
half_t
                     Double
0 Double
pEps Double
vEps (Bool -> Bool
not Bool
revCurves) forall a. [a] -> [a] -> [a]
++
                     CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip CubicBezier Double
q CubicBezier Double
pr Double
umin Double
umax Double
half_t Double
new_tmax
                     Double
0 Double
pEps Double
vEps (Bool -> Bool
not Bool
revCurves)
                | Bool
otherwise ->
                    let (CubicBezier Double
ql, CubicBezier Double
qr) = forall a (b :: * -> *).
(Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> (b a, b a)
splitBezier CubicBezier Double
q Double
0.5
                        half_u :: Double
half_u = Double
umin forall a. Num a => a -> a -> a
+ (Double
umax forall a. Num a => a -> a -> a
- Double
umin) forall a. Fractional a => a -> a -> a
/ Double
2
                    in forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip CubicBezier Double
ql CubicBezier Double
newP Double
umin Double
half_u
                       Double
new_tmin Double
new_tmax Double
newClip Double
pEps Double
vEps (Bool -> Bool
not Bool
revCurves) forall a. [a] -> [a] -> [a]
++
                       CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip CubicBezier Double
qr CubicBezier Double
newP Double
half_u Double
umax Double
new_tmin Double
new_tmax
                       Double
newClip Double
pEps Double
vEps (Bool -> Bool
not Bool
revCurves)
      -- iterate with the curves swapped.
     | Bool
otherwise ->
        forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip CubicBezier Double
q CubicBezier Double
newP Double
umin Double
umax Double
new_tmin
        Double
new_tmax Double
newClip Double
pEps Double
vEps (Bool -> Bool
not Bool
revCurves)

minEps :: Double
minEps :: Double
minEps = Double
1e-8

-- | Find the intersections between two Bezier curves, using the
-- Bezier Clip algorithm. Returns the parameters for both curves.
bezierIntersection :: CubicBezier Double -> CubicBezier Double -> Double -> [(Double, Double)]
bezierIntersection :: CubicBezier Double
-> CubicBezier Double -> Double -> [(Double, Double)]
bezierIntersection CubicBezier Double
p CubicBezier Double
q Double
vEps = CubicBezier Double
-> CubicBezier Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> [(Double, Double)]
bezierClip CubicBezier Double
p CubicBezier Double
q Double
0 Double
1 Double
0 Double
1 Double
0 Double
eps2 Double
vEps Bool
False
  where eps2 :: Double
eps2 = forall a. Ord a => a -> a -> a
max (forall a. Ord a => a -> a -> a
min (forall (b :: * -> *).
GenericBezier b =>
b Double -> Double -> Double
bezierParamTolerance CubicBezier Double
p Double
vEps) (forall (b :: * -> *).
GenericBezier b =>
b Double -> Double -> Double
bezierParamTolerance CubicBezier Double
q Double
vEps)) Double
minEps

-- TODO:
-- following curve generate very large list of intersections
-- let b1 =  CubicBezier {cubicC0 = Point 365.70000000000005 477.40000000000003, cubicC1 = Point 373.3 476.70000000000005, cubicC2 = Point 381.1 476.3, cubicC3 = Point 389.20000000000005 476.3};
--     b2 = CubicBezier {cubicC0 = Point 365.70000000000005 477.40000000000003, cubicC1 = Point 365.70000000000005 476.6, cubicC2 = Point 365.70000000000005 475.8, cubicC3 = Point 365.70000000000005 475.0}

------------------------ Line intersection -------------------------------------
-- Clipping a line uses a simplified version of the Bezier Clip algorithm,
-- and uses the (thin) line itself instead of the fat line.

-- | Find the zero of a 1D bezier 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 :: BernsteinPoly Double -- ^ the bernstein coefficients of the polynomial
               -> Double  -- ^ The lower bound of the interval 
               -> Double  -- ^ The upper bound of the interval
               -> Double  -- ^ The accuracy
               -> [Double] -- ^ The roots found
bezierFindRoot :: BernsteinPoly Double -> Double -> Double -> Double -> [Double]
bezierFindRoot BernsteinPoly Double
p Double
tmin Double
tmax Double
eps
  -- no intersection
  | forall a. Maybe a -> Bool
isNothing Maybe (Double, Double)
chop_interval = []

  -- not enough reduction, so split the curve in case we have
  -- multiple intersections
  | Double
clip forall a. Ord a => a -> a -> Bool
> Double
0.8 =
    let (BernsteinPoly Double
p1, BernsteinPoly Double
p2) = forall a.
(Unbox a, Num a) =>
BernsteinPoly a -> a -> (BernsteinPoly a, BernsteinPoly a)
bernsteinSplit BernsteinPoly Double
newP Double
0.5
        half_t :: Double
half_t = Double
new_tmin forall a. Num a => a -> a -> a
+ (Double
new_tmax forall a. Num a => a -> a -> a
- Double
new_tmin) forall a. Fractional a => a -> a -> a
/ Double
2
    in BernsteinPoly Double -> Double -> Double -> Double -> [Double]
bezierFindRoot BernsteinPoly Double
p1 Double
new_tmin Double
half_t Double
eps forall a. [a] -> [a] -> [a]
++
       BernsteinPoly Double -> Double -> Double -> Double -> [Double]
bezierFindRoot BernsteinPoly Double
p2 Double
half_t Double
new_tmax Double
eps

  -- within tolerance
  | Double
new_tmax forall a. Num a => a -> a -> a
- Double
new_tmin forall a. Ord a => a -> a -> Bool
< Double
eps =
      [Double
new_tmin forall a. Num a => a -> a -> a
+ (Double
new_tmaxforall a. Num a => a -> a -> a
-Double
new_tmin)forall a. Fractional a => a -> a -> a
/Double
2]

      -- iterate
  | Bool
otherwise =
        BernsteinPoly Double -> Double -> Double -> Double -> [Double]
bezierFindRoot BernsteinPoly Double
newP Double
new_tmin Double
new_tmax Double
eps

  where
    chop_interval :: Maybe (Double, Double)
chop_interval = Double -> Double -> [Double] -> Maybe (Double, Double)
chopHull Double
0 Double
0 (forall a. Unbox a => Vector a -> [a]
V.toList forall a b. (a -> b) -> a -> b
$ forall a. BernsteinPoly a -> Vector a
bernsteinCoeffs BernsteinPoly Double
p)
    Just (Double
chop_tmin, Double
chop_tmax) = Maybe (Double, Double)
chop_interval
    newP :: BernsteinPoly Double
newP = forall a.
(Unbox a, Ord a, Fractional a) =>
BernsteinPoly a -> a -> a -> BernsteinPoly a
bernsteinSubsegment BernsteinPoly Double
p Double
chop_tmin Double
chop_tmax
    clip :: Double
clip = Double
chop_tmax forall a. Num a => a -> a -> a
- Double
chop_tmin
    new_tmin :: Double
new_tmin = Double
tmax forall a. Num a => a -> a -> a
* Double
chop_tmin forall a. Num a => a -> a -> a
+ Double
tmin forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
chop_tmin)
    new_tmax :: Double
new_tmax = Double
tmax forall a. Num a => a -> a -> a
* Double
chop_tmax forall a. Num a => a -> a -> a
+ Double
tmin forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
chop_tmax)

-- | Find the intersections of the curve with a line.

-- Apply a transformation to the bezier that maps the line onto the
-- X-axis.  Then we only need to test the Y-values for a zero.
bezierLineIntersections :: CubicBezier Double -> Line Double -> Double -> [Double]
bezierLineIntersections :: CubicBezier Double -> Line Double -> Double -> [Double]
bezierLineIntersections CubicBezier Double
b (Line Point Double
p Point Double
q) Double
eps =
  forall a. (a -> Bool) -> [a] -> [a]
filter (\Double
x -> Double
x forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
x forall a. Ord a => a -> a -> Bool
< Double
1) forall a b. (a -> b) -> a -> b
$
  Double -> Double -> Double -> Double -> [Double]
cubicRoot (Double
p3 forall a. Num a => a -> a -> a
- Double
3forall a. Num a => a -> a -> a
*Double
p2 forall a. Num a => a -> a -> a
+ Double
3forall a. Num a => a -> a -> a
*Double
p1 forall a. Num a => a -> a -> a
- Double
p0) (Double
3forall a. Num a => a -> a -> a
*Double
p2 forall a. Num a => a -> a -> a
- Double
6forall a. Num a => a -> a -> a
*Double
p1 forall a. Num a => a -> a -> a
+ Double
3forall a. Num a => a -> a -> a
*Double
p0) (Double
3forall a. Num a => a -> a -> a
*Double
p1 forall a. Num a => a -> a -> a
- Double
3forall a. Num a => a -> a -> a
*Double
p0) Double
p0
  where (CubicBezier (Point Double
p0 Double
_) (Point Double
p1 Double
_) (Point Double
p2 Double
_) (Point Double
p3 Double
_)) = 
          forall a. HasCallStack => Maybe a -> a
fromJust (forall a.
(Eq a, Fractional a) =>
Transform a -> Maybe (Transform a)
inverse forall a b. (a -> b) -> a -> b
$ forall a. Num a => Point a -> Transform a
translate Point Double
p forall a b. AffineTransform a b => Transform b -> a -> a
$* forall a. Floating a => Point a -> Transform a
rotateVec (Point Double
q forall v. AdditiveGroup v => v -> v -> v
^-^ Point Double
p)) forall a b. AffineTransform a b => Transform b -> a -> a
$* CubicBezier Double
b

-- let cb = (CubicBezier (Point 0 0) (Point 3 4) (Point 10 4) (Point 31 2)); cb1 = fst (splitBezier cb 0.83242); cb2 = CubicBezier {bezierC0 = Point 4.542593123258268 2.7028033902052537, bezierC1 = Point 9.036628467934 3.788306467438, bezierC2 = Point 16.832161 3.4493180000000002, bezierC3 = Point 31.0 2.0}
-- bezierIntersection (CubicBezier (Point 0 0) (Point 3 4) (Point 10 4) (Point 31 2)) (CubicBezier (Point 0 0) (Point 6 8) (Point 2 42) (Point 4 15)) 1e-10