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