Re: [Haskell-cafe] Wildly off-topic: de Boor's algorithm

Andrew Coppin wrote:

Given a suitable definition for Vector2 (i.e., a 2D vector with the appropriate classes), it is delightfully trivial to implement de Casteljau's algorithm:

de_Casteljau :: Scalar -> [Vector2] -> [[Vector2]] de_Casteljau t [p] = [[p]] de_Casteljau t ps = ps : de_Casteljau t (zipWith (line t) ps (tail ps))

line :: Scalar -> Vector2 -> Vector2 -> Vector2 line t a b = (1-t) *| a + t *| b

Now if you want to compute a given point along a Bezier spline, you can do

bezier :: Scalar -> [Vector2] -> Vector2 bezier t ps = head $ last $ de_Casteljau t ps

You can chop one in half with

split :: Scalar -> [Vector2] -> ([Vector2], [Vector2]) split t ps = let pss = de_Casteljau t ps in (map head pss, map last pss)

And any other beautiful incantations.

Now, de Boor's algorithm is a generalisation of de Casteljau's algorithm. It draws B-splines instead of Bezier-splines (since B-splines generalise Bezier-splines). But I think I may have ACTUALLY MELTED MY BRAIN attempting to comprehend it. Can anybody supply a straightforward Haskell implementation?

It took me a while to get around to it, and another while to work it out, but here's what I've come up with. First, here's a restatement of your code with a concrete choice of types (using Data.VectorSpace from the vector-space package) and a few minor stylistic changes just so things will line up with the generalized version better:

import Data.VectorSpace

interp a x y = lerp x y a

deCasteljau [] t = [] deCasteljau ps t = ps : deCasteljau (zipWith (interp t) ps (tail ps)) t

bezier ps = head . last . deCasteljau ps

split ps t = (map head pss, reverse (map last pss)) where pss = deCasteljau ps t

To generalize to De Boor's algorithm, the primary change is the interpolation operation. In De Casteljau's algorithm, every interpolation is over the same fixed interval 0 <= x <= 1. For De Boor's algorithm we need a more general linear interpolation on the arbitrary interval [x0,x1], because all the interpolations in De Boor's recurrence are between pairs of knots, which have arbitrary values instead of just 0 or 1. Because of Haskell's laziness, we can also take care of searching the result table for the correct set of control points at the same time, just by clamping the input to the desired interval and pre-emptively returning the corresponding 'y' if the input is outside that interval. This way, to find the final interpolant we only need to go to the end of the table, as in 'deCasteljau'. Unlike 'deCasteljau', only a portion of the table is actually computed (a triangular portion with the active control points as base and a path from the vertex of the triangle to the final entry in the table).

interp x (x0,x1) (y0,y1) | x < x0 = y0 | x >= x1 = y1 | otherwise = lerp y0 y1 a where a = (x - x0) / (x1 - x0)

Computing the table is now nearly as straightforward as in De Casteljau's algorithm:

deBoor p _ [] x = [] deBoor p (_:us) ds x = ds : deBoor (p-1) us ds' x where ds' = zipWith (interp x) (spans p us) (spans 1 ds)

Making use of a simple list function to select the spans:

spans n xs = zip xs (drop n xs)

Note that the algorithm does not make use of @us!!0@ at all. I believe this is correct, based both on the Wikipedia description of the algorithm and the implementations I've seen. De Boor's recurrence seems to require an irrelevant choice of extra control point that would be, notionally, @ds!!(-1)@. This control point has no actual influence inside the domain of the spline, although it /can/ affect values outside the domain (The domain being taken as the central portion of the knot vector where there are @p+1@ non-zero basis functions, forming a complete basis for the degree @p@ polynomials on that interval). This implementation makes the arbitrary but justifiable choice that the extra control point be identical to the first, so that the position of the first knot is irrelevant. It could alternatively be written to take the extra control point as an argument or as a part of @ds@, in which case the caller would be required to supply an additional control point that does not actually influence the spline. I initially found this result difficult to convince myself of even though it seems very plausible mathematically, because it seems to indicate that in general the position of the first and last knots are utterly irrelevant and I never saw any remarks to that effect in any of my reading on B-splines. Empirically, though, it seems to hold. Moving an internal knot at one end of a basis function does not alter the shape of that function in the segment furthest opposite, which is basically the same effect the first knot should have on the first basis function (and the opposite segment is the only one that falls inside the domain of the spline). It still may be that I'm wrong, and if anyone knows I'd love to hear about it, but I'm presently inclined to believe that this is correct. Finally, the (very simple) driver to evaluate a B-spline:

bspline p us ds = head . last . deBoor p us ds

And a nearly-as-simple driver to evaluate a NURBS curve (type signature included not because it couldn't be inferred but because it takes a bit of sorting out to understand, so I wanted to present it in a slightly more digestible form):

nurbs :: (VectorSpace v, Scalar v ~ s, VectorSpace s, Scalar s ~ Scalar v, Fractional s, Ord s) => Int -> [s] -> [(v, s)] -> s -> v nurbs n us ds = project . bspline n us (map homogenize ds) where project (p,w) = recip w *^ p homogenize (d,w) = (w *^ d, w)

For example, a NURBS circle (0 <= x <= 1):

circle :: Double -> (Double, Double) circle = nurbs 2 us ds where us = [0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1] ds = [ (( 1, 0),1) , (( 1,-1),w) , (( 0,-1),1) , ((-1,-1),w) , ((-1, 0),1) , ((-1, 1),w) , (( 0, 1),1) , (( 1, 1),w) , (( 1, 0),1) ] w = sqrt 0.5

You can also split your B-spline just like a Bezier curve:

simpleSplit p us ds t = (map head dss, reverse (map last dss)) where dss = deBoor p us ds t

The resulting B-splines use the same knot vector as the original and, in general, will have many duplicated control points. The ends of the splines can be cut off as long as you leave @p+1@ repeated control points, where @p@ is the degree of the spline. Alternatively, you can do a bit more work in the split function and cut the knot vector at the split point and insert @p+1@ copies of the split point into the resulting knot vector:

split p us ds t = ((us0, ds0), (us1, ds1)) where us0 = takeWhile (

us1 = replicate (p+1) t ++ dropWhile (<=t) us ds1 = reverse (trimTo (drop (p+1) us1) (map last dss))

dss = deBoor p us ds t

trimTo list xs = zipWith const xs list

-- James