GeographicLib 2.5
PolygonArea.cpp
Go to the documentation of this file.
1/**
2 * \file PolygonArea.cpp
3 * \brief Implementation for GeographicLib::PolygonAreaT class
4 *
5 * Copyright (c) Charles Karney (2010-2023) <karney@alum.mit.edu> and licensed
6 * under the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 **********************************************************************/
9
11
12namespace GeographicLib {
13
14 using namespace std;
15
16 template<class GeodType>
17 int PolygonAreaT<GeodType>::transit(real lon1, real lon2) {
18 // Return 1 or -1 if crossing prime meridian in east or west direction.
19 // Otherwise return zero. longitude = +/-0 considered to be positive.
20 // This is (should be?) compatible with transitdirect which computes
21 // exactly the parity of
22 // int(floor((lon1 + lon12) / 360)) - int(floor(lon1 / 360)))
23 real lon12 = Math::AngDiff(lon1, lon2);
24 lon1 = Math::AngNormalize(lon1);
25 lon2 = Math::AngNormalize(lon2);
26 // N.B. lon12 == 0 gives cross = 0
27 return
28 // edge case lon1 = 180, lon2 = 360->0, lon12 = 180 to give 1
29 lon12 > 0 && ((lon1 < 0 && lon2 >= 0) ||
30 // lon12 > 0 && lon1 > 0 && lon2 == 0 implies lon1 == 180
31 (lon1 > 0 && lon2 == 0)) ? 1 :
32 // non edge case lon1 = -180, lon2 = -360->-0, lon12 = -180
33 (lon12 < 0 && lon1 >= 0 && lon2 < 0 ? -1 : 0);
34 // This was the old method (treating +/- 0 as negative). However, with the
35 // new scheme for handling longitude differences this fails on:
36 // lon1 = -180, lon2 = -360->-0, lon12 = -180 gives 0 not -1.
37 // return
38 // lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
39 // (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
40 }
41
42 // an alternate version of transit to deal with longitudes in the direct
43 // problem.
44 template<class GeodType>
45 int PolygonAreaT<GeodType>::transitdirect(real lon1, real lon2) {
46 // Compute exactly the parity of
47 // int(floor(lon2 / 360)) - int(floor(lon1 / 360))
48 // C++ C remainder -> [-360, 360]
49 // Java % -> (-720, 720) switch to IEEEremainder -> [-360, 360]
50 // JS % -> (-720, 720)
51 // Python fmod -> (-720, 720) swith to Math.remainder
52 // Fortran, Octave skip
53 // If mod function gives result in [-360, 360]
54 // [0, 360) -> 0; [-360, 0) or 360 -> 1
55 // If mod function gives result in (-720, 720)
56 // [0, 360) or [-inf, -360) -> 0; [-360, 0) or [360, inf) -> 1
57 lon1 = remainder(lon1, real(2 * Math::td));
58 lon2 = remainder(lon2, real(2 * Math::td));
59 return ( (lon2 >= 0 && lon2 < Math::td ? 0 : 1) -
60 (lon1 >= 0 && lon1 < Math::td ? 0 : 1) );
61 }
62
63 template<class GeodType>
64 void PolygonAreaT<GeodType>::AddPoint(real lat, real lon) {
65 if (_num == 0) {
66 _lat0 = _lat1 = lat;
67 _lon0 = _lon1 = lon;
68 } else {
69 real s12, S12, t;
70 _earth.GenInverse(_lat1, _lon1, lat, lon, _mask,
71 s12, t, t, t, t, t, S12);
72 _perimetersum += s12;
73 if (!_polyline) {
74 _areasum += S12;
75 _crossings += transit(_lon1, lon);
76 }
77 _lat1 = lat; _lon1 = lon;
78 }
79 ++_num;
80 }
81
82 template<class GeodType>
83 void PolygonAreaT<GeodType>::AddEdge(real azi, real s) {
84 if (_num) { // Do nothing if _num is zero
85 real lat, lon, S12, t;
86 _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
87 lat, lon, t, t, t, t, t, S12);
88 _perimetersum += s;
89 if (!_polyline) {
90 _areasum += S12;
91 _crossings += transitdirect(_lon1, lon);
92 }
93 _lat1 = lat; _lon1 = lon;
94 ++_num;
95 }
96 }
97
98 template<class GeodType>
99 unsigned PolygonAreaT<GeodType>::Compute(bool reverse, bool sign,
100 real& perimeter, real& area) const
101 {
102 real s12, S12, t;
103 if (_num < 2) {
104 perimeter = 0;
105 if (!_polyline)
106 area = 0;
107 return _num;
108 }
109 if (_polyline) {
110 perimeter = _perimetersum();
111 return _num;
112 }
113 _earth.GenInverse(_lat1, _lon1, _lat0, _lon0, _mask,
114 s12, t, t, t, t, t, S12);
115 perimeter = _perimetersum(s12);
116 Accumulator<> tempsum(_areasum);
117 tempsum += S12;
118 int crossings = _crossings + transit(_lon1, _lon0);
119 AreaReduce(tempsum, crossings, reverse, sign);
120 area = real(0) + tempsum();
121 return _num;
122 }
123
124 template<class GeodType>
125 unsigned PolygonAreaT<GeodType>::TestPoint(real lat, real lon,
126 bool reverse, bool sign,
127 real& perimeter, real& area) const
128 {
129 if (_num == 0) {
130 perimeter = 0;
131 if (!_polyline)
132 area = 0;
133 return 1;
134 }
135 perimeter = _perimetersum();
136 real tempsum = _polyline ? 0 : _areasum();
137 int crossings = _crossings;
138 unsigned num = _num + 1;
139 for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
140 real s12, S12, t;
141 _earth.GenInverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
142 i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
143 _mask, s12, t, t, t, t, t, S12);
144 perimeter += s12;
145 if (!_polyline) {
146 tempsum += S12;
147 crossings += transit(i == 0 ? _lon1 : lon,
148 i != 0 ? _lon0 : lon);
149 }
150 }
151
152 if (_polyline)
153 return num;
154
155 AreaReduce(tempsum, crossings, reverse, sign);
156 area = real(0) + tempsum;
157 return num;
158 }
159
160 template<class GeodType>
161 unsigned PolygonAreaT<GeodType>::TestEdge(real azi, real s,
162 bool reverse, bool sign,
163 real& perimeter, real& area) const
164 {
165 if (_num == 0) { // we don't have a starting point!
166 perimeter = Math::NaN();
167 if (!_polyline)
168 area = Math::NaN();
169 return 0;
170 }
171 unsigned num = _num + 1;
172 perimeter = _perimetersum() + s;
173 if (_polyline)
174 return num;
175
176 real tempsum = _areasum();
177 int crossings = _crossings;
178 {
179 real lat, lon, s12, S12, t;
180 _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
181 lat, lon, t, t, t, t, t, S12);
182 tempsum += S12;
183 crossings += transitdirect(_lon1, lon);
184 _earth.GenInverse(lat, lon, _lat0, _lon0, _mask,
185 s12, t, t, t, t, t, S12);
186 perimeter += s12;
187 tempsum += S12;
188 crossings += transit(lon, _lon0);
189 }
190
191 AreaReduce(tempsum, crossings, reverse, sign);
192 area = real(0) + tempsum;
193 return num;
194 }
195
196 template<class GeodType>
197 template<typename T>
198 void PolygonAreaT<GeodType>::AreaReduce(T& area, int crossings,
199 bool reverse, bool sign) const {
200 Remainder(area);
201 if (crossings & 1) area += (area < 0 ? 1 : -1) * _area0/2;
202 // area is with the clockwise sense. If !reverse convert to
203 // counter-clockwise convention.
204 if (!reverse) area *= -1;
205 // If sign put area in (-_area0/2, _area0/2], else put area in [0, _area0)
206 if (sign) {
207 if (area > _area0/2)
208 area -= _area0;
209 else if (area <= -_area0/2)
210 area += _area0;
211 } else {
212 if (area >= _area0)
213 area -= _area0;
214 else if (area < 0)
215 area += _area0;
216 }
217 }
218
219 template class GEOGRAPHICLIB_EXPORT PolygonAreaT<Geodesic>;
220 template class GEOGRAPHICLIB_EXPORT PolygonAreaT<GeodesicExact>;
221 template class GEOGRAPHICLIB_EXPORT PolygonAreaT<Rhumb>;
222
223} // namespace GeographicLib
#define GEOGRAPHICLIB_EXPORT
Definition Constants.hpp:67
GeographicLib::Math::real real
Definition GeodSolve.cpp:28
Header for GeographicLib::PolygonAreaT class.
An accumulator for sums.
static T AngNormalize(T x)
Definition Math.cpp:66
static constexpr int td
degrees per turn
Definition Math.hpp:149
static T NaN()
Definition Math.cpp:277
static T AngDiff(T x, T y, T &e)
Definition Math.cpp:77
void AddPoint(real lat, real lon)
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
void AddEdge(real azi, real s)
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
Namespace for GeographicLib.