1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#ifndef UTIL_GEOMETRY_S2LOOP_H__
4#define UTIL_GEOMETRY_S2LOOP_H__
5
6#include <map>
7using std::map;
8using std::multimap;
9
10#include <vector>
11using std::vector;
12
13
14#include "base/logging.h"
15#include "base/macros.h"
16#include "s2edgeindex.h"
17#include "s2region.h"
18#include "s2latlngrect.h"
19#include "s2edgeutil.h"
20
21class S2Loop;
22// Defined in the cc file. A helper class for AreBoundariesCrossing.
23class WedgeProcessor;
24
25// Indexing structure to efficiently compute intersections.
26class S2LoopIndex: public S2EdgeIndex {
27 public:
28 explicit S2LoopIndex(S2Loop const* loop): loop_(loop) {}
29 virtual ~S2LoopIndex() {}
30
31 // There is no need to overwrite Reset(), as the only data that's
32 // used to implement this class is all contained in the loop data.
33 // void Reset();
34
35 virtual S2Point const* edge_from(int index) const;
36 virtual S2Point const* edge_to(int index) const;
37 virtual int num_edges() const;
38
39 private:
40 S2Loop const* loop_;
41};
42
43// An S2Loop represents a simple spherical polygon. It consists of a single
44// chain of vertices where the first vertex is implicitly connected to the
45// last. All loops are defined to have a CCW orientation, i.e. the interior
46// of the polygon is on the left side of the edges. This implies that a
47// clockwise loop enclosing a small area is interpreted to be a CCW loop
48// enclosing a very large area.
49//
50// Loops are not allowed to have any duplicate vertices (whether adjacent or
51// not), and non-adjacent edges are not allowed to intersect. Loops must have
52// at least 3 vertices. Although these restrictions are not enforced in
53// optimized code, you may get unexpected results if they are violated.
54//
55// Point containment is defined such that if the sphere is subdivided into
56// faces (loops), every point is contained by exactly one face. This implies
57// that loops do not necessarily contain all (or any) of their vertices.
58//
59// TODO(user): When doing operations on two loops, always create the
60// edgeindex for the bigger of the two. Same for polygons.
61class S2Loop : public S2Region {
62 public:
63 // Create an empty S2Loop that should be initialized by calling Init() or
64 // Decode().
65 S2Loop();
66
67 // Convenience constructor that calls Init() with the given vertices.
68 explicit S2Loop(vector<S2Point> const& vertices);
69
70 // Initialize a loop connecting the given vertices. The last vertex is
71 // implicitly connected to the first. All points should be unit length.
72 // Loops must have at least 3 vertices.
73 void Init(vector<S2Point> const& vertices);
74
75 // This parameter should be removed as soon as people stop using the
76 // deprecated version of IsValid.
77 static int const kDefaultMaxAdjacent = 0;
78
79 // Check whether this loop is valid. Note that in debug mode, validity
80 // is checked at loop creation time, so IsValid()
81 // should always return true.
82 bool IsValid() const;
83
84
85 // These two versions are deprecated and ignore max_adjacent.
86 // DEPRECATED.
87 static bool IsValid(vector<S2Point> const& vertices, int max_adjacent);
88 // DEPRECATED.
89 bool IsValid(int max_adjacent) const;
90
91 // Initialize a loop corresponding to the given cell.
92 explicit S2Loop(S2Cell const& cell);
93
94 ~S2Loop();
95
96 // The depth of a loop is defined as its nesting level within its containing
97 // polygon. "Outer shell" loops have depth 0, holes within those loops have
98 // depth 1, shells within those holes have depth 2, etc. This field is only
99 // used by the S2Polygon implementation.
100 int depth() const { return depth_; }
101 void set_depth(int depth) { depth_ = depth; }
102
103 // Return true if this loop represents a hole in its containing polygon.
104 bool is_hole() const { return (depth_ & 1) != 0; }
105
106 // The sign of a loop is -1 if the loop represents a hole in its containing
107 // polygon, and +1 otherwise.
108 int sign() const { return is_hole() ? -1 : 1; }
109
110 int num_vertices() const { return num_vertices_; }
111
112 // For convenience, we make two entire copies of the vertex list available:
113 // vertex(n..2*n-1) is mapped to vertex(0..n-1), where n == num_vertices().
114 S2Point const& vertex(int i) const {
115 DCHECK_GE(i, 0);
116 DCHECK_LT(i, (2 * num_vertices_));
117 int j = i - num_vertices();
118 return vertices_[j >= 0 ? j : i];
119 }
120
121 // Return true if the loop area is at most 2*Pi. Degenerate loops are
122 // handled consistently with S2::RobustCCW(), i.e., if a loop can be
123 // expressed as the union of degenerate or nearly-degenerate CCW triangles,
124 // then it will always be considered normalized.
125 bool IsNormalized() const;
126
127 // Invert the loop if necessary so that the area enclosed by the loop is at
128 // most 2*Pi.
129 void Normalize();
130
131 // Reverse the order of the loop vertices, effectively complementing
132 // the region represented by the loop.
133 void Invert();
134
135 // Return the area of the loop interior, i.e. the region on the left side of
136 // the loop. The return value is between 0 and 4*Pi. (Note that the return
137 // value is not affected by whether this loop is a "hole" or a "shell".)
138 double GetArea() const;
139
140 // Return the true centroid of the loop multiplied by the area of the loop
141 // (see s2.h for details on centroids). The result is not unit length, so
142 // you may want to normalize it. Also note that in general, the centroid
143 // may not be contained by the loop.
144 //
145 // We prescale by the loop area for two reasons: (1) it is cheaper to
146 // compute this way, and (2) it makes it easier to compute the centroid of
147 // more complicated shapes (by splitting them into disjoint regions and
148 // adding their centroids).
149 //
150 // Note that the return value is not affected by whether this loop is a
151 // "hole" or a "shell".
152 S2Point GetCentroid() const;
153
154 // Return the sum of the turning angles at each vertex. The return value is
155 // positive if the loop is counter-clockwise, negative if the loop is
156 // clockwise, and zero if the loop is a great circle. Degenerate and
157 // nearly-degenerate loops are handled consistently with S2::RobustCCW().
158 // So for example, if a loop has zero area (i.e., it is a very small CCW
159 // loop) then the turning angle will always be negative.
160 //
161 // This quantity is also called the "geodesic curvature" of the loop.
162 double GetTurningAngle() const;
163
164 // Return true if the region contained by this loop is a superset of the
165 // region contained by the given other loop.
166 bool Contains(S2Loop const* b) const;
167
168 // Return true if the region contained by this loop intersects the region
169 // contained by the given other loop.
170 bool Intersects(S2Loop const* b) const;
171
172 // Given two loops of a polygon (see s2polygon.h for requirements), return
173 // true if A contains B. This version of Contains() is much cheaper since
174 // it does not need to check whether the boundaries of the two loops cross.
175 bool ContainsNested(S2Loop const* b) const;
176
177 // Return +1 if A contains B (i.e. the interior of B is a subset of the
178 // interior of A), -1 if the boundaries of A and B cross, and 0 otherwise.
179 // Requires that A does not properly contain the complement of B, i.e.
180 // A and B do not contain each other's boundaries. This method is used
181 // for testing whether multi-loop polygons contain each other.
182 //
183 // WARNING: there is a bug in this function - it does not detect loop
184 // crossings in certain situations involving shared edges. CL 2926852 works
185 // around this bug for polygon intersection, but it probably effects other
186 // tests. TODO: fix ContainsOrCrosses() and revert CL 2926852.
187 int ContainsOrCrosses(S2Loop const* b) const;
188
189 // Return true if two loops have the same boundary. This is true if and
190 // only if the loops have the same vertices in the same cyclic order.
191 // (For testing purposes.)
192 bool BoundaryEquals(S2Loop const* b) const;
193
194 // Return true if two loops have the same boundary except for vertex
195 // perturbations. More precisely, the vertices in the two loops must be in
196 // the same cyclic order, and corresponding vertex pairs must be separated
197 // by no more than "max_error". (For testing purposes.)
198 bool BoundaryApproxEquals(S2Loop const* b, double max_error = 1e-15) const;
199
200 // Return true if the two loop boundaries are within "max_error" of each
201 // other along their entire lengths. The two loops may have different
202 // numbers of vertices. More precisely, this method returns true if the two
203 // loops have parameterizations a:[0,1] -> S^2, b:[0,1] -> S^2 such that
204 // distance(a(t), b(t)) <= max_error for all t. You can think of this as
205 // testing whether it is possible to drive two cars all the way around the
206 // two loops such that no car ever goes backward and the cars are always
207 // within "max_error" of each other. (For testing purposes.)
208 bool BoundaryNear(S2Loop const* b, double max_error = 1e-15) const;
209
210 // This method computes the oriented surface integral of some quantity f(x)
211 // over the loop interior, given a function f_tri(A,B,C) that returns the
212 // corresponding integral over the spherical triangle ABC. Here "oriented
213 // surface integral" means:
214 //
215 // (1) f_tri(A,B,C) must be the integral of f if ABC is counterclockwise,
216 // and the integral of -f if ABC is clockwise.
217 //
218 // (2) The result of this function is *either* the integral of f over the
219 // loop interior, or the integral of (-f) over the loop exterior.
220 //
221 // Note that there are at least two common situations where it easy to work
222 // around property (2) above:
223 //
224 // - If the integral of f over the entire sphere is zero, then it doesn't
225 // matter which case is returned because they are always equal.
226 //
227 // - If f is non-negative, then it is easy to detect when the integral over
228 // the loop exterior has been returned, and the integral over the loop
229 // interior can be obtained by adding the integral of f over the entire
230 // unit sphere (a constant) to the result.
231 //
232 // Also requires that the default constructor for T must initialize the
233 // value to zero. (This is true for built-in types such as "double".)
234 template <class T>
235 T GetSurfaceIntegral(T f_tri(S2Point const&, S2Point const&, S2Point const&))
236 const;
237
238 ////////////////////////////////////////////////////////////////////////
239 // S2Region interface (see s2region.h for details):
240
241 // GetRectBound() is guaranteed to return exact results, while GetCapBound()
242 // is conservative.
243 virtual S2Loop* Clone() const;
244 virtual S2Cap GetCapBound() const;
245 virtual S2LatLngRect GetRectBound() const { return bound_; }
246
247 virtual bool Contains(S2Cell const& cell) const;
248 virtual bool MayIntersect(S2Cell const& cell) const;
249 virtual bool VirtualContainsPoint(S2Point const& p) const {
250 return Contains(p); // The same as Contains() below, just virtual.
251 }
252
253 // The point 'p' does not need to be normalized.
254 bool Contains(S2Point const& p) const;
255
256 virtual void Encode(Encoder* const encoder) const;
257 virtual bool Decode(Decoder* const decoder);
258 virtual bool DecodeWithinScope(Decoder* const decoder);
259
260 private:
261 // Internal constructor used only by Clone() that makes a deep copy of
262 // its argument.
263 explicit S2Loop(S2Loop const* src);
264
265 void InitOrigin();
266 void InitBound();
267
268 // Internal implementation of the Decode and DecodeWithinScope methods above.
269 // If within_scope is true, memory is allocated for vertices_ and data
270 // is copied from the decoder using memcpy. If it is false, vertices_
271 // will point to the memory area inside the decoder, and the field
272 // owns_vertices_ is set to false.
273 bool DecodeInternal(Decoder* const decoder,
274 bool within_scope);
275
276 // Internal implementation of the Intersects() method above.
277 bool IntersectsInternal(S2Loop const* b) const;
278
279 // Return an index "first" and a direction "dir" (either +1 or -1) such that
280 // the vertex sequence (first, first+dir, ..., first+(n-1)*dir) does not
281 // change when the loop vertex order is rotated or inverted. This allows
282 // the loop vertices to be traversed in a canonical order. The return
283 // values are chosen such that (first, ..., first+n*dir) are in the range
284 // [0, 2*n-1] as expected by the vertex() method.
285 int GetCanonicalFirstVertex(int* dir) const;
286
287 // Return the index of a vertex at point "p", or -1 if not found.
288 // The return value is in the range 1..num_vertices_ if found.
289 int FindVertex(S2Point const& p) const;
290
291 // This method checks all edges of this loop (A) for intersection
292 // against all edges of B. If there is any shared vertex , the
293 // wedges centered at this vertex are sent to wedge_processor.
294 //
295 // Returns true only when the edges intersect in the sense of
296 // S2EdgeUtil::RobustCrossing, returns false immediately when the
297 // wedge_processor returns true: this means the wedge processor
298 // knows the value of the property that the caller wants to compute,
299 // and no further inspection is needed. For instance, if the
300 // property is "loops intersect", then a wedge intersection is all
301 // it takes to return true.
302 //
303 // See Contains(), Intersects() and ContainsOrCrosses() for the
304 // three uses of this function.
305 bool AreBoundariesCrossing(
306 S2Loop const* b, WedgeProcessor* wedge_processor) const;
307
308 // When the loop is modified (the only cae being Invert() at this time),
309 // the indexing structures need to be deleted as they become invalid.
310 void ResetMutableFields();
311
312 // We store the vertices in an array rather than a vector because we don't
313 // need any STL methods, and computing the number of vertices using size()
314 // would be relatively expensive (due to division by sizeof(S2Point) == 24).
315 // When DecodeWithinScope is used to initialize the loop, we do not
316 // take ownership of the memory for vertices_, and the owns_vertices_ field
317 // is used to prevent deallocation and overwriting.
318 int num_vertices_;
319 S2Point* vertices_;
320 bool owns_vertices_;
321
322 S2LatLngRect bound_;
323 bool origin_inside_;
324 int depth_;
325
326 // Quadtree index structure of this loop's edges.
327 mutable S2LoopIndex index_;
328
329 // Map for speeding up FindVertex: We will compute a map from vertex to
330 // index in the vertex array as soon as there has been enough calls.
331 mutable int num_find_vertex_calls_;
332 mutable map<S2Point, int> vertex_to_index_;
333
334 DISALLOW_EVIL_CONSTRUCTORS(S2Loop);
335};
336
337//////////////////// Implementation Details Follow ////////////////////////
338
339// Since this method is templatized and public, the implementation needs to be
340// in the .h file.
341
342template <class T>
343T S2Loop::GetSurfaceIntegral(T f_tri(S2Point const&, S2Point const&,
344 S2Point const&)) const {
345 // We sum "f_tri" over a collection T of oriented triangles, possibly
346 // overlapping. Let the sign of a triangle be +1 if it is CCW and -1
347 // otherwise, and let the sign of a point "x" be the sum of the signs of the
348 // triangles containing "x". Then the collection of triangles T is chosen
349 // such that either:
350 //
351 // (1) Each point in the loop interior has sign +1, and sign 0 otherwise; or
352 // (2) Each point in the loop exterior has sign -1, and sign 0 otherwise.
353 //
354 // The triangles basically consist of a "fan" from vertex 0 to every loop
355 // edge that does not include vertex 0. These triangles will always satisfy
356 // either (1) or (2). However, what makes this a bit tricky is that
357 // spherical edges become numerically unstable as their length approaches
358 // 180 degrees. Of course there is not much we can do if the loop itself
359 // contains such edges, but we would like to make sure that all the triangle
360 // edges under our control (i.e., the non-loop edges) are stable. For
361 // example, consider a loop around the equator consisting of four equally
362 // spaced points. This is a well-defined loop, but we cannot just split it
363 // into two triangles by connecting vertex 0 to vertex 2.
364 //
365 // We handle this type of situation by moving the origin of the triangle fan
366 // whenever we are about to create an unstable edge. We choose a new
367 // location for the origin such that all relevant edges are stable. We also
368 // create extra triangles with the appropriate orientation so that the sum
369 // of the triangle signs is still correct at every point.
370
371 // The maximum length of an edge for it to be considered numerically stable.
372 // The exact value is fairly arbitrary since it depends on the stability of
373 // the "f_tri" function. The value below is quite conservative but could be
374 // reduced further if desired.
375 static double const kMaxLength = M_PI - 1e-5;
376
377 // The default constructor for T must initialize the value to zero.
378 // (This is true for built-in types such as "double".)
379 T sum = T();
380 S2Point origin = vertex(0);
381 for (int i = 1; i + 1 < num_vertices(); ++i) {
382 // Let V_i be vertex(i), let O be the current origin, and let length(A,B)
383 // be the length of edge (A,B). At the start of each loop iteration, the
384 // "leading edge" of the triangle fan is (O,V_i), and we want to extend
385 // the triangle fan so that the leading edge is (O,V_i+1).
386 //
387 // Invariants:
388 // 1. length(O,V_i) < kMaxLength for all (i > 1).
389 // 2. Either O == V_0, or O is approximately perpendicular to V_0.
390 // 3. "sum" is the oriented integral of f over the area defined by
391 // (O, V_0, V_1, ..., V_i).
392 DCHECK(i == 1 || origin.Angle(vertex(i)) < kMaxLength);
393 DCHECK(origin == vertex(0) || fabs(origin.DotProd(vertex(0))) < 1e-15);
394
395 if (vertex(i+1).Angle(origin) > kMaxLength) {
396 // We are about to create an unstable edge, so choose a new origin O'
397 // for the triangle fan.
398 S2Point old_origin = origin;
399 if (origin == vertex(0)) {
400 // The following point is well-separated from V_i and V_0 (and
401 // therefore V_i+1 as well).
402 origin = S2::RobustCrossProd(vertex(0), vertex(i)).Normalize();
403 } else if (vertex(i).Angle(vertex(0)) < kMaxLength) {
404 // All edges of the triangle (O, V_0, V_i) are stable, so we can
405 // revert to using V_0 as the origin.
406 origin = vertex(0);
407 } else {
408 // (O, V_i+1) and (V_0, V_i) are antipodal pairs, and O and V_0 are
409 // perpendicular. Therefore V_0.CrossProd(O) is approximately
410 // perpendicular to all of {O, V_0, V_i, V_i+1}, and we can choose
411 // this point O' as the new origin.
412 origin = vertex(0).CrossProd(old_origin);
413
414 // Advance the edge (V_0,O) to (V_0,O').
415 sum += f_tri(vertex(0), old_origin, origin);
416 }
417 // Advance the edge (O,V_i) to (O',V_i).
418 sum += f_tri(old_origin, vertex(i), origin);
419 }
420 // Advance the edge (O,V_i) to (O,V_i+1).
421 sum += f_tri(origin, vertex(i), vertex(i+1));
422 }
423 // If the origin is not V_0, we need to sum one more triangle.
424 if (origin != vertex(0)) {
425 // Advance the edge (O,V_n-1) to (O,V_0).
426 sum += f_tri(origin, vertex(num_vertices() - 1), vertex(0));
427 }
428 return sum;
429}
430
431#endif // UTIL_GEOMETRY_S2LOOP_H__
432