Modify

Opened 22 months ago

Closed 8 months ago

#11930 closed Bugs (fixed)

huiller method is inaccurate

Reported by: Charles Karney <charles@…> Owned by: Barend Gehrels
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 8 months ago by awulkiew

Milestone: To Be DeterminedBoost 1.64.0
Resolution: fixed
Status: newclosed

Thanks!

Fixed by this PR: https://github.com/boostorg/geometry/pull/359

The name is changed to area::spherical.

The new results are:

  area        relerr
1.44588e-05 2.10295e-11
1.44588e-06 2.10165e-13
1.44588e-07 2.01378e-15
1.44588e-08 0

Modify Ticket

Change Properties
Set your email in Preferences
Action
as closed The owner will remain Barend Gehrels.
The resolution will be deleted.

Add Comment


E-mail address and name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.