# Exercise for the Reader

## December 5, 2009

### B-Splines: Better Math

Filed under: 3D Graphics, Numerical methods — Seth Porter @ 9:13 pm

In the previous three parts of this series, we’d sorted out why B-splines are interesting curves, attempted to sort out B-splines terminology, and decided that it would be really nice to be able to evaluation them with uniform functions. I’d tried some quick and dirty approaches to force it to be, and had gotten little in response. It was time to actually do some algebra and figure out what was going on.

## A Note About Tools

That means there are going to be a lot of LaTex-created images in this section. Sorry, but it’s a lot less painful than trying to write out math in raw HTML. (Okay, writing this inspired me to search for LaTex to MathML translators, and of course the W3 itself maintains quite a nice list. Oh well, I promise I’ll look into that for next time.) In the meantime, I’m using the wiki built into Trac to write this post, with the TracMath plugin to take inline LaTeX blocks and render them as PNGs. This means that I just type

#!latex
$$N_{i,k}(v)={{\left( T_{k+i} - v\right)\,N_{i+1,k-1}(v)}\over{T_{k+i}-T_{i+1}}} + {{\left(v-T_{i}\right)\,N_{i,k-1}(v)}\over{T_{k+i-1}-T_{i} }}$$


and I get

Well, okay, I don’t write that, because I can’t remember LaTeX math syntax to save my life. (I got pretty good at it in college, but it’s been a while.) Anyway, it’s hard to parse by eye, and I’m very prone to swap a u for a v, which can be unfortunate when you’ve just carefully drawn a distinction between them.

Instead, I’m using Maxima, a free Computer Algebra System. It’s stranger to pick up than Mathematica, but the price is right (and it’s quite powerful once you get used to its ways). It seems to be really flaky on recent Ubuntu releases, which is unfortunate, but I’m primarily using 5.17.0 for Windows and it works like a champ. It has a very nice tex function, which will spit out the TeX formatting for any given expression.

So in practice, I work the math in Maxima, have it convert to TeX, then paste that into the post. (Sometimes I tweak it slightly, for example to pull out a leading “1/2” rather than leaving the whole expression over 2 as Maxima seems to prefer. But for the most part, the expressions in the text are directly out of Maxima.) I have some ideas of writing a Trac plugin to close the loop and directly evaluate and typeset Maxima expressions in the text, but I haven’t got it all worked out yet: the equations in a piece generally build on each other, so you need enough lifecycle awareness to start and maintain a separate Maxima state for each page in flight. Also you probably need some non-displayed helper definitions and intermediate values, so you need a separate bit of syntax to indicate what should be rendered into the page. I think the fact of the existence of the PageOutlineMacro means it’s possible, but the fact that the PageOutlineMacro actually re-parses the whole page means it’s not necessarily easy. (I’d also still have to solve creating and maintaining enough state that the various individual blocks could fetch their expression / formatted results from the initial gather-and-run macro.) Anyway, for the moment I’ll do it the hard way.

## Back to the Math

I’m going to run through a concrete k=3 (quadratic) example with 6 control points (so n=5), and claim that this demonstrates a general solution for k=3. I’m pretty sure it does, but it’s a lot easier to keep track of the indexes with a concrete number. Or so I claim.

So, the knot vector T should be the sequence [0, 0, 0, 1, 2, 3, 4, 4, 4]. Recall that indices are zero-based, and that Ni,1 is defined as

Given this, the Ni,1 terms will start being interesting at i=2, where N2,1= 1 for 0 <= u < 1. From then on, each Ni,1 is non-zero on Ni,1= 1 for i-2 <= u < i-2+1, up through the last value we need at N5,1. Using Maxima to evaluate the raw B-spline equation

for our knot vector T (again, this is the k=3 case), we get a whomping huge expression in terms of the control points P, the Ni,1 “switches”, and of course v. I won’t reproduce the whole thing here (fully expanded, there are 33 terms, and even my new wider columns can’t handle that), but here’s the coefficient of the second control point P1 for flavor:

### Explicit Equations for Each Segment

The good thing is that we don’t ever really need the whole equation at once. Only one of the Ni,1 terms will be 1 at any given time — this is what produces the curve segments. By setting N2,1 through N5,1 to 1 in turn (and the other Ni,1 terms to 0) we can produce the four curve segments. In turn, these are

### Maxima Script

For anyone who’s following along at home, here’s the Maxima to get our results so far:

N(N1, T, i, k, v) := if k = 1 then N1[i] else
(if T[i+k-1] = T[i] then 0 else (v- T[i])* N(N1, T, i, k-1, v) / (T[i+k-1] - T[i]))
+
(if T[i+k] = T[i+1] then 0 else (T[i+k] - v) * N(N1, T, i+1, k-1, v) / (T[i+k] - T[i+1]));

BSpline(P, N1, T, n, k, v) := sum(P[i] * N(N1, T, i, k, v), i, 0, n);
T5Vals : [0, 0, 0, 1, 2, 3, 4, 4, 4];
T5[i] := T5Vals[i+1];

BSpline5_3(v) = BSpline(P, N1, T5, 5, 3, v);

makeN1Helper(curveSeg, n, k) := makelist(if j - k + 3 = curveSeg then 1 else 0, j, 0, n-1);
BSpline5_3_Seg(seg, P, u) := BSpline(P, makeN1Helper(seg, 5, 3), T5, 5, 3, u);
for seg : 1 step 1 thru 4 do disp(
segment[seg] = facsum(ratsimp(
BSpline5_3_Seg(seg, P, v)), P[seg-1], P[seg], P[seg+1]));


Note that Maxima lists are one-based, hence the helper function T5 to preserve the zero-based indexing in the original definitions. Since Ni,0 is never referenced, I can get away with a list starting at index 1 (though I admit this is made more confusing by the choice of a zero-based loop variable. Oops.) Also note that BSpline5_3 wasn’t actually defined as a function there — that was just an equality, to make Maxima print it out nicely. Other than that, this is pretty much what it looks like. Oh, and notice that the recursive N definition has to avoid a couple of divide by zero cases to implement the “where 0/0 is defined to be 0” escape hatch.

### Reparamaterization

Okay, those equations work but they look pretty ugly (and they’re different in each segment, but we’ll work on that in a minute). It took me a while to sort this out, but I sort of tipped my hand back in the “Definitions” post where I introduce w and j, so I’ll just repeat the equation from there: v = u (n-k+2) = j+w-1. Remember that j is the segment number ranging from 1 to n-k+2, which is exactly the subscript of segment in the explicit equations above. So we substitute w-1 in the first equation, 1+w-1 in the second and so forth. Using the Maxima definitions above, this is just evaluating the BSpline5_3_Seg function above at seg+w-1 rather than v:

for seg : 1 step 1 thru 4 do disp(
segment[seg] = facsum(ratsimp(
BSpline5_3_Seg(seg, P, seg+w-1)), P[seg-1], P[seg], P[seg+1]));


(the facsum bit just rearranges the terms) yielding

Now we’re getting somewhere. The middle two equations are identical except for shifted subscripts on P. The first two and last two are different, but we can work that out in a minute. First let’s generalize the uniform version and check our claim that

works for those middle terms. This will also give us a peek at the what’s going on with the end segments.

B3[uniform](P, j, w) := (-P[j]*(2*w^2-2*w-1)+P[j+1]*w^2+P[j-1]*(w-1)^2)/2;
for seg : 1 step 1 thru 4 do disp(
difference[seg] = ratsimp(
BSpline5_3_Seg(seg, P, seg+w-1) - B3[uniform](P, seg, w))
);


yielding

What we’re computing here is the difference between the “real” equation for each of the four segments and what we get from using the uniform equation for all of the segments (ignoring the segment number).

Doesn’t look immediately promising, but there are two good things:

• the two middle segments are as predicted by the uniform function above
• the two end segments’ differences only involve two control points

We know that P1 and P4 are used in the middle two segments as well as the end segments, so we can’t do anything to them. However, the two end points (P0 and P5) are only used in the end segments. So what it we feed some different value (not the current P0, but some combination of P0 and P1) into the uniform equation? Can we make the result be the same as evaluating the “real” B-spline function on the original control points?

In the uniform version of segment 1, we replace P0 with a new variable x, then set it equal to the real first segment of the k=3 B-spline:

We then solve for x hoping we can eliminate all the w terms:

rightAnswer : ratsimp(BSpline5_3_Seg(1, P, 1+w-1));
uniformWithFreeP0 : ratsubst(x, P[0], B3[uniform](P, 1, w));
solve(uniformWithFreeP0 = rightAnswer, x);


and we get (drumroll, please)

(Okay, that didn’t really need a LaTeX block, but it really did take me quite a while to figure out, and I was inexplicably happy when I found it (literally: one very really reason for this whole series is to explain to my wife why I was grinning so much that week). If I’d happened to have run into that formula up front I would have spared myself a great deal of frustrating math trying to figure out how to make sense of what I was seeing. Of course, I also wouldn’t have learned a lot about splines, so maybe it was better this way.)

Repeating the process for the last control point (in this case P5), we get x = 2 P5-P4:

rightAnswer2 : ratsimp(BSpline5_3_Seg(4, P, 4+w-1));
uniform2WithFreeP0 : ratsubst(x, P[5], B3[uniform](P, 4, w));
solve(uniform2WithFreeP0 = rightAnswer2, x);


This is a reassuring result — the equations are all symmetric, so it should be possible to reverse the control points and produce the same curve. Therefore, it’s nice to see that this formulation for the end-points is also symmetric.

So the practical upshot of this is that we now have a way to evaluate all the segments of a B-spline, even the “weird” end segments, with the same uniform equation. To do this, we have to massage the control points a little. Given the input sequence

( P0 , P1 , … , Pn-1 , Pn )

we need to make a copy, so we can feed the uniform function the sequence

( 2 P0 – P1 , P1 , … , Pn-1 , 2 Pn – Pn-1 )

instead. That doesn’t seem like much of a price to pay for freedom (or at least for being able to evaluate all segments of all k=3 B-splines curves with the same straight-line evaluation code).

The same general approach works for k=4 curves (cubic B-splines). Because the cubic blending functions involve more points in each segment, we have to replace both P0 and P1 with free variables, and solve for them simultaneously (using the known equations for the first two segments, instead of just the first segment), but it’s the same logic and it works quite nicely. I haven’t tried to generalize a formula for any value of k, mostly because quadratic and cubic B-splines are far more common than other types.

This took much longer to post than I’d planned, but I hope it was worth the wait for someone. I can absolutely guarantee that this isn’t an original result, but I haven’t actually seen it in Wikipedia or MathWorld or the usual suspects, so I feel entitled to be slightly smug.