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