-- | Bezier curves.
--
module Geometry.Bezier
(
    Bezier(..)

    -- * General properties
    , bezierLength

    -- * Indexing
    , bezierParametric
    , bezierArcParametric

    -- * Subdividing
    , bezierSubdivideEquiparametric
    , bezierSubdivideEquidistant
    , bezierSubdivideCasteljau

    -- * Interpolation
    , bezierSmoothen

    -- * References
    -- $references
)
where



import           Control.DeepSeq
import           Data.Vector     (Vector, (!))
import qualified Data.Vector     as V

import Geometry.Core
import Geometry.LookupTable.Lookup1
import Numerics.ConvergentRecursion
import Numerics.DifferentialEquation
import Numerics.Integrate
import Numerics.LinearEquationSystem



-- $setup
-- >>> import Draw
-- >>> import qualified Graphics.Rendering.Cairo as C



-- $references
--
-- == Arc length parameterization
--
-- * Moving Along a Curve with Specified Speed (2019)
--   by David Eberly
--   https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
--
-- == Smoothening
--
-- * Cubic Bézier Splines
--   by Michael Joost
--   https://www.michael-joost.de/bezierfit.pdf
--
-- * Building Smooth Paths Using Bézier Curves
--   by Stuart Kent
--   https://www.stkent.com/2015/07/03/building-smooth-paths-using-bezier-curves.html



-- | Cubic Bezier curve, defined by start, first/second control points, and end.
data Bezier = Bezier !Vec2 !Vec2 !Vec2 !Vec2 deriving (Bezier -> Bezier -> Bool
(Bezier -> Bezier -> Bool)
-> (Bezier -> Bezier -> Bool) -> Eq Bezier
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Bezier -> Bezier -> Bool
== :: Bezier -> Bezier -> Bool
$c/= :: Bezier -> Bezier -> Bool
/= :: Bezier -> Bezier -> Bool
Eq, Eq Bezier
Eq Bezier
-> (Bezier -> Bezier -> Ordering)
-> (Bezier -> Bezier -> Bool)
-> (Bezier -> Bezier -> Bool)
-> (Bezier -> Bezier -> Bool)
-> (Bezier -> Bezier -> Bool)
-> (Bezier -> Bezier -> Bezier)
-> (Bezier -> Bezier -> Bezier)
-> Ord Bezier
Bezier -> Bezier -> Bool
Bezier -> Bezier -> Ordering
Bezier -> Bezier -> Bezier
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Bezier -> Bezier -> Ordering
compare :: Bezier -> Bezier -> Ordering
$c< :: Bezier -> Bezier -> Bool
< :: Bezier -> Bezier -> Bool
$c<= :: Bezier -> Bezier -> Bool
<= :: Bezier -> Bezier -> Bool
$c> :: Bezier -> Bezier -> Bool
> :: Bezier -> Bezier -> Bool
$c>= :: Bezier -> Bezier -> Bool
>= :: Bezier -> Bezier -> Bool
$cmax :: Bezier -> Bezier -> Bezier
max :: Bezier -> Bezier -> Bezier
$cmin :: Bezier -> Bezier -> Bezier
min :: Bezier -> Bezier -> Bezier
Ord, Int -> Bezier -> ShowS
[Bezier] -> ShowS
Bezier -> String
(Int -> Bezier -> ShowS)
-> (Bezier -> String) -> ([Bezier] -> ShowS) -> Show Bezier
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Bezier -> ShowS
showsPrec :: Int -> Bezier -> ShowS
$cshow :: Bezier -> String
show :: Bezier -> String
$cshowList :: [Bezier] -> ShowS
showList :: [Bezier] -> ShowS
Show)

instance NFData Bezier where rnf :: Bezier -> ()
rnf Bezier
_ = ()

instance Transform Bezier where
    transform :: Transformation -> Bezier -> Bezier
transform Transformation
t (Bezier Vec2
a Vec2
b Vec2
c Vec2
d) = Vec2 -> Vec2 -> Vec2 -> Vec2 -> Bezier
Bezier
        (Transformation -> Vec2 -> Vec2
forall geo. Transform geo => Transformation -> geo -> geo
transform Transformation
t Vec2
a)
        (Transformation -> Vec2 -> Vec2
forall geo. Transform geo => Transformation -> geo -> geo
transform Transformation
t Vec2
b)
        (Transformation -> Vec2 -> Vec2
forall geo. Transform geo => Transformation -> geo -> geo
transform Transformation
t Vec2
c)
        (Transformation -> Vec2 -> Vec2
forall geo. Transform geo => Transformation -> geo -> geo
transform Transformation
t Vec2
d)

-- | <<docs/haddock/Geometry/Bezier/bounding_box_bezier.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Geometry/Bezier/bounding_box_bezier.svg" 150 100 $ \_ -> do
--     let curve = let curveRaw = transform (rotate (deg (-30))) (Bezier (Vec2 0 0) (Vec2 1 5) (Vec2 2.5 (-1)) (Vec2 3 3))
--                     fitToBox = transform (transformBoundingBox curveRaw (shrinkBoundingBox 10 [zero, Vec2 150 100]) (TransformBBSettings FitWidthHeight IgnoreAspect FitAlignCenter))
--                 in fitToBox curveRaw
--     cairoScope $ do
--         setColor (mma 1)
--         C.setDash [1.5,3] 0
--         sketch (boundingBox curve)
--         C.stroke
--     cairoScope $ do
--         C.setLineWidth 2
--         sketch curve
--         C.stroke
-- :}
-- Generated file: size 2KB, crc32: 0x73983106
instance HasBoundingBox Bezier where
    boundingBox :: Bezier -> BoundingBox
boundingBox bezier :: Bezier
bezier@(Bezier Vec2
start Vec2
_ Vec2
_ Vec2
end) = [Vec2] -> BoundingBox
forall a. HasBoundingBox a => a -> BoundingBox
boundingBox ([Vec2
start, Vec2
end] [Vec2] -> [Vec2] -> [Vec2]
forall a. [a] -> [a] -> [a]
++ [Bezier -> Double -> Vec2
bezierParametric Bezier
bezier Double
t | Double
t <- [Double]
extremalTs])
      where

        -- Alg idea: find the roots of the Bezier’s first derivative, collect
        -- points where x/y components are extremal. Hastily written ugly code

        extremalTs :: [Double]
extremalTs = (Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Double
t -> Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>=Double
0 Bool -> Bool -> Bool
&& Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) (Bezier -> [Double]
extremalT Bezier
bezier)

        extremalT :: Bezier -> [Double]
extremalT (Bezier Vec2
a Vec2
b Vec2
c Vec2
d) = Double -> Double -> Double -> [Double]
forall {a}. (Floating a, Ord a) => a -> a -> a -> [a]
solveQuadratic (Vec2 -> Double
x Vec2
a') (Vec2 -> Double
x Vec2
b') (Vec2 -> Double
x Vec2
c') [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Double -> Double -> Double -> [Double]
forall {a}. (Floating a, Ord a) => a -> a -> a -> [a]
solveQuadratic (Vec2 -> Double
y Vec2
a') (Vec2 -> Double
y Vec2
b') (Vec2 -> Double
y Vec2
c')
          where
            -- Coefficients of the first derivative polynomial of the Bezier curve
            a' :: Vec2
a' = ((-Double
3) Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
a) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (Double
9Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
b) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Double
9Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
c Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double
3Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
d
            b' :: Vec2
b' = Double
6Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
a Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Double
12Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
b Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double
6Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
c
            c' :: Vec2
c' = (-Double
3)Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
a Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double
3Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*.Vec2
b

            x :: Vec2 -> Double
x (Vec2 Double
x' Double
_) = Double
x'
            y :: Vec2 -> Double
y (Vec2 Double
_ Double
y') = Double
y'

        solveQuadratic :: a -> a -> a -> [a]
solveQuadratic a
a a
b a
c = case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
discriminant a
0 of
            Ordering
LT -> []
            Ordering
EQ -> [-a
ba -> a -> a
forall a. Fractional a => a -> a -> a
/(a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
a)]
            Ordering
GT -> [ (-a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
sqrt a
discriminant) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
a), (-a
b a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a. Floating a => a -> a
sqrt a
discriminant) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
a)]
          where
            discriminant :: a
discriminant = a
ba -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 a -> a -> a
forall a. Num a => a -> a -> a
- a
4 a -> a -> a
forall a. Num a => a -> a -> a
* a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
c

-- | Point on a 'Bezier' curve.
--
-- \[
-- \text{bezierParametric}_{\mathbf a,\mathbf b,\mathbf c,\mathbf d}(t) = (1-t)^3\,\mathbf a + 3 (1-t)^2 t\,\mathbf b + 3 (1-t) t^2\,\mathbf c + t^3\,\mathbf d
-- \]
bezierParametric
  :: Bezier
  -> Double -- ^ \([0\ldots 1] = [\mathrm{start}\ldots\mathrm{end}]\)
  -> Vec2
bezierParametric :: Bezier -> Double -> Vec2
bezierParametric (Bezier Vec2
a Vec2
b Vec2
c Vec2
d) Double
t
  =      (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t)Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
3     Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
a
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t)Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t   Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
b
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t)  Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
c
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+.           Double
tDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
3 Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
d

-- | First derivative of a 'Bezier' curve, i.e. its velocity vector.
--
-- \[
-- \text{bezierParametric}'_{\mathbf a,\mathbf b,\mathbf c,\mathbf d}(t) = -3(1-t)^2\,\mathbf a + (3+t(-12+9t))\,\mathbf b + (6-9t)t\,\mathbf c + 3t^2\,\mathbf d
-- \]
bezierParametric'
  :: Bezier
  -> Double -- ^ \([0\ldots 1] = [\mathrm{start}\ldots\mathrm{end}]\)
  -> Vec2
bezierParametric' :: Bezier -> Double -> Vec2
bezierParametric' (Bezier Vec2
a Vec2
b Vec2
c Vec2
d) Double
t
  =    (-Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t)Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)    Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
a
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(-Double
12Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
9Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)) Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
b
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. ((Double
6Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
9Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)     Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
c
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)         Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
d

-- | Second derivative of a 'Bezier' curve, i.e. its acceleration vector.
_bezierParametric''
  :: Bezier
  -> Double -- ^ \([0\ldots 1] = [\mathrm{start}\ldots\mathrm{end}]\)
  -> Vec2
_bezierParametric'' :: Bezier -> Double -> Vec2
_bezierParametric'' (Bezier Vec2
a Vec2
b Vec2
c Vec2
d) Double
t
  =    (Double
6Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
6Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)         Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
a
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (-Double
12Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
18Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)      Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
b
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (Double
6Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
18Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)        Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
c
    Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (Double
6Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)           Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
d

-- | Length of a Bézier curve by integration using Simpson’s Rule. It’s comparably
-- slow, but very precise.
--
-- If speed is more important than quality, simply sum up a segmented subdivision.
-- 16 subdivisions is often enough for a reasonably good result:
--
-- >>> curve = Bezier (Vec2 0 0) (Vec2 1 5) (Vec2 2.5 (-1)) (Vec2 3 3)
-- >>> bezierLength curve
-- 5.645...
-- >>> polylineLength . Polyline . bezierSubdivideEquiparametric 16 $ curve
-- 5.612...
bezierLength
    :: Bezier
    -> Double
bezierLength :: Bezier -> Double
bezierLength Bezier
bezier = (Int -> Double) -> Double -> Double
retryExponentiallyUntilPrecision ((Double -> Double) -> Double -> Double -> Int -> Double
forall a. Fractional a => (a -> a) -> a -> a -> Int -> a
integrateSimpson13 Double -> Double
f Double
0 Double
1) Double
1e-6
  where
    f :: Double -> Double
f Double
t = Vec2 -> Double
norm (Bezier -> Double -> Vec2
bezierParametric' Bezier
bezier Double
t)

-- | Trace a 'Bezier' curve with a number of points, using the polynomial curve
-- parameterization. This is very fast, but leads to unevenly spaced points.
--
-- For subdivision by arc length, use 'bezierSubdivideEquidistant'.
bezierSubdivideEquiparametric
    :: Int -- ^ Number of segments
    -> Bezier
    -> [Vec2]
bezierSubdivideEquiparametric :: Int -> Bezier -> [Vec2]
bezierSubdivideEquiparametric Int
n Bezier
bz = (Double -> Vec2) -> [Double] -> [Vec2]
forall a b. (a -> b) -> [a] -> [b]
map (Bezier -> Double -> Vec2
bezierParametric Bezier
bz) [Double]
points
  where
    points :: [Double]
points = (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
x -> Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) [Int
0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

-- | Trace a 'Bezier' curve with a number of evenly spaced points by arc length.
-- This is much more expensive than 'bezierSubdivideEquiparametric', but may be desirable for
-- aesthetic purposes.
--
-- Here it is alongside 'bezierSubdivideEquiparametric':
--
-- <<docs/haddock/Geometry/Bezier/subdivide_s_t_comparison.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Geometry/Bezier/subdivide_s_t_comparison.svg" 300 150 $ \_ -> do
--     let curve = let curveRaw = transform (rotate (deg (-30))) (Bezier (Vec2 0 0) (Vec2 1 5) (Vec2 2.5 (-1)) (Vec2 3 3))
--                     fitToBox = transform (transformBoundingBox curveRaw (Vec2 10 10, Vec2 290 90) (TransformBBSettings FitWidthHeight IgnoreAspect FitAlignCenter))
--                 in fitToBox curveRaw
--         evenlySpaced = bezierSubdivideEquidistant 16 curve
--         unevenlySpaced = bezierSubdivideEquiparametric 16 curve
--         offsetBelow :: Transform geo => geo -> geo
--         offsetBelow = transform (translate (Vec2 0 50))
--     cairoScope $ do
--         setColor $ mma 1
--         sketch curve
--         C.stroke
--         sketch (offsetBelow curve)
--         C.stroke
--     for_ (zip evenlySpaced unevenlySpaced) $ \(e, u') -> do
--         let u = offsetBelow u'
--         let circle p = sketch (Circle p 3) >> C.stroke
--             connect p q = do
--                 let line = resizeLineSymmetric (*0.8) (Line p q)
--                 sketch line
--                 C.stroke
--         cairoScope (setColor (mma 0) >> circle e)
--         cairoScope (setColor (mma 3) >> circle u)
--         cairoScope (setColor (black `withOpacity` 0.2) >> connect e u)
-- :}
-- Generated file: size 17KB, crc32: 0x7c147951
bezierSubdivideEquidistant
    :: Int -- ^ Number of segments
    -> Bezier
    -> [Vec2]
bezierSubdivideEquidistant :: Int -> Bezier -> [Vec2]
bezierSubdivideEquidistant Int
n Bezier
bz = (Double -> Vec2) -> [Double] -> [Vec2]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Vec2
bezier [Double]
distances
  where

    -- The step width should correlate with the length of the curve to get a decent
    -- RK estimator. This allows both large and tiny curves to be subdivided well.
    -- Increasing this to beyond 2^10 shows only pixel-wide changes, if even.
    bezier :: Double -> Vec2
bezier = Bezier -> Double -> Double -> Vec2
bezierArcParametric_ode Bezier
bz (Double
len Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
10)
    len :: Double
len = Bezier -> Double
bezierLength Bezier
bz

    distances :: [Double]
    distances :: [Double]
distances = [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lenDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) | Int
i <- [Int
0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]]

-- | Get the position on a Bezier curve as a fraction of its length, via solving a
-- differential equation. This is /much/ more expensive to compute than 'bezierParametric'.
--
-- This caches the internal LUT when partially applied, so that the following will
-- only compute it once for repeated lookups:
--
-- @
-- let walkOnBezier = 'bezierArcParametric' bezier 0.01
-- 'print' [walkOnBezier d | d <- [0, 0.1 .. 5]]
-- @
bezierArcParametric
    :: Bezier
    -> Double -- ^ Precision parameter (smaller is more precise but slower).
    -> Double -- ^ Distance to walk on the curve. Clips (stops at the end) when asked to »walk too far«.
    -> Vec2   -- ^ Point at that distance
bezierArcParametric :: Bezier -> Double -> Double -> Vec2
bezierArcParametric = Bezier -> Double -> Double -> Vec2
bezierArcParametric_ode

-- There’s also another method to do this using Newton’s method, detialed in the
-- paper (see references on bezier subdivision).
bezierArcParametric_ode
    :: Bezier
    -> Double -- ^ Precision parameter (smaller is more precise but slower).
    -> Double -- ^ Distance to walk on the curve. Clips (stops at the end) when asked to »walk too far«.
    -> Vec2   -- ^ Point at that distance
bezierArcParametric_ode :: Bezier -> Double -> Double -> Vec2
bezierArcParametric_ode Bezier
bz Double
ds
  = let lut :: LookupTable1 Double Double
lut = Bezier -> Double -> LookupTable1 Double Double
s_to_t_lut_ode Bezier
bz Double
ds
    in \Double
s -> let t :: Double
t = LookupTable1 Double Double -> Double -> Double
lookupInterpolated LookupTable1 Double Double
lut Double
s
             in Bezier -> Double -> Vec2
bezierParametric Bezier
bz Double
t

-- | S⇆T lookup table for a Bezier curve
--
-- We do not explicitly add the end point, so too big of a step width will
-- overshoot the end and be cut off. This can be remedied by a small enough step size ;-)
s_to_t_lut_ode
    :: Bezier
    -> Double -- ^ ODE solver step width. Correlates with result precision/length.
    -> LookupTable1 Double Double -- ^ Lookup table
s_to_t_lut_ode :: Bezier -> Double -> LookupTable1 Double Double
s_to_t_lut_ode Bezier
bz Double
ds = Vector (Double, Double) -> LookupTable1 Double Double
forall a b. Vector (a, b) -> LookupTable1 a b
LookupTable1 ([(Double, Double)] -> Vector (Double, Double)
forall {a}. [(a, Double)] -> Vector (a, Double)
sol_to_vec [(Double, Double)]
sol)
  where
    sol_to_vec :: [(a, Double)] -> Vector (a, Double)
sol_to_vec = [(a, Double)] -> Vector (a, Double)
forall a. [a] -> Vector a
V.fromList ([(a, Double)] -> Vector (a, Double))
-> ([(a, Double)] -> [(a, Double)])
-> [(a, Double)]
-> Vector (a, Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((a, Double) -> Bool) -> [(a, Double)] -> [(a, Double)]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\(a
_s, Double
t) -> Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1)

    sol :: [(Double, Double)]
sol = (Double -> Double -> Double)
-> Double -> Double -> Double -> [(Double, Double)]
forall vec.
VectorSpace vec =>
(Double -> vec -> vec)
-> vec -> Double -> Double -> [(Double, vec)]
rungeKuttaConstantStep Double -> Double -> Double
forall {p}. p -> Double -> Double
dt_ds Double
t0 Double
s0 Double
ds

    dt_ds :: p -> Double -> Double
dt_ds p
_s Double
t = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Vec2 -> Double
norm (Bezier -> Double -> Vec2
bezierParametric' Bezier
bz Double
t)

    t0 :: Double
t0 = Double
0
    s0 :: Double
s0 = Double
0

-- | Approximage a Bezier curve with line segments up to a certain precision, using
-- relatively few points.
--
-- The idea behind Casteljau subdivision is that each Bézier curve can be exactly
-- subdivided into two Bézier curves (of same degree). This is done recursively, in
-- this implementation (and commonly) in the middle of the curve. Once a curve
-- segment is flat enough (given by the tolerance parameter), it is simply rendered
-- as a line.
--
-- The algorithm is based on
-- [an excellent blogpost](https://web.archive.org/web/20180307160123/http://antigrain.com/research/adaptive_bezier/index.html)
-- that is sadly only available on the internet archive these days.
--
-- <<docs/haddock/Geometry/Bezier/bezierSubdivideCasteljau.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Geometry/Bezier/bezierSubdivideCasteljau.svg" 500 330 $ \_ -> do
--     let curve = let curveRaw = transform (rotate (deg (-30))) (Bezier (Vec2 0 0) (Vec2 1 5) (Vec2 2.5 (-1)) (Vec2 3 3))
--                     fitToBox = transform (transformBoundingBox curveRaw (shrinkBoundingBox 10 [zero, Vec2 500 200]) (TransformBBSettings FitWidthHeight IgnoreAspect FitAlignCenter))
--                 in fitToBox curveRaw
--         paintOffset = Vec2 0 30
--     for_ (zip [0..] [50,25,10,2]) $ \(i, tolerance) -> cairoScope $ do
--         let points = bezierSubdivideCasteljau tolerance (transform (translate (fromIntegral i *. paintOffset)) curve)
--         setColor (mma i)
--         C.setLineWidth 2
--         sketch (Polyline points) >> C.stroke
--         for_ points $ \p -> sketch (Circle p 3) >> C.fill
--     cairoScope $ do
--         C.setLineWidth 3
--         setColor black
--         sketch (transform (translate (4*.paintOffset)) curve)
--         C.stroke
-- :}
-- Generated file: size 20KB, crc32: 0x679b311c
bezierSubdivideCasteljau :: Double -> Bezier -> [Vec2]
bezierSubdivideCasteljau :: Double -> Bezier -> [Vec2]
bezierSubdivideCasteljau Double
tolerance curve :: Bezier
curve@(Bezier Vec2
pFirst Vec2
_ Vec2
_ Vec2
_) = Vec2
pFirst Vec2 -> [Vec2] -> [Vec2]
forall a. a -> [a] -> [a]
: Bezier -> [Vec2]
go Bezier
curve
  where
    go :: Bezier -> [Vec2]
go (Bezier Vec2
p1 p2 :: Vec2
p2@(Vec2 Double
x2 Double
y2) p3 :: Vec2
p3@(Vec2 Double
x3 Double
y3) p4 :: Vec2
p4@(Vec2 Double
x4 Double
y4)) =
        let
            p12 :: Vec2
p12   = (Vec2
p1   Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p2  ) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2
            p23 :: Vec2
p23   = (Vec2
p2   Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p3  ) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2
            p34 :: Vec2
p34   = (Vec2
p3   Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p4  ) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2
            p123 :: Vec2
p123  = (Vec2
p12  Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p23 ) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2
            p234 :: Vec2
p234  = (Vec2
p23  Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p34 ) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2
            p1234 :: Vec2
p1234 = (Vec2
p123 Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
p234) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2

            dp :: Vec2
dp@(Vec2 Double
dx Double
dy) = Vec2
p4 Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
p1

            -- d2, d3 are the distance from p2, p3 from the line
            -- connecting p1 and p4. A curve is flat when those
            -- two are short together.
            d2 :: Double
d2 = Double -> Double
forall a. Num a => a -> a
abs ((Double
x2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
dy Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
y2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
dx)
            d3 :: Double
d3 = Double -> Double
forall a. Num a => a -> a
abs ((Double
x3Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
dy Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
y3Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y4)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
dx)
            curveIsFlat :: Bool
curveIsFlat = (Double
d2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d3)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
d2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d3) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
toleranceDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Vec2 -> Double
normSquare Vec2
dp
        in if Bool
curveIsFlat
            then -- We return only the last point so we don’t get duplicate
                 -- points for each start+end of adjacent curves.
                 -- The very first point is forgotten by the
                [Vec2
p4]
            else Bezier -> [Vec2]
go (Vec2 -> Vec2 -> Vec2 -> Vec2 -> Bezier
Bezier Vec2
p1 Vec2
p12 Vec2
p123 Vec2
p1234)
                 [Vec2] -> [Vec2] -> [Vec2]
forall a. [a] -> [a] -> [a]
++
                 Bezier -> [Vec2]
go (Vec2 -> Vec2 -> Vec2 -> Vec2 -> Bezier
Bezier Vec2
p1234 Vec2
p234 Vec2
p34 Vec2
p4)

-- | Smoothen a number of points by putting a Bezier curve between each pair.
-- Useful to e.g. make a sketch nicer, or interpolate between points of a crude
-- solution to a differential equation.
--
-- If the first and last point are identical, assume the trajectory is closed, and
-- smoothly interpolate between beginning and end as well, yielding a nice loop.
--
-- <<docs/haddock/Geometry/Bezier/bezierSmoothen.svg>>
--
-- This works even for input loops with at least three points (plus one duplicate for the looping condition start=end). (For smaller inputs the
-- scale is arbitrary so these solutions are excluded here.)
--
-- <<docs/haddock/Geometry/Bezier/bezierSmoothen_tiny_loop.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Geometry/Bezier/bezierSmoothen.svg" 400 300 $ \_ -> do
--     let points = [ Vec2 100 275, Vec2 150 125, Vec2 300 275, Vec2 350 75
--                  , Vec2 250 50, Vec2 75 75, Vec2 75 50, Vec2 225 100
--                  , Vec2 100 275 ]
--         prettyBezier bezier@(Bezier p0 p1 p2 p3) = do
--             cairoScope $ do
--                 setColor black
--                 sketch bezier
--                 C.stroke
--             cairoScope $ do
--                 setColor (mma 0)
--                 sketch (Circle p1 4) >> C.fill
--                 sketch (Line p0 p1) >> C.stroke
--             cairoScope $ do
--                 setColor (mma 1)
--                 sketch (Circle p2 4) >> C.fill
--                 sketch (Line p3 p2) >> C.stroke
--             cairoScope $ do
--                 sketch (Circle p0 5)
--                 setColor (mma 3)
--                 C.fillPreserve
--                 setColor black
--                 C.stroke
--     C.setLineWidth 2
--     for_ (bezierSmoothen points) prettyBezier
-- :}
-- Generated file: size 15KB, crc32: 0xb9f3566d
--
-- >>> :{
-- haddockRender "Geometry/Bezier/bezierSmoothen_tiny_loop.svg" 200 200 $ \_ -> do
--     let points = [ Vec2 20 60
--                  , Vec2 150 50
--                  , Vec2 140 180
--                  , Vec2 20 60 ]
--     let prettyBezier bezier@(Bezier p0 p1 p2 p3) = do
--             cairoScope $ do
--                 setColor black
--                 sketch bezier
--                 C.stroke
--             cairoScope $ do
--                 setColor (mma 0)
--                 sketch (Circle p1 4) >> C.fill
--                 sketch (Line p0 p1) >> C.stroke
--             cairoScope $ do
--                 setColor (mma 1)
--                 sketch (Circle p2 4) >> C.fill
--                 sketch (Line p3 p2) >> C.stroke
--             cairoScope $ do
--                 sketch (Circle p0 5)
--                 setColor (mma 3)
--                 C.fillPreserve
--                 setColor black
--                 C.stroke
--     C.setLineWidth 2
--     for_ (bezierSmoothen points) prettyBezier
-- :}
-- Generated file: size 7KB, crc32: 0xca727498
bezierSmoothen :: Sequential vector => vector Vec2 -> Vector Bezier
bezierSmoothen :: forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> Vector Bezier
bezierSmoothen vector Vec2
vecSequence
    | Vector Vec2 -> Vec2
forall a. Vector a -> a
V.head Vector Vec2
vec Vec2 -> Vec2 -> Bool
forall a. Eq a => a -> a -> Bool
== Vector Vec2 -> Vec2
forall a. Vector a -> a
V.last Vector Vec2
vec = Vector Vec2 -> Vector Bezier
forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> Vector Bezier
bezierSmoothenLoop (Vector Vec2 -> Vector Vec2
forall a. Vector a -> Vector a
V.tail Vector Vec2
vec)
    | Bool
otherwise = Vector Vec2 -> Vector Bezier
forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> Vector Bezier
bezierSmoothenOpen Vector Vec2
vec
    where vec :: Vector Vec2
vec = vector Vec2 -> Vector Vec2
forall a. vector a -> Vector a
forall (f :: * -> *) a. Sequential f => f a -> Vector a
toVector vector Vec2
vecSequence

-- | Smoothen a number of points by putting a Bezier curve between each pair,
-- assuming the curve is open (i.e. the last and first point have nothing to do
-- with each other).
bezierSmoothenOpen :: Sequential vector => vector Vec2 -> Vector Bezier
bezierSmoothenOpen :: forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> Vector Bezier
bezierSmoothenOpen vector Vec2
pointsSequence = (Vec2 -> Vec2 -> Vec2 -> Vec2 -> Bezier)
-> Vector Vec2
-> Vector Vec2
-> Vector Vec2
-> Vector Vec2
-> Vector Bezier
forall a b c d e.
(a -> b -> c -> d -> e)
-> Vector a -> Vector b -> Vector c -> Vector d -> Vector e
V.zipWith4 Vec2 -> Vec2 -> Vec2 -> Vec2 -> Bezier
Bezier Vector Vec2
points Vector Vec2
controlPointsStart Vector Vec2
controlPointsEnd (Vector Vec2 -> Vector Vec2
forall a. Vector a -> Vector a
V.tail Vector Vec2
points)
  where
    points :: Vector Vec2
points = vector Vec2 -> Vector Vec2
forall a. vector a -> Vector a
forall (f :: * -> *) a. Sequential f => f a -> Vector a
toVector vector Vec2
pointsSequence
    n :: Int
n = Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1

    controlPointsStart :: Vector Vec2
controlPointsStart =
        let low :: Vector Double
low   = Int -> Vector Double
lowerDiagonal (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
            diag :: Vector Double
diag  = Int -> Vector Double
diagonal      Int
n
            upper :: Vector Double
upper = Int -> Vector Double
upperDiagonal (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
            rhs :: Vector Vec2
rhs   = Int -> Vector Vec2 -> Vector Vec2
forall vec. VectorSpace vec => Int -> Vector vec -> Vector vec
target        Int
n Vector Vec2
points
        in Vector Double
-> Vector Double -> Vector Double -> Vector Vec2 -> Vector Vec2
forall vec.
VectorSpace vec =>
Vector Double
-> Vector Double -> Vector Double -> Vector vec -> Vector vec
solveTridiagonal Vector Double
low Vector Double
diag Vector Double
upper Vector Vec2
rhs
    controlPointsEnd :: Vector Vec2
controlPointsEnd = Int -> (Int -> Vec2) -> Vector Vec2
forall a. Int -> (Int -> a) -> Vector a
V.generate (Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
controlPointsStart) ((Int -> Vec2) -> Vector Vec2) -> (Int -> Vec2) -> Vector Vec2
forall a b. (a -> b) -> a -> b
$ \Int
i -> case () of
        ()
_ | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 -> (Vector Vec2
points Vector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
! Int
n Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vector Vec2
controlPointsStart Vector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
! (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
2
          | Bool
otherwise -> Double
2 Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. (Vector Vec2
points Vector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
! (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vector Vec2
controlPointsStart Vector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
! (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)

upperDiagonal :: Int -> Vector Double
upperDiagonal :: Int -> Vector Double
upperDiagonal Int
len = Int -> Double -> Vector Double
forall a. Int -> a -> Vector a
V.replicate Int
len Double
1

lowerDiagonal :: Int -> Vector Double
lowerDiagonal :: Int -> Vector Double
lowerDiagonal Int
len = Int -> (Int -> Double) -> Vector Double
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
len ((Int -> Double) -> Vector Double)
-> (Int -> Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ \Int
i -> case () of
    ()
_ | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
lenInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 -> Double
2
      | Bool
otherwise -> Double
1

diagonal :: Int -> Vector Double
diagonal :: Int -> Vector Double
diagonal Int
len = Int -> (Int -> Double) -> Vector Double
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
len ((Int -> Double) -> Vector Double)
-> (Int -> Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ \Int
i -> case () of
    ()
_ | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0     -> Double
2
      | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
lenInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 -> Double
7
      | Bool
otherwise  -> Double
4

target :: VectorSpace vec => Int -> Vector vec -> Vector vec
target :: forall vec. VectorSpace vec => Int -> Vector vec -> Vector vec
target Int
n Vector vec
vertices = Int -> (Int -> vec) -> Vector vec
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
n ((Int -> vec) -> Vector vec) -> (Int -> vec) -> Vector vec
forall a b. (a -> b) -> a -> b
$ \Int
i -> case () of
    ()
_ | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0    ->      Vector vec
vertices Vector vec -> Int -> vec
forall a. Vector a -> Int -> a
! Int
0     vec -> vec -> vec
forall v. VectorSpace v => v -> v -> v
+. Double
2 Double -> vec -> vec
forall v. VectorSpace v => Double -> v -> v
*. Vector vec
vertices Vector vec -> Int -> vec
forall a. Vector a -> Int -> a
! Int
1
      | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1  -> Double
8 Double -> vec -> vec
forall v. VectorSpace v => Double -> v -> v
*. Vector vec
vertices Vector vec -> Int -> vec
forall a. Vector a -> Int -> a
! (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) vec -> vec -> vec
forall v. VectorSpace v => v -> v -> v
+.      Vector vec
vertices Vector vec -> Int -> vec
forall a. Vector a -> Int -> a
! Int
n
      | Bool
otherwise -> Double
4 Double -> vec -> vec
forall v. VectorSpace v => Double -> v -> v
*. Vector vec
vertices Vector vec -> Int -> vec
forall a. Vector a -> Int -> a
! Int
i     vec -> vec -> vec
forall v. VectorSpace v => v -> v -> v
+. Double
2 Double -> vec -> vec
forall v. VectorSpace v => Double -> v -> v
*. Vector vec
vertices Vector vec -> Int -> vec
forall a. Vector a -> Int -> a
! (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)

-- | Like 'bezierSmoothen', but will smoothly connect the start and the end of the
-- given trajectory as well.
bezierSmoothenLoop
    :: Sequential vector
    => vector Vec2  -- ^ [0,1,2,3 ... 99]
    -> Vector Bezier
bezierSmoothenLoop :: forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> Vector Bezier
bezierSmoothenLoop vector Vec2
pointsSequence
    | Vector Vec2 -> Int
forall a. Vector a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Vector Vec2
points Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
2 = Vector Bezier
forall a. Vector a
V.empty -- Is possible, but not unique. The »radius«/size is arbitrary.
    | Bool
otherwise =
        -- The idea is this: we can artificially lengthen a closed trajectory by
        -- wrapping it onto itself. We then interpolate it as if it was open, and later
        -- forget the end parts again. In the overlap, we get a smooth transition.
        -- The overlap has to be 3 points:
        --   1. 0-th derivative so we umm continue in the right direction
        --   2. 1st derivative, so the corner is not pointy
        --   3. 2nd derivative, so the corner is smooth.
        let pointsWrapped :: Vector Vec2
pointsWrapped = Int -> Vector Vec2 -> Vector Vec2
forall {a}. Int -> Vector a -> Vector a
takeLast Int
3 Vector Vec2
points Vector Vec2 -> Vector Vec2 -> Vector Vec2
forall a. Semigroup a => a -> a -> a
<> Vector Vec2
points Vector Vec2 -> Vector Vec2 -> Vector Vec2
forall a. Semigroup a => a -> a -> a
<> Int -> Vector Vec2 -> Vector Vec2
forall {a}. Int -> Vector a -> Vector a
V.take Int
3 Vector Vec2
points
                          -- NB <=2 points is handled in the guard at the top

        in Int -> Int -> Vector Bezier -> Vector Bezier
forall a. Int -> Int -> Vector a -> Vector a
V.slice Int
3 (Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points) (Vector Vec2 -> Vector Bezier
forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> Vector Bezier
bezierSmoothenOpen Vector Vec2
pointsWrapped)
            -- Notes:
            --   1. The slice starts at point 2 because points 0+1 are responsible for continuity of the 1st+2nd derivative.
            --   2. (implicit) we also discard the 1st+2nd last elements for the same reason in the slice.
    where
      points :: Vector Vec2
points = vector Vec2 -> Vector Vec2
forall a. vector a -> Vector a
forall (f :: * -> *) a. Sequential f => f a -> Vector a
toVector vector Vec2
pointsSequence
      takeLast :: Int -> Vector a -> Vector a
takeLast Int
n Vector a
vec = Int -> Vector a -> Vector a
forall {a}. Int -> Vector a -> Vector a
V.drop (Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
vec Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n) Vector a
vec