00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/PolygonArea.hpp>
00011
00012 namespace GeographicLib {
00013
00014 using namespace std;
00015
00016 template <class GeodType>
00017 void PolygonAreaT<GeodType>::AddPoint(real lat, real lon) {
00018 lon = Math::AngNormalize(lon);
00019 if (_num == 0) {
00020 _lat0 = _lat1 = lat;
00021 _lon0 = _lon1 = lon;
00022 } else {
00023 real s12, S12, t;
00024 _earth.GenInverse(_lat1, _lon1, lat, lon, _mask, s12, t, t, t, t, t, S12);
00025 _perimetersum += s12;
00026 if (!_polyline) {
00027 _areasum += S12;
00028 _crossings += transit(_lon1, lon);
00029 }
00030 _lat1 = lat; _lon1 = lon;
00031 }
00032 ++_num;
00033 }
00034
00035 template <class GeodType>
00036 void PolygonAreaT<GeodType>::AddEdge(real azi, real s) {
00037 if (_num) {
00038 real lat, lon, S12, t;
00039 _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
00040 lat, lon, t, t, t, t, t, S12);
00041 _perimetersum += s;
00042 if (!_polyline) {
00043 _areasum += S12;
00044 _crossings += transit(_lon1, lon);
00045 }
00046 _lat1 = lat; _lon1 = lon;
00047 ++_num;
00048 }
00049 }
00050
00051 template <class GeodType>
00052 unsigned PolygonAreaT<GeodType>::Compute(bool reverse, bool sign,
00053 real& perimeter, real& area) const {
00054 real s12, S12, t;
00055 if (_num < 2) {
00056 perimeter = 0;
00057 if (!_polyline)
00058 area = 0;
00059 return _num;
00060 }
00061 if (_polyline) {
00062 perimeter = _perimetersum();
00063 return _num;
00064 }
00065 _earth.GenInverse(_lat1, _lon1, _lat0, _lon0, _mask,
00066 s12, t, t, t, t, t, S12);
00067 perimeter = _perimetersum(s12);
00068 Accumulator<> tempsum(_areasum);
00069 tempsum += S12;
00070 int crossings = _crossings + transit(_lon1, _lon0);
00071 if (crossings & 1)
00072 tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
00073
00074
00075 if (!reverse)
00076 tempsum *= -1;
00077
00078 if (sign) {
00079 if (tempsum > _area0/2)
00080 tempsum -= _area0;
00081 else if (tempsum <= -_area0/2)
00082 tempsum += _area0;
00083 } else {
00084 if (tempsum >= _area0)
00085 tempsum -= _area0;
00086 else if (tempsum < 0)
00087 tempsum += _area0;
00088 }
00089 area = 0 + tempsum();
00090 return _num;
00091 }
00092
00093 template <class GeodType>
00094 unsigned PolygonAreaT<GeodType>::TestPoint(real lat, real lon,
00095 bool reverse, bool sign,
00096 real& perimeter, real& area) const
00097 {
00098 if (_num == 0) {
00099 perimeter = 0;
00100 if (!_polyline)
00101 area = 0;
00102 return 1;
00103 }
00104 perimeter = _perimetersum();
00105 real tempsum = _polyline ? 0 : _areasum();
00106 int crossings = _crossings;
00107 unsigned num = _num + 1;
00108 for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
00109 real s12, S12, t;
00110 _earth.GenInverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
00111 i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
00112 _mask, s12, t, t, t, t, t, S12);
00113 perimeter += s12;
00114 if (!_polyline) {
00115 tempsum += S12;
00116 crossings += transit(i == 0 ? _lon1 : lon,
00117 i != 0 ? _lon0 : lon);
00118 }
00119 }
00120
00121 if (_polyline)
00122 return num;
00123
00124 if (crossings & 1)
00125 tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
00126
00127
00128 if (!reverse)
00129 tempsum *= -1;
00130
00131 if (sign) {
00132 if (tempsum > _area0/2)
00133 tempsum -= _area0;
00134 else if (tempsum <= -_area0/2)
00135 tempsum += _area0;
00136 } else {
00137 if (tempsum >= _area0)
00138 tempsum -= _area0;
00139 else if (tempsum < 0)
00140 tempsum += _area0;
00141 }
00142 area = 0 + tempsum;
00143 return num;
00144 }
00145
00146 template <class GeodType>
00147 unsigned PolygonAreaT<GeodType>::TestEdge(real azi, real s,
00148 bool reverse, bool sign,
00149 real& perimeter, real& area) const {
00150 if (_num == 0) {
00151 perimeter = Math::NaN();
00152 if (!_polyline)
00153 area = Math::NaN();
00154 return 0;
00155 }
00156 unsigned num = _num + 1;
00157 perimeter = _perimetersum() + s;
00158 if (_polyline)
00159 return num;
00160
00161 real tempsum = _areasum();
00162 int crossings = _crossings;
00163 {
00164 real lat, lon, s12, S12, t;
00165 _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
00166 lat, lon, t, t, t, t, t, S12);
00167 tempsum += S12;
00168 crossings += transit(_lon1, lon);
00169 _earth.GenInverse(lat, lon, _lat0, _lon0, _mask, s12, t, t, t, t, t, S12);
00170 perimeter += s12;
00171 tempsum += S12;
00172 crossings += transit(lon, _lon0);
00173 }
00174
00175 if (crossings & 1)
00176 tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
00177
00178
00179 if (!reverse)
00180 tempsum *= -1;
00181
00182 if (sign) {
00183 if (tempsum > _area0/2)
00184 tempsum -= _area0;
00185 else if (tempsum <= -_area0/2)
00186 tempsum += _area0;
00187 } else {
00188 if (tempsum >= _area0)
00189 tempsum -= _area0;
00190 else if (tempsum < 0)
00191 tempsum += _area0;
00192 }
00193 area = 0 + tempsum;
00194 return num;
00195 }
00196
00197 template class GEOGRAPHICLIB_EXPORT PolygonAreaT<Geodesic>;
00198 template class GEOGRAPHICLIB_EXPORT PolygonAreaT<GeodesicExact>;
00199
00200 }