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@…


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);
      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


  • 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


Fixed by this PR:

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

Add Comment

Modify Ticket

Change Properties
Set your email in Preferences
as closed The owner will remain barendgehrels.
The resolution will be deleted. Next status will be 'reopened'.

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

Note: See TracTickets for help on using tickets.