GeographicLib  1.21
PolygonArea.hpp
Go to the documentation of this file.
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