PolygonArea.java
Go to the documentation of this file.
1 
8 package net.sf.geographiclib;
9 
60 public class PolygonArea {
61 
62  private Geodesic _earth;
63  private double _area0; // Full ellipsoid area
64  private boolean _polyline; // Assume polyline (don't close and skip area)
65  private int _mask;
66  private int _num;
67  private int _crossings;
68  private Accumulator _areasum, _perimetersum;
69  private double _lat0, _lon0, _lat1, _lon1;
70  private static int transit(double lon1, double lon2) {
71  // Return 1 or -1 if crossing prime meridian in east or west direction.
72  // Otherwise return zero.
73  // Compute lon12 the same way as Geodesic.Inverse.
74  lon1 = GeoMath.AngNormalize(lon1);
75  lon2 = GeoMath.AngNormalize(lon2);
76  double lon12 = GeoMath.AngDiff(lon1, lon2).first;
77  int cross =
78  lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
79  (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
80  return cross;
81  }
82  // an alternate version of transit to deal with longitudes in the direct
83  // problem.
84  private static int transitdirect(double lon1, double lon2) {
85  // We want to compute exactly
86  // int(floor(lon2 / 360)) - int(floor(lon1 / 360))
87  // Since we only need the parity of the result we can use std::remquo but
88  // this is buggy with g++ 4.8.3 and requires C++11. So instead we do
89  lon1 = lon1 % 720.0; lon2 = lon2 % 720.0;
90  return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
91  ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
92  }
93 
101  public PolygonArea(Geodesic earth, boolean polyline) {
102  _earth = earth;
103  _area0 = _earth.EllipsoidArea();
104  _polyline = polyline;
107  (_polyline ? GeodesicMask.NONE :
109  _perimetersum = new Accumulator(0);
110  if (!_polyline)
111  _areasum = new Accumulator(0);
112  Clear();
113  }
114 
118  public void Clear() {
119  _num = 0;
120  _crossings = 0;
121  _perimetersum.Set(0);
122  if (!_polyline) _areasum.Set(0);
123  _lat0 = _lon0 = _lat1 = _lon1 = Double.NaN;
124  }
125 
134  public void AddPoint(double lat, double lon) {
135  lon = GeoMath.AngNormalize(lon);
136  if (_num == 0) {
137  _lat0 = _lat1 = lat;
138  _lon0 = _lon1 = lon;
139  } else {
140  GeodesicData g = _earth.Inverse(_lat1, _lon1, lat, lon, _mask);
141  _perimetersum.Add(g.s12);
142  if (!_polyline) {
143  _areasum.Add(g.S12);
144  _crossings += transit(_lon1, lon);
145  }
146  _lat1 = lat; _lon1 = lon;
147  }
148  ++_num;
149  }
150 
160  public void AddEdge(double azi, double s) {
161  if (_num > 0) { // Do nothing if _num is zero
162  GeodesicData g = _earth.Direct(_lat1, _lon1, azi, s, _mask);
163  _perimetersum.Add(g.s12);
164  if (!_polyline) {
165  _areasum.Add(g.S12);
166  _crossings += transitdirect(_lon1, g.lon2);
167  }
168  _lat1 = g.lat2; _lon1 = g.lon2;
169  ++_num;
170  }
171  }
172 
184  public PolygonResult Compute() { return Compute(false, true); }
201  public PolygonResult Compute(boolean reverse, boolean sign) {
202  if (_num < 2)
203  return new PolygonResult(_num, 0, _polyline ? Double.NaN : 0);
204  if (_polyline)
205  return new PolygonResult(_num, _perimetersum.Sum(), Double.NaN);
206 
207  GeodesicData g = _earth.Inverse(_lat1, _lon1, _lat0, _lon0, _mask);
208  Accumulator tempsum = new Accumulator(_areasum);
209  tempsum.Add(g.S12);
210  int crossings = _crossings + transit(_lon1, _lon0);
211  if ((crossings & 1) != 0)
212  tempsum.Add((tempsum.Sum() < 0 ? 1 : -1) * _area0/2);
213  // area is with the clockwise sense. If !reverse convert to
214  // counter-clockwise convention.
215  if (!reverse)
216  tempsum.Negate();
217  // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
218  if (sign) {
219  if (tempsum.Sum() > _area0/2)
220  tempsum.Add(-_area0);
221  else if (tempsum.Sum() <= -_area0/2)
222  tempsum.Add(+_area0);
223  } else {
224  if (tempsum.Sum() >= _area0)
225  tempsum.Add(-_area0);
226  else if (tempsum.Sum() < 0)
227  tempsum.Add(+_area0);
228  }
229  return
230  new PolygonResult(_num, _perimetersum.Sum(g.s12), 0 + tempsum.Sum());
231  }
232 
256  public PolygonResult TestPoint(double lat, double lon,
257  boolean reverse, boolean sign) {
258  if (_num == 0)
259  return new PolygonResult(1, 0, _polyline ? Double.NaN : 0);
260 
261  double perimeter = _perimetersum.Sum();
262  double tempsum = _polyline ? 0 : _areasum.Sum();
263  int crossings = _crossings;
264  int num = _num + 1;
265  for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
266  GeodesicData g =
267  _earth.Inverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
268  i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
269  _mask);
270  perimeter += g.s12;
271  if (!_polyline) {
272  tempsum += g.S12;
273  crossings += transit(i == 0 ? _lon1 : lon,
274  i != 0 ? _lon0 : lon);
275  }
276  }
277 
278  if (_polyline)
279  return new PolygonResult(num, perimeter, Double.NaN);
280 
281  if ((crossings & 1) != 0)
282  tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
283  // area is with the clockwise sense. If !reverse convert to
284  // counter-clockwise convention.
285  if (!reverse)
286  tempsum *= -1;
287  // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
288  if (sign) {
289  if (tempsum > _area0/2)
290  tempsum -= _area0;
291  else if (tempsum <= -_area0/2)
292  tempsum += _area0;
293  } else {
294  if (tempsum >= _area0)
295  tempsum -= _area0;
296  else if (tempsum < 0)
297  tempsum += _area0;
298  }
299  return new PolygonResult(num, perimeter, 0 + tempsum);
300  }
301 
323  public PolygonResult TestEdge(double azi, double s,
324  boolean reverse, boolean sign) {
325  if (_num == 0) // we don't have a starting point!
326  return new PolygonResult(0, Double.NaN, Double.NaN);
327 
328  int num = _num + 1;
329  double perimeter = _perimetersum.Sum() + s;
330  if (_polyline)
331  return new PolygonResult(num, perimeter, Double.NaN);
332 
333  double tempsum = _areasum.Sum();
334  int crossings = _crossings;
335  {
336  double lat, lon, s12, S12, t;
337  GeodesicData g =
338  _earth.Direct(_lat1, _lon1, azi, false, s, _mask);
339  tempsum += g.S12;
340  crossings += transitdirect(_lon1, g.lon2);
341  g = _earth.Inverse(g.lat2, g.lon2, _lat0, _lon0, _mask);
342  perimeter += g.s12;
343  tempsum += g.S12;
344  crossings += transit(g.lon2, _lon0);
345  }
346 
347  if ((crossings & 1) != 0)
348  tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
349  // area is with the clockwise sense. If !reverse convert to
350  // counter-clockwise convention.
351  if (!reverse)
352  tempsum *= -1;
353  // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
354  if (sign) {
355  if (tempsum > _area0/2)
356  tempsum -= _area0;
357  else if (tempsum <= -_area0/2)
358  tempsum += _area0;
359  } else {
360  if (tempsum >= _area0)
361  tempsum -= _area0;
362  else if (tempsum < 0)
363  tempsum += _area0;
364  }
365 
366  return new PolygonResult(num, perimeter, 0 + tempsum);
367  }
368 
374  public double MajorRadius() { return _earth.MajorRadius(); }
375 
380  public double Flattening() { return _earth.Flattening(); }
381 
390  public Pair CurrentPoint() { return new Pair(_lat1, _lon1); }
391 }
static const double lat
static int transit(double lon1, double lon2)
PolygonResult TestPoint(double lat, double lon, boolean reverse, boolean sign)
void AddPoint(double lat, double lon)
void g(const string &key, int i)
Definition: testBTree.cpp:41
GeodesicData Direct(double lat1, double lon1, double azi1, double s12)
Definition: Geodesic.java:313
PolygonResult Compute(boolean reverse, boolean sign)
GeodesicData Inverse(double lat1, double lon1, double lat2, double lon2)
Definition: Geodesic.java:615
void AddEdge(double azi, double s)
Point3 cross(const Point3 &p, const Point3 &q, OptionalJacobian< 3, 3 > H1, OptionalJacobian< 3, 3 > H2)
cross product
Definition: Point3.cpp:64
PolygonResult TestEdge(double azi, double s, boolean reverse, boolean sign)
EIGEN_DEVICE_FUNC const SignReturnType sign() const
RealScalar s
void reverse(const MatrixType &m)
static int transitdirect(double lon1, double lon2)
static const double lon
static double AngNormalize(double x)
Definition: GeoMath.java:191
Point2 t(10, 10)
static Pair AngDiff(double x, double y)
Definition: GeoMath.java:221
PolygonArea(Geodesic earth, boolean polyline)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:14