Cubic bezier curves - get Y for given X - special case where X of control points is increasing
Asked Answered
M

2

9

I've read a few discussions regarding finding Y at X for a cubic Bezier curve, and have also read this article regarding the matter.

My case is more constrained than the general one, and I wonder if there's a better solution than the general ones mentioned in the above discussions.

My case:

  • The X value of the different control points is increasing. Ie: X3 > X2 > X1 > X0.
  • Also, as a result of above, X(t) is strictly monotonically increasing as well.

Is there any efficient algorithm that takes such constraints into account?

Miramirabeau answered 16/8, 2018 at 14:47 Comment(0)
L
17

First off: this answer only works because your control point constraint means we're always dealing with a parametric equivalent of a normal function. Which is obviously what you want in this case, but anyone who finds this answer in the future should be aware that this answer is based on the assumption that there is only one y value for any given x value...

This is absolutely not true for Bezier curves in general

With that said, we know that—even though we've expressed this curve as a parametric curve in two dimensions—we're dealing with a curve that for all intents and purposes must have some unknown function of the form y = f(x). We also know that as long as we know the "t" value that uniquely belongs to a specific x (which is only the case because of your strictly monotonically increasing coefficients property), we can compute y as y = By(t), so the question is: can we compute the t value that we need to plug into By(t), given some known x value?

To which the answer is: yes, we can.

First, any x value we start off has its own x = Bx(t), so given that we know x, we should be able to solve that function and find the corresponding value(s) of t that leads to that x.

let's look at the function for x(t):

x(t) = a(1-t)³ + 3b(1-t)²t + 3c(1-t)t² + dt³

We can rewrite this to a plain polynomial form as:

x(t) = (-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + a

This is a standard cubic polynomial, using only known constants as coefficients, and we can trivially rewrite this to:

x = (-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + a
0 = (-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + (a-x)

So now we just... solve this equation: we know every value except for t, we just need some mathematical insight to tell us how to do this.

...Of course "just" is not the right qualifier here, there is nothing "just" about finding the roots of a cubic function, but thankfully, Gerolano Cardano laid the ground works to determining the roots back in the 16th century, using complex numbers. Before anyone had even invented complex numbers. Quite a feat! But this is a programming answer, not a history lesson, so let's get implementing:

Given some known value for x, and a knowledge of our coordinates a, b, c, and d, we can implement our root-finding as follows:

// Find the roots for a cubic polynomial with bernstein coefficients
// {pa, pb, pc, pd}. The function will first convert those to the
// standard polynomial coefficients, and then run through Cardano's
// formula for finding the roots of a depressed cubic curve.
double[] findRoots(double x, double[] coordinates) {
  double
    pa = coordinates[0],
    pb = coordinates[1],
    pc = coordinates[2],
    pd = coordinates[3],
    pa3 = 3 * pa,
    pb3 = 3 * pb,
    pc3 = 3 * pc,
    a = -pa  +   pb3 - pc3 + pd,     
    b =  pa3 - 2*pb3 + pc3, 
    c = -pa3 +   pb3, 
    d =  pa  -     x;

  // Fun fact: any Bezier curve may (accidentally or on purpose)
  // perfectly model any lower order curve, so we want to test 
  // for that: lower order curves are much easier to root-find.
  if (approximately(a, 0)) {
    // this is not a cubic curve.
    if (approximately(b, 0)) {
      // in fact, this is not a quadratic curve either.
      if (approximately(c, 0)) {
        // in fact in fact, there are no solutions.
        return new double[]{};
      }
      // linear solution:
      return new double[]{-d / c};
    }
    // quadratic solution:
    double
      q = sqrt(c * c - 4 * b * d), 
      b2 = 2 * b;
    return new double[]{
      (q - c) / b2, 
      (-c - q) / b2
    };
  }

  // At this point, we know we need a cubic solution,
  // and the above a/b/c/d values were technically
  // a pre-optimized set because a might be zero and
  // that would cause the following divisions to error.

  b /= a;
  c /= a;
  d /= a;

  double
    b3 = b / 3,
    p = (3 * c - b*b) / 3, 
    p3 = p / 3, 
    q = (2 * b*b*b - 9 * b * c + 27 * d) / 27, 
    q2 = q / 2, 
    discriminant = q2*q2 + p3*p3*p3, 
    u1, v1;

  // case 1: three real roots, but finding them involves complex
  // maths. Since we don't have a complex data type, we use trig
  // instead, because complex numbers have nice geometric properties.
  if (discriminant < 0) {
    double
      mp3 = -p/3,
      r = sqrt(mp3*mp3*mp3), 
      t = -q / (2 * r), 
      cosphi = t < -1 ? -1 : t > 1 ? 1 : t, 
      phi = acos(cosphi), 
      crtr = crt(r), 
      t1 = 2 * crtr;
    return new double[]{
      t1 * cos(phi / 3) - b3,
      t1 * cos((phi + TAU) / 3) - b3,
      t1 * cos((phi + 2 * TAU) / 3) - b3
    };
  }

  // case 2: three real roots, but two form a "double root",
  // and so will have the same resultant value. We only need
  // to return two values in this case.
  else if (discriminant == 0) {
    u1 = q2 < 0 ? crt(-q2) : -crt(q2);
    return new double[]{
      2 * u1 - b3,
      -u1 - b3
    };
  }

  // case 3: one real root, 2 complex roots. We don't care about
  // complex results so we just ignore those and directly compute
  // that single real root.
  else {
    double sd = sqrt(discriminant);
    u1 = crt(-q2 + sd);
    v1 = crt(q2 + sd);
    return new double[]{u1 - v1 - b3};
  }
}

Okay, that's quite the slab of code, with quite a few additionals:

  • crt() is the cuberoot function. We actually don't care about complex numbers in this case so the easier way to implement this is with a def, or macro, or ternary, or whatever shorthand your language of choice offers: crt(x) = x < 0 ? -pow(-x, 1f/3f) : pow(x, 1f/3f);.
  • tau is just 2π. It's useful to have around when you're doing geometry programming.
  • approximately is a function that compares a value to a very small interval around the target because IEEE floating point numerals are jerks. Basically we're talking about approximately(a,b) = return abs(a-b) < 0.000001 or something.

The rest should be fairly self-explanatory, if a little java-esque (I'm using Processing for these kind of things).

With this implementation, we can write our implementation to find y, given x. Now, the maths may yield more than one root, but our very specific initial conditions means that we know that there is only one root in the [0,1] interval that is used for Bezier curves, so we can simply check which of the three values fits that criterium, and then compute y(t) with that root.

double xCoordinates = [...];
double yCoordinates = [...];

...

double x = some value we know!
double[] roots = findRoots(x, xCoordinates);

double t;
if (roots.length > 0) {
  for (double r: roots) {
    if (r < 0 || r > 1) continue;
    t = r;
    break; }}

double y = compute(t, yCoordinates);

And that's it, we're done: we now have the "t" value that we can use to get the associated "y" value.

Life answered 16/8, 2018 at 18:32 Comment(11)
Thanks! I'll work out the code, and see if it returns the expected values.Miramirabeau
Hmm.. I'm getting negative values for t. Here's the code I used (I named the variables according to the graphics. It's C#, Vector has 2 double components). Do you spot any error?Miramirabeau
aaaaallright so I wrote out a working solution (in Processing) so let's edit this answer. I needed to do this work anyway for a new section on pomax.github.io/bezierinfo so at least we both win =)Valuator
I didn't realize you wrote that primer :) Did you add the section already to the online site? what number is it? I'll translate the above code to C# and see what I get.Miramirabeau
Hmm.. I'm getting slightly worse performance than the binary search, probably due to the abundant use of trigonometric & sqrt functions. Couldn't fully validate the code above; I get mixed results, but didn't dig to see is the error are on my end or the above. I'll do with the binary search for now. I'll mark this as an answer, though I'm not able to verify it.Miramirabeau
I'm starting the writeup for it today - whether it'll end up done remains to be seen. The trig shouldn't be a real problem, C# should have the same numerical precision as Java. I tested it by sweeping a curve, if you have same coordinates for curves where things seem to go wrong, I can see what that does in my tests.Valuator
(Also yes, if your object is speed: always, always, always build a LUT and binary search your way through that. It'll be imprecise, but for showing graphics we don't care about precision, we only care about discrete precision: if we find the right pixel we're done, no need to waste any more more cpu cycles. If you need a result for "doing more work with", however, use this)Valuator
I noticed some undefined behaviour too, so I fixed the answer: the discriminant should be a true <0, ==0, >0 test, not using approximately(). I also added a few value caching vars, and this behaves both correctly for me now, and pretty performant, too.Valuator
@Mike'Pomax'Kamermans Btw in my part of world the math notation for intervals is <>,(),(>,<) so I edited it back but I left the rest of your edit untouched.Serge
@Serge but remember that SO is explicitly centered around English as main language of conduct, so for answers with maths that includes using the English conventions for symbols and notation. While using <> may be your local convention, for an answer on SO it's better to stick with the English convention, namely [a,b] for closed intervals, and (a,b) for open intervals, using (a,b] and [a,b) for half-open intervals.Valuator
If you're wondering how to use this on control points, like a css cubic-bezier(0.4, 0.0, 0.2, 1): That curve is defined by (0, 0), (0.4, 0), (0.2, 1), (1, 1). So in order to use findRoots, call it like findRoots(x, 0f, controlA.x, controlB.x, 1f) Run the resulting roots through that roots.length check in the answer. I had to modify it to remove/add the float tolerance of 0.00001. This gets you the t corresponding to your x. Run that through the polynomial bezier equation to get y: y = controlA.y * 3 * (1 - t) * (1 - t) * t + controlB.y * 3 * (1 - t) * t * t + t * t * tNibble
S
2

If a binary search is too complex, there is still an O(1) approach but its fairly limited. I assume you are using a 4 control point (p0(x0,y0),p1(x1,y1),p2(x2,y2),p3(x3,y3)) cubic Bezier parametrized by some t in the interval [0.0 , 1.0] so:

t = 0.0 -> x(t) = x0, y(t) = y0;
t = 1.0 -> x(t) = x3, y(t) = y3;

First lets forget about Beziers for a while and use a catmull-rom curve instead, which is just an alternative way to represent the same curve. To convert between the 2 cubics use these:

// BEzier to Catmull-Rom
const double m=6.0;
X0 = x3+(x0-x1)*m; Y0 = y3+(y0-y1)*m;
X1 = x0;           Y1 = y0;
X2 = x3;           Y2 = y3;
X3 = x0+(x3-x2)*m; Y3 = y0+(y3-y2)*m;

// Catmull-Rom to Bezier
const double m=1.0/6.0;
x0 = X1;           y0 = Y1;
x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
x3 = X2;           y3 = Y2;

where (xi,yi) are Bezier control points and (Xi,Yi) are Catmull-Rom points. Now if the X distance between all control points have the same distance:

(X3-X2) == (X2-X1) == (X1-X0)

then the X coordinate is linear with t. That means we can compute t directly from X:

t = (X-X1)/(X2-X1);

Now we can compute the Y for any X directly. So if you can chose the control points then chose them so they meet the X distance condition.

If the condition is not met you can try to change control points so it is satisfied (by binary search, by subdividing the cubics into more patches etc...) but beware changing control points can change the shape of the resulting curve if not careful.

Serge answered 16/8, 2018 at 16:33 Comment(1)
Thanks! I can't choose the X values at will, as they're dependent on other factors. But Mike's answer seems to offer a path.Miramirabeau

© 2022 - 2024 — McMap. All rights reserved.