s2edgeutil.h

// Copyright 2005 Google Inc. All Rights Reserved.

#ifndef UTIL_GEOMETRY_S2EDGEUTIL_H__
#define UTIL_GEOMETRY_S2EDGEUTIL_H__

#include "base/logging.h"
#include "base/macros.h"
#include "s2.h"
#include "s2latlngrect.h"

class S2LatLngRect;

// This class contains various utility functions related to edges.  It
// collects together common code that is needed to implement polygonal
// geometry such as polylines, loops, and general polygons.
class S2EdgeUtil {
 public:
  // This class allows a vertex chain v0, v1, v2, ... to be efficiently
  // tested for intersection with a given fixed edge AB.
  class EdgeCrosser {
   public:
    // AB is the given fixed edge, and C is the first vertex of the vertex
    // chain.  All parameters must point to fixed storage that persists for
    // the lifetime of the EdgeCrosser object.
    inline EdgeCrosser(S2Point const* a, S2Point const* b, S2Point const* c);

    // Call this function when your chain 'jumps' to a new place.
    inline void RestartAt(S2Point const* c);

    // This method is equivalent to calling the S2EdgeUtil::RobustCrossing()
    // function (defined below) on the edges AB and CD.  It returns +1 if
    // there is a crossing, -1 if there is no crossing, and 0 if two points
    // from different edges are the same.  Returns 0 or -1 if either edge is
    // degenerate.  As a side effect, it saves vertex D to be used as the next
    // vertex C.
    inline int RobustCrossing(S2Point const* d);

    // This method is equivalent to the S2EdgeUtil::EdgeOrVertexCrossing()
    // method defined below.  It is similar to RobustCrossing, but handles
    // cases where two vertices are identical in a way that makes it easy to
    // implement point-in-polygon containment tests.
    inline bool EdgeOrVertexCrossing(S2Point const* d);

   private:
    // This function handles the "slow path" of RobustCrossing(), which does
    // not need to be inlined.
    int RobustCrossingInternal(S2Point const* d);

    // The fields below are all constant.
    S2Point const* const a_;
    S2Point const* const b_;
    S2Point const a_cross_b_;

    // The fields below are updated for each vertex in the chain.
    S2Point const* c_;     // Previous vertex in the vertex chain.
    int acb_;              // The orientation of the triangle ACB.
  };

  // This class computes a bounding rectangle that contains all edges
  // defined by a vertex chain v0, v1, v2, ...  All vertices must be unit
  // length.  Note that the bounding rectangle of an edge can be larger than
  // the bounding rectangle of its endpoints, e.g. consider an edge that
  // passes through the north pole.
  class RectBounder {
   public:
    RectBounder() : bound_(S2LatLngRect::Empty()) {}

    // This method is called to add each vertex to the chain.  'b'
    // must point to fixed storage that persists through the next call
    // to AddPoint.  This means that if you don't store all of your
    // points for the lifetime of the bounder, you must at least store
    // the last two points and alternate which one you use for the
    // next point.
    void AddPoint(S2Point const* b);

    // Return the bounding rectangle of the edge chain that connects the
    // vertices defined so far.
    S2LatLngRect GetBound() const { return bound_; }

   private:
    S2Point const* a_;      // The previous vertex in the chain.
    S2LatLng a_latlng_;     // The corresponding latitude-longitude.
    S2LatLngRect bound_;    // The current bounding rectangle.
  };

  // The purpose of this class is to find edges that intersect a given
  // longitude interval.  It can be used as an efficient rejection test when
  // attempting to find edges that intersect a given region.  It accepts a
  // vertex chain v0, v1, v2, ...  and returns a boolean value indicating
  // whether each edge intersects the specified longitude interval.
  class LongitudePruner {
   public:
    // 'interval' is the longitude interval to be tested against, and
    // 'v0' is the first vertex of edge chain.
    LongitudePruner(S1Interval const& interval, S2Point const& v0);

    // Returns true if the edge (v0, v1) intersects the given longitude
    // interval, and then saves 'v1' to be used as the next 'v0'.
    inline bool Intersects(S2Point const& v1);

   private:
    S1Interval interval_;    // The interval to be tested against.
    double lng0_;            // The longitude of the next v0.
  };

  // Return true if edge AB crosses CD at a point that is interior
  // to both edges.  Properties:
  //
  //  (1) SimpleCrossing(b,a,c,d) == SimpleCrossing(a,b,c,d)
  //  (2) SimpleCrossing(c,d,a,b) == SimpleCrossing(a,b,c,d)
  static bool SimpleCrossing(S2Point const& a, S2Point const& b,
                             S2Point const& c, S2Point const& d);

  // Like SimpleCrossing, except that points that lie exactly on a line are
  // arbitrarily classified as being on one side or the other (according to
  // the rules of S2::RobustCCW).  It returns +1 if there is a crossing, -1
  // if there is no crossing, and 0 if any two vertices from different edges
  // are the same.  Returns 0 or -1 if either edge is degenerate.
  // Properties of RobustCrossing:
  //
  //  (1) RobustCrossing(b,a,c,d) == RobustCrossing(a,b,c,d)
  //  (2) RobustCrossing(c,d,a,b) == RobustCrossing(a,b,c,d)
  //  (3) RobustCrossing(a,b,c,d) == 0 if a==c, a==d, b==c, b==d
  //  (3) RobustCrossing(a,b,c,d) <= 0 if a==b or c==d
  //
  // Note that if you want to check an edge against a *chain* of other
  // edges, it is much more efficient to use an EdgeCrosser (above).
  static int RobustCrossing(S2Point const& a, S2Point const& b,
                            S2Point const& c, S2Point const& d);

  // Given two edges AB and CD where at least two vertices are identical
  // (i.e. RobustCrossing(a,b,c,d) == 0), this function defines whether the
  // two edges "cross" in a such a way that point-in-polygon containment tests
  // can be implemented by counting the number of edge crossings.  The basic
  // rule is that a "crossing" occurs if AB is encountered after CD during a
  // CCW sweep around the shared vertex starting from a fixed reference point.
  //
  // Note that according to this rule, if AB crosses CD then in general CD
  // does not cross AB.  However, this leads to the correct result when
  // counting polygon edge crossings.  For example, suppose that A,B,C are
  // three consecutive vertices of a CCW polygon.  If we now consider the edge
  // crossings of a segment BP as P sweeps around B, the crossing number
  // changes parity exactly when BP crosses BA or BC.
  //
  // Useful properties of VertexCrossing (VC):
  //
  //  (1) VC(a,a,c,d) == VC(a,b,c,c) == false
  //  (2) VC(a,b,a,b) == VC(a,b,b,a) == true
  //  (3) VC(a,b,c,d) == VC(a,b,d,c) == VC(b,a,c,d) == VC(b,a,d,c)
  //  (3) If exactly one of a,b equals one of c,d, then exactly one of
  //      VC(a,b,c,d) and VC(c,d,a,b) is true
  //
  // It is an error to call this method with 4 distinct vertices.
  static bool VertexCrossing(S2Point const& a, S2Point const& b,
                             S2Point const& c, S2Point const& d);

  // A convenience function that calls RobustCrossing() to handle cases
  // where all four vertices are distinct, and VertexCrossing() to handle
  // cases where two or more vertices are the same.  This defines a crossing
  // function such that point-in-polygon containment tests can be implemented
  // by simply counting edge crossings.
  static bool EdgeOrVertexCrossing(S2Point const& a, S2Point const& b,
                                   S2Point const& c, S2Point const& d);

  // Given two edges AB and CD such that RobustCrossing() is true, return
  // their intersection point.  Useful properties of GetIntersection (GI):
  //
  //  (1) GI(b,a,c,d) == GI(a,b,d,c) == GI(a,b,c,d)
  //  (2) GI(c,d,a,b) == GI(a,b,c,d)
  //
  // The returned intersection point X is guaranteed to be close to the edges
  // AB and CD, but if the edges intersect at a very small angle then X may
  // not be close to the true mathematical intersection point P.  See the
  // description of "kIntersectionTolerance" below for details.
  static S2Point GetIntersection(S2Point const& a, S2Point const& b,
                                 S2Point const& c, S2Point const& d);

  // This distance is an upper bound on the distance from the intersection
  // point returned by GetIntersection() to either of the two edges that were
  // intersected.  In particular, if "x" is the intersection point, then
  // GetDistance(x, a, b) and GetDistance(x, c, d) will both be smaller than
  // this value.  The intersection tolerance is also large enough such if it
  // is passed as the "vertex_merge_radius" of an S2PolygonBuilder, then the
  // intersection point will be spliced into the edges AB and/or CD if they
  // are also supplied to the S2PolygonBuilder.
  static S1Angle const kIntersectionTolerance;

  // Given a point X and an edge AB, return the distance ratio AX / (AX + BX).
  // If X happens to be on the line segment AB, this is the fraction "t" such
  // that X == Interpolate(A, B, t).  Requires that A and B are distinct.
  static double GetDistanceFraction(S2Point const& x,
                                    S2Point const& a, S2Point const& b);

  // Return the point X along the line segment AB whose distance from A is the
  // given fraction "t" of the distance AB.  Does NOT require that "t" be
  // between 0 and 1.  Note that all distances are measured on the surface of
  // the sphere, so this is more complicated than just computing (1-t)*a + t*b
  // and normalizing the result.
  static S2Point Interpolate(double t, S2Point const& a, S2Point const& b);

  // Like Interpolate(), except that the parameter "ax" represents the desired
  // distance from A to the result X rather than a fraction between 0 and 1.
  static S2Point InterpolateAtDistance(S1Angle const& ax,
                                       S2Point const& a, S2Point const& b);

  // A slightly more efficient version of InterpolateAtDistance() that can be
  // used when the distance AB is already known.
  static S2Point InterpolateAtDistance(S1Angle const& ax,
                                       S2Point const& a, S2Point const& b,
                                       S1Angle const& ab);

  // Return the minimum distance from X to any point on the edge AB.  All
  // arguments should be unit length.  The result is very accurate for small
  // distances but may have some numerical error if the distance is large
  // (approximately Pi/2 or greater).  The case A == B is handled correctly.
  static S1Angle GetDistance(S2Point const& x,
                             S2Point const& a, S2Point const& b);

  // A slightly more efficient version of GetDistance() where the cross
  // product of the two endpoints has been precomputed.  The cross product
  // does not need to be normalized, but should be computed using
  // S2::RobustCrossProd() for the most accurate results.
  static S1Angle GetDistance(S2Point const& x,
                             S2Point const& a, S2Point const& b,
                             S2Point const& a_cross_b);

  // Return the point along the edge AB that is closest to the point X.
  // The fractional distance of this point along the edge AB can be obtained
  // using GetDistanceFraction() above.
  static S2Point GetClosestPoint(S2Point const& x,
                                 S2Point const& a, S2Point const& b);

  // A slightly more efficient version of GetClosestPoint() where the cross
  // product of the two endpoints has been precomputed.  The cross product
  // does not need to be normalized, but should be computed using
  // S2::RobustCrossProd() for the most accurate results.
  static S2Point GetClosestPoint(S2Point const& x,
                                 S2Point const& a, S2Point const& b,
                                 S2Point const& a_cross_b);

  // Return true if every point on edge B=b0b1 is no further than "tolerance"
  // from some point on edge A=a0a1.
  // Requires that tolerance is less than 90 degrees.
  static bool IsEdgeBNearEdgeA(S2Point const& a0, S2Point const& a1,
                               S2Point const& b0, S2Point const& b1,
                               S1Angle const& tolerance);

  // For an edge chain (x0, x1, x2), a wedge is the region to the left
  // of the edges. More precisely, it is the union of all the rays
  // from x1x0 to x1x2, clockwise.
  // The following are Wedge comparison functions for two wedges A =
  // (a0, ab1, a2) and B = (b0, a12, b2). These are used in S2Loops.

  // Returns true if wedge A fully contains or is equal to wedge B.
  static bool WedgeContains(S2Point const& a0, S2Point const& ab1,
                            S2Point const& a2, S2Point const& b0,
                            S2Point const& b2);

  // Returns true if the intersection of the two wedges is not empty.
  static bool WedgeIntersects(S2Point const& a0, S2Point const& ab1,
                              S2Point const& a2, S2Point const& b0,
                              S2Point const& b2);

  // Detailed relation from wedges A to wedge B.
  enum WedgeRelation {
    WEDGE_EQUALS,
    WEDGE_PROPERLY_CONTAINS,  // A is a strict superset of B.
    WEDGE_IS_PROPERLY_CONTAINED,  // A is a strict subset of B.
    WEDGE_PROPERLY_OVERLAPS,  // All of A intsect B, A-B and B-A are non-empty.
    WEDGE_IS_DISJOINT,  // A is disjoint from B
  };

  // Return the relation from wedge A to B.
  static WedgeRelation GetWedgeRelation(
      S2Point const& a0, S2Point const& ab1, S2Point const& a2,
      S2Point const& b0, S2Point const& b2);

  DISALLOW_IMPLICIT_CONSTRUCTORS(S2EdgeUtil);  // Contains only static methods.
};

inline S2EdgeUtil::EdgeCrosser::EdgeCrosser(
    S2Point const* a, S2Point const* b, S2Point const* c)
    : a_(a), b_(b), a_cross_b_(a_->CrossProd(*b_)) {
  RestartAt(c);
}

inline void S2EdgeUtil::EdgeCrosser::RestartAt(S2Point const* c) {
  c_ = c;
  acb_ = -S2::RobustCCW(*a_, *b_, *c_, a_cross_b_);
}

inline int S2EdgeUtil::EdgeCrosser::RobustCrossing(S2Point const* d) {
  // For there to be an edge crossing, the triangles ACB, CBD, BDA, DAC must
  // all be oriented the same way (CW or CCW).  We keep the orientation of ACB
  // as part of our state.  When each new point D arrives, we compute the
  // orientation of BDA and check whether it matches ACB.  This checks whether
  // the points C and D are on opposite sides of the great circle through AB.

  // Recall that RobustCCW is invariant with respect to rotating its
  // arguments, i.e. ABC has the same orientation as BDA.
  int bda = S2::RobustCCW(*a_, *b_, *d, a_cross_b_);
  int result;
  if (bda == -acb_ && bda != 0) {
    result = -1;  // Most common case -- triangles have opposite orientations.
  } else if ((bda & acb_) == 0) {
    result = 0;  // At least one value is zero -- two vertices are identical.
  } else {  // Slow path.
    DCHECK_EQ(acb_, bda);
    DCHECK_NE(0, bda);
    result = RobustCrossingInternal(d);
  }
  // Now save the current vertex D as the next vertex C, and also save the
  // orientation of the new triangle ACB (which is opposite to the current
  // triangle BDA).
  c_ = d;
  acb_ = -bda;
  return result;
}

inline bool S2EdgeUtil::EdgeCrosser::EdgeOrVertexCrossing(S2Point const* d) {
  // We need to copy c_ since it is clobbered by RobustCrossing().
  S2Point const* c = c_;
  int crossing = RobustCrossing(d);
  if (crossing < 0) return false;
  if (crossing > 0) return true;
  return VertexCrossing(*a_, *b_, *c, *d);
}

inline bool S2EdgeUtil::LongitudePruner::Intersects(S2Point const& v1) {
  double lng1 = S2LatLng::Longitude(v1).radians();
  bool result = interval_.Intersects(S1Interval::FromPointPair(lng0_, lng1));
  lng0_ = lng1;
  return result;
}

#endif  // UTIL_GEOMETRY_S2EDGEUTIL_H__