<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div><div>On Aug 5, 2010, at 4:31 PM, Andrew Coppin wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite"><div><a href="mailto:mokus@deepbondi.net">mokus@deepbondi.net</a> wrote:<br><blockquote type="cite">Andrew Coppin wrote:<br></blockquote><blockquote type="cite"> <br></blockquote><blockquote type="cite"><blockquote type="cite">Given a suitable definition for Vector2 (i.e., a 2D vector with the<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">appropriate classes), it is delightfully trivial to implement de<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">Casteljau's algorithm:<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">de_Casteljau :: Scalar -> [Vector2] -> [[Vector2]]<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">de_Casteljau t [p] = [[p]]<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">de_Casteljau t ps = ps : de_Casteljau t (zipWith (line t) ps (tail ps))<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">line :: Scalar -> Vector2 -> Vector2 -> Vector2<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">line t a b = (1-t) *| a + t *| b<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">Now, de Boor's algorithm is a generalisation of de Casteljau's<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">algorithm. It draws B-splines instead of Bezier-splines (since B-splines<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">generalise Bezier-splines). But I think I may have ACTUALLY MELTED MY<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">BRAIN attempting to comprehend it. Can anybody supply a straightforward<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">Haskell implementation?<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> <br></blockquote></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">It took me a while to get around to it, and another while to work it out,<br></blockquote><blockquote type="cite">but here's what I've come up with.<br></blockquote><br>OK, cool.<br><br><blockquote type="cite">First, here's a restatement of your<br></blockquote><blockquote type="cite">code with a concrete choice of types (using Data.VectorSpace from the<br></blockquote><blockquote type="cite">vector-space package)<br></blockquote><br>My types *are* concrete. They're from AC-Vector. ;-)<br></div></blockquote><div><br></div><div>Nope, they were abstract until you told me where I can get an implementation ;)</div><br><blockquote type="cite"><div><blockquote type="cite">To generalize to De Boor's algorithm, the primary change is the<br></blockquote><blockquote type="cite">interpolation operation. In De Casteljau's algorithm, every interpolation<br></blockquote><blockquote type="cite">is over the same fixed interval 0 <= x <= 1. For De Boor's algorithm we<br></blockquote><blockquote type="cite">need a more general linear interpolation on the arbitrary interval<br></blockquote><blockquote type="cite">[x0,x1], because all the interpolations in De Boor's recurrence are<br></blockquote><blockquote type="cite">between pairs of knots, which have arbitrary values instead of just 0 or<br></blockquote><blockquote type="cite">1.<br></blockquote><blockquote type="cite"> <br></blockquote><br>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.<font class="Apple-style-span" color="#000000"><font class="Apple-style-span" color="#144FAE"><br></font></font></div></blockquote><div><br></div><div>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.</div><div><br></div><div>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.</div><br><blockquote type="cite"><div><blockquote type="cite">Computing the table is now nearly as straightforward as in De Casteljau's<br></blockquote><blockquote type="cite">algorithm:<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite"> <br></blockquote><blockquote type="cite"><blockquote type="cite">deBoor p _ [] x = []<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">deBoor p (_:us) ds x = ds : deBoor (p-1) us ds' x<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> where<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> ds' = zipWith (interp x) (spans p us) (spans 1 ds)<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> <br></blockquote></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Making use of a simple list function to select the spans:<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite"> <br></blockquote><blockquote type="cite"><blockquote type="cite">spans n xs = zip xs (drop n xs)<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> <br></blockquote></blockquote><br>Don't you need an uncurry in there? Or am I missing something?<br></div></blockquote><div><br></div>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.</div><div><br><blockquote type="cite"><div><blockquote type="cite">Note that the algorithm does not make use of @us!!0@ at all.<br></blockquote><br>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...)<br><br><blockquote type="cite">I believe<br></blockquote><blockquote type="cite">this is correct, based both on the Wikipedia description of the algorithm<br></blockquote><blockquote type="cite">and the implementations I've seen.<br></blockquote><br>Heh, prove it. ;-)<br><br>(I guess just go draw some splines and see if they look... spliney.)</div></blockquote><div><br></div><div>Oh, believe me, I've done a lot of the informal stuff like that to check my sanity.</div><div><br></div><div>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. ;)</div><div><br></div><blockquote type="cite"><div><blockquote type="cite">Finally, the (very simple) driver to evaluate a B-spline:<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite"> <br></blockquote><blockquote type="cite"><blockquote type="cite">bspline p us ds = head . last . deBoor p us ds<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> <br></blockquote></blockquote><br>Yeah, figures.<br></div></blockquote><div><br></div>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.</div><div><br><blockquote type="cite"><div><br><blockquote type="cite">And a nearly-as-simple driver to evaluate a NURBS curve (type signature<br></blockquote><blockquote type="cite">included not because it couldn't be inferred but because it takes a bit of<br></blockquote><blockquote type="cite">sorting out to understand, so I wanted to present it in a slightly more<br></blockquote><blockquote type="cite">digestible form):<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite"> <br></blockquote><blockquote type="cite"><blockquote type="cite">nurbs :: (VectorSpace v, Scalar v ~ s,<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> VectorSpace s, Scalar s ~ Scalar v,<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> Fractional s, Ord s) =><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> Int -> [s] -> [(v, s)] -> s -> v<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">nurbs n us ds = project . bspline n us (map homogenize ds)<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> where<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> project (p,w) = recip w *^ p<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> homogenize (d,w) = (w *^ d, w)<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"> <br></blockquote></blockquote><br>Homogenous coordinate systems? Ouch.<br></div></blockquote><div><br></div><div>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.</div><br><blockquote type="cite"><div><br><blockquote type="cite">You can also split your B-spline just like a Bezier curve:</blockquote></div></blockquote>... (snip) ...<br><blockquote type="cite"><div>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?)<br></div></blockquote><div><br></div><div>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.</div><div><br></div><div>There are a fair number of other interesting properties and algorithms relating to B-splines at:</div><div><br></div><div><a href="http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/">http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/</a></div><div><br></div><div>I found these course notes to be a very useful resource.</div></div><br><div>-- James</div></body></html>