I have implemented code in Python that attempts to find the perpendicular points of a higher order bezier curve of order $n$ (where $n$ is around 25 control points) parametrized in standard form on $t \in [0, 1]$.
First, I convert the curve to its polynomial (power basis) representation. Then, I solve the equation $(\vec Q(t) - P) \cdot \vec Q'(t) = 0$ using numpy's np.roots, where $\vec Q(t)$ is the parametric polynomial curve in 2D represented by its 2D coefficients, $P$ is the input point, and $\vec Q'(t)$ is the derivative of the curve.
For low order curves (<10) the solution seems mostly stable. For very low orders (<5) it seems entirely stable. However, for the higher order curves, I am finding that it is only stable near low values of t (around $t < 0.5$ to $t < 0.8$ depending on the curve's shape) but as soon as $t$ approaches $1$, the solution becomes very unstable.
I have determined this is not due to the shape of the curve, since if I simply reverse the ordering of the control points, essentially reparameterizing the curve so that the new parameter $s = 1 -t$, I find that the end of the curve near $t = 1$ (now $s = 0$) has a stable solution. Practically, I can therefore solve for both instances of the curve (flipped and not flipped) and combine the roots with some additional checking.
I would like to know, however why this is happening, and if there are potentially better solutions to solve the instability than reversing the control points and solving both systems? I have noticed that the resultant polynomial (of order $2n - 2$ i.e. around order 48) has coefficients of varying orders of magnitude (from $10^{21}$ to $10^6$). I have tried normalizing these coefficients but it doesn't seem to help.
Any suggestions? Is this a known issue? And why does it occur for t near 1 but not 0? (my intuition says its because the higher terms are made smaller since a number near zero multiplied by itself many times becomes small). Could it have to do with how the original polynomial is built from the Bezier? Help would be much appreciated
Converting from the Bernstein basis (Bézier form) to the power basis is almost always a bad idea, from the point of view of numerical stability.
So, yes, what you’re encountering is a known problem. There are a couple of papers by Farouki and Rajan that discuss the numerical character of the Bernstein (and other) bases. I suggest starting here:
R.T. Farouki, V.T. Rajan.
Algorithms for polynomials in Bernstein form.
Computer Aided Geometric Design.
Volume 5, Issue 1, June 1988, Pages 1-26
Yes, probably. If you look at the math that you’re doing to get the polynomial coefficients from the Bézier control points, I expect you’ll find places where you’re subtracting very large numbers, so getting subtractive cancellation. If you did all your calculations in high precision, I’d expect it to work OK.
Your curve-reversing idea is a good workaround.