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

mokus@deepbondi.net wrote:

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, 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.

OK, cool.

First, here's a restatement of your code with a concrete choice of types (using Data.VectorSpace from the vector-space package)

My types *are* concrete. They're from AC-Vector. ;-)

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.

Right. That's basically what I figured. I'm having trouble wrapping my brain about the exact relationship between the number of control points, the number of knots, the degree of the spline, which knots and control points are active at a given parameter value, etc. There seems to be an utterly *huge* avenue for off-by-one errors here.

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).

Right. Because only a subset of the control points affect a given region of spline. (That's what makes B-splines superior, after all...)

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)

Don't you need an uncurry in there? Or am I missing something?

Note that the algorithm does not make use of @us!!0@ at all.

Yeah, that's the killer, isn't it? Figuring out exactly which combination of list function you need to pipe the correct values to the right places. (You might remember be asking about "expression dye" for this exact reason...)

I believe this is correct, based both on the Wikipedia description of the algorithm and the implementations I've seen.

Heh, prove it. ;-) (I guess just go draw some splines and see if they look... spliney.)

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.

Um... OK, I'm confused.

Finally, the (very simple) driver to evaluate a B-spline:

...

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

Yeah, figures.

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)

Homogenous coordinate systems? Ouch.

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

Now this I did not know... (Wikipedia doesn't list *all* properties that B-splines have. For example... how do you compute an "offset" curve for an arbitrary B-spline?)