GeographicLib  2.0
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 namespace 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  using std::remainder;
49  // C++ C remainder -> [-360, 360]
50  // Java % -> (-720, 720) switch to IEEEremainder -> [-360, 360]
51  // JS % -> (-720, 720)
52  // Python fmod -> (-720, 720) swith to Math.remainder
53  // Fortran, Octave skip
54  // If mod function gives result in [-360, 360]
55  // [0, 360) -> 0; [-360, 0) or 360 -> 1
56  // If mod function gives result in (-720, 720)
57  // [0, 360) or [-inf, -360) -> 0; [-360, 0) or [360, inf) -> 1
58  lon1 = remainder(lon1, real(2 * Math::td));
59  lon2 = remainder(lon2, real(2 * Math::td));
60  return ( (lon2 >= 0 && lon2 < Math::td ? 0 : 1) -
61  (lon1 >= 0 && lon1 < Math::td ? 0 : 1) );
62  }
63 
64  template <class GeodType>
65  void PolygonAreaT<GeodType>::AddPoint(real lat, real lon) {
66  if (_num == 0) {
67  _lat0 = _lat1 = lat;
68  _lon0 = _lon1 = lon;
69  } else {
70  real s12, S12, t;
71  _earth.GenInverse(_lat1, _lon1, lat, lon, _mask,
72  s12, t, t, t, t, t, S12);
73  _perimetersum += s12;
74  if (!_polyline) {
75  _areasum += S12;
76  _crossings += transit(_lon1, lon);
77  }
78  _lat1 = lat; _lon1 = lon;
79  }
80  ++_num;
81  }
82 
83  template <class GeodType>
84  void PolygonAreaT<GeodType>::AddEdge(real azi, real s) {
85  if (_num) { // Do nothing if _num is zero
86  real lat, lon, S12, t;
87  _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
88  lat, lon, t, t, t, t, t, S12);
89  _perimetersum += s;
90  if (!_polyline) {
91  _areasum += S12;
92  _crossings += transitdirect(_lon1, lon);
93  }
94  _lat1 = lat; _lon1 = lon;
95  ++_num;
96  }
97  }
98 
99  template <class GeodType>
100  unsigned PolygonAreaT<GeodType>::Compute(bool reverse, bool sign,
101  real& perimeter, real& area) const
102  {
103  real s12, S12, t;
104  if (_num < 2) {
105  perimeter = 0;
106  if (!_polyline)
107  area = 0;
108  return _num;
109  }
110  if (_polyline) {
111  perimeter = _perimetersum();
112  return _num;
113  }
114  _earth.GenInverse(_lat1, _lon1, _lat0, _lon0, _mask,
115  s12, t, t, t, t, t, S12);
116  perimeter = _perimetersum(s12);
117  Accumulator<> tempsum(_areasum);
118  tempsum += S12;
119  int crossings = _crossings + transit(_lon1, _lon0);
120  AreaReduce(tempsum, crossings, reverse, sign);
121  area = real(0) + tempsum();
122  return _num;
123  }
124 
125  template <class GeodType>
126  unsigned PolygonAreaT<GeodType>::TestPoint(real lat, real lon,
127  bool reverse, bool sign,
128  real& perimeter, real& area) const
129  {
130  if (_num == 0) {
131  perimeter = 0;
132  if (!_polyline)
133  area = 0;
134  return 1;
135  }
136  perimeter = _perimetersum();
137  real tempsum = _polyline ? 0 : _areasum();
138  int crossings = _crossings;
139  unsigned num = _num + 1;
140  for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
141  real s12, S12, t;
142  _earth.GenInverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
143  i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
144  _mask, s12, t, t, t, t, t, S12);
145  perimeter += s12;
146  if (!_polyline) {
147  tempsum += S12;
148  crossings += transit(i == 0 ? _lon1 : lon,
149  i != 0 ? _lon0 : lon);
150  }
151  }
152 
153  if (_polyline)
154  return num;
155 
156  AreaReduce(tempsum, crossings, reverse, sign);
157  area = real(0) + tempsum;
158  return num;
159  }
160 
161  template <class GeodType>
162  unsigned PolygonAreaT<GeodType>::TestEdge(real azi, real s,
163  bool reverse, bool sign,
164  real& perimeter, real& area) const
165  {
166  if (_num == 0) { // we don't have a starting point!
167  perimeter = Math::NaN();
168  if (!_polyline)
169  area = Math::NaN();
170  return 0;
171  }
172  unsigned num = _num + 1;
173  perimeter = _perimetersum() + s;
174  if (_polyline)
175  return num;
176 
177  real tempsum = _areasum();
178  int crossings = _crossings;
179  {
180  real lat, lon, s12, S12, t;
181  _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
182  lat, lon, t, t, t, t, t, S12);
183  tempsum += S12;
184  crossings += transitdirect(_lon1, lon);
185  _earth.GenInverse(lat, lon, _lat0, _lon0, _mask,
186  s12, t, t, t, t, t, S12);
187  perimeter += s12;
188  tempsum += S12;
189  crossings += transit(lon, _lon0);
190  }
191 
192  AreaReduce(tempsum, crossings, reverse, sign);
193  area = real(0) + tempsum;
194  return num;
195  }
196 
197  template <class GeodType>
198  template <typename T>
199  void PolygonAreaT<GeodType>::AreaReduce(T& area, int crossings,
200  bool reverse, bool sign) const {
201  Remainder(area);
202  if (crossings & 1) area += (area < 0 ? 1 : -1) * _area0/2;
203  // area is with the clockwise sense. If !reverse convert to
204  // counter-clockwise convention.
205  if (!reverse) area *= -1;
206  // If sign put area in (-_area0/2, _area0/2], else put area in [0, _area0)
207  if (sign) {
208  if (area > _area0/2)
209  area -= _area0;
210  else if (area <= -_area0/2)
211  area += _area0;
212  } else {
213  if (area >= _area0)
214  area -= _area0;
215  else if (area < 0)
216  area += _area0;
217  }
218  }
219 
223 
224 } // 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:65
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:84
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12