{-# LANGUAGE BangPatterns, MultiWayIf #-}
module Geom2D.CubicBezier.Approximate
       (approximatePath, approximateQuadPath, approximatePathMax, approximateQuadPathMax,
        approximateCubic)
       where
import Geom2D
import Geom2D.CubicBezier.Basic
import Geom2D.CubicBezier.Numeric
import Data.Maybe
import Data.List
import qualified Data.Vector.Unboxed as V
import qualified Data.Map as M
import Data.Function

interpolate :: (Num a) => a -> a -> a -> a
interpolate :: forall a. Num a => a -> a -> a -> a
interpolate a
a a
b a
x = (a
1forall a. Num a => a -> a -> a
-a
x)forall a. Num a => a -> a -> a
*a
a forall a. Num a => a -> a -> a
+ a
xforall a. Num a => a -> a -> a
*a
b

-- | Approximate a function with piecewise cubic bezier splines using
-- a least-squares fit, within the given tolerance.  Each subcurve is
-- approximated by using a finite number of samples.  It is recommended
-- to avoid changes in direction by subdividing the original function
-- at points of inflection.
approximatePath :: (V.Unbox a, Ord a, Floating a) =>
                   (a -> (Point a, Point a)) -- ^ The function to approximate and it's derivative
                -> Int
                   -- ^ The number of discrete samples taken to
                   -- approximate each subcurve.  More samples are
                   -- more precise but take more time to calculate.
                   -- For good precision 16 is a good candidate.
                -> a                         -- ^ The tolerance
                -> a                         -- ^ The lower parameter of the function      
                -> a                         -- ^ The upper parameter of the function
                -> Bool
                -- ^ Calculate the result faster, but with more
                -- subcurves.  Runs typically 10 times faster, but
                -- generates 50% more subcurves.   Useful for interactive use.
                -> [CubicBezier a]
approximatePath :: forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> Int -> a -> a -> a -> Bool -> [CubicBezier a]
approximatePath a -> (Point a, Point a)
f Int
n a
tol a
tmin a
tmax Bool
fast
  | a
err forall a. Ord a => a -> a -> Bool
< a
tol = [CubicBezier a
curve]
  | Bool
otherwise = forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> Int -> a -> a -> a -> Bool -> [CubicBezier a]
approximatePath' a -> (Point a, Point a)
f Int
n a
tol a
tmin a
tmax Bool
fast
  where
    (CubicBezier a
curve, a
err) = forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> (a -> (Point a, Point a)) -> a -> a -> Int -> (CubicBezier a, a)
approx1cubic Int
n a -> (Point a, Point a)
f a
tmin a
tmax (if Bool
fast then Int
0 else Int
5)
{-# SPECIALIZE approximatePath :: (Double -> (DPoint, DPoint)) -> Int -> Double
                               -> Double -> Double -> Bool -> [CubicBezier Double]  #-}

-- | Approximate a function with piecewise quadratic bezier splines
-- using a least-squares fit, within the given tolerance.  It is
-- recommended to avoid changes in direction by subdividing the
-- original function at points of inflection.
approximateQuadPath :: (Show a, V.Unbox a, Ord a, Floating a) =>
                   (a -> (Point a, Point a)) -- ^ The function to approximate and it's derivative
                -> a                         -- ^ The tolerance
                -> a                         -- ^ The lower parameter of the function      
                -> a                         -- ^ The upper parameter of the function
                -> Bool                      
                -- ^ Calculate the result faster, but with more
                -- subcurves.
                -> [QuadBezier a]
approximateQuadPath :: forall a.
(Show a, Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> a -> Bool -> [QuadBezier a]
approximateQuadPath a -> (Point a, Point a)
f a
tol a
tmin a
tmax Bool
fast
  | a
err forall a. Ord a => a -> a -> Bool
< a
tol = [QuadBezier a
curve]
  | Bool
otherwise = forall a.
(Show a, Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> a -> Bool -> [QuadBezier a]
approximateQuad' a -> (Point a, Point a)
f a
tol a
tmin a
tmax Bool
fast
  where
    curve :: QuadBezier a
curve = forall a.
(Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> QuadBezier a
approx1quad a -> (Point a, Point a)
f a
tmin a
tmax
    err :: a
err = forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> QuadBezier a -> a -> a -> Int -> a
maxDist a -> (Point a, Point a)
f QuadBezier a
curve a
tmin a
tmax (if Bool
fast then Int
0 else Int
5)
{-# SPECIALIZE approximateQuadPath :: (Double -> (DPoint, DPoint)) -> Double ->
    Double -> Double -> Bool -> [QuadBezier Double] #-}
      
-- find the distance between the function at t and the quadratic bezier.
-- calculate the value and derivative at t, and improve the closeness of t.
quadDist :: (V.Unbox a, Ord a, Floating a) =>
            (a -> (Point a, Point a)) -> QuadBezier a -> a
         -> a -> Int -> a -> a
quadDist :: forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> QuadBezier a -> a -> a -> Int -> a -> a
quadDist a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax Int
maxiter a
t =
  forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> QuadBezier a
-> a
-> a
-> a
-> Point a
-> Int
-> (Point a, Point a)
-> a
quadDist' a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax a
t (forall a b. (a, b) -> a
fst (a -> (Point a, Point a)
f forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> a -> a -> a
interpolate a
tmin a
tmax a
t)) Int
maxiter (forall a (b :: * -> *).
(Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> (Point a, Point a)
evalBezierDeriv QuadBezier a
qb a
t)

quadDist' :: (V.Unbox a, Ord a, Floating a) =>
             (a -> (Point a, Point a))
          -> QuadBezier a -> a -> a -> a -> Point a -> Int -> (Point a, Point a) -> a
quadDist' :: forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> QuadBezier a
-> a
-> a
-> a
-> Point a
-> Int
-> (Point a, Point a)
-> a
quadDist' a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax a
t Point a
p Int
maxiter (Point a
b, Point a
b')
  | Int
maxiter forall a. Ord a => a -> a -> Bool
<= Int
1 Bool -> Bool -> Bool
|| forall a. Num a => a -> a
abs (a
err2forall a. Num a => a -> a -> a
-a
err1) forall a. Ord a => a -> a -> Bool
<= a
err1 forall a. Num a => a -> a -> a
* (a
1forall a. Fractional a => a -> a -> a
/a
8) = a
err2
  | Bool
otherwise = forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> QuadBezier a
-> a
-> a
-> a
-> Point a
-> Int
-> (Point a, Point a)
-> a
quadDist' a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax (a
tforall a. Num a => a -> a -> a
+a
ndist) Point a
p (Int
maxiterforall a. Num a => a -> a -> a
-Int
1) (Point a
b2, Point a
b2')
  where dp :: Point a
dp = Point a
p forall v. AdditiveGroup v => v -> v -> v
^-^ Point a
b
        err1 :: a
err1 = forall a. Floating a => Point a -> a
vectorMag Point a
dp
        -- distance from p to the normal at b(t) / velocity at t
        ndist :: a
ndist = (Point a
dp forall a. Num a => Point a -> Point a -> a
^.^ Point a
b') forall a. Fractional a => a -> a -> a
/ (Point a
b' forall a. Num a => Point a -> Point a -> a
^.^ Point a
b')
        (Point a
b2, Point a
b2') = forall a (b :: * -> *).
(Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> (Point a, Point a)
evalBezierDeriv QuadBezier a
qb (a
tforall a. Num a => a -> a -> a
+a
ndist)
        err2 :: a
err2 = forall a. Floating a => Point a -> Point a -> a
vectorDistance Point a
p Point a
b2
        
-- find maximum distance using golden section search
maxDist :: (V.Unbox a, Ord a, Floating a) =>
           (a -> (Point a, Point a)) ->
           QuadBezier a -> a -> a -> Int -> a
maxDist :: forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> QuadBezier a -> a -> a -> Int -> a
maxDist a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax Int
maxiter =
  forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> QuadBezier a -> a -> a -> Int -> a -> a
quadDist a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax Int
maxiter forall a b. (a -> b) -> a -> b
$
  forall a. (Ord a, Floating a) => (a -> a) -> Int -> a
goldSearch (forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> QuadBezier a -> a -> a -> Int -> a -> a
quadDist a -> (Point a, Point a)
f QuadBezier a
qb a
tmin a
tmax Int
maxiter) Int
4

approxquad :: (Ord a, Floating a) =>
              Point a -> Point a -> Point a -> Point a -> QuadBezier a
approxquad :: forall a.
(Ord a, Floating a) =>
Point a -> Point a -> Point a -> Point a -> QuadBezier a
approxquad Point a
p0 Point a
p0' Point a
p1' Point a
p1
  | forall a. Num a => a -> a
abs (forall a. Point a -> a
pointY Point a
q') forall a. Ord a => a -> a -> Bool
< forall a. Num a => a -> a
abs (forall a. Point a -> a
pointX Point a
q'forall a. Num a => a -> a -> a
*a
1e-3) = 
    forall a. Point a -> Point a -> Point a -> QuadBezier a
QuadBezier Point a
p0 (forall a. Num a => Point a -> Point a -> a -> Point a
interpolateVector Point a
p0 Point a
p1 a
0.5) Point a
p1
  | Bool
otherwise = forall a. Point a -> Point a -> Point a -> QuadBezier a
QuadBezier Point a
p0 (Point a
p1forall v. AdditiveGroup v => v -> v -> v
^+^Point a
p1'forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*a
t) Point a
p1
  where
    q :: Point a
q = forall a. Floating a => Point a -> Transform a
rotateVec (forall a. Num a => Point a -> Point a
flipVector Point a
p0') forall a b. AffineTransform a b => Transform b -> a -> a
$* Point a
p1forall v. AdditiveGroup v => v -> v -> v
^-^Point a
p0
    q' :: Point a
q' = forall a. Floating a => Point a -> Transform a
rotateVec (forall a. Num a => Point a -> Point a
flipVector Point a
p0') forall a b. AffineTransform a b => Transform b -> a -> a
$* Point a
p1'
    t :: a
t = - forall a. Point a -> a
pointY Point a
q forall a. Fractional a => a -> a -> a
/ forall a. Point a -> a
pointY Point a
q'

approx1quad :: (Ord a, Floating a) =>
               (a -> (Point a, Point a)) -> a -> a -> QuadBezier a
approx1quad :: forall a.
(Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> QuadBezier a
approx1quad a -> (Point a, Point a)
f a
tmin a
tmax =
  forall a.
(Ord a, Floating a) =>
Point a -> Point a -> Point a -> Point a -> QuadBezier a
approxquad Point a
p0 Point a
p0' Point a
p1' Point a
p1
  where (Point a
p0, Point a
p0') = a -> (Point a, Point a)
f a
tmin
        (Point a
p1, Point a
p1') = a -> (Point a, Point a)
f a
tmax

splitQuad :: (Show a, V.Unbox a, Ord a, Floating a) =>
                a -> a -> (a -> (Point a, Point a))
                -> a -> a -> Int -> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad :: forall a.
(Show a, Unbox a, Ord a, Floating a) =>
a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad a
node a
offset a -> (Point a, Point a)
f a
tmin a
tmax Int
maxiter =
  forall a.
(Show a, Unbox a, Ord a, Floating a) =>
a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> Int
-> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad' a
node a
offset a -> (Point a, Point a)
f a
tmin a
tmax Int
maxiter Int
maxiter
  
splitQuad' :: (Show a, V.Unbox a, Ord a, Floating a) =>
                a -> a -> (a -> (Point a, Point a))
                -> a -> a -> Int -> Int -> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad' :: forall a.
(Show a, Unbox a, Ord a, Floating a) =>
a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> Int
-> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad' a
node a
offset a -> (Point a, Point a)
f a
tmin a
tmax Int
maxiter Int
maxiter2
  | Int
maxiter forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| (a
err0 forall a. Ord a => a -> a -> Bool
< a
2forall a. Num a => a -> a -> a
*a
err1 Bool -> Bool -> Bool
&& a
err0 forall a. Ord a => a -> a -> Bool
> a
err1forall a. Fractional a => a -> a -> a
/a
2) =
      (a
tmid, a
err0, QuadBezier a
curve0, a
err1, QuadBezier a
curve1)
  | Bool
otherwise =
    forall a.
(Show a, Unbox a, Ord a, Floating a) =>
a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> Int
-> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad' (if a
err0 forall a. Ord a => a -> a -> Bool
< a
err1 then a
nodeforall a. Num a => a -> a -> a
+a
offset else a
nodeforall a. Num a => a -> a -> a
-a
offset)
    (a
offsetforall a. Fractional a => a -> a -> a
/a
2) a -> (Point a, Point a)
f a
tmin a
tmax (Int
maxiterforall a. Num a => a -> a -> a
-Int
1) Int
maxiter2
  where
    tmid :: a
tmid = forall a. Num a => a -> a -> a -> a
interpolate a
tmin a
tmax a
node
    curve0 :: QuadBezier a
curve0 = forall a.
(Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> QuadBezier a
approx1quad a -> (Point a, Point a)
f a
tmin a
tmid 
    err0 :: a
err0 = forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> QuadBezier a -> a -> a -> Int -> a
maxDist a -> (Point a, Point a)
f QuadBezier a
curve0 a
tmin a
tmid Int
maxiter2
    curve1 :: QuadBezier a
curve1 = forall a.
(Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> QuadBezier a
approx1quad a -> (Point a, Point a)
f a
tmid a
tmax 
    err1 :: a
err1 = forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> QuadBezier a -> a -> a -> Int -> a
maxDist a -> (Point a, Point a)
f QuadBezier a
curve1 a
tmid a
tmax Int
maxiter2

approximateQuad' :: (Show a, V.Unbox a, Ord a, Floating a) =>
                    (a -> (Point a, Point a)) -> 
                    a -> a -> a -> Bool ->
                    [QuadBezier a]
approximateQuad' :: forall a.
(Show a, Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> a -> Bool -> [QuadBezier a]
approximateQuad' a -> (Point a, Point a)
f a
tol a
tmin a
tmax Bool
fast =
  (if a
err0 forall a. Ord a => a -> a -> Bool
<= a
tol
   then [QuadBezier a
curve0]
   else forall a.
(Show a, Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> a -> Bool -> [QuadBezier a]
approximateQuad' a -> (Point a, Point a)
f a
tol a
tmin a
tmid Bool
fast) forall a. [a] -> [a] -> [a]
++
  (if a
err1 forall a. Ord a => a -> a -> Bool
<= a
tol
   then [QuadBezier a
curve1]
   else forall a.
(Show a, Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> a -> a -> a -> Bool -> [QuadBezier a]
approximateQuad' a -> (Point a, Point a)
f a
tol a
tmid a
tmax Bool
fast)
  where
    (a
tmid, a
err0, QuadBezier a
curve0, a
err1, QuadBezier a
curve1) =
     forall a.
(Show a, Unbox a, Ord a, Floating a) =>
a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad a
0.5 a
0.25 a -> (Point a, Point a)
f a
tmin a
tmax (if Bool
fast then Int
0 else Int
5)

approximatePath' :: (V.Unbox a, Ord a, Floating a) =>
                    (a -> (Point a, Point a)) -> Int ->
                    a -> a -> a -> Bool ->
                    [CubicBezier a]
approximatePath' :: forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> Int -> a -> a -> a -> Bool -> [CubicBezier a]
approximatePath' a -> (Point a, Point a)
f Int
n a
tol a
tmin a
tmax Bool
fast =
  (if a
err0 forall a. Ord a => a -> a -> Bool
<= a
tol
   then [CubicBezier a
curve0]
   else forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> Int -> a -> a -> a -> Bool -> [CubicBezier a]
approximatePath' a -> (Point a, Point a)
f Int
n a
tol a
tmin a
tmid Bool
fast) forall a. [a] -> [a] -> [a]
++
  (if a
err1 forall a. Ord a => a -> a -> Bool
<= a
tol
   then [CubicBezier a
curve1]
   else forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> Int -> a -> a -> a -> Bool -> [CubicBezier a]
approximatePath' a -> (Point a, Point a)
f Int
n a
tol a
tmid a
tmax Bool
fast)
  where
    (a
tmid, a
err0, CubicBezier a
curve0, a
err1, CubicBezier a
curve1) =
      forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, CubicBezier a, a, CubicBezier a)
splitCubic Int
n a
0.5 a
0.25 a -> (Point a, Point a)
f a
tmin a
tmax (if Bool
fast then Int
0 else Int
5)
--{-# SPECIALIZE approximatePath' :: (Double -> (Point Double, Point Double)) -> Int -> Double -> Double -> Double -> [CubicBezier Double]  #-}      

-- | Like approximatePath, but limit the number of subcurves.
approximatePathMax :: (V.Unbox a, Floating a, Ord a) =>
                      Int                        -- ^ The maximum number of subcurves
                   -> (a -> (Point a, Point a))    -- ^ The function to approximate and it's derivative
                   -> Int
                   -- ^ The number of discrete samples taken to
                   -- approximate each subcurve.  More samples are
                   -- more precise but take more time to calculate.
                   -- For good precision 16 is a good candidate.
                   -> a                          -- ^ The tolerance
                   -> a                          -- ^ The lower parameter of the function      
                   -> a                          -- ^ The upper parameter of the function
                   -> Bool
                   -- ^ Calculate the result faster, but with more
                   -- subcurves.  Runs faster (typically 10 times),
                   -- but generates more subcurves (about 50%).
                   -- Useful for interactive use.
                   -> [CubicBezier a]
approximatePathMax :: forall a.
(Unbox a, Floating a, Ord a) =>
Int
-> (a -> (Point a, Point a))
-> Int
-> a
-> a
-> a
-> Bool
-> [CubicBezier a]
approximatePathMax Int
m a -> (Point a, Point a)
f Int
samples a
tol a
tmin a
tmax Bool
fast =
  forall a b.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> a
-> Int
-> Bool
-> (a
    -> a
    -> (a -> (Point a, Point a))
    -> a
    -> a
    -> Int
    -> (a, a, b, a, b))
-> Map a (FunctionSegment b a)
-> [b]
approxMax a -> (Point a, Point a)
f a
tol Int
m Bool
fast (forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, CubicBezier a, a, CubicBezier a)
splitCubic Int
samples) Map a (FunctionSegment (CubicBezier a) a)
segments
  where segments :: Map a (FunctionSegment (CubicBezier a) a)
segments = forall k a. k -> a -> Map k a
M.singleton a
err (forall b a. a -> a -> b -> FunctionSegment b a
FunctionSegment a
tmin a
tmax CubicBezier a
outline)
        (Point a
p0, Point a
p0') = a -> (Point a, Point a)
f a
tmin
        (Point a
p1, Point a
p1') = a -> (Point a, Point a)
f a
tmax
        ts :: Vector a
ts = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\Int
i -> forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
iforall a. Fractional a => a -> a -> a
/(forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
samplesforall a. Num a => a -> a -> a
+a
1) forall a. a -> a -> a
`asTypeOf` a
tmin) forall a b. (a -> b) -> a -> b
$
             forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN (Int
1::Int) Int
samples
        points :: Vector (Point a)
points = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> (Point a, Point a)
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a -> a -> a
interpolate a
tmin a
tmax) Vector a
ts
        curveCb :: CubicBezier a
curveCb = forall a. Point a -> Point a -> Point a -> Point a -> CubicBezier a
CubicBezier Point a
p0 (Point a
p0forall v. AdditiveGroup v => v -> v -> v
^+^Point a
p0') (Point a
p1forall v. AdditiveGroup v => v -> v -> v
^-^Point a
p1') Point a
p1
        (CubicBezier a
outline, a
err) =
          forall a.
(Unbox a, Ord a, Floating a) =>
CubicBezier a
-> Vector (Point a)
-> Maybe (Vector a)
-> Int
-> (CubicBezier a, a)
approximateCubic CubicBezier a
curveCb Vector (Point a)
points (forall a. a -> Maybe a
Just Vector a
ts) (if Bool
fast then Int
0 else Int
5)
{-# SPECIALIZE approximatePathMax ::
    Int -> (Double -> (Point Double, Point Double)) -> Int                      
    -> Double -> Double -> Double -> Bool -> [CubicBezier Double] #-}
data FunctionSegment b a = FunctionSegment {
  forall b a. FunctionSegment b a -> a
fsTmin :: !a,  -- the least t param of the segment in the original curve
  forall b a. FunctionSegment b a -> a
_fsTmax :: !a,  -- the max t param of the segment in the original curve
  forall b a. FunctionSegment b a -> b
fsCurve :: b -- the curve segment
  }

-- | Like approximateQuadPath, but limit the number of subcurves.
approximateQuadPathMax :: (V.Unbox a, Show a, Floating a, Ord a) =>
                      Int                        -- ^ The maximum number of subcurves
                   -> (a -> (Point a, Point a))    -- ^ The function to approximate and it's derivative
                   -> a                          -- ^ The tolerance
                   -> a                          -- ^ The lower parameter of the function      
                   -> a                          -- ^ The upper parameter of the function
                   -> Bool
                   -- ^ Calculate the result faster, but with more
                   -- subcurves.  Runs faster, but generates more
                   -- subcurves.  Useful for interactive use.
                   -> [QuadBezier a]
approximateQuadPathMax :: forall a.
(Unbox a, Show a, Floating a, Ord a) =>
Int
-> (a -> (Point a, Point a))
-> a
-> a
-> a
-> Bool
-> [QuadBezier a]
approximateQuadPathMax Int
m a -> (Point a, Point a)
f a
tol a
tmin a
tmax Bool
fast =
  forall a b.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> a
-> Int
-> Bool
-> (a
    -> a
    -> (a -> (Point a, Point a))
    -> a
    -> a
    -> Int
    -> (a, a, b, a, b))
-> Map a (FunctionSegment b a)
-> [b]
approxMax a -> (Point a, Point a)
f a
tol Int
m Bool
fast forall a.
(Show a, Unbox a, Ord a, Floating a) =>
a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, QuadBezier a, a, QuadBezier a)
splitQuad Map a (FunctionSegment (QuadBezier a) a)
segments
  where segments :: Map a (FunctionSegment (QuadBezier a) a)
segments = forall k a. k -> a -> Map k a
M.singleton a
err (forall b a. a -> a -> b -> FunctionSegment b a
FunctionSegment a
tmin a
tmax QuadBezier a
curveQd)
        (Point a
p0, Point a
p0') = a -> (Point a, Point a)
f a
tmin
        (Point a
p1, Point a
p1') = a -> (Point a, Point a)
f a
tmax
        curveQd :: QuadBezier a
curveQd = forall a.
(Ord a, Floating a) =>
Point a -> Point a -> Point a -> Point a -> QuadBezier a
approxquad Point a
p0 Point a
p0' Point a
p1' Point a
p1
        err :: a
err = forall a.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a)) -> QuadBezier a -> a -> a -> Int -> a
maxDist a -> (Point a, Point a)
f QuadBezier a
curveQd a
tmin a
tmax (if Bool
fast then Int
0 else Int
5)
{-# SPECIALIZE approximateQuadPathMax ::
    Int -> (Double -> (Point Double, Point Double))
    -> Double -> Double -> Double -> Bool -> [QuadBezier Double] #-}
-- Keep a map from maxError to FunctionSegment for each subsegment to keep
-- track of the segment with the maximum error.  This ensures a n
-- log(n) execution time, rather than n^2 when a list is used.
approxMax :: (V.Unbox a, Ord a, Floating a) =>
             (a -> (Point a, Point a)) -> a -> Int -> Bool
          -> (a -> a -> (a -> (Point a, Point a)) -> a -> a -> Int -> (a, a, b, a, b))
          -> M.Map a (FunctionSegment b a)
          -> [b]
approxMax :: forall a b.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> a
-> Int
-> Bool
-> (a
    -> a
    -> (a -> (Point a, Point a))
    -> a
    -> a
    -> Int
    -> (a, a, b, a, b))
-> Map a (FunctionSegment b a)
-> [b]
approxMax a -> (Point a, Point a)
f a
tol Int
maxCurves Bool
fast a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, b, a, b)
splitBez Map a (FunctionSegment b a)
segments
  | (Int
maxCurves forall a. Ord a => a -> a -> Bool
<= Int
1) Bool -> Bool -> Bool
|| (a
err forall a. Ord a => a -> a -> Bool
< a
tol) =
    forall a b. (a -> b) -> [a] -> [b]
map forall b a. FunctionSegment b a -> b
fsCurve forall a b. (a -> b) -> a -> b
$ forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (forall a. Ord a => a -> a -> Ordering
compare forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` forall b a. FunctionSegment b a -> a
fsTmin) forall a b. (a -> b) -> a -> b
$
    forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall k a. Map k a -> [(k, a)]
M.toList Map a (FunctionSegment b a)
segments
  | Bool
otherwise = forall a b.
(Unbox a, Ord a, Floating a) =>
(a -> (Point a, Point a))
-> a
-> Int
-> Bool
-> (a
    -> a
    -> (a -> (Point a, Point a))
    -> a
    -> a
    -> Int
    -> (a, a, b, a, b))
-> Map a (FunctionSegment b a)
-> [b]
approxMax a -> (Point a, Point a)
f a
tol (Int
maxCurvesforall a. Num a => a -> a -> a
-Int
1) Bool
fast a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, b, a, b)
splitBez forall a b. (a -> b) -> a -> b
$
                forall k a. Ord k => k -> a -> Map k a -> Map k a
M.insert a
err_l (forall b a. a -> a -> b -> FunctionSegment b a
FunctionSegment a
t_min a
t_mid b
curve_l) forall a b. (a -> b) -> a -> b
$
                forall k a. Ord k => k -> a -> Map k a -> Map k a
M.insert a
err_r (forall b a. a -> a -> b -> FunctionSegment b a
FunctionSegment a
t_mid a
t_max b
curve_r)
                Map a (FunctionSegment b a)
newSegments
  where
    ((a
err, FunctionSegment a
t_min a
t_max b
_), Map a (FunctionSegment b a)
newSegments) =
      forall k a. Map k a -> ((k, a), Map k a)
M.deleteFindMax Map a (FunctionSegment b a)
segments
    (a
t_mid, a
err_l, b
curve_l, a
err_r, b
curve_r) =
      a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, b, a, b)
splitBez a
0.5 a
0.25 a -> (Point a, Point a)
f a
t_min a
t_max (if Bool
fast then Int
0 else Int
5)
      
splitCubic :: (V.Unbox a, Ord a, Floating a) =>
                Int -> a -> a -> (a -> (Point a, Point a))
                -> a -> a -> Int -> (a, a, CubicBezier a, a, CubicBezier a)
splitCubic :: forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, CubicBezier a, a, CubicBezier a)
splitCubic Int
n a
node a
offset a -> (Point a, Point a)
f a
tmin a
tmax Int
maxiter
  | Int
maxiter forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| (a
err0 forall a. Ord a => a -> a -> Bool
< a
2forall a. Num a => a -> a -> a
*a
err1 Bool -> Bool -> Bool
&& a
err0 forall a. Ord a => a -> a -> Bool
> a
err1forall a. Fractional a => a -> a -> a
/a
2) =
      (a
tmid, a
err0, CubicBezier a
curve0, a
err1, CubicBezier a
curve1)
  | Bool
otherwise = 
      forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> a
-> a
-> (a -> (Point a, Point a))
-> a
-> a
-> Int
-> (a, a, CubicBezier a, a, CubicBezier a)
splitCubic Int
n (if a
err0 forall a. Ord a => a -> a -> Bool
< a
err1 then a
nodeforall a. Num a => a -> a -> a
+a
offset else a
nodeforall a. Num a => a -> a -> a
-a
offset)
      (a
offsetforall a. Fractional a => a -> a -> a
/a
2) a -> (Point a, Point a)
f a
tmin a
tmax (Int
maxiterforall a. Num a => a -> a -> a
-Int
1)
  where
    tmid :: a
tmid = forall a. Num a => a -> a -> a -> a
interpolate a
tmin a
tmax a
node
    (CubicBezier a
curve0, a
err0) = forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> (a -> (Point a, Point a)) -> a -> a -> Int -> (CubicBezier a, a)
approx1cubic Int
n a -> (Point a, Point a)
f a
tmin a
tmid Int
maxiter
    (CubicBezier a
curve1, a
err1) = forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> (a -> (Point a, Point a)) -> a -> a -> Int -> (CubicBezier a, a)
approx1cubic Int
n a -> (Point a, Point a)
f a
tmid a
tmax Int
maxiter
    
approx1cubic :: (V.Unbox a, Ord a, Floating a) =>
           Int -> (a -> (Point a, Point a)) -> a -> a ->
           Int -> (CubicBezier a, a)
approx1cubic :: forall a.
(Unbox a, Ord a, Floating a) =>
Int
-> (a -> (Point a, Point a)) -> a -> a -> Int -> (CubicBezier a, a)
approx1cubic Int
n a -> (Point a, Point a)
f a
t0 a
t1 Int
maxiter =
  forall a.
(Unbox a, Ord a, Floating a) =>
CubicBezier a
-> Vector (Point a)
-> Maybe (Vector a)
-> Int
-> (CubicBezier a, a)
approximateCubic CubicBezier a
curveCb Vector (Point a)
points (forall a. a -> Maybe a
Just Vector a
ts) Int
maxiter
  where (Point a
p0, Point a
p0') = a -> (Point a, Point a)
f a
t0
        (Point a
p1, Point a
p1') = a -> (Point a, Point a)
f a
t1
        ts :: Vector a
ts = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\Int
i -> forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
iforall a. Fractional a => a -> a -> a
/(forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nforall a. Num a => a -> a -> a
+a
1))
             (forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Int
1 Int
n :: V.Vector Int)
        points :: Vector (Point a)
points = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> (Point a, Point a)
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a -> a -> a
interpolate a
t0 a
t1) Vector a
ts
        curveCb :: CubicBezier a
curveCb = forall a. Point a -> Point a -> Point a -> Point a -> CubicBezier a
CubicBezier Point a
p0 (Point a
p0forall v. AdditiveGroup v => v -> v -> v
^+^Point a
p0') (Point a
p1forall v. AdditiveGroup v => v -> v -> v
^+^Point a
p1') Point a
p1
{-# SPECIALIZE approx1cubic ::  Int -> (Double -> (Point Double, Point Double)) -> Double -> Double -> Int -> (CubicBezier Double, Double) #-}

-- | @approximateCubic b pts maxiter@ finds the least squares fit of a bezier
-- curve to the points @pts@.  The resulting bezier has the same first
-- and last control point as the curve @b@, and have tangents colinear with @b@.
approximateCubic :: (V.Unbox a, Ord a, Floating a) =>
                    CubicBezier a         -- ^ Curve
                    -> V.Vector (Point a) -- ^ Points
                    -> Maybe (V.Vector a) -- ^ Params.  Approximate if Nothing
                    -> Int                -- ^ Maximum iterations
                    -> (CubicBezier a, a) -- ^ result curve and maximum error
approximateCubic :: forall a.
(Unbox a, Ord a, Floating a) =>
CubicBezier a
-> Vector (Point a)
-> Maybe (Vector a)
-> Int
-> (CubicBezier a, a)
approximateCubic CubicBezier a
curve Vector (Point a)
pts Maybe (Vector a)
mbTs Int
maxiter =
  let ts :: Vector a
ts = forall a. a -> Maybe a -> a
fromMaybe (forall a.
(Unbox a, Floating a) =>
Point a -> Point a -> Vector (Point a) -> Vector a
approximateParams (forall a. CubicBezier a -> Point a
cubicC0 CubicBezier a
curve) (forall a. CubicBezier a -> Point a
cubicC3 CubicBezier a
curve) Vector (Point a)
pts) Maybe (Vector a)
mbTs
      curve2 :: CubicBezier a
curve2 = forall a. a -> Maybe a -> a
fromMaybe CubicBezier a
curve forall a b. (a -> b) -> a -> b
$ forall a.
(Unbox a, Fractional a, Eq a) =>
CubicBezier a
-> Vector (Point a) -> Vector a -> Maybe (CubicBezier a)
lsqDist CubicBezier a
curve Vector (Point a)
pts Vector a
ts
      (Vector (Point a)
bt, Vector (Point a)
bt') = forall a b.
(Unbox a, Unbox b) =>
Vector (a, b) -> (Vector a, Vector b)
V.unzip forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a (b :: * -> *).
(Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> (Point a, Point a)
evalBezierDeriv CubicBezier a
curve2) Vector a
ts
      err :: a
err = forall a. (Unbox a, Ord a) => Vector a -> a
V.maximum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Floating a => Point a -> Point a -> a
vectorDistance Vector (Point a)
pts Vector (Point a)
bt
      (CubicBezier a
c, Vector a
_, Vector a
_, a
err2, Vector (Point a)
_) =
        forall a. a -> Maybe a -> a
fromMaybe (CubicBezier a
curve2, Vector a
ts, forall a. HasCallStack => a
undefined, a
err, forall a. HasCallStack => a
undefined) forall a b. (a -> b) -> a -> b
$
        forall a.
(Unbox a, Ord a, Floating a) =>
CubicBezier a
-> Vector (Point a)
-> Vector a
-> Int
-> a
-> Vector (Point a)
-> Vector (Point a)
-> Maybe (CubicBezier a, Vector a, Vector a, a, Vector (Point a))
approximateCubic' CubicBezier a
curve2 Vector (Point a)
pts Vector a
ts Int
maxiter a
err Vector (Point a)
bt Vector (Point a)
bt'
  in (CubicBezier a
c, a
err2)
{-# SPECIALIZE approximateCubic :: CubicBezier Double -> V.Vector (Point Double)
  -> Maybe (V.Vector Double) -> Int -> (CubicBezier Double, Double) #-}

-- find (a, b) which minimises ∑ᵢ(a*aᵢ + b*bᵢ + epsᵢ)²
leastSquares :: (V.Unbox a, Fractional a, Eq a) =>
                V.Vector a -> V.Vector a -> V.Vector a -> Maybe (a, a)
leastSquares :: forall a.
(Unbox a, Fractional a, Eq a) =>
Vector a -> Vector a -> Vector a -> Maybe (a, a)
leastSquares Vector a
as Vector a
bs Vector a
epses = forall a.
(Eq a, Fractional a) =>
a -> a -> a -> a -> a -> a -> Maybe (a, a)
solveLinear2x2 a
a a
b a
c a
d a
e a
f
  where
    square :: a -> a
square a
x = a
xforall a. Num a => a -> a -> a
*a
x
    a :: a
a = forall a. (Unbox a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map forall a. Num a => a -> a
square Vector a
as
    b :: a
b = forall a. (Unbox a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(*) Vector a
as Vector a
bs
    c :: a
c = forall a. (Unbox a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(*) Vector a
as Vector a
epses
    d :: a
d = a
b
    e :: a
e = forall a. (Unbox a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map forall a. Num a => a -> a
square Vector a
bs
    f :: a
f = forall a. (Unbox a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(*) Vector a
bs Vector a
epses
{-# SPECIALIZE leastSquares ::V.Vector Double -> V.Vector Double -> V.Vector Double -> Maybe (Double, Double) #-}

-- find the least squares between the points pᵢ and B(tᵢ) for
-- bezier curve B, where pts contains the points pᵢ and ts
-- the values of tᵢ .
-- The tangent at the beginning and end is maintained.
-- Since the start and end point remains the same,
-- we need to find the new value of p2' = p1 + α₁ * (p2 - p1)
-- and p₃' = p4 + α2 * (p3 - p4)
-- minimizing (∑|B(tᵢ) - pᵢ|²) gives a linear equation
-- with two unknown values (α₁ and α₂)
lsqDist :: (V.Unbox a, Fractional a, Eq a) =>
           CubicBezier a
           -> V.Vector (Point a) -> V.Vector a -> Maybe (CubicBezier a)
lsqDist :: forall a.
(Unbox a, Fractional a, Eq a) =>
CubicBezier a
-> Vector (Point a) -> Vector a -> Maybe (CubicBezier a)
lsqDist (CubicBezier (Point !a
p1x !a
p1y) (Point !a
p2x !a
p2y) (Point !a
p3x !a
p3y) (Point !a
p4x !a
p4y)) Vector (Point a)
pts Vector a
ts = let
  calcParams :: a -> Point a -> (a, a, a, a, a, a)
calcParams a
t (Point a
px a
py) = let
    t2 :: a
t2 = a
t forall a. Num a => a -> a -> a
* a
t; t3 :: a
t3 = a
t2 forall a. Num a => a -> a -> a
* a
t
    ax :: a
ax = a
3 forall a. Num a => a -> a -> a
* (a
p2x forall a. Num a => a -> a -> a
- a
p1x) forall a. Num a => a -> a -> a
* (a
t3 forall a. Num a => a -> a -> a
- a
2 forall a. Num a => a -> a -> a
* a
t2 forall a. Num a => a -> a -> a
+ a
t)
    ay :: a
ay = a
3 forall a. Num a => a -> a -> a
* (a
p2y forall a. Num a => a -> a -> a
- a
p1y) forall a. Num a => a -> a -> a
* (a
t3 forall a. Num a => a -> a -> a
- a
2 forall a. Num a => a -> a -> a
* a
t2 forall a. Num a => a -> a -> a
+ a
t)
    bx :: a
bx = a
3 forall a. Num a => a -> a -> a
* (a
p3x forall a. Num a => a -> a -> a
- a
p4x) forall a. Num a => a -> a -> a
* (a
t2 forall a. Num a => a -> a -> a
- a
t3)
    by :: a
by = a
3 forall a. Num a => a -> a -> a
* (a
p3y forall a. Num a => a -> a -> a
- a
p4y) forall a. Num a => a -> a -> a
* (a
t2 forall a. Num a => a -> a -> a
- a
t3)
    cx :: a
cx = (a
p4x forall a. Num a => a -> a -> a
- a
p1x) forall a. Num a => a -> a -> a
* (a
3 forall a. Num a => a -> a -> a
* a
t2 forall a. Num a => a -> a -> a
- a
2 forall a. Num a => a -> a -> a
* a
t3) forall a. Num a => a -> a -> a
+ a
p1x forall a. Num a => a -> a -> a
- a
px
    cy :: a
cy = (a
p4y forall a. Num a => a -> a -> a
- a
p1y) forall a. Num a => a -> a -> a
* (a
3 forall a. Num a => a -> a -> a
* a
t2 forall a. Num a => a -> a -> a
- a
2 forall a. Num a => a -> a -> a
* a
t3) forall a. Num a => a -> a -> a
+ a
p1y forall a. Num a => a -> a -> a
- a
py
    in (a
ax forall a. Num a => a -> a -> a
* a
ax forall a. Num a => a -> a -> a
+ a
ay forall a. Num a => a -> a -> a
* a
ay,
        a
ax forall a. Num a => a -> a -> a
* a
bx forall a. Num a => a -> a -> a
+ a
ay forall a. Num a => a -> a -> a
* a
by,
        a
ax forall a. Num a => a -> a -> a
* a
cx forall a. Num a => a -> a -> a
+ a
ay forall a. Num a => a -> a -> a
* a
cy,
        a
bx forall a. Num a => a -> a -> a
* a
ax forall a. Num a => a -> a -> a
+ a
by forall a. Num a => a -> a -> a
* a
ay,
        a
bx forall a. Num a => a -> a -> a
* a
bx forall a. Num a => a -> a -> a
+ a
by forall a. Num a => a -> a -> a
* a
by,
        a
bx forall a. Num a => a -> a -> a
* a
cx forall a. Num a => a -> a -> a
+ a
by forall a. Num a => a -> a -> a
* a
cy)
  add6 :: (a, b, c, d, e, f) -> (a, b, c, d, e, f) -> (a, b, c, d, e, f)
add6 (!a
a,!b
b,!c
c,!d
d,!e
e,!f
f) (!a
a',!b
b',!c
c',!d
d',!e
e',!f
f') =
    (a
aforall a. Num a => a -> a -> a
+a
a',b
bforall a. Num a => a -> a -> a
+b
b',c
cforall a. Num a => a -> a -> a
+c
c',d
dforall a. Num a => a -> a -> a
+d
d',e
eforall a. Num a => a -> a -> a
+e
e',f
fforall a. Num a => a -> a -> a
+f
f')
  ( a
as, a
bs, a
cs, a
ds, a
es, a
fs ) = forall a. Unbox a => (a -> a -> a) -> Vector a -> a
V.foldl1' forall {a} {b} {c} {d} {e} {f}.
(Num a, Num b, Num c, Num d, Num e, Num f) =>
(a, b, c, d, e, f) -> (a, b, c, d, e, f) -> (a, b, c, d, e, f)
add6 forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith a -> Point a -> (a, a, a, a, a, a)
calcParams Vector a
ts Vector (Point a)
pts
  in do (a
alpha1, a
alpha2) <- forall a.
(Eq a, Fractional a) =>
a -> a -> a -> a -> a -> a -> Maybe (a, a)
solveLinear2x2 a
as a
bs a
cs a
ds a
es a
fs 
        let cp1 :: Point a
cp1 = forall a. a -> a -> Point a
Point (a
alpha1 forall a. Num a => a -> a -> a
* (a
p2x forall a. Num a => a -> a -> a
- a
p1x) forall a. Num a => a -> a -> a
+ a
p1x) (a
alpha1 forall a. Num a => a -> a -> a
* (a
p2y forall a. Num a => a -> a -> a
- a
p1y) forall a. Num a => a -> a -> a
+ a
p1y)
            cp2 :: Point a
cp2 = forall a. a -> a -> Point a
Point (a
alpha2 forall a. Num a => a -> a -> a
* (a
p3x forall a. Num a => a -> a -> a
- a
p4x) forall a. Num a => a -> a -> a
+ a
p4x) (a
alpha2 forall a. Num a => a -> a -> a
* (a
p3y forall a. Num a => a -> a -> a
- a
p4y) forall a. Num a => a -> a -> a
+ a
p4y)
        forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall a. Point a -> Point a -> Point a -> Point a -> CubicBezier a
CubicBezier (forall a. a -> a -> Point a
Point a
p1x a
p1y) Point a
cp1 Point a
cp2 (forall a. a -> a -> Point a
Point a
p4x a
p4y)
{-# SPECIALIZE lsqDist :: CubicBezier Double
           -> V.Vector (Point Double) -> V.Vector Double -> Maybe (CubicBezier Double) #-}

-- calculate the least Squares bezier curve by choosing approximate values
-- of t, and iterating again with an improved estimate of t, by taking the
-- the values of t for which the points are closest to the curve
approximateCubic' :: (V.Unbox a, Ord a, Floating a) =>
                     CubicBezier a
                  -> V.Vector (Point a) -> V.Vector a
                  -> Int -> a -> V.Vector (Point a)
                  -> V.Vector (Point a)
                  -> Maybe (CubicBezier a, V.Vector a, V.Vector a, a, V.Vector (Point a))
approximateCubic' :: forall a.
(Unbox a, Ord a, Floating a) =>
CubicBezier a
-> Vector (Point a)
-> Vector a
-> Int
-> a
-> Vector (Point a)
-> Vector (Point a)
-> Maybe (CubicBezier a, Vector a, Vector a, a, Vector (Point a))
approximateCubic' (CubicBezier Point a
p1 Point a
p2 Point a
p3 Point a
p4) Vector (Point a)
pts Vector a
ts Int
maxiter a
err Vector (Point a)
bt Vector (Point a)
bt' = do
  let dir1 :: Vector (Point a)
dir1 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map ((forall a b. AffineTransform a b => Transform b -> a -> a
$* (Point a
p2forall v. AdditiveGroup v => v -> v -> v
^-^Point a
p1)) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Floating a => Point a -> Transform a
rotateVec forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => Point a -> Point a
flipVector) Vector (Point a)
bt'
      dir2 :: Vector (Point a)
dir2 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map ((forall a b. AffineTransform a b => Transform b -> a -> a
$* (Point a
p3forall v. AdditiveGroup v => v -> v -> v
^-^Point a
p4)) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Floating a => Point a -> Transform a
rotateVec forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => Point a -> Point a
flipVector) Vector (Point a)
bt'
      ps :: Vector (Point a)
ps = forall a b c d.
(Unbox a, Unbox b, Unbox c, Unbox d) =>
(a -> b -> c -> d) -> Vector a -> Vector b -> Vector c -> Vector d
V.zipWith3 (\Point a
b Point a
b' Point a
p ->
                        forall a. Floating a => Point a -> Transform a
rotateVec (forall a. Num a => Point a -> Point a
flipVector Point a
b') forall a b. AffineTransform a b => Transform b -> a -> a
$*
                        (Point a
pforall v. AdditiveGroup v => v -> v -> v
^-^Point a
b)) Vector (Point a)
bt Vector (Point a)
bt' Vector (Point a)
pts
      errs :: Vector a
errs = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a. Num a => a -> a
negateforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Point a -> a
pointY) Vector (Point a)
ps
      as :: Vector a
as = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (\Point a
d a
t -> a
3forall a. Num a => a -> a -> a
*forall a. Point a -> a
pointY Point a
dforall a. Num a => a -> a -> a
*(a
1forall a. Num a => a -> a -> a
-a
t)forall a. Num a => a -> a -> a
*(a
1forall a. Num a => a -> a -> a
-a
t)forall a. Num a => a -> a -> a
*a
t)
           Vector (Point a)
dir1 Vector a
ts
      bs :: Vector a
bs = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (\Point a
d a
t -> a
3forall a. Num a => a -> a -> a
*forall a. Point a -> a
pointY Point a
dforall a. Num a => a -> a -> a
*(a
1forall a. Num a => a -> a -> a
-a
t)forall a. Num a => a -> a -> a
*a
tforall a. Num a => a -> a -> a
*a
t)
           Vector (Point a)
dir2 Vector a
ts
  (a
a,a
b) <- forall a.
(Unbox a, Fractional a, Eq a) =>
Vector a -> Vector a -> Vector a -> Maybe (a, a)
leastSquares Vector a
as Vector a
bs Vector a
errs
  let newTs :: Vector a
newTs = forall a b c d e f.
(Unbox a, Unbox b, Unbox c, Unbox d, Unbox e, Unbox f) =>
(a -> b -> c -> d -> e -> f)
-> Vector a
-> Vector b
-> Vector c
-> Vector d
-> Vector e
-> Vector f
V.zipWith5 (\a
t Point a
p Point a
d1 Point a
d2 Point a
b' ->
                           forall a. Ord a => a -> a -> a
max a
0 forall a b. (a -> b) -> a -> b
$ forall a. Ord a => a -> a -> a
min a
1 forall a b. (a -> b) -> a -> b
$
                           a
t forall a. Num a => a -> a -> a
+ (forall a. Point a -> a
pointX Point a
p forall a. Num a => a -> a -> a
- a
3forall a. Num a => a -> a -> a
*(a
1forall a. Num a => a -> a -> a
-a
t)forall a. Num a => a -> a -> a
*a
tforall a. Num a => a -> a -> a
*(a
aforall a. Num a => a -> a -> a
*forall a. Point a -> a
pointX Point a
d1forall a. Num a => a -> a -> a
*(a
1forall a. Num a => a -> a -> a
-a
t) forall a. Num a => a -> a -> a
+
                                                      a
bforall a. Num a => a -> a -> a
*forall a. Point a -> a
pointX Point a
d2forall a. Num a => a -> a -> a
*a
t)) forall a. Fractional a => a -> a -> a
/
                           forall a. Floating a => Point a -> a
vectorMag Point a
b')
              Vector a
ts Vector (Point a)
ps Vector (Point a)
dir1 Vector (Point a)
dir2 Vector (Point a)
bt'
      newCurve :: CubicBezier a
newCurve = forall a. Point a -> Point a -> Point a -> Point a -> CubicBezier a
CubicBezier Point a
p1 (Point a
p2 forall v. AdditiveGroup v => v -> v -> v
^+^ a
aforall v. VectorSpace v => Scalar v -> v -> v
*^(Point a
p2forall v. AdditiveGroup v => v -> v -> v
^-^Point a
p1)) (Point a
p3 forall v. AdditiveGroup v => v -> v -> v
^+^ a
bforall v. VectorSpace v => Scalar v -> v -> v
*^(Point a
p3forall v. AdditiveGroup v => v -> v -> v
^-^Point a
p4)) Point a
p4
      (Vector (Point a)
bt2,Vector (Point a)
bt2') = forall a b.
(Unbox a, Unbox b) =>
Vector (a, b) -> (Vector a, Vector b)
V.unzip forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a (b :: * -> *).
(Unbox a, Fractional a, GenericBezier b) =>
b a -> a -> (Point a, Point a)
evalBezierDeriv CubicBezier a
newCurve) Vector a
newTs
      err2 :: Vector a
err2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Floating a => Point a -> Point a -> a
vectorDistance Vector (Point a)
pts Vector (Point a)
bt2
      maxErr :: a
maxErr = forall a. (Unbox a, Ord a) => Vector a -> a
V.maximum Vector a
err2
      -- alternative method for finding the t values:
      -- newTs = V.zipWith (-) ts (V.zipWith (calcDeltaT newCurve) pts ts)
  if Int
maxiter forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| forall a. Num a => a -> a
abs(a
err forall a. Num a => a -> a -> a
- a
maxErr) forall a. Ord a => a -> a -> Bool
<= a
errforall a. Fractional a => a -> a -> a
/a
8
    then forall (m :: * -> *) a. Monad m => a -> m a
return (CubicBezier a
newCurve, Vector a
newTs, Vector a
err2, a
maxErr, Vector (Point a)
bt2)
    else forall a.
(Unbox a, Ord a, Floating a) =>
CubicBezier a
-> Vector (Point a)
-> Vector a
-> Int
-> a
-> Vector (Point a)
-> Vector (Point a)
-> Maybe (CubicBezier a, Vector a, Vector a, a, Vector (Point a))
approximateCubic' CubicBezier a
newCurve Vector (Point a)
pts Vector a
newTs (Int
maxiterforall a. Num a => a -> a -> a
-Int
1) a
maxErr Vector (Point a)
bt2 Vector (Point a)
bt2'
{-# SPECIALIZE approximateCubic' ::
  CubicBezier Double -> V.Vector (Point Double) -> V.Vector Double
  -> Int -> Double -> V.Vector (Point Double)
  -> V.Vector (Point Double)
  -> Maybe (CubicBezier Double, V.Vector Double, V.Vector Double, Double, V.Vector (Point Double)) #-}


-- approximate t by calculating the distances between all points
-- and dividing by the total sum
approximateParams :: (V.Unbox a, Floating a) =>
                     Point a -> Point a -> V.Vector (Point a) -> V.Vector a
approximateParams :: forall a.
(Unbox a, Floating a) =>
Point a -> Point a -> Vector (Point a) -> Vector a
approximateParams Point a
start Point a
end Vector (Point a)
pts 
  | forall a. Unbox a => Vector a -> Bool
V.null Vector (Point a)
pts = forall a. Unbox a => Vector a
V.empty
  | Bool
otherwise =
      let dists :: Vector a
dists = forall a. Unbox a => Int -> (Int -> a) -> Vector a
V.generate (forall a. Unbox a => Vector a -> Int
V.length Vector (Point a)
pts)
                  (\Int
i -> if Int
i forall a. Eq a => a -> a -> Bool
== Int
0
                         then forall a. Floating a => Point a -> Point a -> a
vectorDistance Point a
start (forall a. Unbox a => Vector a -> Int -> a
V.unsafeIndex Vector (Point a)
pts Int
0)
                         else forall a. Floating a => Point a -> Point a -> a
vectorDistance (forall a. Unbox a => Vector a -> Int -> a
V.unsafeIndex Vector (Point a)
pts (Int
iforall a. Num a => a -> a -> a
-Int
1)) (forall a. Unbox a => Vector a -> Int -> a
V.unsafeIndex Vector (Point a)
pts Int
i))
          total :: a
total = forall a. (Unbox a, Num a) => Vector a -> a
V.sum Vector a
dists forall a. Num a => a -> a -> a
+ forall a. Floating a => Point a -> Point a -> a
vectorDistance (forall a. Unbox a => Vector a -> a
V.last Vector (Point a)
pts) Point a
end
      in forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a. Fractional a => a -> a -> a
/ a
total) forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => (a -> a -> a) -> Vector a -> Vector a
V.scanl1 forall a. Num a => a -> a -> a
(+) Vector a
dists
{-# SPECIALIZE approximateParams ::
   Point Double -> Point Double -> V.Vector (Point Double) -> V.Vector Double #-}

-- Alternative method for finding the next t values, using
-- Newton-Rafphson.  There is no noticable difference in speed or
-- efficiency.

-- calcDeltaT :: CubicBezier -> Point -> Double -> Double
-- calcDeltaT curve (Point !ptx !pty) t = let
--   (Point bezx bezy, Point dbezx dbezy, Point ddbezx ddbezy, _) = evalBezierDerivs curve t
--   in ((bezx - ptx) * dbezx + (bezy - pty) * dbezy) /
--      (dbezx * dbezx + dbezy * dbezy + (bezx - ptx) * ddbezx + (bezy - pty) * ddbezy)