{-# 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          (sort, tails)
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 = [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
ta -> a -> a
forall a. Num a => a -> a -> a
-a
a)a -> a -> a
forall a. Fractional a => a -> a -> a
/(a
ba -> a -> a
forall a. Num a => a -> a -> a
-a
a))

-- | @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] -> [[a]]
windows :: forall a. Int -> [a] -> [[a]]
windows Int
k = ([a] -> Bool) -> [[a]] -> [[a]]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
k) (Int -> Bool) -> ([a] -> Int) -> [a] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) ([[a]] -> [[a]]) -> ([a] -> [[a]]) -> [a] -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
k) ([[a]] -> [[a]]) -> ([a] -> [[a]]) -> [a] -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [[a]]
forall a. [a] -> [[a]]
tails

-- | @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
k [a]
xs = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k ([a] -> a
forall a. HasCallStack => [a] -> a
head [a]
xs) [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 ([a] -> a
forall a. HasCallStack => [a] -> a
last [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 -> [n]
_knots :: [n] }

mkPolarPt :: Ord n => Point v n -> [n] -> PolarPt v n
mkPolarPt :: forall n (v :: * -> *). Ord n => Point v n -> [n] -> PolarPt v n
mkPolarPt Point v n
pt [n]
kts = Point v n -> [n] -> PolarPt v n
forall (v :: * -> *) n. Point v n -> [n] -> PolarPt v n
PP Point v n
pt ([n] -> [n]
forall a. Ord a => [a] -> [a]
sort [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 [n]
kts1) (PP Point v n
pt2 [n]
kts2)
  = Point v n -> [n] -> PolarPt v n
forall n (v :: * -> *). Ord n => Point v n -> [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 ([n] -> n
forall a. HasCallStack => [a] -> a
head [n]
kts1) ([n] -> n
forall a. HasCallStack => [a] -> a
last [n]
kts2) n
newKt Point v n
pt1 Point v n
pt2)
      (n
newKt n -> [n] -> [n]
forall a. a -> [a] -> [a]
: Int -> [n] -> [n]
forall a. Int -> [a] -> [a]
drop Int
1 [n]
kts1)
  where
    newKt :: n
newKt = [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. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
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
1n -> 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 -> [n] -> PolarPt v n)
-> BSpline v n -> [[n]] -> [PolarPt v n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Point v n -> [n] -> PolarPt v n
forall n (v :: * -> *). Ord n => Point v n -> [n] -> PolarPt v n
mkPolarPt (Int -> BSpline v n -> BSpline v n
forall a. Int -> [a] -> [a]
extend Int
2 BSpline v n
controls) (Int -> [n] -> [[n]]
forall a. Int -> [a] -> [[a]]
windows Int
3 ([n] -> [[n]]) -> [n] -> [[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        = ([PolarPt v n] -> (PolarPt v n, PolarPt v n))
-> [[PolarPt v n]] -> [(PolarPt v n, PolarPt v n)]
forall a b. (a -> b) -> [a] -> [b]
map [PolarPt v n] -> (PolarPt v n, PolarPt v n)
forall {v :: * -> *} {n}.
(Additive v, Fractional n, Ord n) =>
[PolarPt v n] -> (PolarPt v n, PolarPt v n)
combineC (Int -> [PolarPt v n] -> [[PolarPt v n]]
forall a. Int -> [a] -> [[a]]
windows Int
2 [PolarPt v n]
controls')
    combineC :: [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 [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                   = ([(PolarPt v n, PolarPt v n)] -> PolarPt v n)
-> [[(PolarPt v n, PolarPt v n)]] -> [PolarPt v n]
forall a b. (a -> b) -> [a] -> [b]
map [(PolarPt v n, PolarPt v n)] -> PolarPt v n
forall {v :: * -> *} {n}.
(Additive v, Fractional n, Ord n) =>
[(PolarPt v n, PolarPt v n)] -> PolarPt v n
combineE (Int
-> [(PolarPt v n, PolarPt v n)] -> [[(PolarPt v n, PolarPt v n)]]
forall a. Int -> [a] -> [[a]]
windows Int
2 [(PolarPt v n, PolarPt v n)]
bezierControls)
    combineE :: [(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 [(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) -> [PolarPt v n] -> FixedSegment v n)
-> [(PolarPt v n, PolarPt v n)]
-> [[PolarPt v n]]
-> [FixedSegment v n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (PolarPt v n, PolarPt v n) -> [PolarPt v n] -> FixedSegment v n
forall {v :: * -> *} {n}.
(PolarPt v n, PolarPt v n) -> [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] -> [[PolarPt v n]]
forall a. Int -> [a] -> [[a]]
windows Int
2 [PolarPt v n]
bezierEnds)
      where
        mkBezier :: (PolarPt v n, PolarPt v n) -> [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)
_ [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