GeographicLib  1.37
PolygonArea.hpp
Go to the documentation of this file.
1 /**
2  * \file PolygonArea.hpp
3  * \brief Header for GeographicLib::PolygonArea class
4  *
5  * Copyright (c) Charles Karney (2010-2014) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP)
11 #define GEOGRAPHICLIB_POLYGONAREA_HPP 1
12 
16 
17 namespace GeographicLib {
18 
19  /**
20  * \brief Polygon areas
21  *
22  * This computes the area of a polygon whose edges are geodesics using the
23  * method given in Section 6 of
24  * - C. F. F. Karney,
25  * <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
26  * Algorithms for geodesics</a>,
27  * J. Geodesy <b>87</b>, 43--55 (2013);
28  * DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
29  * 10.1007/s00190-012-0578-z</a>;
30  * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
31  * geod-addenda.html</a>.
32  *
33  * This class lets you add vertices and edges one at a time to the polygon.
34  * The sequence must start with a vertex and thereafter vertices and edges
35  * can be added in any order. Any vertex after the first creates a new edge
36  * which is the ''shortest'' geodesic from the previous vertex. In some
37  * cases there may be two or many such shortest geodesics and the area is
38  * then not uniquely defined. In this case, either add an intermediate
39  * vertex or add the edge ''as'' an edge (by defining its direction and
40  * length).
41  *
42  * The area and perimeter are accumulated in two times the standard floating
43  * point precision to guard against the loss of accuracy with many-sided
44  * polygons. At any point you can ask for the perimeter and area so far.
45  * There's an option to treat the points as defining a polyline instead of a
46  * polygon; in that case, only the perimeter is computed.
47  *
48  * This is a templated class to allow it to be used with either Geodesic and
49  * GeodesicExact. GeographicLib::PolygonArea and
50  * GeographicLib::PolygonAreaExact are typedefs for these two cases.
51  *
52  * @tparam GeodType the geodesic class to use.
53  *
54  * Example of use:
55  * \include example-PolygonArea.cpp
56  *
57  * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility
58  * providing access to the functionality of PolygonArea.
59  **********************************************************************/
60 
61  template <class GeodType = Geodesic>
62  class PolygonAreaT {
63  private:
64  typedef Math::real real;
65  GeodType _earth;
66  real _area0; // Full ellipsoid area
67  bool _polyline; // Assume polyline (don't close and skip area)
68  unsigned _mask;
69  unsigned _num;
70  int _crossings;
71  Accumulator<> _areasum, _perimetersum;
72  real _lat0, _lon0, _lat1, _lon1;
73  static inline int transit(real lon1, real lon2) {
74  // Return 1 or -1 if crossing prime meridian in east or west direction.
75  // Otherwise return zero.
76  // Compute lon12 the same way as Geodesic::Inverse.
77  lon1 = Math::AngNormalize(lon1);
78  lon2 = Math::AngNormalize(lon2);
79  real lon12 = Math::AngDiff(lon1, lon2);
80  int cross =
81  lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
82  (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
83  return cross;
84  }
85  public:
86 
87  /**
88  * Constructor for PolygonAreaT.
89  *
90  * @param[in] earth the Geodesic object to use for geodesic calculations.
91  * @param[in] polyline if true that treat the points as defining a polyline
92  * instead of a polygon (default = false).
93  **********************************************************************/
94  PolygonAreaT(const GeodType& earth, bool polyline = false)
95  : _earth(earth)
96  , _area0(_earth.EllipsoidArea())
97  , _polyline(polyline)
98  , _mask(GeodType::LATITUDE | GeodType::LONGITUDE | GeodType::DISTANCE |
99  (_polyline ? GeodType::NONE : GeodType::AREA))
100  { Clear(); }
101 
102  /**
103  * Clear PolygonAreaT, allowing a new polygon to be started.
104  **********************************************************************/
105  void Clear() {
106  _num = 0;
107  _crossings = 0;
108  _areasum = 0;
109  _perimetersum = 0;
110  _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN();
111  }
112 
113  /**
114  * Add a point to the polygon or polyline.
115  *
116  * @param[in] lat the latitude of the point (degrees).
117  * @param[in] lon the longitude of the point (degrees).
118  *
119  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
120  * lon should be in the range [&minus;540&deg;, 540&deg;).
121  **********************************************************************/
122  void AddPoint(real lat, real lon);
123 
124  /**
125  * Add an edge to the polygon or polyline.
126  *
127  * @param[in] azi azimuth at current point (degrees).
128  * @param[in] s distance from current point to next point (meters).
129  *
130  * \e azi should be in the range [&minus;540&deg;, 540&deg;). This does
131  * nothing if no points have been added yet. Use PolygonAreaT::CurrentPoint
132  * to determine the position of the new vertex.
133  **********************************************************************/
134  void AddEdge(real azi, real s);
135 
136  /**
137  * Return the results so far.
138  *
139  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
140  * traversal counts as a positive area.
141  * @param[in] sign if true then return a signed result for the area if
142  * the polygon is traversed in the "wrong" direction instead of returning
143  * the area for the rest of the earth.
144  * @param[out] perimeter the perimeter of the polygon or length of the
145  * polyline (meters).
146  * @param[out] area the area of the polygon (meters<sup>2</sup>); only set
147  * if \e polyline is false in the constructor.
148  * @return the number of points.
149  **********************************************************************/
150  unsigned Compute(bool reverse, bool sign,
151  real& perimeter, real& area) const;
152 
153  /**
154  * Return the results assuming a tentative final test point is added;
155  * however, the data for the test point is not saved. This lets you report
156  * a running result for the perimeter and area as the user moves the mouse
157  * cursor. Ordinary floating point arithmetic is used to accumulate the
158  * data for the test point; thus the area and perimeter returned are less
159  * accurate than if PolygonAreaT::AddPoint and PolygonAreaT::Compute are
160  * used.
161  *
162  * @param[in] lat the latitude of the test point (degrees).
163  * @param[in] lon the longitude of the test point (degrees).
164  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
165  * traversal counts as a positive area.
166  * @param[in] sign if true then return a signed result for the area if
167  * the polygon is traversed in the "wrong" direction instead of returning
168  * the area for the rest of the earth.
169  * @param[out] perimeter the approximate perimeter of the polygon or length
170  * of the polyline (meters).
171  * @param[out] area the approximate area of the polygon
172  * (meters<sup>2</sup>); only set if polyline is false in the
173  * constructor.
174  * @return the number of points.
175  *
176  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
177  * lon should be in the range [&minus;540&deg;, 540&deg;).
178  **********************************************************************/
179  unsigned TestPoint(real lat, real lon, bool reverse, bool sign,
180  real& perimeter, real& area) const;
181 
182  /**
183  * Return the results assuming a tentative final test point is added via an
184  * azimuth and distance; however, the data for the test point is not saved.
185  * This lets you report a running result for the perimeter and area as the
186  * user moves the mouse cursor. Ordinary floating point arithmetic is used
187  * to accumulate the data for the test point; thus the area and perimeter
188  * returned are less accurate than if PolygonAreaT::AddEdge and
189  * PolygonAreaT::Compute are used.
190  *
191  * @param[in] azi azimuth at current point (degrees).
192  * @param[in] s distance from current point to final test point (meters).
193  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
194  * traversal counts as a positive area.
195  * @param[in] sign if true then return a signed result for the area if
196  * the polygon is traversed in the "wrong" direction instead of returning
197  * the area for the rest of the earth.
198  * @param[out] perimeter the approximate perimeter of the polygon or length
199  * of the polyline (meters).
200  * @param[out] area the approximate area of the polygon
201  * (meters<sup>2</sup>); only set if polyline is false in the
202  * constructor.
203  * @return the number of points.
204  *
205  * \e azi should be in the range [&minus;540&deg;, 540&deg;).
206  **********************************************************************/
207  unsigned TestEdge(real azi, real s, bool reverse, bool sign,
208  real& perimeter, real& area) const;
209 
210  /// \cond SKIP
211  /**
212  * <b>DEPRECATED</b>
213  * The old name for PolygonAreaT::TestPoint.
214  **********************************************************************/
215  unsigned TestCompute(real lat, real lon, bool reverse, bool sign,
216  real& perimeter, real& area) const {
217  return TestPoint(lat, lon, reverse, sign, perimeter, area);
218  }
219  /// \endcond
220 
221  /** \name Inspector functions
222  **********************************************************************/
223  ///@{
224  /**
225  * @return \e a the equatorial radius of the ellipsoid (meters). This is
226  * the value inherited from the Geodesic object used in the constructor.
227  **********************************************************************/
228 
229  Math::real MajorRadius() const { return _earth.MajorRadius(); }
230 
231  /**
232  * @return \e f the flattening of the ellipsoid. This is the value
233  * inherited from the Geodesic object used in the constructor.
234  **********************************************************************/
235  Math::real Flattening() const { return _earth.Flattening(); }
236 
237  /**
238  * Report the previous vertex added to the polygon or polyline.
239  *
240  * @param[out] lat the latitude of the point (degrees).
241  * @param[out] lon the longitude of the point (degrees).
242  *
243  * If no points have been added, then NaNs are returned. Otherwise, \e lon
244  * will be in the range [&minus;180&deg;, 180&deg;).
245  **********************************************************************/
246  void CurrentPoint(real& lat, real& lon) const
247  { lat = _lat1; lon = _lon1; }
248  ///@}
249  };
250 
251  /**
252  * @relates PolygonAreaT
253  *
254  * Polygon areas using Geodesic. This should be used if the flattening is
255  * small.
256  **********************************************************************/
258 
259  /**
260  * @relates PolygonAreaT
261  *
262  * Polygon areas using GeodesicExact. (But note that the implementation of
263  * areas in GeodesicExact uses a high order series and this is only accurate
264  * for modest flattenings.)
265  **********************************************************************/
267 
268 } // namespace GeographicLib
269 
270 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP
static T AngNormalize(T x)
Definition: Math.hpp:399
void CurrentPoint(real &lat, real &lon) const
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
static T NaN()
Definition: Math.hpp:460
unsigned TestPoint(real lat, real lon, bool reverse, bool sign, real &perimeter, real &area) const
Definition: PolygonArea.cpp:94
void AddEdge(real azi, real s)
Definition: PolygonArea.cpp:36
Math::real Flattening() const
PolygonAreaT(const GeodType &earth, bool polyline=false)
Definition: PolygonArea.hpp:94
An accumulator for sums.
Definition: Accumulator.hpp:40
Header for GeographicLib::Geodesic class.
PolygonAreaT< GeodesicExact > PolygonAreaExact
Header for GeographicLib::Accumulator class.
unsigned Compute(bool reverse, bool sign, real &perimeter, real &area) const
Definition: PolygonArea.cpp:52
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T AngDiff(T x, T y)
Definition: Math.hpp:429
void AddPoint(real lat, real lon)
Definition: PolygonArea.cpp:17
Header for GeographicLib::GeodesicExact class.
Math::real MajorRadius() const
PolygonAreaT< Geodesic > PolygonArea