GeographicLib
1.21
|
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