Getting y from x co-ord for cubic bezier curve, fast Newton-Raphson method

Question

Given the points of a Bezier curve (P0, P1, P2, P3) in 2D, I would like to find the y co-ordinate for a given x co-ordinate. The problem is well defined because of the following restrictions:

• P0 = (0,0), P3 = (1,1)
• P1 = (t, 1-t) for t between 0, 1
• P2 = 1 - P1 (x and y)

I have the following function to calculate the answer, having put in all the restrictions above into the Bezier curve formula here CubicBezier.html. I am using Newton-Raphson to work out the parameter of the point I want, and I know this works as I don't let the loop finish until it has (within defined tolerance).

I am using this function to apply a contrast filter to an image. For this 0.5 gives the same image back, 0.0 does max reduction in contrast and 1.0 does max increase.

EDIT The following function has been corrected and now works perfectly.

/*
 * Parameters: p - x co-ord of P1, 
 *             x - x we want to find y for
 *
 * This method is unstable for p ~= 0.5, maybe needs refinement.
 */

#include <iostream>
#include <math.h>

#define ITER_TOL  0.00001

float maths::bezier(float p, float x)
{
    if(p < 0.f || p > 1.f || x < 0.f || x > 1.f)
    {
        std::cerr << "Both parameters must be between 0 and 1, returning dummy value" << std::endl;
        return 0.f;
    }
    //First guess for u
    float u = x;
    //Coefficients of curve (for x and y co-ord)
    float x3 = 6 * p - 2;
    float x2 = 3 - 9 * p;
    float x1 = 3 * p;
    float x0 = -x;

    float y3 = 6 * (1-p) - 2;
    float y2 = 3 - 9 * (1-p);
    float y1 = 3 * (1-p);

    //Newton-Raphson refinement
    for(int i=0; fabs(x3*u*u*u + x2*u*u + x1*u + x0) > ITER_TOL && i<1000; i++)
    {
        u = u - (x3*u*u*u + x2*u*u + x1*u + x0) /
                (3*x3*u*u + 2*x2*u + x1);
        //std::cout << i << ": " << u << std::endl;
        //Deal with non-convergence
        if(i==999)
        {
            std::cerr << "Warning, Newton-Raphson method did not converge in Maths.cpp, returning dummy" << std::endl;
            return 0.f;
        }
    }
    //Calculate y co-ord
    return y3*u*u*u + y2*u*u + y1*u;
}

If we set p = 0.5, we should get a straight line, but when I do this for a linspace and plot the points, I get a bend between 0.5 and 1.0. Can anyone see why this is happening, and if there is anything I can do about it?

Answer 1

I compiled your code and noticed that the loop runs only 0 or 1 iterations. Probably it's because a fabs is missing somewhere in the convergence check?