On Aug 5, 2010, at 4:31 PM, Andrew Coppin wrote:
> mo...@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
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
