I'm working on a software program which currently utilizes cubic Bézier splines for generating continuous well-defined function output (only 1 output value per input) in the time domain with varying sample rates. Initially I utilized some root finding code using Cardano's method, in order to convert time domain (X axis) positions into Bezier curve t position values [0, 1]. However, this is proving to have been a poor choice of algorithms for my purpose due to numeric instability. Also the current solution I am using utilizes trigonometry functions, which can be costly on many computer platforms, especially embedded ones.
I am now looking into other algorithms, but it also got me thinking as to whether cubic Bézier curves are even the optimal solution, when the application is not constrained to any particular spline type.
Ideal spline solution properties:
- Intuitive user interface control of curve segments
- Very stable "well defined" interpolation sampled output
- Software optimization friendly floating point algorithm for CPU and possibly FPGA or DSP embedded solutions
Cubic Bézier curve splines and their 2 non-intersecting control points are familiar in many software programs for animation and illustration. So ideally I would continue to use this type of spline (at least on the front end) and choose the optimal solution algorithm which meets these requirements. One solution I am currently working with utilizes the Newton–Raphson method in combination with a binary search to find a t [0, 1] value given an X position, such as with this project: https://github.com/gre/bezier-easing However, this was designed for rendering curves to screens and so I'm not yet sure if it will meet my numeric stability requirements. It does seem that the algorithm utilizes fairly basic math operations though, but if there is a better spline type to use in this type of situation, without having to go through this iterative approach of finding the t value for a given X, that may be a better option. Any helpful pointers would be greatly appreciated.
I ended up using the solution described at the end of my question, namely the Newton–Raphson method in combination with a binary search to find a t [0, 1] value given an X position. Code for this can be found with this project: https://github.com/gre/bezier-easing
I was also able to optimize it to the point where it is comparable in performance to the original Cardano method I was using, which was much more complex, difficult to work with, and quite numerically unstable. The performance/accuracy trade off the Newton-Raphson method can easily be tuned as well, just by adding more iterations to the loops. One downside is that it is not as conducive to SIMD optimizations, because each iteration depends on results from the previous one.