There are few alternative solutions I'm aware of, but which also have some drawbacks:
- Iterate the spline with constant time step from the beginning of the spline (or previous position) while evaluating the length until you find the segment in the spline that's given distance away from the beginning. Then linearly interpolate between the endpoints of the found segment to get C0 continuously interpolated value. The drawback is that you have to incrementally iterate over the spline which with larger steps can be time consuming. This isn't really a big issue with a single camera track, but if you have to interpolate many splines in a frame (e.g. for character animation) it can be. It's also necessary to maintain state information about the interpolation to gain from temporal coherency in interpolation, which complicates the operation. The biggest drawback of this method is thus unpredictable performance characteristics, i.e. from one call to the next there can be several orders of magnitude difference in CPU usage depending on the interpolation precision requirements and the length of the step, which is bad particularly for games which should aim for steady performance.
- Pre-calculate look-up-table to map distance on the spline to a time value. While interpolating over the spline get segment from the table for the required distance and linearly interpolate the time value between previous and next values in the table to get C0 continuous interpolation. Then use this value to interpolate the spline. Main drawback is the memory requirements for storing the table, whose size depends on the desired precision.
- Interpolate over the spline by dividing the step length with spline speed, which you get by evaluating derivate of the spline at given time value. This technique doesn't survive through larger time steps, isn't frame rate independent and has issues with splines which have 0-derivates.
Basic Idea
While I was Googling for solutions, I stumbled upon a solution proposed by Nils Pipenbrinck in DevMaster forums, where you would create inverse cubic spline length function to approximate the mapping from (normalized) distance along on the spline to t. In other words, instead of trying to directly interpolate the spline with t, you first interpolate another spline which maps distance along the spline to the t value, which is then used to interpolate the original path spline. However, approximating the "distance to t" function with cubic polynomial showed quite a large error in some of my test situations so I ended up binning this exact solution.
Instead, I tried a variation of directly approximating the spline distance function (i.e. function which maps given t value to the distance on the spline) with cubic polynomial and upon interpolation calculating inverse of the polynomial to get proper t value for evaluating the spline. Because the distance function is approximated with a cubic polynomial there exists a closed form inverse function which is reasonably efficient to evaluate. The precision in this approach was clearly superior in my tests, so I ended up using this approach which has the same memory overhead as Nils' solution (i.e. two floats for tangent values and one for the spline length) even though the run-time CPU overhead is slightly bigger. For performance improvements it would be also possible to cache some of the coefficients needed for inverse calculation.
Calculating Spline Distance Function Parameters
To calculate the length of the spline in brute-force way is simple to implement by iterating over the spline and summing up the length of linear spline segments. Calculation of the tangents is slightly more involved though, which I implemented by applying cubic least-squares fitting to the spline distance function, i.e. approximating the function with cubic polynomial f(t)=at3+bt2+ct+d. The spline distance function maps given t value to the distance along the spline from the beginning of the spline, so this function can be easily created by stepping through the spline with constant delta t values and summing up the linear spline fragment lengths. The coefficient calculations for least-squares fitting can be done while evaluating the spline length, so there's no need to store the evaluated values or spline fragment lengths anywhere.
Below is the function for the Hermite spline tangent calculation in Spin-X Engine (If you are not familiar with templates, just think of "typename math<T>::scalar_t" as a verbose way of writing "float"
). This function can be used for any dimensional cubic spline as long as the type implements basic + and - operators, multiplication by scalar with * operator, norm() function and math<T> template class specialization, so it's possible for example to use it for 1d curve using float/double types or 2d & 3d curves using vec2f/vec2d & vec3f/vec3d/simd_vec3f types respectively defined in Spin-X Engine:Code cpp:
template<typename T> typename math<T>::scalar_t spline_distance_func(const cubic_spline<T> &s_, vec2<typename math<T>::scalar_t> &hermite_tangents_, unsigned num_steps_) { // calculate coefficients for least-squares cubic fit of the spline distance function (monotonic in range [0, 1]) PFC_ASSERT(num_steps_); typedef typename math<T>::scalar_t scalar_t; const scalar_t step=scalar_t(1)/scalar_t(num_steps_); T v0=s_.d; vec4<scalar_t> coeffs(0, 0, 0, 0); scalar_t len=0; scalar_t x0=0, x1=0, x2=0, x3=0, x4=0, x5=0, x6=0; for(unsigned i=1; i<=num_steps_; ++i) { // calculate cumulative distance and least-squares fit coefficients scalar_t si=scalar_t(i)*step, si2=si*si, si3=si2*si; T v1=evaluate(s_, si); len+=norm(v1-v0); v0=v1; coeffs+=vec4<scalar_t>(len*si3, len*si2, len*si, len); x0+=scalar_t(1.0); x1+=si; x2+=si2; x3+=si3; x4+=si2*si2; x5+=si3*si2; x6+=si3*si3; } // calculate cubic spline coefficients for normalized distance function coeffs*=rcp_z(len); mat44<scalar_t> m(x6, x5, x4, x3, x5, x4, x3, x2, x4, x3, x2, x1, x3, x2, x1, x0); coeffs*=inv(m); // calculate Hermite tangents for the spline distance approximation hermite_tangents_.x=coeffs.z; hermite_tangents_.y=coeffs.x*scalar_t(3.0)+coeffs.y*scalar_t(2.0)+coeffs.z; return len; }
Quick Introduction to Least-Squares Fit Algorithm
If you are not familiar what least-squares fit does or how it works, here's a short version: It takes a set of points and finds coefficients for a polynomial (that is, the values for a, b, c and d in the cubic polynomial mentioned above) that defines an Nth degree spline (3rd degree in our case, thus there are 4 coefficients) so that the sum of squared distances of the points from the spline is minimized (more precisely, vertical not perpendicular distances). It may sound complicated to implement, but the algorithm is quite simple once you get it.
To explain the algorithm shortly, you evaluate the function (calculate spline length up to a given t-value in our case), multiply that value with a coefficient vector whose length depends on the order of the polynomial we are fitting (e.g. for 3rd degree there's 4 coefficients, thus it's 4d vector), where each vector component holds the sample position raised to increasing power (e.g. when evaluating a point at t=0.2, then the vector holds values [0.23, 0.22, 0.21, 1] = [0.008, 0.04, 0.2, 1]) and add all these vectors together. You also add up another vector whose length for an Nth degree spline is 2N+1 and similarly to the other vector contains sample position rised to increasing power (e.g. for 3rd degree spline and t=0.2, the 7d vector contains values [0.26, 0.25, 0.24, 0.23, 0.22, 0.21, 1]).
Once all samples are calculated and added up to these two vectors, the larger vector is split up into a (N+1)x(N+1) diagonally symmetric matrix like shown below and the inverse of the matrix is calculated, which is multiplied by the coefficient vector to get the polynomial coefficients. Note that if you use fixed number of samples and positions, this matrix can be pre-calculated, so for higher degree polynomials you could use e.g. Matlab to pre-calculate the matrix instead of having to write code for N-dimensional matrix inverse. If you are still confused, you can have a look at the spline_distance_func() function for implementation details.
,where pi = sample point position and f(pi) = value at the sample position.
Cubic polynomial coefficients for f(t)=at3+bt2+ct+d are then calculated as follows:
Interpolation
Once we have pre-calculated Hermite spline tangents for the distance function, we need to generate cubic polynomial for it so that at f(0)=0 and f(1)=1. This can be done for example by using make_hermite_spline() function and setting p0=0, p1=1 and t0 & t1 arguments to the tangent values returned from spline_distance_func() function. While interpolating over the spline, the inverse of the cubic spline polynomial is calculated for given normalized distance value (i.e. divide the distance on the spline by its length, which is also returned by spline_distance_func()). The inverse calculation of cubic polynomial can be a bit heavy on math since it requires some understanding of complex number exponents, but you can check inv_evaluate() function for details below:
Code cpp:
template<typename T> T inv_evaluate(const cubic_spline<T> &ds_, T t_) { PFC_ASSERT_PEDANTIC(t_>=T(0.0) && t_<=T(1.0)); if(ds_.a) { // calculate inverse of the cubic spline distance function T r3a=rcp(T(3.0)*ds_.a), br3a=ds_.b*r3a; T q=ds_.c*r3a-br3a*br3a, q3=q*q*q; T r=T(1.5)*r3a*(t_-ds_.d+ds_.c*br3a)-br3a*br3a*br3a; T d=q3+r*r; if(d>=0) { // single root T s=sqrt(d); return cubicrt(r+s)+cubicrt(r-s)-br3a; } // triple root (get root in range [0, 1]) T sth, cth, sq=sqrt(-q); sincos(sth, cth, acos(r/(sq*sq*sq))*T(1.0/3.0)); float z=sq*(math<T>::sqrt3*sth-cth)-br3a; static const T epsilon=T(0.00001); if(z>=-epsilon && z<=T(1.0)+epsilon) return sat(z); z=-sq*(math<T>::sqrt3*sth+cth)-br3a; if(z>=-epsilon && z<=T(1.0)+epsilon) return sat(z); return T(2.0)*sq*cth-br3a; } if(ds_.b) { // calculate inverse of the quadratic spline distance function T sq=ds_.c*ds_.c-T(4.0)*ds_.b*(ds_.d-t_); if(sq<=0) return 0; return (ds_.c-sqrt(sq))/(T(2.0)*ds_.b); } // for the rest use linear f(t)=t return t_; }
The Results
Below are results of the approximation for couple of cubic splines. The red knots visualize how the values are distributed when interpolating the spline with constant delta t increments, while blue knots show more even distribution of the knots on the spline with the constant speed interpolation approximation presented here.
These two graphs below demonstrate the respective spline distance changes of the above cubic splines. The green line represents the ideal constant speed line, and more closely a curve follows it the more constant speed there is in the spline interpolation.
Menu
Categories
Upcoming Events
vBulletin Message