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