Interpolating values between interval, interpolation as per Bezier curve

Question

To implement a 2D animation I am looking for interpolating values between two key frames with the velocity of change defined by a Bezier curve. The problem is Bezier curve is represented in parametric form whereas requirement is to be able to evaluate the value for a particular time.

To elaborate, lets say the value of 10 and 40 is to be interpolated across 4 seconds with the value changing not constantly but as defined by a bezier curve represented as 0,0 0.2,0.3 0.5,0.5 1,1. Now if I am drawing at 24 frames per second, I need to evaluate the value for every frame. How can I do this ? I looked at De Casteljau algorithm and thought that dividing the curve into 24*4 pieces for 4 seconds would solve my problem but that sounds erroneous as time is along the "x" axis and not along the curve.

To further simplify If I draw the curve in a plane, the x axis represents the time and the y axis the value I am looking for. What I actually require is to to be able to find out "y" corresponding to "x". Then I can divide x in 24 divisions and know the value for each frame

Answer 1

I was facing the same problem: Every animation package out there seems to use Bézier curves to control values over time, but there is no information out there on how to implement a Bézier curve as a y(x) function. So here is what I came up with.

A standard cubic Bézier curve in 2D space can be defined by the four points P₀=(x₀, y₀) .. P₃=(x₃, y₃).
P₀ and P₃ are the end points of the curve, while P₁ and P₂ are the handles affecting its shape. Using a parameter t ϵ [0, 1], the x and y coordinates for any given point along the curve can then be determined using the equations
A) x = (1-t)³x₀ + 3t(1-t)²x₁ + 3t²(1-t)x₂ + t³x₃ and
B) y = (1-t)³y₀ + 3t(1-t)²y₁ + 3t²(1-t)y₂ + t³y₃.

What we want is a function y(x) that, given an x coordinate, will return the corresponding y coordinate of the curve. For this to work, the curve must move monotonically from left to right, so that it doesn't occupy the same x coordinate more than once on different y positions. The easiest way to ensure this is to restrict the input points so that x₀ < x₃ and x₁, x₂ ϵ [x₀, x₃]. In other words, P₀ must be to the left of P₃ with the two handles between them.

In order to calculate y for a given x, we must first determine t from x. Getting y from t is then a simple matter of applying t to equation B.

I see two ways of determining t for a given y.

First, you might try a binary search for t. Start with a lower bound of 0 and an upper bound of 1 and calculate x for these values for t via equation A. Keep bisecting the interval until you get a reasonably close approximation. While this should work fine, it will neither be particularly fast nor very precise (at least not both at once).

The second approach is to actually solve equation A for t. That's a bit tough to implement because the equation is cubic. On the other hand, calculation becomes really fast and yields precise results.

Equation A can be rewritten as
(-x₀+3x₁-3x₂+x₃)t³ + (3x₀-6x₁+3x₂)t² + (-3x₀+3x₁)t + (x₀-x) = 0.
Inserting your actual values for x₀..x₃, we get a cubic equation of the form at³ + bt² + c*t + d = 0 for which we know there is only one solution within [0, 1]. We can now solve this equation using an algorithm like the one posted in this Stack Overflow answer.

The following is a little C# class demonstrating this approach. It should be simple enough to convert it to a language of your choice.

using System;

public class Point {
    public Point(double x, double y) {
        X = x;
        Y = y;
    }
    public double X { get; private set; }
    public double Y { get; private set; }
}

public class BezierCurve {
    public BezierCurve(Point p0, Point p1, Point p2, Point p3) {
        P0 = p0;
        P1 = p1;
        P2 = p2;
        P3 = p3;
    }

    public Point P0 { get; private set; }
    public Point P1 { get; private set; }
    public Point P2 { get; private set; }
    public Point P3 { get; private set; }

    public double? GetY(double x) {
        // Determine t
        double t;
        if (x == P0.X) {
            // Handle corner cases explicitly to prevent rounding errors
            t = 0;
        } else if (x == P3.X) {
            t = 1;
        } else {
            // Calculate t
            double a = -P0.X + 3 * P1.X - 3 * P2.X + P3.X;
            double b = 3 * P0.X - 6 * P1.X + 3 * P2.X;
            double c = -3 * P0.X + 3 * P1.X;
            double d = P0.X - x;
            double? tTemp = SolveCubic(a, b, c, d);
            if (tTemp == null) return null;
            t = tTemp.Value;
        }

        // Calculate y from t
        return Cubed(1 - t) * P0.Y
            + 3 * t * Squared(1 - t) * P1.Y
            + 3 * Squared(t) * (1 - t) * P2.Y
            + Cubed(t) * P3.Y;
    }

    // Solves the equation ax³+bx²+cx+d = 0 for x ϵ ℝ
    // and returns the first result in [0, 1] or null.
    private static double? SolveCubic(double a, double b, double c, double d) {
        if (a == 0) return SolveQuadratic(b, c, d);
        if (d == 0) return 0;

        b /= a;
        c /= a;
        d /= a;
        double q = (3.0 * c - Squared(b)) / 9.0;
        double r = (-27.0 * d + b * (9.0 * c - 2.0 * Squared(b))) / 54.0;
        double disc = Cubed(q) + Squared(r);
        double term1 = b / 3.0;

        if (disc > 0) {
            double s = r + Math.Sqrt(disc);
            s = (s < 0) ? -CubicRoot(-s) : CubicRoot(s);
            double t = r - Math.Sqrt(disc);
            t = (t < 0) ? -CubicRoot(-t) : CubicRoot(t);

            double result = -term1 + s + t;
            if (result >= 0 && result <= 1) return result;
        } else if (disc == 0) {
            double r13 = (r < 0) ? -CubicRoot(-r) : CubicRoot(r);

            double result = -term1 + 2.0 * r13;
            if (result >= 0 && result <= 1) return result;

            result = -(r13 + term1);
            if (result >= 0 && result <= 1) return result;
        } else {
            q = -q;
            double dum1 = q * q * q;
            dum1 = Math.Acos(r / Math.Sqrt(dum1));
            double r13 = 2.0 * Math.Sqrt(q);

            double result = -term1 + r13 * Math.Cos(dum1 / 3.0);
            if (result >= 0 && result <= 1) return result;

            result = -term1 + r13 * Math.Cos((dum1 + 2.0 * Math.PI) / 3.0);
            if (result >= 0 && result <= 1) return result;

            result = -term1 + r13 * Math.Cos((dum1 + 4.0 * Math.PI) / 3.0);
            if (result >= 0 && result <= 1) return result;
        }

        return null;
    }

    // Solves the equation ax² + bx + c = 0 for x ϵ ℝ
    // and returns the first result in [0, 1] or null.
    private static double? SolveQuadratic(double a, double b, double c) {
        double result = (-b + Math.Sqrt(Squared(b) - 4 * a * c)) / (2 * a);
        if (result >= 0 && result <= 1) return result;

        result = (-b - Math.Sqrt(Squared(b) - 4 * a * c)) / (2 * a);
        if (result >= 0 && result <= 1) return result;

        return null;
    }

    private static double Squared(double f) { return f * f; }

    private static double Cubed(double f) { return f * f * f; }

    private static double CubicRoot(double f) { return Math.Pow(f, 1.0 / 3.0); }
}