{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

-- |
-- Module      :  Diagrams.CubicSpline.Boehm
-- Copyright   :  (c) 2015 diagrams-lib team (see LICENSE)
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  diagrams-discuss@googlegroups.com
--
-- Boehm's algorithm for converting a cubic B-spline into a sequence
-- of cubic Bezier curves.
--
-- See
--
--   * Thomas W. Sederberg, /An Introduction to B-Spline Curves/,
--     <http://web.archive.org/web/20120227050519/http://tom.cs.byu.edu/~455/bs.pdf>
--
--   * Lyle Ramshaw, /Blossoming: A Connect-the-Dots Approach to Splines/,
--     <http://www.hpl.hp.com/techreports/Compaq-DEC/SRC-RR-19.pdf>
module Diagrams.CubicSpline.Boehm (
  BSpline,
  bsplineToBeziers,
  bspline,
) where

import Data.List.NonEmpty (NonEmpty (..))
import qualified Data.List.NonEmpty as NE
import Diagrams.Core (N, Point, V, origin)
import Diagrams.Located (at, loc, unLoc)
import Diagrams.Segment (FixedSegment (..), fromFixedSeg)
import Diagrams.TrailLike (TrailLike, fromLocSegments)
import Diagrams.Util (iterateN)
import Linear.Vector (Additive)
import Linear.Vector.Compat (lerp)

type BSpline v n = NonEmpty (Point v n)

-- | @affineCombo a b t x y@ computes an affine combination of x and y
--   which lies at parameter t, if x has parameter a and y has parameter b.
--   The usual @lerp@ arises by giving x parameter 0 and y parameter 1.
affineCombo :: (Additive f, Fractional a) => a -> a -> a -> f a -> f a -> f a
affineCombo :: forall (f :: * -> *) a.
(Additive f, Fractional a) =>
a -> a -> a -> f a -> f a -> f a
affineCombo a
a a
b a
t = a -> f a -> f a -> f a
forall (f :: * -> *) a.
(Additive f, Num a) =>
a -> f a -> f a -> f a
lerp ((a
t a -> a -> a
forall a. Num a => a -> a -> a
- a
a) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
b a -> a -> a
forall a. Num a => a -> a -> a
- a
a))

takeMaybe :: Int -> [a] -> Maybe [a]
takeMaybe :: forall a. Int -> [a] -> Maybe [a]
takeMaybe Int
0 [a]
_ = [a] -> Maybe [a]
forall a. a -> Maybe a
Just []
takeMaybe Int
_ [] = Maybe [a]
forall a. Maybe a
Nothing
takeMaybe Int
k (a
a : [a]
as) = (a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [a]) -> Maybe [a] -> Maybe [a]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> [a] -> Maybe [a]
forall a. Int -> [a] -> Maybe [a]
takeMaybe (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [a]
as

-- | @windows k xs@ yields all the length-@k@ windows from @xs@, e.g.
--   @windows 3 [a,b,c,d,e] == [[a,b,c], [b,c,d], [c,d,e]]@.
windows :: Int -> [a] -> [NonEmpty a]
windows :: forall a. Int -> [a] -> [NonEmpty a]
windows Int
0 [a]
_ = []
windows Int
_ [] = []
windows Int
k (a
a : [a]
as) = case Int -> [a] -> Maybe [a]
forall a. Int -> [a] -> Maybe [a]
takeMaybe (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [a]
as of
  Just [a]
as' -> (a
a a -> [a] -> NonEmpty a
forall a. a -> [a] -> NonEmpty a
:| [a]
as') NonEmpty a -> [NonEmpty a] -> [NonEmpty a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [NonEmpty a]
forall a. Int -> [a] -> [NonEmpty a]
windows Int
k [a]
as
  Maybe [a]
Nothing -> []

-- | @extend k xs@ extends @xs@ on both ends by prepending @k@ copies
-- of its head and appending @k@ copies of its last element.  For example,
-- @extend 2 [1..5] == [1,1,1,2,3,4,5,5,5]@.
extend :: Int -> [a] -> [a]
extend :: forall a. Int -> [a] -> [a]
extend Int
_ [] = []
extend Int
k (a
x : [a]
xs) = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
x [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xs) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k (NonEmpty a -> a
forall a. NonEmpty a -> a
NE.last (a
x a -> [a] -> NonEmpty a
forall a. a -> [a] -> NonEmpty a
:| [a]
xs))

-- | A "polar point" is a point along with three knot values.
--   We consider the "blossom" of a cubic spline, a 3-ary symmetric
--   polynomial; a polar point consists of 3 values paired with the
--   output of the blossom at those input values.  Blossoms have nice
--   affine properties so this makes it easy to keep track of how
--   points may be combined to yield other points of interest.
--
--   Invariant: knot values are in nondecreasing order.
data PolarPt v n = PP {forall (v :: * -> *) n. PolarPt v n -> Point v n
unPP :: Point v n, forall (v :: * -> *) n. PolarPt v n -> NonEmpty n
_knots :: NonEmpty n}

mkPolarPt :: Ord n => Point v n -> NonEmpty n -> PolarPt v n
mkPolarPt :: forall n (v :: * -> *).
Ord n =>
Point v n -> NonEmpty n -> PolarPt v n
mkPolarPt Point v n
pt NonEmpty n
kts = Point v n -> NonEmpty n -> PolarPt v n
forall (v :: * -> *) n. Point v n -> NonEmpty n -> PolarPt v n
PP Point v n
pt (NonEmpty n -> NonEmpty n
forall a. Ord a => NonEmpty a -> NonEmpty a
NE.sort NonEmpty n
kts)

-- | Precondition: the knots of the two polar points overlap, like abc
--   and bcd.  The @Int@ should be 0 or 1, indicating which knot to
--   replicate (0 means to replicate b, yielding bbc, 1 means to
--   replicate c, yielding bcc).
combine ::
  (Additive v, Fractional n, Ord n) =>
  Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
combine :: forall (v :: * -> *) n.
(Additive v, Fractional n, Ord n) =>
Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
combine Int
k (PP Point v n
pt1 NonEmpty n
kts1) (PP Point v n
pt2 NonEmpty n
kts2) =
  Point v n -> NonEmpty n -> PolarPt v n
forall n (v :: * -> *).
Ord n =>
Point v n -> NonEmpty n -> PolarPt v n
mkPolarPt
    (n -> n -> n -> Point v n -> Point v n -> Point v n
forall (f :: * -> *) a.
(Additive f, Fractional a) =>
a -> a -> a -> f a -> f a -> f a
affineCombo (NonEmpty n -> n
forall a. NonEmpty a -> a
NE.head NonEmpty n
kts1) (NonEmpty n -> n
forall a. NonEmpty a -> a
NE.last NonEmpty n
kts2) n
newKt Point v n
pt1 Point v n
pt2)
    (n
newKt n -> [n] -> NonEmpty n
forall a. a -> [a] -> NonEmpty a
:| Int -> NonEmpty n -> [n]
forall a. Int -> NonEmpty a -> [a]
NE.drop Int
1 NonEmpty n
kts1)
 where
  newKt :: n
newKt = NonEmpty n -> [n]
forall a. NonEmpty a -> [a]
NE.toList NonEmpty n
kts2 [n] -> Int -> n
forall a. HasCallStack => [a] -> Int -> a
!! Int
k

-- | Convert a uniform cubic B-spline to a sequence of cubic beziers.
--   (/Uniform/ refers to the fact that the knots are assumed to be
--   evenly spaced, with no duplicates.)  The knots at the end are
--   replicated so the cubic spline begins and ends at the first and
--   last control points, tangent to the line from the end control
--   point to the next.
bsplineToBeziers ::
  (Additive v, Fractional n, Num n, Ord n) =>
  BSpline v n ->
  [FixedSegment v n]
bsplineToBeziers :: forall (v :: * -> *) n.
(Additive v, Fractional n, Num n, Ord n) =>
BSpline v n -> [FixedSegment v n]
bsplineToBeziers BSpline v n
controls = [FixedSegment v n]
beziers
 where
  n :: Int
n = BSpline v n -> Int
forall a. NonEmpty a -> Int
NE.length BSpline v n
controls
  numKnots :: Int
numKnots = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2
  knots :: [n]
knots = Int -> (n -> n) -> n -> [n]
forall a. Int -> (a -> a) -> a -> [a]
iterateN Int
numKnots (n -> n -> n
forall a. Num a => a -> a -> a
+ n
1 n -> n -> n
forall a. Fractional a => a -> a -> a
/ (Int -> n
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
numKnots n -> n -> n
forall a. Num a => a -> a -> a
- n
1)) n
0

  -- The control points are P(a,b,c), P(b,c,d), P(c,d,e), and so on.
  controls' :: [PolarPt v n]
controls' = (Point v n -> NonEmpty n -> PolarPt v n)
-> [Point v n] -> [NonEmpty n] -> [PolarPt v n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Point v n -> NonEmpty n -> PolarPt v n
forall n (v :: * -> *).
Ord n =>
Point v n -> NonEmpty n -> PolarPt v n
mkPolarPt (Int -> [Point v n] -> [Point v n]
forall a. Int -> [a] -> [a]
extend Int
2 (BSpline v n -> [Point v n]
forall a. NonEmpty a -> [a]
NE.toList BSpline v n
controls)) (Int -> [n] -> [NonEmpty n]
forall a. Int -> [a] -> [NonEmpty a]
windows Int
3 ([n] -> [NonEmpty n]) -> [n] -> [NonEmpty n]
forall a b. (a -> b) -> a -> b
$ Int -> [n] -> [n]
forall a. Int -> [a] -> [a]
extend Int
2 [n]
knots)

  -- The bezier internal control points are affine combinations of
  -- the spline control points.
  bezierControls :: [(PolarPt v n, PolarPt v n)]
bezierControls = (NonEmpty (PolarPt v n) -> (PolarPt v n, PolarPt v n))
-> [NonEmpty (PolarPt v n)] -> [(PolarPt v n, PolarPt v n)]
forall a b. (a -> b) -> [a] -> [b]
map NonEmpty (PolarPt v n) -> (PolarPt v n, PolarPt v n)
forall {v :: * -> *} {n}.
(Additive v, Fractional n, Ord n) =>
NonEmpty (PolarPt v n) -> (PolarPt v n, PolarPt v n)
combineC (Int -> [PolarPt v n] -> [NonEmpty (PolarPt v n)]
forall a. Int -> [a] -> [NonEmpty a]
windows Int
2 [PolarPt v n]
controls')
  combineC :: NonEmpty (PolarPt v n) -> (PolarPt v n, PolarPt v n)
combineC (PolarPt v n
pabc :| [PolarPt v n
pbcd]) = (Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
forall (v :: * -> *) n.
(Additive v, Fractional n, Ord n) =>
Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
combine Int
0 PolarPt v n
pabc PolarPt v n
pbcd, Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
forall (v :: * -> *) n.
(Additive v, Fractional n, Ord n) =>
Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
combine Int
1 PolarPt v n
pabc PolarPt v n
pbcd)
  combineC NonEmpty (PolarPt v n)
_ = [Char] -> (PolarPt v n, PolarPt v n)
forall a. HasCallStack => [Char] -> a
error [Char]
"combineC must be called on a list of length 2"

  -- The bezier end points are affine combinations of the bezier
  -- control points.
  bezierEnds :: [PolarPt v n]
bezierEnds = (NonEmpty (PolarPt v n, PolarPt v n) -> PolarPt v n)
-> [NonEmpty (PolarPt v n, PolarPt v n)] -> [PolarPt v n]
forall a b. (a -> b) -> [a] -> [b]
map NonEmpty (PolarPt v n, PolarPt v n) -> PolarPt v n
forall {v :: * -> *} {n}.
(Additive v, Fractional n, Ord n) =>
NonEmpty (PolarPt v n, PolarPt v n) -> PolarPt v n
combineE (Int
-> [(PolarPt v n, PolarPt v n)]
-> [NonEmpty (PolarPt v n, PolarPt v n)]
forall a. Int -> [a] -> [NonEmpty a]
windows Int
2 [(PolarPt v n, PolarPt v n)]
bezierControls)
  combineE :: NonEmpty (PolarPt v n, PolarPt v n) -> PolarPt v n
combineE ((PolarPt v n
_, PolarPt v n
pabb) :| [(PolarPt v n
pbbc, PolarPt v n
_)]) = Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
forall (v :: * -> *) n.
(Additive v, Fractional n, Ord n) =>
Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
combine Int
0 PolarPt v n
pabb PolarPt v n
pbbc
  combineE NonEmpty (PolarPt v n, PolarPt v n)
_ = [Char] -> PolarPt v n
forall a. HasCallStack => [Char] -> a
error [Char]
"combineE must be called on a list of length 2"

  -- Finally, we actually put together the generated bezier segments.
  beziers :: [FixedSegment v n]
beziers = ((PolarPt v n, PolarPt v n)
 -> NonEmpty (PolarPt v n) -> FixedSegment v n)
-> [(PolarPt v n, PolarPt v n)]
-> [NonEmpty (PolarPt v n)]
-> [FixedSegment v n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (PolarPt v n, PolarPt v n)
-> NonEmpty (PolarPt v n) -> FixedSegment v n
forall {v :: * -> *} {n}.
(PolarPt v n, PolarPt v n)
-> NonEmpty (PolarPt v n) -> FixedSegment v n
mkBezier (Int -> [(PolarPt v n, PolarPt v n)] -> [(PolarPt v n, PolarPt v n)]
forall a. Int -> [a] -> [a]
drop Int
1 [(PolarPt v n, PolarPt v n)]
bezierControls) (Int -> [PolarPt v n] -> [NonEmpty (PolarPt v n)]
forall a. Int -> [a] -> [NonEmpty a]
windows Int
2 [PolarPt v n]
bezierEnds)
   where
    mkBezier :: (PolarPt v n, PolarPt v n)
-> NonEmpty (PolarPt v n) -> FixedSegment v n
mkBezier (PolarPt v n
paab, PolarPt v n
pabb) (PolarPt v n
paaa :| [PolarPt v n
pbbb]) =
      Point v n
-> Point v n -> Point v n -> Point v n -> FixedSegment v n
forall (v :: * -> *) n.
Point v n
-> Point v n -> Point v n -> Point v n -> FixedSegment v n
FCubic (PolarPt v n -> Point v n
forall (v :: * -> *) n. PolarPt v n -> Point v n
unPP PolarPt v n
paaa) (PolarPt v n -> Point v n
forall (v :: * -> *) n. PolarPt v n -> Point v n
unPP PolarPt v n
paab) (PolarPt v n -> Point v n
forall (v :: * -> *) n. PolarPt v n -> Point v n
unPP PolarPt v n
pabb) (PolarPt v n -> Point v n
forall (v :: * -> *) n. PolarPt v n -> Point v n
unPP PolarPt v n
pbbb)
    mkBezier (PolarPt v n, PolarPt v n)
_ NonEmpty (PolarPt v n)
_ = [Char] -> FixedSegment v n
forall a. HasCallStack => [Char] -> a
error [Char]
"mkBezier must be called on a list of length 2"

-- Note that the above algorithm works in any dimension but is
-- very specific to *cubic* splines.  This can of course be
-- generalized to higher degree splines but keeping track of
-- everything gets a bit more complicated; to be honest I am not
-- quite sure how to do it.

-- | Generate a uniform cubic B-spline from the given control points.
--   The spline starts and ends at the first and last control points,
--   and is tangent to the line to the second(-to-last) control point.
--   It does not necessarily pass through any of the other control
--   points.
--
--   <<diagrams/src_Diagrams_CubicSpline_Boehm_bsplineEx.svg#diagram=bsplineEx&width=300>>
--
--   > pts = map p2 [(0,0), (2,3), (5,-2), (-4,1), (0,3)]
--   > spot = circle 0.2 # fc blue # lw none
--   > bsplineEx = mconcat
--   >   [ position (zip pts (repeat spot))
--   >   , bspline pts
--   >   ]
--   >   # frame 0.5
bspline :: (TrailLike t, V t ~ v, N t ~ n) => BSpline v n -> t
bspline :: forall t (v :: * -> *) n.
(TrailLike t, V t ~ v, N t ~ n) =>
BSpline v n -> t
bspline = Located [Segment Closed v n] -> t
Located [Segment Closed (V t) (N t)] -> t
forall t. TrailLike t => Located [Segment Closed (V t) (N t)] -> t
fromLocSegments (Located [Segment Closed v n] -> t)
-> (BSpline v n -> Located [Segment Closed v n])
-> BSpline v n
-> t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Located (Segment Closed v n)] -> Located [Segment Closed v n]
forall {a}.
(Additive (V a), Num (N a)) =>
[Located a] -> Located [a]
fixup ([Located (Segment Closed v n)] -> Located [Segment Closed v n])
-> (BSpline v n -> [Located (Segment Closed v n)])
-> BSpline v n
-> Located [Segment Closed v n]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (FixedSegment v n -> Located (Segment Closed v n))
-> [FixedSegment v n] -> [Located (Segment Closed v n)]
forall a b. (a -> b) -> [a] -> [b]
map FixedSegment v n -> Located (Segment Closed v n)
forall n (v :: * -> *).
(Num n, Additive v) =>
FixedSegment v n -> Located (Segment Closed v n)
fromFixedSeg ([FixedSegment v n] -> [Located (Segment Closed v n)])
-> (BSpline v n -> [FixedSegment v n])
-> BSpline v n
-> [Located (Segment Closed v n)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. BSpline v n -> [FixedSegment v n]
forall (v :: * -> *) n.
(Additive v, Fractional n, Num n, Ord n) =>
BSpline v n -> [FixedSegment v n]
bsplineToBeziers
 where
  fixup :: [Located a] -> Located [a]
fixup [] = [] [a] -> Point (V [a]) (N [a]) -> Located [a]
forall a. a -> Point (V a) (N a) -> Located a
`at` Point (V a) (N a)
Point (V [a]) (N [a])
forall (f :: * -> *) a. (Additive f, Num a) => Point f a
origin
  fixup (Located a
b1 : [Located a]
rest) = (Located a -> a
forall a. Located a -> a
unLoc Located a
b1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (Located a -> a) -> [Located a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Located a -> a
forall a. Located a -> a
unLoc [Located a]
rest) [a] -> Point (V [a]) (N [a]) -> Located [a]
forall a. a -> Point (V a) (N a) -> Located a
`at` Located a -> Point (V a) (N a)
forall a. Located a -> Point (V a) (N a)
loc Located a
b1