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 with can be said to come from x = Bx(t), so given that we know x, we should be able to 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)3 + 3b(1-t)2t + 3c(1-t)t2 + dt3
We can rewrite this to a plain polynomial form as:
x(t) = (-a + 3b- 3c + d)t3 + (3a - 6b + 3c)t2 + (-3a + 3b)t + a
This is a standard cubic polynomial, with only known constants as coefficients, and we can trivially rewrite this to:
(-a + 3b- 3c + d)t3 + (3a - 6b + 3c)t2 + (-3a + 3b)t + (a-x) = 0
You might be wondering "where did all the -x go for the other values a, b, c, and d?" and the answer there is that they all cancel out, so the only one we actually end up needing to subtract is the one all the way at the end. Handy!
So now we just... solve this equation: we know everything except 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 pa, double pb, double pc, double pd) {
double
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. Which is a little more involved than just calling the function because cubic roots are complicated things. You can get up to three roots back. But only one of them will be applicable to our "t interval" [0,1]:
double x = some value we know!
double[] roots = getTforX(x);
double t;
if (roots.length > 0) {
for (double _t: roots) {
if (_t<0 || _t>1) continue;
t = _t;
break;
}
}
And that's it, we're done: we now have the "t" value that we can use to get the associated "y" value.
