▲ 164 points
back
18 comments
Why don’t you try the analytic method for cubics? The Abel-Ruffini only says that the anti-derivative is not elementary, not that it doesn’t exist. The eliptic integrals are non-elementary, but calculating them on a computer is not really different in principle than e.g. computing natural logarithm.
If you plug the desired integral into Mathematica or Wolfram Alpha, it should just spit out the relevant anti derivative expressed in terms of some standard elliptic function, which can then be calculated using some series expansion (or some library of math functions).
You're really trying to nerd snipe, aren't you? I spent some quality time with elliptic functions as the solutions to the elastica, and also with the Fresnel integral which of course generates the Euler spiral. So yes, it wouldn't be surprising if there were a closed form solution, and that it might even be more efficient to compute.
So, this would be a good topic for someone else to show off their math skills :)
Haha, of course it wouldn’t really be much of an improvement over what you already have, and clearly I’m too lazy to do it myself — I would just enjoy looking at it :)
For those unfamiliar with the term 'nerd sniping' (as was I): https://xkcd.com/356/
Given that the desired integral will be a sixth order polynomial, it's not even that "we can't find the intergral" (we can), it's that there is no closed form solution once we've found it.
Integrate[Sqrt[a + b u + c u^2 + d u^3], u]
No, that’s precisely my point: this expression can be evaluated faster and more accurately than by applying generic numeric integration method. Computing this EllipticF function is not different than computing natural logarithm, and if you don’t mind latter, there is not much reason to shy from former.
To be specific
r1,r2,r3 are roots of a+bu+cu^2+u^3,
Cleaning up the output a bit, you get
(2 (b c EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)]-9 a d EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)]+(c+3 d u) (a+u (b+u (c+d u))) Sqrt[((-u+r2) (u-r3))/(r2-r3)^2] (r2-r3)+2 c^2 (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)] (EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] r1+EllipticE[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (-r1+r3))-6 b d (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)] (EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] r1+EllipticE[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (-r1+r3))))/(15 d Sqrt[a+u (b+u (c+d u))] Sqrt[-(((u-r2) (u-r3))/(r2-r3)^2)] (r2-r3)).
But I guess if EllipticF is fast to compute, finding the roots can also be done pretty quickly, so yes its probably not that bad.
You need to integrate the square root of a quartic polynomial – the derivative of a cubic is a quadratic, then you square the norm of that, which gives a quartic. I spent a little time in Mathematica and it there is a solution, but it seems very unlikely to me it's going to be more efficient than numerical integration - it involves finding the roots of a quartic and has nine elliptic function evaluations. From a slide I found[0] it looks like it takes on the order of 100ns to compute an elliptic function, so at best it's still going to be slower than the adaptive quadrature even with an accuracy of 1e-12 or so. Plus there are a lot of divisions by differences of roots, and as each one of those approaches zero there will be numerical stability issues.
Long story short, it looks doable, but is a lot of work to get right and has almost no chance of beating out numerical integration in performance.
If you're concerned about numerical inaccuracy of an expression, it might be worthwhile trying Herbie (https://herbie.uwplse.org/), a system that will try to rearrange expressions to minimize numerical error even in corner cases. It usually makes expressions bigger and slower, but since it works by brute force sometimes it even finds simplifications.
> Or, to put it another way, it’s [sampling] astonishingly slow when high accuracy is desired
Raph, did you try adaptive sampling? How does the performance compare if you sample with an upper bound on change in angle, versus regular sampling?
Also curious if you wouldn’t mind expanding on use cases that need this much accuracy. The main things I’ve had to use Bezier arc length for are graphics related, e.g. hair rendering, and in that domain you don’t need much accuracy at all.
These quadrature techniques are strictly more powerful than adaptive sampling. Part of the motivation for this blog post was to raise awareness of quadrature, as it's an excellent general method when you want to numerically approximate an integral.
One reason to compute with excess accuracy is so that you can compute partial derivatives of arclength with respect to some parameter in the inner loop of an optimizer, to drive Newton or related solving. If the estimate has significant errors, then it can throw off the outer optimization loops. Of course, this only makes sense if it's cheap to compute the inner terms; if you have to pay dearly for accuracy, then use a more robust optimizer.
I've got nice analytical solution for the quadratic bezier:
https://astiopin.github.io/2019/01/04/qbez-parabola.html
Cool. Does it include your newly designed spline? https://raphlinus.github.io/curves/2018/12/21/new-spline.htm...
No, but I intend to use it as the foundational layer that I build the production implementation on top of.
what's with all the bézier stuff, lately? :D
The simple answer is that these days I go where the muse takes me. I'm basically on a year-long sabbatical after leaving Google this summer, where I'm not particularly worrying about making money. I did my PhD on tools for interactive curve design, and I was frustrated while I was at Google because IP issues made working on curves complicated. So when I had some new ideas, it was like the dam bursting, and I've been enjoying lots of motivation and creativity. I started the codebase for kurbo on Christmas and wrote 1300 lines of code that day. I'll keep doing it while it's fun.
Clearly, people like reading about it (it's also really easy to miss that there's a surprisingly high number of Bezier-related posts getting on the front page year round. They usually get there and then fall off again in 8-12 hours =)