00001 /** 00002 * \file PolygonArea.hpp 00003 * \brief Header for GeographicLib::PolygonArea class 00004 * 00005 * Copyright (c) Charles Karney (2010-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP) 00011 #define GEOGRAPHICLIB_POLYGONAREA_HPP 1 00012 00013 #include <GeographicLib/Geodesic.hpp> 00014 #include <GeographicLib/GeodesicExact.hpp> 00015 #include <GeographicLib/Accumulator.hpp> 00016 00017 namespace GeographicLib { 00018 00019 /** 00020 * \brief Polygon areas 00021 * 00022 * This computes the area of a polygon whose edges are geodesics using the 00023 * method given in Section 6 of 00024 * - C. F. F. Karney, 00025 * <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> 00026 * Algorithms for geodesics</a>, 00027 * J. Geodesy <b>87</b>, 43--55 (2013); 00028 * DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> 00029 * 10.1007/s00190-012-0578-z</a>; 00030 * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html"> 00031 * geod-addenda.html</a>. 00032 * 00033 * This class lets you add vertices and edges one at a time to the polygon. 00034 * The sequence must start with a vertex and thereafter vertices and edges 00035 * can be added in any order. Any vertex after the first creates a new edge 00036 * which is the ''shortest'' geodesic from the previous vertex. In some 00037 * cases there may be two or many such shortest geodesics and the area is 00038 * then not uniquely defined. In this case, either add an intermediate 00039 * vertex or add the edge ''as'' an edge (by defining its direction and 00040 * length). 00041 * 00042 * The area and perimeter are accumulated in two times the standard floating 00043 * point precision to guard against the loss of accuracy with many-sided 00044 * polygons. At any point you can ask for the perimeter and area so far. 00045 * There's an option to treat the points as defining a polyline instead of a 00046 * polygon; in that case, only the perimeter is computed. 00047 * 00048 * This is a templated class to allow it to be used with either Geodesic and 00049 * GeodesicExact. GeographicLib::PolygonArea and 00050 * GeographicLib::PolygonAreaExact are typedefs for these two cases. 00051 * 00052 * @tparam GeodType the geodesic class to use. 00053 * 00054 * Example of use: 00055 * \include example-PolygonArea.cpp 00056 * 00057 * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility 00058 * providing access to the functionality of PolygonArea. 00059 **********************************************************************/ 00060 00061 template <class GeodType = Geodesic> 00062 class PolygonAreaT { 00063 private: 00064 typedef Math::real real; 00065 GeodType _earth; 00066 real _area0; // Full ellipsoid area 00067 bool _polyline; // Assume polyline (don't close and skip area) 00068 unsigned _mask; 00069 unsigned _num; 00070 int _crossings; 00071 Accumulator<> _areasum, _perimetersum; 00072 real _lat0, _lon0, _lat1, _lon1; 00073 static inline int transit(real lon1, real lon2) { 00074 // Return 1 or -1 if crossing prime meridian in east or west direction. 00075 // Otherwise return zero. 00076 // Compute lon12 the same way as Geodesic::Inverse. 00077 lon1 = Math::AngNormalize(lon1); 00078 lon2 = Math::AngNormalize(lon2); 00079 real lon12 = Math::AngDiff(lon1, lon2); 00080 int cross = 00081 lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 : 00082 (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0); 00083 return cross; 00084 } 00085 public: 00086 00087 /** 00088 * Constructor for PolygonAreaT. 00089 * 00090 * @param[in] earth the Geodesic object to use for geodesic calculations. 00091 * @param[in] polyline if true that treat the points as defining a polyline 00092 * instead of a polygon (default = false). 00093 **********************************************************************/ 00094 PolygonAreaT(const GeodType& earth, bool polyline = false) 00095 : _earth(earth) 00096 , _area0(_earth.EllipsoidArea()) 00097 , _polyline(polyline) 00098 , _mask(GeodType::LATITUDE | GeodType::LONGITUDE | GeodType::DISTANCE | 00099 (_polyline ? GeodType::NONE : GeodType::AREA)) 00100 { Clear(); } 00101 00102 /** 00103 * Clear PolygonAreaT, allowing a new polygon to be started. 00104 **********************************************************************/ 00105 void Clear() { 00106 _num = 0; 00107 _crossings = 0; 00108 _areasum = 0; 00109 _perimetersum = 0; 00110 _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN(); 00111 } 00112 00113 /** 00114 * Add a point to the polygon or polyline. 00115 * 00116 * @param[in] lat the latitude of the point (degrees). 00117 * @param[in] lon the longitude of the point (degrees). 00118 * 00119 * \e lat should be in the range [−90°, 90°] and \e 00120 * lon should be in the range [−540°, 540°). 00121 **********************************************************************/ 00122 void AddPoint(real lat, real lon); 00123 00124 /** 00125 * Add an edge to the polygon or polyline. 00126 * 00127 * @param[in] azi azimuth at current point (degrees). 00128 * @param[in] s distance from current point to next point (meters). 00129 * 00130 * \e azi should be in the range [−540°, 540°). This does 00131 * nothing if no points have been added yet. Use PolygonAreaT::CurrentPoint 00132 * to determine the position of the new vertex. 00133 **********************************************************************/ 00134 void AddEdge(real azi, real s); 00135 00136 /** 00137 * Return the results so far. 00138 * 00139 * @param[in] reverse if true then clockwise (instead of counter-clockwise) 00140 * traversal counts as a positive area. 00141 * @param[in] sign if true then return a signed result for the area if 00142 * the polygon is traversed in the "wrong" direction instead of returning 00143 * the area for the rest of the earth. 00144 * @param[out] perimeter the perimeter of the polygon or length of the 00145 * polyline (meters). 00146 * @param[out] area the area of the polygon (meters<sup>2</sup>); only set 00147 * if \e polyline is false in the constructor. 00148 * @return the number of points. 00149 **********************************************************************/ 00150 unsigned Compute(bool reverse, bool sign, 00151 real& perimeter, real& area) const; 00152 00153 /** 00154 * Return the results assuming a tentative final test point is added; 00155 * however, the data for the test point is not saved. This lets you report 00156 * a running result for the perimeter and area as the user moves the mouse 00157 * cursor. Ordinary floating point arithmetic is used to accumulate the 00158 * data for the test point; thus the area and perimeter returned are less 00159 * accurate than if PolygonAreaT::AddPoint and PolygonAreaT::Compute are 00160 * used. 00161 * 00162 * @param[in] lat the latitude of the test point (degrees). 00163 * @param[in] lon the longitude of the test point (degrees). 00164 * @param[in] reverse if true then clockwise (instead of counter-clockwise) 00165 * traversal counts as a positive area. 00166 * @param[in] sign if true then return a signed result for the area if 00167 * the polygon is traversed in the "wrong" direction instead of returning 00168 * the area for the rest of the earth. 00169 * @param[out] perimeter the approximate perimeter of the polygon or length 00170 * of the polyline (meters). 00171 * @param[out] area the approximate area of the polygon 00172 * (meters<sup>2</sup>); only set if polyline is false in the 00173 * constructor. 00174 * @return the number of points. 00175 * 00176 * \e lat should be in the range [−90°, 90°] and \e 00177 * lon should be in the range [−540°, 540°). 00178 **********************************************************************/ 00179 unsigned TestPoint(real lat, real lon, bool reverse, bool sign, 00180 real& perimeter, real& area) const; 00181 00182 /** 00183 * Return the results assuming a tentative final test point is added via an 00184 * azimuth and distance; however, the data for the test point is not saved. 00185 * This lets you report a running result for the perimeter and area as the 00186 * user moves the mouse cursor. Ordinary floating point arithmetic is used 00187 * to accumulate the data for the test point; thus the area and perimeter 00188 * returned are less accurate than if PolygonAreaT::AddEdge and 00189 * PolygonAreaT::Compute are used. 00190 * 00191 * @param[in] azi azimuth at current point (degrees). 00192 * @param[in] s distance from current point to final test point (meters). 00193 * @param[in] reverse if true then clockwise (instead of counter-clockwise) 00194 * traversal counts as a positive area. 00195 * @param[in] sign if true then return a signed result for the area if 00196 * the polygon is traversed in the "wrong" direction instead of returning 00197 * the area for the rest of the earth. 00198 * @param[out] perimeter the approximate perimeter of the polygon or length 00199 * of the polyline (meters). 00200 * @param[out] area the approximate area of the polygon 00201 * (meters<sup>2</sup>); only set if polyline is false in the 00202 * constructor. 00203 * @return the number of points. 00204 * 00205 * \e azi should be in the range [−540°, 540°). 00206 **********************************************************************/ 00207 unsigned TestEdge(real azi, real s, bool reverse, bool sign, 00208 real& perimeter, real& area) const; 00209 00210 /// \cond SKIP 00211 /** 00212 * <b>DEPRECATED</b> 00213 * The old name for PolygonAreaT::TestPoint. 00214 **********************************************************************/ 00215 unsigned TestCompute(real lat, real lon, bool reverse, bool sign, 00216 real& perimeter, real& area) const { 00217 return TestPoint(lat, lon, reverse, sign, perimeter, area); 00218 } 00219 /// \endcond 00220 00221 /** \name Inspector functions 00222 **********************************************************************/ 00223 ///@{ 00224 /** 00225 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00226 * the value inherited from the Geodesic object used in the constructor. 00227 **********************************************************************/ 00228 00229 Math::real MajorRadius() const { return _earth.MajorRadius(); } 00230 00231 /** 00232 * @return \e f the flattening of the ellipsoid. This is the value 00233 * inherited from the Geodesic object used in the constructor. 00234 **********************************************************************/ 00235 Math::real Flattening() const { return _earth.Flattening(); } 00236 00237 /** 00238 * Report the previous vertex added to the polygon or polyline. 00239 * 00240 * @param[out] lat the latitude of the point (degrees). 00241 * @param[out] lon the longitude of the point (degrees). 00242 * 00243 * If no points have been added, then NaNs are returned. Otherwise, \e lon 00244 * will be in the range [−180°, 180°). 00245 **********************************************************************/ 00246 void CurrentPoint(real& lat, real& lon) const 00247 { lat = _lat1; lon = _lon1; } 00248 ///@} 00249 }; 00250 00251 /** 00252 * @relates PolygonAreaT 00253 * 00254 * Polygon areas using Geodesic. This should be used if the flattening is 00255 * small. 00256 **********************************************************************/ 00257 typedef PolygonAreaT<Geodesic> PolygonArea; 00258 00259 /** 00260 * @relates PolygonAreaT 00261 * 00262 * Polygon areas using GeodesicExact. (But note that the implementation of 00263 * areas in GeodesicExact uses a high order series and this is only accurate 00264 * for modest flattenings.) 00265 **********************************************************************/ 00266 typedef PolygonAreaT<GeodesicExact> PolygonAreaExact; 00267 00268 } // namespace GeographicLib 00269 00270 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP