GeographicLib  1.21
PolygonArea.cpp
Go to the documentation of this file.
00001 /**
00002  * \file PolygonArea.cpp
00003  * \brief Implementation 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 #include <GeographicLib/PolygonArea.hpp>
00011 
00012 #define GEOGRAPHICLIB_POLYGONAREA_CPP \
00013   "$Id: ae2fce0b24653309ca8835d962b1a3e047a6768a $"
00014 
00015 RCSID_DECL(GEOGRAPHICLIB_POLYGONAREA_CPP)
00016 RCSID_DECL(GEOGRAPHICLIB_POLYGONAREA_HPP)
00017 RCSID_DECL(GEOGRAPHICLIB_ACCUMULATOR_HPP)
00018 
00019 namespace GeographicLib {
00020 
00021   using namespace std;
00022 
00023   void PolygonArea::AddPoint(real lat, real lon) throw() {
00024     if (_num == 0) {
00025       _lat0 = _lat1 = lat;
00026       _lon0 = _lon1 = lon;
00027     } else {
00028       real s12, S12, t;
00029       _earth.GenInverse(_lat1, _lon1, lat, lon, _mask, s12, t, t, t, t, t, S12);
00030       _perimetersum += s12;
00031       if (!_polyline) {
00032         _areasum += S12;
00033         _crossings += transit(_lon1, lon);
00034       }
00035       _lat1 = lat;
00036       _lon1 = lon;
00037     }
00038     ++_num;
00039   }
00040 
00041   unsigned PolygonArea::Compute(bool reverse, bool sign,
00042                                 real& perimeter, real& area) const throw() {
00043     real s12, S12, t;
00044     if (_num < 2) {
00045       perimeter = 0;
00046       if (!_polyline)
00047         area = 0;
00048       return _num;
00049     }
00050     if (_polyline) {
00051       perimeter = _perimetersum();
00052       return _num;
00053     }
00054     _earth.GenInverse(_lat1, _lon1, _lat0, _lon0, _mask,
00055                       s12, t, t, t, t, t, S12);
00056     perimeter = _perimetersum(s12);
00057     Accumulator<real> tempsum(_areasum);
00058     tempsum += S12;
00059     int crossings = _crossings + transit(_lon1, _lon0);
00060     if (crossings & 1)
00061       tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
00062     // area is with the clockwise sense.  If !reverse convert to
00063     // counter-clockwise convention.
00064     if (!reverse)
00065       tempsum *= -1;
00066     // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
00067     if (sign) {
00068       if (tempsum > _area0/2)
00069         tempsum -= _area0;
00070       else if (tempsum <= -_area0/2)
00071         tempsum += _area0;
00072     } else {
00073       if (tempsum >= _area0)
00074         tempsum -= _area0;
00075       else if (tempsum < 0)
00076         tempsum += _area0;
00077     }
00078     area = 0 + tempsum();
00079     return _num;
00080   }
00081 
00082   unsigned PolygonArea::TestCompute(real lat, real lon, bool reverse, bool sign,
00083                                     real& perimeter, real& area) const throw() {
00084     if (_num == 0) {
00085       perimeter = 0;
00086       if (!_polyline)
00087         area = 0;
00088       return 1;
00089     }
00090     perimeter = _perimetersum();
00091     real tempsum = _polyline ? 0 : _areasum();
00092     int crossings = _crossings, num = _num + 1;
00093     for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
00094       real s12, S12, t;
00095       _earth.GenInverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
00096                         i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
00097                         _mask, s12, t, t, t, t, t, S12);
00098       perimeter += s12;
00099       if (!_polyline) {
00100         tempsum += S12;
00101         crossings += transit(i == 0 ? _lon1 : lon,
00102                              i != 0 ? _lon0 : lon);
00103       }
00104     }
00105 
00106     if (_polyline)
00107       return num;
00108 
00109     if (crossings & 1)
00110       tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
00111     // area is with the clockwise sense.  If !reverse convert to
00112     // counter-clockwise convention.
00113     if (!reverse)
00114       tempsum *= -1;
00115     // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
00116     if (sign) {
00117       if (tempsum > _area0/2)
00118         tempsum -= _area0;
00119       else if (tempsum <= -_area0/2)
00120         tempsum += _area0;
00121     } else {
00122       if (tempsum >= _area0)
00123         tempsum -= _area0;
00124       else if (tempsum < 0)
00125         tempsum += _area0;
00126     }
00127     area = 0 + tempsum;
00128     return num;
00129   }
00130 
00131 } // namespace GeographicLib