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

[Andrew Coppin]

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

> 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 (<t) us ++ replicate (p+1) t
>>         ds0 = trimTo (drop (p+1) us0) (map head dss)
>>
>>         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?)