From: Casey Muratori Discontinuous curve report [LONG] 2003-06-06 15:25 I recently updated my curve fitter in Granny to support animations with 0th order and 1st order discontinuities, since there actually have been a fair number of customers using this sort of thing (first order is obvious - bouncing balls, etc. - but zeroth order was used to do teleportation of stuff, which people like to use for tricks and such). I rewrote the curve fitter entirely and have some things to report that I think might be interesting. I never see this stuff written in papers, so I think it's all "trade knowledge", sadly. Or maybe I am just not reading the right papers (usually this is where Peter-Pike jumps in with some good references :) First, some background. The original curve fitter in Granny is designed to work by oversampling, by at least 2x, some set of data from an art package like MAX or Maya. So, if the artist makes a 30 frame animation, the exporters sample it at 60 frames and then give it to the curve fitter. This way, it's never ambiguous what was supposed to happen _in between_ the frames - you know, nyquist and all that. You can oversample by more than that, but it's not strictly necessary. It produces a least-squares-optimally-fit bspline curve (NUNRBS - Non Uniform Non Rational B-spline) to the samples for some prescribed error tolerance. This way, the engine only ever sees nice, easy-to-evaluate b-splines, no matter what random-ass plug-ins or special sauce the artists may have used to generate their motion, AND it always guaranteed to be as continuous as you want (ie., pick the degree, and you have one less than that many continuous derivatives), so you can slow down time to any timestep you want and still non-popping, non-jerky motion. A plus for bullet-time fans, which apparently is everybody. Now, the problem with all this smoothness is, of course, if the animation wasn't supposed to be smooth. So, what I wanted to do, for this particular update, was modify the b-spline processing so that it could be given certain user-specified tolerance levels, and any time these tolerances were violated, it would preserve the violation as a discontinuity in the curve. Since I work on a licensable product, an ancillary goal here was that if I could make all this work without changing the run-time side of things at all, then people could continue to freely mix-and-match versions of our exporter with versions of our run-time library, which is always a plus for development. Luckily for me, it just so happens that with b-spline curves, if you insert two knots at the same place in the curve, you reduce the continuity of the curve by one degree. So, a quadratic b-spline with two knots in one place is a piecewise linear curve. Insert a third, and its a step curve. How convenient. You would ordinarily think that you would have to change the way you evaluate the curve for this to work, but it just so happens (conveniently again), that you don't. Because you always search for a knot range in the curve that is greater-than-or-equal-to the current time, as soon as the current time is greater than the double (or tripled) knot, it skips _all_ those knot values and moves on to the next knot range, automatically. So you are never evaluating the region _in between_ two knots with the same time value, you never divide by zero (which would otherwise happen), and you are a happy puppy. If you don't intuitively know why this is, write out the NUNRBS blend pyramid and see what happens for yourself. You always dodge the bullet at each rung. So, really, what this all boiled down to was the fact that I needed to modify the solver to be able to take curves that had doubled or tripled (etc.) knots in them. Since it is a least-squares solver, this means it's trying to solve the equation A^T A x = A^T b, where A is a matrix of b-spline coefficients (it's samples times knots in size), b are the samples, and x is the resulting control points for the knot times. Note that I am NOT solving for knot times, because that would be non-linear. I use heuristics to pick knot times, and just run the solver multiple times (since it's O(n) in the number of samples, and extremely fast). Under normal (NON-discontinuous) solves, you have an A matrix where each row has degree + 1 entries (degree being the degree of the spline), since that's how many knots are involved in each b-spline evaluation (and hence affect each sample). The degree + 1 blocks span as many rows as there are samples that fall within the same knot range, and then the next block starts one column to the right, and so on (this is obvious from the way b-splines work). So you are dealing with matrices that look like this (here I show quadratic b-splines): * * * - - - * * * - - - * * * - - - - * * * - - - - * * * - - - * * * - - - * * * - - - * * * - - - - * * * - - - * * * Where the *'s are non-zero entries, and the - are zeros. When you square a matrix of this form, you get this: * * * - - + * * * - + + * * * - + + * * - - + + * for obvious reasons. It's symmetric (the +'s are not new informatiom, just transposes of the stuff to the right of the diagonal), and the degree + 1 bands of information starting from the diagonal on each row are all you need to know about it to know it fully. Now, what happens when you double or triple a knot is that you skip that many additional columns in A: * * * - - - * * * - - - * * * - - - - - * * * - <- ha HA! KNOT DOUBLED! - - * * * - - - * * * - - - * * * - - - * * * - - - - * * * - - - * * * Which then shrinks your bandwidth in A^T A around that knot: * * * - - + * * - - + + * * * - - + * * - - + + * Now, there are two problems here. The first problem is that you will get cancellation as you lose bandwidth, due to the nature of the b-splines. The second problem is that if you have two doubled knots _in a row_, you will get actual zeroes in the matrix along the diagonal, also a problem. But is it really a problem? No. It is not a problem. The reason it is not a problem is because of the nature of the b-splines. When you have such agressive knot doubling going on, the shrinking of the matrix produces _banded_ shrinking, so the worse you will get is this kind of thing: * * * - - - - + * * - - - - + + * - - - - - - - - - - - - - - - * * * - - - - + * * - - - - + + * where the missing row could be any number of rows (ie., it will still be square, you'll just have that many empty columns and rows around that point, depending on the degree of over-knottage). But that's a fine looking matrix, if I do say so myself. And I do. Because all you have to do to handle this kind of matrix is modify your factoring algorithm to ignore the rows that have no data. This is essentially equivalent to pivoting them around: * * * - - - - + * * - - - - + + * - - - - - - - * * * - - - - + * * - - - - + + * - - - - - - - - But why bother pivoting when you don't have to pivot? That's what I say. During solving, you simply substitute in _reasonable_ knot values any time you come to a row with no data. This is perfectly happy and valid, because you are in the null-space of the matrix, so you're welcome to do whatever you please. Using the same value as the next knot (which you've already solved if you're back-substituting), for example, is a fine choice, and what I picked. This nice part of doing it like this, with no actual pivoting and such, is that all you have to do is introduce a simple epsilon check when you use the diagonal, and substitute in a suitable alternative path at that time, rather than having to deal with the complexities normally involved in pivotting. Now, I know what you are all saying at this point. You are saying, "But Casey - WHY have you not simply broken the curve into continuous sections, solved THOSE, and then reassembled them into a single curve?" I don't blame you for saying that, because that would be a very reasonable strategy. However, what that fails to take into consideration is that the extra DOFs introduced by the double knots _can be used by the solver_ to improve the accuracy of the curve. So, allowing the solver to solve for _as many values as it can_, always guarantees that it will produce the best curve for our money, where by "best" I mean "least-squares closest to the samples", and by "money" I mean "number of knots". Because remember, just because a knot is doubled or tripled _doesn't mean the control points are the same_... it just means the knot value is. So you want to let the solver pick _multiple control points_ at that knot value, so that it can control the left-side or right-side derivatives, since hey, it's got an extra knot to do that with, so why let it go to waste? I would also like to mention in passing, that whilst doing this, I took the time to rewrite the entire curve solver (which dated back to 2000), and I am very happy to report that _this_ time, since I knew so much about the problem after having worked with it exclusively for a whole month back then (when I originally wrote it), that I was able to write it completely deterministically sparse, meaning that no matrix is ever actually constructed, and no ranges are ever remembered. It is an entirely fixed-size, perfectly packed algorithm that does the minimal amount of work and completes all operations at the soonest possible time, even though it is _still_ completely general (ie., you can ask it to solve for _any_ degree of spline - although presumably you would not ask for degrees > 3). It was a lot of fun and runs over ten times faster than my original one, which used a more general "sparse-matrix class" kind of approach (which is ass). So I am very proud of that and filled with the glee that one would expect to be filled with if one became filled with glee at such things. Now, what I am NOT particularly happy with, on the other hand, is that I do not know of a guarantably stable way to _automatically_ determine the tolerances for determining what is and is not a discontinuity. It seems like maybe there is no way, because depending on what the artist intended, two sets of samples with nearly identical jumps might be meant to be interpreted totally differently. So, it may be that you always want the user to specify tolerances, or perhaps even for them to mark frames where they wish the discontinuities to be (which I could imagine doing in the future by having special-purpose queries into the art tool to look for purposefully-set discontinuous places in the f-curves). The solver doesn't care where you tell it to double up on the knots, so it'd be easy to put something more sophisticated in along these lines. Cheers, - Casey