consider using Eriksson’s formula in preference to L’Huilier’s for the area of a spherical triangle
- Dominant language
- C++
- Stars
- 2.7k
- Forks
- 357
- Avg merge
- 11h 4m
- Merged PRs (30d)
- 1
Description
[Currently S2 calculates](https://github.com/google/s2geometry/blob/master/src/s2/s2measures.cc#L56-L110) the signed area *E* of an oriented spherical triangle with vertices at unit vectors [*a*, *b*, *c*] by first finding the measures of the central angles between the vectors, then applying an angle-measure based formula (the spherical analog of Heron’s formula in the plane). This involves a bunch of expensive circular functions, and has many corners for rounding errors to sneak into the calculation.
If we expand the current computation it might be written:
> *α* = atan2( √[(*b* × *c*) · (*b* × *c*)], *b* · *c* )
> *β* = atan2( √[(*c* × *a*) · (*c* × *a*)], *c* · *a* )
> *γ* = atan2( √[(*a* × *b*) · (*a* × *b*)], *a* · *b* )
> *σ* = ½(*α* + *β* + *γ*)
>
> *s* = sgn(*a* · (*b* × *c*))
> *E* = 4*s* atan( √[ tan½*σ* tan½(*σ* − *α*) tan½(*σ* − *β*) tan½(*σ* − *γ*) ] )
Instead, consider using a formula from
Euler (1781) [“De mensura angulorum solidorum”](https://scholarlycommons.pacific.edu/euler-works/514/) (non-vector variant)
Van Oosterom & Strackee (1983) [“The solid angle of a plane triangle”](https://doi.org/10.1109/TBME.1983.325207). *IEEE transactions on Biomedical Engineering* 2: 125–126.
Eriksson (1990) [“On the Measure of Solid Angles”](https://www.jstor.org/stable/2691141). *Mathematics Magazine*, 63(3), 184–187
> *E* = 2 atan2(*a* · (*b* × *c*), 1 + *a*·*b* + *b*·*c* + *c*·*a*)
> (or)
> *E* = 2 atan2(*a* · ((*b* − *a*) × (*c* − *a*)), (*b* + *a*) · (*c* + *a*))
Notice that the numerator of the tangent needs to be calculated anyway to get the sign in the current version (it is the standard orientation predicate for a triple of unit vectors). The denominator takes 2 vector-vector additions and 1 dot product.
Overall, we cut the number of circular function evaluations from 8 down to 1, cut the dot products from 7 to 2, cut the cross products from 4 to 1, and cut square roots from 4 to 0.
The second line above avoids loss of significance in the typical case where triangle vertices are close together by computing the differences between vectors up front. If even more worried about possible precision loss for small skinny triangles, we could even replace (*b* − *a*) × (*c* − *a*) with (½(*b* + *c*) − *a*) × (*c* − *b*) or similar, as is done in `RobustCrossProduct`.
* * *
I don’t know that anyone has carefully worked through the error analysis of this formula, or tried to figure out if it needs special touch-ups in extreme situations (e.g. *a* almost antipodes to *b*), but I would expect it to be a significantly easier thing to study than the L’Huilier formula, as it sticks to basic floating point arithmetic.
Contributor guide
Assessment
This issue has not been assessed yet.