I've been working on finding the inflection points of a cubic Bezier curve using the method described in a paper Hain, Venkat, Racherla, and Langan - Fast, Precise Flattening of Cubic Bézier Segment Offset Curves I recently read. The curve is defined by four control points, and the method involves calculating certain coefficients from these control points and then using them to determine the inflection points.
Given a Bézier curve defined by control points: $$ P_0 = (x_0, y_0), \quad P_1 = (x_1, y_1), \quad P_2 = (x_2, y_2), \quad P_3 = (x_3, y_3). $$ Using the Bézier basis matrix, the coefficients in terms of the control points are: $$ \begin{align*} a_x &= -x_0 + 3x_1 - 3x_2 + x_3 \\ b_x &= 3x_0 - 6x_1 + 3x_2 \\ c_x &= -3x_0 + 3x_1 \\ d_x &= x_0 \\ \end{align*} $$ \begin{align*} a_y &= -y_0 + 3y_1 - 3y_2 + y_3 \\ b_y &= 3y_0 - 6y_1 + 3y_2 \\ c_y &= -3y_0 + 3y_1 \\ d_y &= y_0. \\ \end{align*} To find the inflection points, we consider the first and second derivatives of the curve with respect to the parameter $(t)$. At inflection points, the component of the acceleration (second derivative) perpendicular to the velocity (first derivative) is zero. This leads to the equation:
$$ \frac{dx}{dt} \cdot \frac{d^2y}{dt^2} - \frac{d^2x}{dt^2} \cdot \frac{dy}{dt} = (3a_xt^2 + 2b_xt + c_x)(6a_yt + 2b_y) - (6a_xt + 2b_x)(3a_yt^2 + 2b_yt + c_y) = 0. $$ Solving this question for $t$ yields: $$ t_\text{cusp} = -\frac{1}{2} \times \frac{aycx - axcy}{aybx - axby} $$ and $$ t_1 = t_\text{cusp} - \sqrt{t_\text{cusp}^2 \times \frac{1}{3} \left( \frac{b_yc_x - b_xc_y}{a_yb_x - a_xb_y} \right)} $$ $$ t_2 = t_\text{cusp} + \sqrt{t_\text{cusp}^2 \times \frac{1}{3} \left( \frac{b_yc_x - b_xc_y}{a_yb_x - a_xb_y} \right)} $$ where $t_1$ and $t_2$ are possible inflection points.
However, I've encountered an issue with certain sets of control points where the denominator $a_yb_x - a_xb_y$ resolves to $0$, causing a division by zero error. For instance, consider the control points: $$ P_0 = (0, 0), \quad P_1 = (1, 1), \quad P_2 = (2, -1), \quad P_3 = (3, 0). $$ Plugging these values in we get: \begin{align*} a_x &= -0 + 3\cdot0 - 3\cdot1 + 0 \\ a_x &= -2\\[1em] b_x &= 3\cdot0 - 6\cdot0 + 3\cdot1 \\ b_x &= 3\\[1em] a_y &= -2 + 3\cdot-1 - 3\cdot3 + 0 \\ a_y &= -14 \\[1em] b_y &= 3\cdot2 - 6\cdot-1 + 3\cdot3 \\ b_y &= 21. \end{align*} Using these values, we can see that $a_yb_x-a_xb_y$ resolves to $0$: \begin{align*} a_yb_x-a_xb_y &= -14\cdot3-(-2\cdot21) = -42-(-42) = 0. \end{align*}
For this set, the method fails due to the division by zero. Visually inspecting the curve I expect there to be a point roughly or perhaps exactly at the halfway point.
However, when I slightly perturb one of the control points, say turning the 3 to 2.9999 in $p_3$: $$ P_0 = (0, 0), \quad P_1 = (1, 1), \quad P_2 = (2, -1), \quad P_3 = (2.9999, 0) $$ the method calculates an inflection point at approximately $t = 0.5000041667371988$.
I've attached two images for clarity. The first image shows the curve defined by the control points in blue. The second image shows the same curve in blue, but with the curvature plotted in red, clearly indicating a zero curvature (inflection point) halfway along the curve.
I suspect that if the inflection point exists in the centre of the curve, this case arises though I have no idea how to prove that.
Can anyone provide insights into why the method fails for the original set of control points and how to handle such cases? Is there something fundamental that I am missing, the paper does not seem to mention this case?
