4
votes

I have four points that make a cubic bezier curve:

P1 = (10, 5)
P2 = (9, 12)
P3 = (24, -2)
P4 = (25, 3)

Now I want to find the inflection point of this curve. I googled and/but everyone is referring to one site: http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html

Unfortunately the applets are not working and I simply fail to put it all together. Could someone please be so kind to show me how to calculate the inflection point of my curve?

enter image description here

2

2 Answers

2
votes

To elaborate on MBo's answer, C(t) = X' * Y'' - X'' * Y' = 0 is pretty much the pseudocode you want already, because the first and second derivatives are really easy to calculate.

Following this explanation of Bezier derivatives, and a generic Bezier function with coordinates (x1,y1)...(x4,y4) we get:

fx(t) = x1 (1-t)³ + 3·x2·(1-t)²·t + 3·x3·(1-t)·t² + x4·t³
fx'(t) = a·(1-t)² + 2·b·(1-t)·t + c·t²

where a = 3(x2-x1), b = 3(x3-x2), and c = 3(x4-x3), and:

fx''(t) = u·(1-t) + v·t

where u = 2(b-a) and v = 2(c-b). And of course the same goes for the y component:

fy(t) = y1 (1-t)³ + 3·y2·(1-t)²·t + 3·y3·(1-t)·t² + y4·t³
fy'(t) = a'·(1-t)² + 2·b'·(1-t)·t + c'·t²
fy''(t) = u'·(1-t) + v'·t

where a' is the same as a but with y values, etc.

Working out the maths for C(t) = fx'(t)·fy''(t) - fx''(t)·fy'(t) is annoying, but that's why we own computers. If you own a raspberry pi, you own a license for Mathematica, so let's put that to use:

enter image description here

That's a massive formula, but finding the inflections for an "arbitrary" curve is a bit silly, because Bezier curves are invariant to linear affine transforms, so the the t value of the inflection point stays the same whether we check "the real curve" or whether we rotate/translate/scale the curve so that it has more convenient coordinates. Like translating it so that (x1,y1) ends up being (0,0), and (x4,y4) lies on the x-axis so that y4 is zero.

If we do that, then we get a vastly simpler formula:

enter image description here

How much simpler is this? well:

18 times:
  -     x3 * y2
  + 3 * x3 * y2 * t
  - 3 * x3 * y2 * t^2
  -     x4 * y2 * t
  + 2 * x4 * y2 * t^2
  +     x2 * y3
  - 3 * x2 * y3 * t
  + 3 * x2 * y3 * t^2
  -     x4 * y3 * t^2

Which, since we're programming, is a ton of cacheable values. Taking:

a = x3 * y2
b = x4 * y2
c = x2 * y3
d = x4 * y3

we can simplify C(t) as:

1/18 * C(t) = -a + 3at - 3at^2 - bt + 2bt^2 + c - 3ct + 3ct^2 - dt^2
= -3at^2 + 2bt^2 + 3ct^2 - dt^2 + 3at - bt - 3ct - a + c
= (-3a + 2b + 3c - d)t^2 + (3a - b - 3c)t + (c - a)

Putting that factor 18 back in, this is just a simple quadratic formula, for which we can find the roots using the quadratic root identity with more simple values:

v1 = (-3a + 2b + 3c - d) * 18
v2 = (3a - b - 3c) * 18
v3 = (c - a) * 18

And, provided 3a +d is not equal to 2b+3c (because there are no roots if that is the case), we get:

  sqr = sqrt(v2^2 - 4 * v1 * v3)
    d = 2 * v1
root1 =  (sqr - v2) / d
root2 = -(sqr + v2) / d

Throw away roots that don't fall in the Bezier interval [0,1], and what you're left with are the t values for which your original curve inflects.

It's just that I am not able to put things together in a way so that I end up having some nice pseudo-code. That is also why I gave some points and coordinates. I'd like to see someone calculating this with real numbers

Especially on Stackoverflow, it's better not to be lazy up front and instead commit to probably having to learn a new thing.

4
votes

For Bezier curve you have two parametric equations X(t) and Y(t). To determine inflection point of parametric curve, you need to find where curve curvature (wiki) changes the sign. So you need to find the first and the second derivatives of above functions and solve equation:

C(t) = X' * Y'' - X'' * Y' = 0

The first derivative is quadratic, the second is linear, so equation is cubic for t and may have up to 3 solutions.

Edit: Have read linked article, equation is simplified to quadratic one and may have up to 2 solutions.

If solution(s) exists in t range 0..1, you have also to check whether it is real inflection point - check that C'(t) <> 0 at this t value.

Example: blue circles are inflection points (two catched)

enter image description here

Pieces of real code (Delphi)
Given:
P is array of control points
Cf is coefficients of Bezier in power basis

P, Cf: array[0..3] of TPoint

//calculate Bezier coefficients
Cf[3].x := p[3].x - 3 * p[2].x + 3 * p[1].x - p[0].x;
Cf[2].x := 3 * (p[0].x - 2 * p[1].x + p[2].x);
Cf[1].x := 3 * (p[1].x - p[0].x);
Cf[0].x := p[0].x;
//the same for Y

//find parameters of quadratic equation
// a*t^2 + b*t + c = 0
a := 3 * (cf[2].X *cf[3].Y - cf[2].Y *cf[3].X);
b := 3 * (cf[1].X *cf[3].Y - cf[1].Y *cf[3].X);
c := cf[1].X *cf[2].Y - cf[1].Y *cf[2].X;

//here solve quadratic equations, find t parameters
//don't forget a lot of special cases like a=0, D<0, D=0, t outside 0..1 range
Discriminant := b * b - 4 * a * c;
....