Sign in to follow this  
yoavhacohen

Monotonic cubic interpolation

Recommended Posts

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]

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this