[Boost-bugs] [Boost C++ Libraries] #11928: surveyor area method is inaccurate

Subject: [Boost-bugs] [Boost C++ Libraries] #11928: surveyor area method is inaccurate
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2016-01-21 19:00:40

#11928: surveyor area method is inaccurate
----------------------------------------+---------------------------
 Reporter:  Charles Karney <charles@…>  |      Owner:  barendgehrels
     Type:  Bugs                        |     Status:  new
Milestone:  To Be Determined            |  Component:  geometry
  Version:  Boost 1.58.0                |   Severity:  Problem
 Keywords:  area                        |
----------------------------------------+---------------------------
 The following code illustrates a problem with the "surveyor" strategy
 for computing planar areas.  It computes the area of a square of area
 2 and the same square translated by 10^8^ in x and y.  bg::area
 incorrectly computes the area in the second case as 0.

 {{{#!c++
 #include <iostream>
 #include <boost/geometry.hpp>

 namespace bg = boost::geometry;
 typedef bg::model::point<double, 2, bg::cs::cartesian> pt;

 int main() {
   bg::model::polygon<pt, false, false> poly;
   // A simple square of area 2
   bg::read_wkt("POLYGON((1 0,0 1,-1 0,0 -1))", poly);
   // Computes the correct area (= 2)
   std::cout << "area " << bg::area(poly) << "\n";
   // The same square shifted by (10^8 10^8)
   bg::read_wkt("POLYGON((100000001 100000000, \
                          100000000 100000001, \
                           99999999 100000000, \
                          100000000  99999999))", poly);
   // Should give 2 but instead gives 0
   std::cout << "area " << bg::area(poly) << "\n";
 }
 }}}

 Running this code gives
 {{{
 area 2
 area 0
 }}}

 Discussion:
 * Surely it's a bad idea to embed the algorithm name "surveyor" into
   the API.  Computing 2d areas is really straightforward.  Users
   shouldn't need to worry about the underlying method.  Maybe
   "planar_area" is a better name.
 * In any case "surveyor" is a rather silly name of the type beloved of
   Wikipedia editors.  There's no evidence that surveyors actually used
   this method.
 * It's more expensive (at 2*n multiplies and 2*n additions) than
   applying the trapezium rule (otherwise known as "doing the
   integral"), i.e., summing (x[i+1] - x[i]) * (y[i+1] + y[i]) / 2,
   which clocks in at n multiplies and 3*n additions.
 * Finally, as this example illustrates, it's subject to catastrophic
   loss of accuracy, while the trapezium rule will give the exact
   result.  (Kahan contrasts evaluating x*x - y*y compared with (x - y) *
   (x + y) with x = 1e15+1, y = 1e15-1.  The former result is wrong
   by 1.5% while the latter is exact.)

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/11928>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.