Modify ↓
Opened 17 months ago
Closed 3 months ago
#11930 closed Bugs (fixed)
huiller method is inaccurate
Reported by: | Charles Karney <charles@…> | Owned by: | barendgehrels |
---|---|---|---|
Milestone: | Boost 1.64.0 | Component: | geometry |
Version: | Boost 1.58.0 | Severity: | Problem |
Keywords: | area | Cc: | charles.karney@… |
Description
The following code illustrates a problem with the "huiller" strategy for computing spherical areas. It computes the area of a sequence of progressively skinny triangles. The relative errors in the computed areas are much larger than they need be.
#include <iostream> #include <cmath> #include <boost/geometry.hpp> namespace bg = boost::geometry; typedef bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> > pt; int main() { bg::model::polygon<pt, false, false> poly; bg::read_wkt("POLYGON((0 45,0 -45,0 0))", poly); double c = 0.014458780939650811; std::cout << " area relerr" << "\n"; for (int i = 3; i <= 6; ++i) { double h = std::pow(10.0, -i); poly.outer()[2].set<0>(h); double area = bg::area(poly), // Relative error in this estimate is 2e-11 for h = 0.001 area0 = h * c; std::cout << area << " " << (area - area0)/(area0) << "\n"; } }
Running this code gives
area relerr 1.44588e-05 -2.0871e-06 1.44594e-06 4.55941e-05 1.43961e-07 -0.00433739 0 -1
Discussion:
- Surely it's a bad idea to embed the algorithm name "huiller" into the API. Computing spherical areas is really straightforward. Users shouldn't need to worry about the underlying method. Maybe "spherical_area" is a better name.
- In any case, you've misspelled the guy's name. It's either "l'Huilier" or "l'Huillier" with two "i"s.
- Furthermore, you're needlessly blackening l'Huilier's name by applying his formula for the area of a spherical triangle in terms of its sides.
- The three-side formula is a terrible choice because if a + b is nearly equal to c (as is the case in my example), the triangle is badly condititioned.
- Much better is to use the formula of the area of a triangle in terms of two sides and the included angle. See the last two formulas in the Wikipedia article on spherical excess. This nicely reduces to the familar trapezium rule in the plane limit.
- This formula comes with a sign attached, so there's no need to kludge a sign on afterwards as there is with the current implementation.
- This formula together with an ellipsoidal correction is used by GeographicLib for computing geodesic areas.
Attachments (0)
Change History (1)
comment:1 Changed 3 months ago by awulkiew
- Milestone changed from To Be Determined to Boost 1.64.0
- Resolution set to fixed
- Status changed from new to closed
Note: See
TracTickets for help on using
tickets.
Thanks!
Fixed by this PR: https://github.com/boostorg/geometry/pull/359
The name is changed to area::spherical.
The new results are: