1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#ifndef UTIL_GEOMETRY_S2EDGEUTIL_H__
4#define UTIL_GEOMETRY_S2EDGEUTIL_H__
5
6#include "base/logging.h"
7#include "base/macros.h"
8#include "s2.h"
9#include "s2latlngrect.h"
10
11class S2LatLngRect;
12
13// This class contains various utility functions related to edges. It
14// collects together common code that is needed to implement polygonal
15// geometry such as polylines, loops, and general polygons.
16class S2EdgeUtil {
17 public:
18 // This class allows a vertex chain v0, v1, v2, ... to be efficiently
19 // tested for intersection with a given fixed edge AB.
20 class EdgeCrosser {
21 public:
22 // AB is the given fixed edge, and C is the first vertex of the vertex
23 // chain. All parameters must point to fixed storage that persists for
24 // the lifetime of the EdgeCrosser object.
25 inline EdgeCrosser(S2Point const* a, S2Point const* b, S2Point const* c);
26
27 // Call this function when your chain 'jumps' to a new place.
28 inline void RestartAt(S2Point const* c);
29
30 // This method is equivalent to calling the S2EdgeUtil::RobustCrossing()
31 // function (defined below) on the edges AB and CD. It returns +1 if
32 // there is a crossing, -1 if there is no crossing, and 0 if two points
33 // from different edges are the same. Returns 0 or -1 if either edge is
34 // degenerate. As a side effect, it saves vertex D to be used as the next
35 // vertex C.
36 inline int RobustCrossing(S2Point const* d);
37
38 // This method is equivalent to the S2EdgeUtil::EdgeOrVertexCrossing()
39 // method defined below. It is similar to RobustCrossing, but handles
40 // cases where two vertices are identical in a way that makes it easy to
41 // implement point-in-polygon containment tests.
42 inline bool EdgeOrVertexCrossing(S2Point const* d);
43
44 private:
45 // This function handles the "slow path" of RobustCrossing(), which does
46 // not need to be inlined.
47 int RobustCrossingInternal(S2Point const* d);
48
49 // The fields below are all constant.
50 S2Point const* const a_;
51 S2Point const* const b_;
52 S2Point const a_cross_b_;
53
54 // The fields below are updated for each vertex in the chain.
55 S2Point const* c_; // Previous vertex in the vertex chain.
56 int acb_; // The orientation of the triangle ACB.
57 };
58
59 // This class computes a bounding rectangle that contains all edges
60 // defined by a vertex chain v0, v1, v2, ... All vertices must be unit
61 // length. Note that the bounding rectangle of an edge can be larger than
62 // the bounding rectangle of its endpoints, e.g. consider an edge that
63 // passes through the north pole.
64 class RectBounder {
65 public:
66 RectBounder() : bound_(S2LatLngRect::Empty()) {}
67
68 // This method is called to add each vertex to the chain. 'b'
69 // must point to fixed storage that persists through the next call
70 // to AddPoint. This means that if you don't store all of your
71 // points for the lifetime of the bounder, you must at least store
72 // the last two points and alternate which one you use for the
73 // next point.
74 void AddPoint(S2Point const* b);
75
76 // Return the bounding rectangle of the edge chain that connects the
77 // vertices defined so far.
78 S2LatLngRect GetBound() const { return bound_; }
79
80 private:
81 S2Point const* a_; // The previous vertex in the chain.
82 S2LatLng a_latlng_; // The corresponding latitude-longitude.
83 S2LatLngRect bound_; // The current bounding rectangle.
84 };
85
86 // The purpose of this class is to find edges that intersect a given
87 // longitude interval. It can be used as an efficient rejection test when
88 // attempting to find edges that intersect a given region. It accepts a
89 // vertex chain v0, v1, v2, ... and returns a boolean value indicating
90 // whether each edge intersects the specified longitude interval.
91 class LongitudePruner {
92 public:
93 // 'interval' is the longitude interval to be tested against, and
94 // 'v0' is the first vertex of edge chain.
95 LongitudePruner(S1Interval const& interval, S2Point const& v0);
96
97 // Returns true if the edge (v0, v1) intersects the given longitude
98 // interval, and then saves 'v1' to be used as the next 'v0'.
99 inline bool Intersects(S2Point const& v1);
100
101 private:
102 S1Interval interval_; // The interval to be tested against.
103 double lng0_; // The longitude of the next v0.
104 };
105
106 // Return true if edge AB crosses CD at a point that is interior
107 // to both edges. Properties:
108 //
109 // (1) SimpleCrossing(b,a,c,d) == SimpleCrossing(a,b,c,d)
110 // (2) SimpleCrossing(c,d,a,b) == SimpleCrossing(a,b,c,d)
111 static bool SimpleCrossing(S2Point const& a, S2Point const& b,
112 S2Point const& c, S2Point const& d);
113
114 // Like SimpleCrossing, except that points that lie exactly on a line are
115 // arbitrarily classified as being on one side or the other (according to
116 // the rules of S2::RobustCCW). It returns +1 if there is a crossing, -1
117 // if there is no crossing, and 0 if any two vertices from different edges
118 // are the same. Returns 0 or -1 if either edge is degenerate.
119 // Properties of RobustCrossing:
120 //
121 // (1) RobustCrossing(b,a,c,d) == RobustCrossing(a,b,c,d)
122 // (2) RobustCrossing(c,d,a,b) == RobustCrossing(a,b,c,d)
123 // (3) RobustCrossing(a,b,c,d) == 0 if a==c, a==d, b==c, b==d
124 // (3) RobustCrossing(a,b,c,d) <= 0 if a==b or c==d
125 //
126 // Note that if you want to check an edge against a *chain* of other
127 // edges, it is much more efficient to use an EdgeCrosser (above).
128 static int RobustCrossing(S2Point const& a, S2Point const& b,
129 S2Point const& c, S2Point const& d);
130
131 // Given two edges AB and CD where at least two vertices are identical
132 // (i.e. RobustCrossing(a,b,c,d) == 0), this function defines whether the
133 // two edges "cross" in a such a way that point-in-polygon containment tests
134 // can be implemented by counting the number of edge crossings. The basic
135 // rule is that a "crossing" occurs if AB is encountered after CD during a
136 // CCW sweep around the shared vertex starting from a fixed reference point.
137 //
138 // Note that according to this rule, if AB crosses CD then in general CD
139 // does not cross AB. However, this leads to the correct result when
140 // counting polygon edge crossings. For example, suppose that A,B,C are
141 // three consecutive vertices of a CCW polygon. If we now consider the edge
142 // crossings of a segment BP as P sweeps around B, the crossing number
143 // changes parity exactly when BP crosses BA or BC.
144 //
145 // Useful properties of VertexCrossing (VC):
146 //
147 // (1) VC(a,a,c,d) == VC(a,b,c,c) == false
148 // (2) VC(a,b,a,b) == VC(a,b,b,a) == true
149 // (3) VC(a,b,c,d) == VC(a,b,d,c) == VC(b,a,c,d) == VC(b,a,d,c)
150 // (3) If exactly one of a,b equals one of c,d, then exactly one of
151 // VC(a,b,c,d) and VC(c,d,a,b) is true
152 //
153 // It is an error to call this method with 4 distinct vertices.
154 static bool VertexCrossing(S2Point const& a, S2Point const& b,
155 S2Point const& c, S2Point const& d);
156
157 // A convenience function that calls RobustCrossing() to handle cases
158 // where all four vertices are distinct, and VertexCrossing() to handle
159 // cases where two or more vertices are the same. This defines a crossing
160 // function such that point-in-polygon containment tests can be implemented
161 // by simply counting edge crossings.
162 static bool EdgeOrVertexCrossing(S2Point const& a, S2Point const& b,
163 S2Point const& c, S2Point const& d);
164
165 // Given two edges AB and CD such that RobustCrossing() is true, return
166 // their intersection point. Useful properties of GetIntersection (GI):
167 //
168 // (1) GI(b,a,c,d) == GI(a,b,d,c) == GI(a,b,c,d)
169 // (2) GI(c,d,a,b) == GI(a,b,c,d)
170 //
171 // The returned intersection point X is guaranteed to be close to the edges
172 // AB and CD, but if the edges intersect at a very small angle then X may
173 // not be close to the true mathematical intersection point P. See the
174 // description of "kIntersectionTolerance" below for details.
175 static S2Point GetIntersection(S2Point const& a, S2Point const& b,
176 S2Point const& c, S2Point const& d);
177
178 // This distance is an upper bound on the distance from the intersection
179 // point returned by GetIntersection() to either of the two edges that were
180 // intersected. In particular, if "x" is the intersection point, then
181 // GetDistance(x, a, b) and GetDistance(x, c, d) will both be smaller than
182 // this value. The intersection tolerance is also large enough such if it
183 // is passed as the "vertex_merge_radius" of an S2PolygonBuilder, then the
184 // intersection point will be spliced into the edges AB and/or CD if they
185 // are also supplied to the S2PolygonBuilder.
186 static S1Angle const kIntersectionTolerance;
187
188 // Given a point X and an edge AB, return the distance ratio AX / (AX + BX).
189 // If X happens to be on the line segment AB, this is the fraction "t" such
190 // that X == Interpolate(A, B, t). Requires that A and B are distinct.
191 static double GetDistanceFraction(S2Point const& x,
192 S2Point const& a, S2Point const& b);
193
194 // Return the point X along the line segment AB whose distance from A is the
195 // given fraction "t" of the distance AB. Does NOT require that "t" be
196 // between 0 and 1. Note that all distances are measured on the surface of
197 // the sphere, so this is more complicated than just computing (1-t)*a + t*b
198 // and normalizing the result.
199 static S2Point Interpolate(double t, S2Point const& a, S2Point const& b);
200
201 // Like Interpolate(), except that the parameter "ax" represents the desired
202 // distance from A to the result X rather than a fraction between 0 and 1.
203 static S2Point InterpolateAtDistance(S1Angle const& ax,
204 S2Point const& a, S2Point const& b);
205
206 // A slightly more efficient version of InterpolateAtDistance() that can be
207 // used when the distance AB is already known.
208 static S2Point InterpolateAtDistance(S1Angle const& ax,
209 S2Point const& a, S2Point const& b,
210 S1Angle const& ab);
211
212 // Return the minimum distance from X to any point on the edge AB. All
213 // arguments should be unit length. The result is very accurate for small
214 // distances but may have some numerical error if the distance is large
215 // (approximately Pi/2 or greater). The case A == B is handled correctly.
216 static S1Angle GetDistance(S2Point const& x,
217 S2Point const& a, S2Point const& b);
218
219 // A slightly more efficient version of GetDistance() where the cross
220 // product of the two endpoints has been precomputed. The cross product
221 // does not need to be normalized, but should be computed using
222 // S2::RobustCrossProd() for the most accurate results.
223 static S1Angle GetDistance(S2Point const& x,
224 S2Point const& a, S2Point const& b,
225 S2Point const& a_cross_b);
226
227 // Return the point along the edge AB that is closest to the point X.
228 // The fractional distance of this point along the edge AB can be obtained
229 // using GetDistanceFraction() above.
230 static S2Point GetClosestPoint(S2Point const& x,
231 S2Point const& a, S2Point const& b);
232
233 // A slightly more efficient version of GetClosestPoint() where the cross
234 // product of the two endpoints has been precomputed. The cross product
235 // does not need to be normalized, but should be computed using
236 // S2::RobustCrossProd() for the most accurate results.
237 static S2Point GetClosestPoint(S2Point const& x,
238 S2Point const& a, S2Point const& b,
239 S2Point const& a_cross_b);
240
241 // Return true if every point on edge B=b0b1 is no further than "tolerance"
242 // from some point on edge A=a0a1.
243 // Requires that tolerance is less than 90 degrees.
244 static bool IsEdgeBNearEdgeA(S2Point const& a0, S2Point const& a1,
245 S2Point const& b0, S2Point const& b1,
246 S1Angle const& tolerance);
247
248 // For an edge chain (x0, x1, x2), a wedge is the region to the left
249 // of the edges. More precisely, it is the union of all the rays
250 // from x1x0 to x1x2, clockwise.
251 // The following are Wedge comparison functions for two wedges A =
252 // (a0, ab1, a2) and B = (b0, a12, b2). These are used in S2Loops.
253
254 // Returns true if wedge A fully contains or is equal to wedge B.
255 static bool WedgeContains(S2Point const& a0, S2Point const& ab1,
256 S2Point const& a2, S2Point const& b0,
257 S2Point const& b2);
258
259 // Returns true if the intersection of the two wedges is not empty.
260 static bool WedgeIntersects(S2Point const& a0, S2Point const& ab1,
261 S2Point const& a2, S2Point const& b0,
262 S2Point const& b2);
263
264 // Detailed relation from wedges A to wedge B.
265 enum WedgeRelation {
266 WEDGE_EQUALS,
267 WEDGE_PROPERLY_CONTAINS, // A is a strict superset of B.
268 WEDGE_IS_PROPERLY_CONTAINED, // A is a strict subset of B.
269 WEDGE_PROPERLY_OVERLAPS, // All of A intsect B, A-B and B-A are non-empty.
270 WEDGE_IS_DISJOINT, // A is disjoint from B
271 };
272
273 // Return the relation from wedge A to B.
274 static WedgeRelation GetWedgeRelation(
275 S2Point const& a0, S2Point const& ab1, S2Point const& a2,
276 S2Point const& b0, S2Point const& b2);
277
278 DISALLOW_IMPLICIT_CONSTRUCTORS(S2EdgeUtil); // Contains only static methods.
279};
280
281inline S2EdgeUtil::EdgeCrosser::EdgeCrosser(
282 S2Point const* a, S2Point const* b, S2Point const* c)
283 : a_(a), b_(b), a_cross_b_(a_->CrossProd(*b_)) {
284 RestartAt(c);
285}
286
287inline void S2EdgeUtil::EdgeCrosser::RestartAt(S2Point const* c) {
288 c_ = c;
289 acb_ = -S2::RobustCCW(*a_, *b_, *c_, a_cross_b_);
290}
291
292inline int S2EdgeUtil::EdgeCrosser::RobustCrossing(S2Point const* d) {
293 // For there to be an edge crossing, the triangles ACB, CBD, BDA, DAC must
294 // all be oriented the same way (CW or CCW). We keep the orientation of ACB
295 // as part of our state. When each new point D arrives, we compute the
296 // orientation of BDA and check whether it matches ACB. This checks whether
297 // the points C and D are on opposite sides of the great circle through AB.
298
299 // Recall that RobustCCW is invariant with respect to rotating its
300 // arguments, i.e. ABC has the same orientation as BDA.
301 int bda = S2::RobustCCW(*a_, *b_, *d, a_cross_b_);
302 int result;
303 if (bda == -acb_ && bda != 0) {
304 result = -1; // Most common case -- triangles have opposite orientations.
305 } else if ((bda & acb_) == 0) {
306 result = 0; // At least one value is zero -- two vertices are identical.
307 } else { // Slow path.
308 DCHECK_EQ(acb_, bda);
309 DCHECK_NE(0, bda);
310 result = RobustCrossingInternal(d);
311 }
312 // Now save the current vertex D as the next vertex C, and also save the
313 // orientation of the new triangle ACB (which is opposite to the current
314 // triangle BDA).
315 c_ = d;
316 acb_ = -bda;
317 return result;
318}
319
320inline bool S2EdgeUtil::EdgeCrosser::EdgeOrVertexCrossing(S2Point const* d) {
321 // We need to copy c_ since it is clobbered by RobustCrossing().
322 S2Point const* c = c_;
323 int crossing = RobustCrossing(d);
324 if (crossing < 0) return false;
325 if (crossing > 0) return true;
326 return VertexCrossing(*a_, *b_, *c, *d);
327}
328
329inline bool S2EdgeUtil::LongitudePruner::Intersects(S2Point const& v1) {
330 double lng1 = S2LatLng::Longitude(v1).radians();
331 bool result = interval_.Intersects(S1Interval::FromPointPair(lng0_, lng1));
332 lng0_ = lng1;
333 return result;
334}
335
336#endif // UTIL_GEOMETRY_S2EDGEUTIL_H__
337