Monotonic cubic interpolation

Started by
2 comments, last by Dave Eberly 14 years, 8 months ago
Hi, In this paper: http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/smoke.pdf In appendix B, there is a description of monotonic cubic interpolation. I have tried to implement it with the code below, but there are sometimes overshoots, although i set the slopes to zeros whenever they have a sign different from delta_k. What am I doing wrong? Thanks in advance, Yoav float MonotonicCubicInterpolation(float t, float fk2, float fk1, float fk, float fk_1) { float d_k0 = (fk1 - fk_1) / 2.0; float d_k1 = (fk2 - fk) / 2.0; float delta_k = (fk1 - fk); // This condition is suppose to eliminate overshoots: if (sign(d_k0) != sign(delta_k) || sign(d_k1) != sign(delta_k)) { d_k0 = 0; d_k1 = 0; } float a3 = d_k0 + d_k1 - 2*delta_k; // EDITED: it was: a3 = d_k0 + d_k1 - delta_k. Thanks to Dave float a2 = 3*delta_k - 2*d_k0 - d_k1; float a1 = d_k0; float a0 = fk; float result = a3*pow(t,3) + a2*t*t + a1*t + a0; if (abs(result - (fk1 + fk) / 2.0) > abs(fk1 - fk) / 2.0) { // overshoot, error! assert(false); } return result; } [Edited by - yoavhacohen on August 16, 2009 2:17:22 PM]
Advertisement
Should you have this instead? float a3 = d_k0 + d_k1 - 2*delta_k;
Thanks, I already saw that error and fix it, but the same problem is still there (the asser(false)).

Any additional suggestions?
Three issues.

First, the PDF has the error for a3 (which you now corrected). The PDF also assumes that t[k+1]-t[k] = 1 for all k. I do not know whether the authors intended this. If your t[k] sequence does not satisfy the difference conditions, all bets are off.

Second, I do not understand the necessary condition for monotonicity that is mentioned in the PDF. Those conditions involve only the signs of terms linear in the f[k], which does not seem right to me. For the cubic curve f(t) to be monotone on [t[k],t[k+1]], say, monotonic increasing, then f'(t) > 0 on this interval. That means the number of real-valued roots of f'(t) on the interval is zero. You can count the real-valued roots using a Sturm sequence of polynomials evaluated at t[k] and t[k+1]. If p0(t), p1(t), and p2(t) are the Sturm sequence for f'(t), the polynomial p2(t) is a constant that involves products of coefficients from f'(t), which are quadratic expressions (not linear).

Third, I do not understand your 'assert' and how it relates to monotonicity. This does not appear in the PDF you mention. Instead, you would need an assertion about the sign of f'(result) and sign of delta_k.

This topic is closed to new replies.

Advertisement