GeographicLib
1.21
|
00001 /** 00002 * \file PolygonArea.hpp 00003 * \brief Header for GeographicLib::PolygonArea class 00004 * 00005 * Copyright (c) Charles Karney (2010, 2011) <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 \ 00012 "$Id: 7a339f312a9c977b9fccad3c0c8bfa9009d863e2 $" 00013 00014 #include <GeographicLib/Geodesic.hpp> 00015 #include <GeographicLib/Constants.hpp> 00016 #include <GeographicLib/Accumulator.hpp> 00017 00018 namespace GeographicLib { 00019 00020 /** 00021 * \brief Polygon Areas. 00022 * 00023 * This computes the area of a geodesic polygon using the method given 00024 * Section 15 of 00025 * - C. F. F. Karney, 00026 * <a href="http://arxiv.org/abs/1102.1215v1">Geodesics 00027 * on an ellipsoid of revolution</a>, 00028 * Feb. 2011; 00029 * preprint 00030 * <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>. 00031 * . 00032 * See also Section 6 of 00033 * - C. F. F. Karney, 00034 * <a href="http://arxiv.org/abs/1109.4448">Algorithms for geodesics</a>, 00035 * Sept. 2011; 00036 * preprint 00037 * <a href="http://arxiv.org/abs/1109.4448">arxiv:1109.4448</a>. 00038 * 00039 * This class lets you add vertices one at a time to the polygon. The area 00040 * and perimeter are accumulated in two times the standard floating point 00041 * precision to guard against the loss of accuracy with many-sided polygons. 00042 * At any point you can ask for the perimeter and area so far. There's an 00043 * option to treat the points as defining a polyline instead of a polygon; in 00044 * that case, only the perimeter is computed. 00045 * 00046 * Example of use: 00047 * \include example-PolygonArea.cpp 00048 * 00049 * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility 00050 * providing access to the functionality of PolygonArea. 00051 **********************************************************************/ 00052 00053 class GEOGRAPHIC_EXPORT PolygonArea { 00054 private: 00055 typedef Math::real real; 00056 Geodesic _earth; 00057 real _area0; // Full ellipsoid area 00058 bool _polyline; // Assume polyline (don't close and skip area) 00059 unsigned _mask; 00060 unsigned _num; 00061 int _crossings; 00062 Accumulator<real> _areasum, _perimetersum; 00063 real _lat0, _lon0, _lat1, _lon1; 00064 // Copied from Geodesic class 00065 static inline real AngNormalize(real x) throw() { 00066 // Place angle in [-180, 180). Assumes x is in [-540, 540). 00067 // 00068 // g++ 4.4.4 holds a temporary in an extended register causing an error 00069 // with the triangle 89,0.1;89,90.1;89,-179.9. The volatile declaration 00070 // fixes this. (The bug probably triggered because transit and 00071 // AngNormalize are inline functions. So don't port this change over to 00072 // Geodesic.hpp.) 00073 volatile real y = x; 00074 return y >= 180 ? y - 360 : (y < -180 ? y + 360 : y); 00075 } 00076 static inline int transit(real lon1, real lon2) { 00077 // Return 1 or -1 if crossing prime meridian in east or west direction. 00078 // Otherwise return zero. 00079 lon1 = AngNormalize(lon1); 00080 lon2 = AngNormalize(lon2); 00081 // treat lon12 = -180 as an eastward geodesic, so convert to 180. 00082 real lon12 = -AngNormalize(lon1 - lon2); // In (-180, 180] 00083 int cross = 00084 lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 : 00085 (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0); 00086 return cross; 00087 } 00088 public: 00089 00090 /** 00091 * Constructor for PolygonArea. 00092 * 00093 * @param[in] earth the Geodesic object to use for geodesic calculations. 00094 * By default this uses the WGS84 ellipsoid. 00095 * @param[in] polyline if true that treat the points as defining a polyline 00096 * instead of a polygon (default = false). 00097 **********************************************************************/ 00098 PolygonArea(const Geodesic& earth, bool polyline = false) throw() 00099 : _earth(earth) 00100 , _area0(_earth.EllipsoidArea()) 00101 , _polyline(polyline) 00102 , _mask(Geodesic::DISTANCE | (_polyline ? 0 : Geodesic::AREA)) 00103 { 00104 Clear(); 00105 } 00106 00107 /** 00108 * Clear PolygonArea, allowing a new polygon to be started. 00109 **********************************************************************/ 00110 void Clear() throw() { 00111 _num = 0; 00112 _crossings = 0; 00113 _areasum = 0; 00114 _perimetersum = 0; 00115 _lat0 = _lon0 = _lat1 = _lon1 = 0; 00116 } 00117 00118 /** 00119 * Add a point to the polygon or polyline. 00120 * 00121 * @param[in] lat the latitude of the point (degrees). 00122 * @param[in] lon the latitude of the point (degrees). 00123 * 00124 * \e lat should be in the range [-90, 90] and \e lon should be in the 00125 * range [-180, 360]. 00126 **********************************************************************/ 00127 void AddPoint(real lat, real lon) throw(); 00128 00129 /** 00130 * Return the results so far. 00131 * 00132 * @param[in] reverse if true then clockwise (instead of counter-clockwise) 00133 * traversal counts as a positive area. 00134 * @param[in] sign if true then return a signed result for the area if 00135 * the polygon is traversed in the "wrong" direction instead of returning 00136 * the area for the rest of the earth. 00137 * @param[out] perimeter the perimeter of the polygon or length of the 00138 * polyline (meters). 00139 * @param[out] area the area of the polygon (meters^2); only set if 00140 * polyline is false in the constructor. 00141 * @return the number of points. 00142 **********************************************************************/ 00143 unsigned Compute(bool reverse, bool sign, 00144 real& perimeter, real& area) const throw(); 00145 00146 /** 00147 * Return the results assuming a tentative final test point is added; 00148 * however, the data for the test point is not saved. This lets you report 00149 * a running result for the perimeter and area as the user moves the mouse 00150 * cursor. Ordinary floating point arithmetic is used to accumulate the 00151 * data for the test point; thus the area and perimeter returned are less 00152 * accurate than if AddPoint and Compute are used. 00153 * 00154 * @param[in] lat the latitude of the test point (degrees). 00155 * @param[in] lon the longitude of the test point (degrees). 00156 * @param[in] reverse if true then clockwise (instead of counter-clockwise) 00157 * traversal counts as a positive area. 00158 * @param[in] sign if true then return a signed result for the area if 00159 * the polygon is traversed in the "wrong" direction instead of returning 00160 * the area for the rest of the earth. 00161 * @param[out] perimeter the approximate perimeter of the polygon or length 00162 * of the polyline (meters). 00163 * @param[out] area the approximate area of the polygon (meters^2); only 00164 * set if polyline is false in the constructor. 00165 * @return the number of points. 00166 * 00167 * \e lat should be in the range [-90, 90] and \e lon should be in the 00168 * range [-180, 360]. 00169 **********************************************************************/ 00170 unsigned TestCompute(real lat, real lon, bool reverse, bool sign, 00171 real& perimeter, real& area) const throw(); 00172 00173 /** \name Inspector functions 00174 **********************************************************************/ 00175 ///@{ 00176 /** 00177 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00178 * the value inherited from the Geodesic object used in the constructor. 00179 **********************************************************************/ 00180 00181 Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } 00182 00183 /** 00184 * @return \e f the flattening of the ellipsoid. This is the value 00185 * inherited from the Geodesic object used in the constructor. 00186 **********************************************************************/ 00187 Math::real Flattening() const throw() { return _earth.Flattening(); } 00188 ///@} 00189 }; 00190 00191 } // namespace GeographicLib 00192 00193 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP