Wildly off-topic: de Boor's algorithm

older
what's the best environment for...

Andrew Coppin

15 Jul 2010 15 Jul '10

1:29 p.m.

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?

Show replies by date

mokus＠deepbondi.net

24 Jul 24 Jul

4:33 p.m.

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

Andrew Coppin

5 Aug 5 Aug

4:31 p.m.

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

James Andrew Cook

6 Aug 6 Aug

8:38 a.m.

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

...

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

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

mokus＠deepbondi.net

10:03 a.m.

mokus@deepbondi.net wrote:

...

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.

How embarrassing, I managed to get this simple math wrong. That's what I get for trying to think in the morning without either my notes or my coffee, I suppose. At least, I have that standard excuse to fall back on, so I'll take advantage of it. "n = m - 1 + p" should read "n + p = m - 1" - I added "p" to the wrong side. It is indeed the property that makes the zip line up, but not by the math I gave. Instantiating n and m as the respective "length" expressions makes the equation:

...

length ds + p = length us - 1

With subsequent derivation in light of the fact that the 1st control point is being ignored (which introduces a "- 1" on the RHS, if I'm not mistaken again), as mentioned earlier, and the fact that the expression in question is supposed to return a list one shorter than ds:

...

length ds + p - 1 = length us - 1 = length (tail us) length ds - 1 = length (tail us) - p length (spans 1 ds) = length (spans p (tail us))

(Under the assumption that all the lists involved are long enough for the tail and drop operations) And all is (hopefully) right again. -- James

Andrew Coppin

10:04 a.m.

mokus@deepbondi.net wrote:

...

How embarrassing, I managed to get this simple math wrong. That's what I get for trying to think in the morning without either my notes or my coffee, I suppose.

I _said_ it was tricky, didn't I? ;-)

5402

Age (days ago)

5424

Last active (days ago)

List overview

Download

5 comments

3 participants

participants (3)

Andrew Coppin
James Andrew Cook
mokus＠deepbondi.net