GeographicLib  1.35
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-2011) <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  * Example of use:
49  * \include example-PolygonArea.cpp
50  *
51  * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility
52  * providing access to the functionality of PolygonArea.
53  **********************************************************************/
54 
56  private:
57  typedef Math::real real;
58  Geodesic _earth;
59  real _area0; // Full ellipsoid area
60  bool _polyline; // Assume polyline (don't close and skip area)
61  unsigned _mask;
62  unsigned _num;
63  int _crossings;
64  Accumulator<real> _areasum, _perimetersum;
65  real _lat0, _lon0, _lat1, _lon1;
66  static inline int transit(real lon1, real lon2) throw() {
67  // Return 1 or -1 if crossing prime meridian in east or west direction.
68  // Otherwise return zero.
69  // Compute lon12 the same way as Geodesic::Inverse.
70  lon1 = Math::AngNormalize(lon1);
71  lon2 = Math::AngNormalize(lon2);
72  real lon12 = Math::AngDiff(lon1, lon2);
73  int cross =
74  lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
75  (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
76  return cross;
77  }
78  public:
79 
80  /**
81  * Constructor for PolygonArea.
82  *
83  * @param[in] earth the Geodesic object to use for geodesic calculations.
84  * By default this uses the WGS84 ellipsoid.
85  * @param[in] polyline if true that treat the points as defining a polyline
86  * instead of a polygon (default = false).
87  **********************************************************************/
88  PolygonArea(const Geodesic& earth, bool polyline = false) throw()
89  : _earth(earth)
90  , _area0(_earth.EllipsoidArea())
91  , _polyline(polyline)
92  , _mask(Geodesic::LATITUDE | Geodesic::LONGITUDE | Geodesic::DISTANCE |
93  (_polyline ? Geodesic::NONE : Geodesic::AREA))
94  { Clear(); }
95 
96  /**
97  * Clear PolygonArea, allowing a new polygon to be started.
98  **********************************************************************/
99  void Clear() throw() {
100  _num = 0;
101  _crossings = 0;
102  _areasum = 0;
103  _perimetersum = 0;
104  _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN<real>();
105  }
106 
107  /**
108  * Add a point to the polygon or polyline.
109  *
110  * @param[in] lat the latitude of the point (degrees).
111  * @param[in] lon the longitude of the point (degrees).
112  *
113  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
114  * lon should be in the range [&minus;540&deg;, 540&deg;).
115  **********************************************************************/
116  void AddPoint(real lat, real lon) throw();
117 
118  /**
119  * Add an edge to the polygon or polyline.
120  *
121  * @param[in] azi azimuth at current point (degrees).
122  * @param[in] s distance from current point to next point (meters).
123  *
124  * \e azi should be in the range [&minus;540&deg;, 540&deg;). This does
125  * nothing if no points have been added yet. Use PolygonArea::CurrentPoint
126  * to determine the position of the new vertex.
127  **********************************************************************/
128  void AddEdge(real azi, real s) throw();
129 
130  /**
131  * Return the results so far.
132  *
133  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
134  * traversal counts as a positive area.
135  * @param[in] sign if true then return a signed result for the area if
136  * the polygon is traversed in the "wrong" direction instead of returning
137  * the area for the rest of the earth.
138  * @param[out] perimeter the perimeter of the polygon or length of the
139  * polyline (meters).
140  * @param[out] area the area of the polygon (meters<sup>2</sup>); only set
141  * if \e polyline is false in the constructor.
142  * @return the number of points.
143  **********************************************************************/
144  unsigned Compute(bool reverse, bool sign,
145  real& perimeter, real& area) const throw();
146 
147  /**
148  * Return the results assuming a tentative final test point is added;
149  * however, the data for the test point is not saved. This lets you report
150  * a running result for the perimeter and area as the user moves the mouse
151  * cursor. Ordinary floating point arithmetic is used to accumulate the
152  * data for the test point; thus the area and perimeter returned are less
153  * accurate than if PolygonArea::AddPoint and PolygonArea::Compute are
154  * used.
155  *
156  * @param[in] lat the latitude of the test point (degrees).
157  * @param[in] lon the longitude of the test point (degrees).
158  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
159  * traversal counts as a positive area.
160  * @param[in] sign if true then return a signed result for the area if
161  * the polygon is traversed in the "wrong" direction instead of returning
162  * the area for the rest of the earth.
163  * @param[out] perimeter the approximate perimeter of the polygon or length
164  * of the polyline (meters).
165  * @param[out] area the approximate area of the polygon
166  * (meters<sup>2</sup>); only set if polyline is false in the
167  * constructor.
168  * @return the number of points.
169  *
170  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
171  * lon should be in the range [&minus;540&deg;, 540&deg;).
172  **********************************************************************/
173  unsigned TestPoint(real lat, real lon, bool reverse, bool sign,
174  real& perimeter, real& area) const throw();
175 
176  /**
177  * Return the results assuming a tentative final test point is added via an
178  * azimuth and distance; however, the data for the test point is not saved.
179  * This lets you report a running result for the perimeter and area as the
180  * user moves the mouse cursor. Ordinary floating point arithmetic is used
181  * to accumulate the data for the test point; thus the area and perimeter
182  * returned are less accurate than if PolygonArea::AddEdge and
183  * PolygonArea::Compute are used.
184  *
185  * @param[in] azi azimuth at current point (degrees).
186  * @param[in] s distance from current point to final test point (meters).
187  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
188  * traversal counts as a positive area.
189  * @param[in] sign if true then return a signed result for the area if
190  * the polygon is traversed in the "wrong" direction instead of returning
191  * the area for the rest of the earth.
192  * @param[out] perimeter the approximate perimeter of the polygon or length
193  * of the polyline (meters).
194  * @param[out] area the approximate area of the polygon
195  * (meters<sup>2</sup>); only set if polyline is false in the
196  * constructor.
197  * @return the number of points.
198  *
199  * \e azi should be in the range [&minus;540&deg;, 540&deg;).
200  **********************************************************************/
201  unsigned TestEdge(real azi, real s, bool reverse, bool sign,
202  real& perimeter, real& area) const throw();
203 
204  /// \cond SKIP
205  /**
206  * <b>DEPRECATED</b>
207  * The old name for PolygonArea::TestPoint.
208  **********************************************************************/
209  unsigned TestCompute(real lat, real lon, bool reverse, bool sign,
210  real& perimeter, real& area) const throw() {
211  return TestPoint(lat, lon, reverse, sign, perimeter, area);
212  }
213  /// \endcond
214 
215  /** \name Inspector functions
216  **********************************************************************/
217  ///@{
218  /**
219  * @return \e a the equatorial radius of the ellipsoid (meters). This is
220  * the value inherited from the Geodesic object used in the constructor.
221  **********************************************************************/
222 
223  Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
224 
225  /**
226  * @return \e f the flattening of the ellipsoid. This is the value
227  * inherited from the Geodesic object used in the constructor.
228  **********************************************************************/
229  Math::real Flattening() const throw() { return _earth.Flattening(); }
230 
231  /**
232  * Report the previous vertex added to the polygon or polyline.
233  *
234  * @param[out] lat the latitude of the point (degrees).
235  * @param[out] lon the longitude of the point (degrees).
236  *
237  * If no points have been added, then NaNs are returned. Otherwise, \e lon
238  * will be in the range [&minus;180&deg;, 180&deg;).
239  **********************************************************************/
240  void CurrentPoint(real& lat, real& lon) const throw()
241  { lat = _lat1; lon = _lon1; }
242  ///@}
243  };
244 
245 } // namespace GeographicLib
246 
247 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP
static T AngNormalize(T x)
Definition: Math.hpp:388
Math::real Flattening() const
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
void CurrentPoint(real &lat, real &lon) const
PolygonArea(const Geodesic &earth, bool polyline=false)
Definition: PolygonArea.hpp:88
Header for GeographicLib::Geodesic class.
Header for GeographicLib::Accumulator class.
static T AngDiff(T x, T y)
Definition: Math.hpp:418
Header for GeographicLib::Constants class.
Math::real MajorRadius() const
Geodesic calculations
Definition: Geodesic.hpp:169