Cubic Bezier curve: given X,Y, solve for t

[abeall]

I'm programming some commands in a graphic application, and once again my extremely limited mathematical knowledge has found me rather stumped.

Here is my cubic bezier function, written in JavaScript syntax(this is not actually JavaScript, though):

function getBezier(percent,p1,cp1,cp2,p2) {
function b1(t) { return t*t*t }
function b2(t) { return 3*t*t*(1-t) }
function b3(t) { return 3*t*(1-t)*(1-t) }
function b4(t) { return (1-t)*(1-t)*(1-t) }
var pos = {x:0,y:0};
pos.x = p1.x*b1(percent) + cp1.x*b2(percent) + cp2.x*b3(percent) + p2.x*b4(percent);
pos.y = p1.y*b1(percent) + cp1.y*b2(percent) + cp2.y*b3(percent) + p2.y*b4(percent);
return pos;
}

You give it a cubic bezier curve defined by four points (two "anchors" which intersect the curve, and two "control points" which affect the curvature), and a time value(written as "percent" above) from 0 to 1, and it returns a point(x,y) at that time.

I'm now in a situation where I need just the opposite. I have a point which I know is on the bezier curve, and I need to find the time value.

Cubic bezier curves can self intersect, in which case I suppose there's a rare possibility that a point actually has two valid time values, but I'm not worried about that, I don't think it can occur in my environment.

After being completely befuddled by the formula for cubic beziers presented on wikipedia (http://upload.wikimedia.org/math/e/2/9/e29b31672536b36c91cb7be0ce8c2ff7.png) I'm left unable to solve for t. Any help in that direction, if it's even possible, would be appreciated.

Thanks!
- aaron

 

[EnumaElish]

In Mathematica 6.0, I define:

pos[t_] = p1*b1[t] + cp1*b2[t] + cp2*b3[t] + p2*b4[t];

where b1[t_] = t*t*t, etc.,

then solve pos[t] = a for t:

Solve[pos[t] == a, t];

which returns 2 complex roots and one real root.

The real root =
-((cp1 - 2 cp2 + p2)/(-3 cp1 + 3 cp2 + p1 -
p2)) - (2^(
1/3) (-9 cp1^2 + 9 cp1 cp2 - 9 cp2^2 + 9 cp2 p1 + 9 cp1 p2 -
9 p1 p2))/(3 (-3 cp1 + 3 cp2 + p1 - p2) (243 a cp1^2 -
54 cp1^3 - 486 a cp1 cp2 + 81 cp1^2 cp2 + 243 a cp2^2 +
81 cp1 cp2^2 - 54 cp2^3 - 162 a cp1 p1 + 162 a cp2 p1 +
81 cp1 cp2 p1 - 162 cp2^2 p1 + 27 a p1^2 + 162 a cp1 p2 -
162 cp1^2 p2 - 162 a cp2 p2 + 81 cp1 cp2 p2 - 54 a p1 p2 +
81 cp1 p1 p2 + 81 cp2 p1 p2 - 27 p1^2 p2 + 27 a p2^2 -
27 p1 p2^2 + Sqrt(4 (-9 cp1^2 + 9 cp1 cp2 - 9 cp2^2 +
9 cp2 p1 + 9 cp1 p2 - 9 p1 p2)^3 + (243 a cp1^2 -
54 cp1^3 - 486 a cp1 cp2 + 81 cp1^2 cp2 + 243 a cp2^2 +
81 cp1 cp2^2 - 54 cp2^3 - 162 a cp1 p1 + 162 a cp2 p1 +
81 cp1 cp2 p1 - 162 cp2^2 p1 + 27 a p1^2 + 162 a cp1 p2 -
162 cp1^2 p2 - 162 a cp2 p2 + 81 cp1 cp2 p2 -
54 a p1 p2 + 81 cp1 p1 p2 + 81 cp2 p1 p2 - 27 p1^2 p2 +
27 a p2^2 - 27 p1 p2^2)^2))^(1/3)) + (243 a cp1^2 -
54 cp1^3 - 486 a cp1 cp2 + 81 cp1^2 cp2 + 243 a cp2^2 +
81 cp1 cp2^2 - 54 cp2^3 - 162 a cp1 p1 + 162 a cp2 p1 +
81 cp1 cp2 p1 - 162 cp2^2 p1 + 27 a p1^2 + 162 a cp1 p2 -
162 cp1^2 p2 - 162 a cp2 p2 + 81 cp1 cp2 p2 - 54 a p1 p2 +
81 cp1 p1 p2 + 81 cp2 p1 p2 - 27 p1^2 p2 + 27 a p2^2 -
27 p1 p2^2 + Sqrt(4 (-9 cp1^2 + 9 cp1 cp2 - 9 cp2^2 +
9 cp2 p1 + 9 cp1 p2 - 9 p1 p2)^3 + (243 a cp1^2 -
54 cp1^3 - 486 a cp1 cp2 + 81 cp1^2 cp2 + 243 a cp2^2 +
81 cp1 cp2^2 - 54 cp2^3 - 162 a cp1 p1 + 162 a cp2 p1 +
81 cp1 cp2 p1 - 162 cp2^2 p1 + 27 a p1^2 + 162 a cp1 p2 -
162 cp1^2 p2 - 162 a cp2 p2 + 81 cp1 cp2 p2 - 54 a p1 p2 +
81 cp1 p1 p2 + 81 cp2 p1 p2 - 27 p1^2 p2 + 27 a p2^2 -
27 p1 p2^2)^2))^(1/3)/(3 2^(1/3) (-3 cp1 + 3 cp2 + p1 - p2))

I hope this is useful. If I misunderstood the question let me know.

 

[dodo]

... the formula for cubic beziers presented on wikipedia (http://upload.wikimedia.org/math/e/2/9/e29b31672536b36c91cb7be0ce8c2ff7.png ...

This formula should be familiar to you, since it is exactly what you have implemented. The B and P in bold typeface are (fixed) vectors with x and y components, so that you actually have two (non-linear) equations on the single unknown t. The same two lines on your code that calculate pos.x and pos.y.

I would try the following:
A) On each of the two equations, distribute all the products and group similar powers of t together, until you finish with a 3-degree polynomial: an expression of the form $$a t^3 + b t^2 + c t + d$$, where each of a, b, c and d will probably be a big expression in parenthesis but, nevertheless, a constant (given an input point, and anchor and control points, all fixed).
B) Then you should be able to multiply each equation by some constant (actually, multiply each equation by the a coefficient of the other), and subtract one equation from the other, to eliminate the $$t^3$$ term, ending with one polynomial of degree 2, $$a' t^2 + b' t + c'$$, which you probably know how to solve.

 

[abeall]

Thank-you both!

Sorry, I should have mentioned I did not create the original getBezier function. After staring it and comparing to the wiki cubic bezier formula for awhile, I somewhat see the translation, though I don't think I could have come up with it on my own.

Enuma,

What is the "real root" and how will it help me find t? Does the expression evaluate to t? It's an immense expression, I've never seen anything like it!

Dodo,

Thanks for the advice. Unfortunately I don't follow. What do you mean by:
A) "all the products" - you mean all the points along the curve, as calculated by getBezier? "group similar powers of t" - what is a power of t? what would constitute 'similar' and in what way are they grouped? "until you finish" - finish what? The grouping, presumably, but I don't really know what grouping means, what exactly is being done at this step. "with a 3-degree polynomial" - meaning something cubed? what should be cubed?
B) "each equation" - at this point, what is each equation? You mean the two equation in the original getBezier function which calculate x,y respectively? "coefficient of the other" - what are the coefficients, or how do I find them? "which you probably know how to solve" - I'm afraid I don't have a clue. I'm not even sure what a and b represent at this point, much less what the little line above them represents.

 

[mda]

A) On each of the two equations, distribute all the products and group similar powers of t together, until you finish with a 3-degree polynomial: an expression of the form $$a t^3 + b t^2 + c t + d$$, where each of a, b, c and d will probably be a big expression in parenthesis but, nevertheless, a constant (given an input point, and anchor and control points, all fixed).

abeall, what dodo means is to rewrite the bezier equation as:
(3*cp2-p2-3*cp1+p1) * t^3 + (3*p2+3*cp1-6*cp2) * t^2 + (-3*p2+3*cp2) * t + p2 = chosen point

 

[abeall]

abeall, what dodo means is to rewrite the bezier equation as:
(3*cp2-p2-3*cp1+p1) * t^3 + (3*p2+3*cp1-6*cp2) * t^2 + (-3*p2+3*cp2) * t + p2 = chosen point

Wow, I have no idea how you did that, but I can see how that is the expression Dodo described. Thank-you very much.

Now for step B... I'm a little lost because t is an unknown, it is what I'm trying to solve, and x and y are known. So I would have thought the equation would need to written somehow as t = ... so that I can find t. Is that not where this is going?

 

[Ben Niehoff]

For a general cubic,

$$at^3 + bt^2 + ct + d = 0$$

the solution is

$$\begin{array}{rcl}t &=& -\frac{b}{3a} + \frac{1}{3a} \left[ K \frac{ \sqrt[3]{36abc -108a^2d-8b^3+12\sqrt{12ac^3-3b^2c^2-54 bcd+81a^2d^2+12ab^3d}}}{2}\right \\ && \left \quad - \bar{K} (3ac-b^2)\frac{2}{\sqrt[3]{36abc-108a^2d-8b^3+12\sqrt{12ac^3 3b^2c^2-54abcd+81a^2d^2+12ab^3d}}} \right] \end{array}$$

where K is given by anyone of the following (where $i = \sqrt{-1}$):

$$K \in \left\{1, \quad -\frac12+i\frac{\sqrt3}2, \quad -\frac12-i\frac{\sqrt3}2 \right\}$$

and $\bar{K}$ is the complex conjugate of K.

As you can see, this is a bit of a pain to try to implement. Usually solutions are found by approximate means instead.

 

[abeall]

As you can see, this is a bit of a pain to try to implement. Usually solutions are found by approximate means instead.

Wow. Yeah. An approximation would be good. Are there any well known approximation methods I could research? The only thing I can think of is to walk the bezier at a small time interval, say 0.01, then sort the points by distance to target point, and refine the closest or closest two until I'm within desired range of error to the target point. Unfortunately, I'm a little hesitant to walk the bezier curve because this command was going to be essentially a faster alternative to another method I'm using that basically does what I just described.

But here's another problem I need to solve in a related situation of this same subject: is it possible to find an X value given a Y value and no time value? Or must should an approximation again be used instead?

 

[hotvette]

As stated in previous posts, it really boils down to solving a cubic equation. If you really think there aren't multiple t's for the x you are looking for, then Newton Raphson should work. I've used it with cubic beziers for the exact same purpose you state and it works well. Assume you are using the Bezier formulation where t=0,1. For the starting point, one possibility is to calculate x for t=0 and t=1, then by simple ratio estimate the t for the desired x, then iterate from there. There are more elaborate methods. Of course, this may or may not be faster than the method you currently use.

The method of walking through the curve can work, but you don't need such a fine increment. Just using a few points (e.g. 4 or 5) should bracket the interval such that you can converge on the solution with very few Newton Raphson iterations (3 or 4 perhaps).

Even better would be to use Halley's method:

http://mathworld.wolfram.com/HalleysMethod.html

 

[EnumaElish]

What is the "real root" and how will it help me find t? Does the expression evaluate to t? It's an immense expression, I've never seen anything like it!

A cubic equation has 3 roots (t's that solve it). In general the roots may be real (a number that does not include the imaginary Sqrt[-1]) or complex (a number that includes the imaginary Sqrt[-1]).

In your case, two roots (t's) are complex, one is real.

Mathematica is pretty good with this type of a problem.

 

[abeall]

Thanks for the explanation, Enuma. That made sense. And thank you hotvette, I am now working on a Newton Raphson approximation.

 

[math_cripple]

For a general cubic,

$$at^3 + bt^2 + ct + d = 0$$

the solution is

$$\begin{array}{rcl}t &=& -\frac{b}{3a} + \frac{1}{3a} \left[ K \frac{ \sqrt[3]{36abc -108a^2d-8b^3+12\sqrt{12ac^3-3b^2c^2-54 bcd+81a^2d^2+12ab^3d}}}{2}\right \\ && \left \quad - \bar{K} (3ac-b^2)\frac{2}{\sqrt[3]{36abc-108a^2d-8b^3+12\sqrt{12ac^3 3b^2c^2-54abcd+81a^2d^2+12ab^3d}}} \right] \end{array}$$

where K is given by anyone of the following (where $i = \sqrt{-1}$):

$$K \in \left\{1, \quad -\frac12+i\frac{\sqrt3}2, \quad -\frac12-i\frac{\sqrt3}2 \right\}$$

and $\bar{K}$ is the complex conjugate of K.

As you can see, this is a bit of a pain to try to implement. Usually solutions are found by approximate means instead.

Can somebody verify this equation? I notice that the terms in the square root are very similar but not identical. I wonder if this is correct or a mistake.

The first term is:
12ac^3-3b^2c^2-54 bcd+81a^2d^2+12ab^3d
The second is:
12ac^3 3b^2c^2-54abcd+81a^2d^2+12ab^3d

Notice the second lacks a '-' in front of the 3b (making it a multiplier) and there is an 'a' in 54abcd where the a is missing from the first.


Also, I'm currently using Newton-Raphson. I have cases that seem to be sensitive to the 'initial guess' and get very bad values if the initial guess is wrong. Is Halley's method less sensitive to this problem?

Thank you