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