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

[James Andrew Cook]

On Aug 5, 2010, at 4:31 PM, Andrew Coppin wrote:

> mokus at 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. ;-)

Nope, they were abstract until you told me where I can get an implementation ;)

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

It took me a while to get the intuition right on those, but here's a quick sketch.  Let n = number of control points, m = number of knots, and p = degree.  For p = 0 (constant segments), each control point corresponds to one span of the knot vector, so n = m - 1.  Each time you move up a degree, the basis functions span one more segment of the knot vector.  Thus, to keep the same number of control points you have to add a knot when you add a degree, so n = m - 1 + p.  Equivalently, n - p = m - 1, which incidentally makes an appearance in the deBoor function below when we zip "spans p us" with "spans 1 ds".  This is just the property that ensures that  the length of these lists is equal.

As for off-by-one errors, I agree.  I am reasonably sure that I at least don't have any major ones, but I would only be slightly surprised if there are corner cases I have handled incorrectly.  After quite a bit of plotting and deliberately introducing errors in the implementation I found that in most of the "obvious" places for off-by-one type errors (such as getting the wrong pair of knots for the pair of control points) the results were spectacularly sensitive to those errors, so at least those ones I think are hard to get wrong if you're actually looking at the splines the implementation draws.

>> 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?

Not unless I pasted the wrong version of the code.  I'm not sure where you mean - perhaps you're thinking of the arity of interp?  The interp function takes 3 arguments, 2 of which are pairs - those pairs are the ones being supplied by zipWith above, and are the results of this function.  I had a version that used "zipWith4 (interp t) us (drop p us) ds (tail ps)" instead but I felt it was overly complicated, so I abstracted out the "spans" function.

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

Oh, believe me, I've done a lot of the informal stuff like that to check my sanity.

The next few paragraphs that were apparently confusing (probably due to my habit of rambling at times) were a gesture in the direction I believe a proof is to be found.  If I feel ambitious I'll probably revisit that line of thinking and attempt to formalize it, now that I've been called out on it. ;)

>> Finally, the (very simple) driver to evaluate a B-spline:
>> 
>>  
>>> bspline p us ds = head . last . deBoor p us ds
>>>    
> 
> Yeah, figures.

Note that, as pretty and nice as it is to be able to make the driver that simple, it could most definitely be improved.  This simple version adds an unnecessary O(n-m) traversal back from the end of the table, for one thing, but more significantly it cannot handle infinite splines (precisely due to that traversal as well).  It also has to create list (:)-cells for many of the table entries it subsequently discards unevaluated.  A more clever implementation would extract precisely the control points (and corresponding set of knots) that contribute to the value of the spline at the point of evaluation, and call deBoor on that sub-spline.

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

They never bothered me that much.  In this case, I actually found them quite natural when I thought about the fact that NURBS were largely motivated (IIUC) by the desire to represent a larger set of conic sections.  In particular, with B-splines we can already represent a parabola, since it's just a quadratic.  If we put that spline in space, _on_ the plane that defines it as a conic section, then project it to a different plane, we have the conic section corresponding to that new plane.  The circle, in the case of homogeneous coordinates where everything is projected to w=1 and the projection is along the w-axis.

> 
>> You can also split your B-spline just like a Bezier curve:
... (snip) ...
> 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?)

I'm assuming here that you mean the curve defined by moving each point some fixed distance along the spline's normal.  In general, I doubt that that curve is even polynomial.  For example, B-splines can have cusps, and at those points the usual interpretation (which I assume is in some what mathematically natural, as a limit of the process or something) fillets those with a circle.   I may be quite wrong here, I am by no means an expert in geometry.

There are a fair number of other interesting properties and algorithms relating to B-splines at:

http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/

I found these course notes to be a very useful resource.

-- James
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.haskell.org/pipermail/haskell-cafe/attachments/20100806/01103c76/attachment-0001.html