PolygonArea.java

Go to the documentation of this file.

package net.sf.geographiclib;

public class PolygonArea {

  private Geodesic _earth;
  private double _area0;        // Full ellipsoid area
  private boolean _polyline;    // Assume polyline (don't close and skip area)
  private int _mask;
  private int _num;
  private int _crossings;
  private Accumulator _areasum, _perimetersum;
  private double _lat0, _lon0, _lat1, _lon1;
  private static int transit(double lon1, double lon2) {
    // Return 1 or -1 if crossing prime meridian in east or west direction.
    // Otherwise return zero.
    // Compute lon12 the same way as Geodesic.Inverse.
    lon1 = GeoMath.AngNormalize(lon1);
    lon2 = GeoMath.AngNormalize(lon2);
    double lon12 = GeoMath.AngDiff(lon1, lon2).first;
    int cross =
      lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
      (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
    return cross;
  }
  // an alternate version of transit to deal with longitudes in the direct
  // problem.
  private static int transitdirect(double lon1, double lon2) {
    // We want to compute exactly
    //   int(floor(lon2 / 360)) - int(floor(lon1 / 360))
    // Since we only need the parity of the result we can use std::remquo but
    // this is buggy with g++ 4.8.3 and requires C++11.  So instead we do
    lon1 = lon1 % 720.0; lon2 = lon2 % 720.0;
    return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
             ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
    }

  public PolygonArea(Geodesic earth, boolean polyline) {
    _earth = earth;
    _area0 = _earth.EllipsoidArea();
    _polyline = polyline;
     _mask = GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE |
       GeodesicMask.DISTANCE |
       (_polyline ? GeodesicMask.NONE :
        GeodesicMask.AREA | GeodesicMask.LONG_UNROLL);
     _perimetersum = new Accumulator(0);
     if (!_polyline)
       _areasum = new Accumulator(0);
     Clear();
  }

  public void Clear() {
    _num = 0;
    _crossings = 0;
    _perimetersum.Set(0);
    if (!_polyline) _areasum.Set(0);
    _lat0 = _lon0 = _lat1 = _lon1 = Double.NaN;
  }

  public void AddPoint(double lat, double lon) {
    lon = GeoMath.AngNormalize(lon);
    if (_num == 0) {
      _lat0 = _lat1 = lat;
      _lon0 = _lon1 = lon;
    } else {
      GeodesicData g = _earth.Inverse(_lat1, _lon1, lat, lon, _mask);
      _perimetersum.Add(g.s12);
      if (!_polyline) {
        _areasum.Add(g.S12);
        _crossings += transit(_lon1, lon);
      }
      _lat1 = lat; _lon1 = lon;
    }
    ++_num;
  }

  public void AddEdge(double azi, double s) {
    if (_num > 0) {             // Do nothing if _num is zero
      GeodesicData g = _earth.Direct(_lat1, _lon1, azi, s, _mask);
      _perimetersum.Add(g.s12);
      if (!_polyline) {
        _areasum.Add(g.S12);
        _crossings += transitdirect(_lon1, g.lon2);
      }
      _lat1 = g.lat2; _lon1 = g.lon2;
      ++_num;
    }
  }

  public PolygonResult Compute() { return Compute(false, true); }
  public PolygonResult Compute(boolean reverse, boolean sign) {
    if (_num < 2)
      return new PolygonResult(_num, 0, _polyline ? Double.NaN : 0);
    if (_polyline)
      return new PolygonResult(_num, _perimetersum.Sum(), Double.NaN);

    GeodesicData g = _earth.Inverse(_lat1, _lon1, _lat0, _lon0, _mask);
    Accumulator tempsum = new Accumulator(_areasum);
    tempsum.Add(g.S12);
    int crossings = _crossings + transit(_lon1, _lon0);
    if ((crossings & 1) != 0)
      tempsum.Add((tempsum.Sum() < 0 ? 1 : -1) * _area0/2);
    // area is with the clockwise sense.  If !reverse convert to
    // counter-clockwise convention.
    if (!reverse)
      tempsum.Negate();
    // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
    if (sign) {
      if (tempsum.Sum() > _area0/2)
        tempsum.Add(-_area0);
      else if (tempsum.Sum() <= -_area0/2)
        tempsum.Add(+_area0);
    } else {
      if (tempsum.Sum() >= _area0)
        tempsum.Add(-_area0);
      else if (tempsum.Sum() < 0)
        tempsum.Add(+_area0);
    }
    return
      new PolygonResult(_num, _perimetersum.Sum(g.s12), 0 + tempsum.Sum());
  }

  public PolygonResult TestPoint(double lat, double lon,
                                 boolean reverse, boolean sign) {
    if (_num == 0)
      return new PolygonResult(1, 0, _polyline ? Double.NaN : 0);

    double perimeter = _perimetersum.Sum();
    double tempsum = _polyline ? 0 : _areasum.Sum();
    int crossings = _crossings;
    int num = _num + 1;
    for (int i = 0; i < (_polyline ? 1 : 2); ++i) {
      GeodesicData g =
        _earth.Inverse(i == 0 ? _lat1 : lat, i == 0 ? _lon1 : lon,
                       i != 0 ? _lat0 : lat, i != 0 ? _lon0 : lon,
                       _mask);
      perimeter += g.s12;
      if (!_polyline) {
        tempsum += g.S12;
        crossings += transit(i == 0 ? _lon1 : lon,
                             i != 0 ? _lon0 : lon);
      }
    }

    if (_polyline)
      return new PolygonResult(num, perimeter, Double.NaN);

    if ((crossings & 1) != 0)
      tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
    // area is with the clockwise sense.  If !reverse convert to
    // counter-clockwise convention.
    if (!reverse)
      tempsum *= -1;
    // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
    if (sign) {
      if (tempsum > _area0/2)
        tempsum -= _area0;
      else if (tempsum <= -_area0/2)
        tempsum += _area0;
    } else {
      if (tempsum >= _area0)
        tempsum -= _area0;
      else if (tempsum < 0)
        tempsum += _area0;
    }
    return new PolygonResult(num, perimeter, 0 + tempsum);
  }

  public PolygonResult TestEdge(double azi, double s,
                                boolean reverse, boolean sign) {
    if (_num == 0)              // we don't have a starting point!
      return new PolygonResult(0, Double.NaN, Double.NaN);

    int num = _num + 1;
    double perimeter = _perimetersum.Sum() + s;
    if (_polyline)
      return new PolygonResult(num, perimeter, Double.NaN);

    double tempsum =  _areasum.Sum();
    int crossings = _crossings;
    {
      double lat, lon, s12, S12, t;
      GeodesicData g =
        _earth.Direct(_lat1, _lon1, azi, false, s, _mask);
      tempsum += g.S12;
      crossings += transitdirect(_lon1, g.lon2);
      g = _earth.Inverse(g.lat2, g.lon2, _lat0, _lon0, _mask);
      perimeter += g.s12;
      tempsum += g.S12;
      crossings += transit(g.lon2, _lon0);
    }

    if ((crossings & 1) != 0)
      tempsum += (tempsum < 0 ? 1 : -1) * _area0/2;
    // area is with the clockwise sense.  If !reverse convert to
    // counter-clockwise convention.
    if (!reverse)
      tempsum *= -1;
    // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
    if (sign) {
      if (tempsum > _area0/2)
        tempsum -= _area0;
      else if (tempsum <= -_area0/2)
        tempsum += _area0;
    } else {
      if (tempsum >= _area0)
        tempsum -= _area0;
      else if (tempsum < 0)
        tempsum += _area0;
    }

    return new PolygonResult(num, perimeter, 0 + tempsum);
  }

  public double MajorRadius() { return _earth.MajorRadius(); }

  public double Flattening() { return _earth.Flattening(); }

  public Pair CurrentPoint() { return new Pair(_lat1, _lon1); }
}