GeographicLib  1.51
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-2019) <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  void PolygonAreaT<GeodType>::AddPoint(real lat, real lon) {
18  lat = Math::LatFix(lat);
19  lon = Math::AngNormalize(lon);
20  if (_num == 0) {
21  _lat0 = _lat1 = lat;
22  _lon0 = _lon1 = lon;
23  } else {
24  real s12, S12, t;
25  _earth.GenInverse(_lat1, _lon1, lat, lon, _mask,
26  s12, t, t, t, t, t, S12);
27  _perimetersum += s12;
28  if (!_polyline) {
29  _areasum += S12;
30  _crossings += transit(_lon1, lon);
31  }
32  _lat1 = lat; _lon1 = lon;
33  }
34  ++_num;
35  }
36 
37  template <class GeodType>
38  void PolygonAreaT<GeodType>::AddEdge(real azi, real s) {
39  if (_num) { // Do nothing if _num is zero
40  real lat, lon, S12, t;
41  _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
42  lat, lon, t, t, t, t, t, S12);
43  _perimetersum += s;
44  if (!_polyline) {
45  _areasum += S12;
46  _crossings += transitdirect(_lon1, lon);
47  lon = Math::AngNormalize(lon);
48  }
49  _lat1 = lat; _lon1 = lon;
50  ++_num;
51  }
52  }
53 
54  template <class GeodType>
55  unsigned PolygonAreaT<GeodType>::Compute(bool reverse, bool sign,
56  real& perimeter, real& area) const
57  {
58  real s12, S12, t;
59  if (_num < 2) {
60  perimeter = 0;
61  if (!_polyline)
62  area = 0;
63  return _num;
64  }
65  if (_polyline) {
66  perimeter = _perimetersum();
67  return _num;
68  }
69  _earth.GenInverse(_lat1, _lon1, _lat0, _lon0, _mask,
70  s12, t, t, t, t, t, S12);
71  perimeter = _perimetersum(s12);
72  Accumulator<> tempsum(_areasum);
73  tempsum += S12;
74  int crossings = _crossings + transit(_lon1, _lon0);
75  AreaReduce(tempsum, crossings, reverse, sign);
76  area = 0 + tempsum();
77  return _num;
78  }
79 
80  template <class GeodType>
81  unsigned PolygonAreaT<GeodType>::TestPoint(real lat, real lon,
82  bool reverse, bool sign,
83  real& perimeter, real& area) const
84  {
85  if (_num == 0) {
86  perimeter = 0;
87  if (!_polyline)
88  area = 0;
89  return 1;
90  }
91  perimeter = _perimetersum();
92  real tempsum = _polyline ? 0 : _areasum();
93  int crossings = _crossings;
94  unsigned num = _num + 1;
95  for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
96  real s12, S12, t;
97  _earth.GenInverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
98  i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
99  _mask, s12, t, t, t, t, t, S12);
100  perimeter += s12;
101  if (!_polyline) {
102  tempsum += S12;
103  crossings += transit(i == 0 ? _lon1 : lon,
104  i != 0 ? _lon0 : lon);
105  }
106  }
107 
108  if (_polyline)
109  return num;
110 
111  AreaReduce(tempsum, crossings, reverse, sign);
112  area = 0 + tempsum;
113  return num;
114  }
115 
116  template <class GeodType>
117  unsigned PolygonAreaT<GeodType>::TestEdge(real azi, real s,
118  bool reverse, bool sign,
119  real& perimeter, real& area) const
120  {
121  if (_num == 0) { // we don't have a starting point!
122  perimeter = Math::NaN();
123  if (!_polyline)
124  area = Math::NaN();
125  return 0;
126  }
127  unsigned num = _num + 1;
128  perimeter = _perimetersum() + s;
129  if (_polyline)
130  return num;
131 
132  real tempsum = _areasum();
133  int crossings = _crossings;
134  {
135  real lat, lon, s12, S12, t;
136  _earth.GenDirect(_lat1, _lon1, azi, false, s, _mask,
137  lat, lon, t, t, t, t, t, S12);
138  tempsum += S12;
139  crossings += transitdirect(_lon1, lon);
140  lon = Math::AngNormalize(lon);
141  _earth.GenInverse(lat, lon, _lat0, _lon0, _mask,
142  s12, t, t, t, t, t, S12);
143  perimeter += s12;
144  tempsum += S12;
145  crossings += transit(lon, _lon0);
146  }
147 
148  AreaReduce(tempsum, crossings, reverse, sign);
149  area = 0 + tempsum;
150  return num;
151  }
152 
153  template <class GeodType>
154  template <typename T>
155  void PolygonAreaT<GeodType>::AreaReduce(T& area, int crossings, bool reverse,
156  bool sign) const {
157  Remainder(area);
158  if (crossings & 1) area += (area < 0 ? 1 : -1) * _area0/2;
159  // area is with the clockwise sense. If !reverse convert to
160  // counter-clockwise convention.
161  if (!reverse) area *= -1;
162  // If sign put area in (-_area0/2, _area0/2], else put area in [0, _area0)
163  if (sign) {
164  if (area > _area0/2)
165  area -= _area0;
166  else if (area <= -_area0/2)
167  area += _area0;
168  } else {
169  if (area >= _area0)
170  area -= _area0;
171  else if (area < 0)
172  area += _area0;
173  }
174  }
175 
179 
180 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:405
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:66
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
static T LatFix(T x)
Definition: Math.hpp:418
void AddEdge(real azi, real s)
Definition: PolygonArea.cpp:38
An accumulator for sums.
Definition: Accumulator.hpp:40
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
unsigned TestPoint(real lat, real lon, bool reverse, bool sign, real &perimeter, real &area) const
Definition: PolygonArea.cpp:81
static T NaN()
Definition: Math.cpp:268
unsigned Compute(bool reverse, bool sign, real &perimeter, real &area) const
Definition: PolygonArea.cpp:55
void AddPoint(real lat, real lon)
Definition: PolygonArea.cpp:17
Header for GeographicLib::PolygonAreaT class.