PolygonArea.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP)
#define GEOGRAPHICLIB_POLYGONAREA_HPP 1

#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicExact.hpp>
#include <GeographicLib/Rhumb.hpp>
#include <GeographicLib/Accumulator.hpp>

namespace GeographicLib {

  template <class GeodType = Geodesic>
  class PolygonAreaT {
  private:
    typedef Math::real real;
    GeodType _earth;
    real _area0;                // Full ellipsoid area
    bool _polyline;             // Assume polyline (don't close and skip area)
    unsigned _mask;
    unsigned _num;
    int _crossings;
    Accumulator<> _areasum, _perimetersum;
    real _lat0, _lon0, _lat1, _lon1;
    static int transit(real lon1, real lon2) {
      // Return 1 or -1 if crossing prime meridian in east or west direction.
      // Otherwise return zero.
      // Compute lon12 the same way as Geodesic::Inverse.
      lon1 = Math::AngNormalize(lon1);
      lon2 = Math::AngNormalize(lon2);
      real lon12 = Math::AngDiff(lon1, lon2);
      // Treat 0 as negative in these tests.  This balances +/- 180 being
      // treated as positive, i.e., +180.
      int cross =
        lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
        (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
      return cross;
    }
    // an alternate version of transit to deal with longitudes in the direct
    // problem.
    static int transitdirect(real lon1, real lon2) {
      // We want to compute exactly
      //   int(floor(lon2 / 360)) - int(floor(lon1 / 360))
      // Since we only need the parity of the result we can use std::remquo;
      // but this is buggy with g++ 4.8.3 (glibc version < 2.22), see
      //   https://sourceware.org/bugzilla/show_bug.cgi?id=17569
      // and requires C++11.  So instead we do
#if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
      using std::remainder;
      lon1 = remainder(lon1, real(720)); lon2 = remainder(lon2, real(720));
      return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) -
               (lon1 >= 0 && lon1 < 360 ? 0 : 1) );
#else
      using std::fmod;
      lon1 = fmod(lon1, real(720)); lon2 = fmod(lon2, real(720));
      return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
               ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
#endif
    }
  public:

    PolygonAreaT(const GeodType& earth, bool polyline = false)
      : _earth(earth)
      , _area0(_earth.EllipsoidArea())
      , _polyline(polyline)
      , _mask(GeodType::LATITUDE | GeodType::LONGITUDE | GeodType::DISTANCE |
              (_polyline ? GeodType::NONE :
               GeodType::AREA | GeodType::LONG_UNROLL))
    { Clear(); }

    void Clear() {
      _num = 0;
      _crossings = 0;
      _areasum = 0;
      _perimetersum = 0;
      _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN();
    }

    void AddPoint(real lat, real lon);

    void AddEdge(real azi, real s);

    unsigned Compute(bool reverse, bool sign,
                     real& perimeter, real& area) const;

    unsigned TestPoint(real lat, real lon, bool reverse, bool sign,
                       real& perimeter, real& area) const;

    unsigned TestEdge(real azi, real s, bool reverse, bool sign,
                      real& perimeter, real& area) const;


    Math::real MajorRadius() const { return _earth.MajorRadius(); }

    Math::real Flattening() const { return _earth.Flattening(); }

    void CurrentPoint(real& lat, real& lon) const
    { lat = _lat1; lon = _lon1; }
  };

  typedef PolygonAreaT<Geodesic> PolygonArea;

  typedef PolygonAreaT<GeodesicExact> PolygonAreaExact;

  typedef PolygonAreaT<Rhumb> PolygonAreaRhumb;

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_POLYGONAREA_HPP