1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#ifndef UTIL_GEOMETRY_S2_H_
4#define UTIL_GEOMETRY_S2_H_
5
6#include <algorithm>
7using std::min;
8using std::max;
9using std::swap;
10using std::reverse;
11
12#include <hash_map>
13using __gnu_cxx::hash_map;
14 // To have template struct hash<T> defined
15#include "base/basictypes.h"
16#include "base/logging.h"
17#include "base/macros.h"
18#include "base/port.h" // for HASH_NAMESPACE_DECLARATION_START
19#include "util/math/vector3-inl.h"
20#include "util/math/matrix3x3.h"
21
22// An S2Point represents a point on the unit sphere as a 3D vector. Usually
23// points are normalized to be unit length, but some methods do not require
24// this. See util/math/vector3-inl.h for the methods available. Among other
25// things, there are overloaded operators that make it convenient to write
26// arithmetic expressions (e.g. (1-x)*p1 + x*p2).
27typedef Vector3_d S2Point;
28
29#include<hash_set>
30namespace __gnu_cxx {
31
32
33template<> struct hash<S2Point> {
34 size_t operator()(S2Point const& p) const;
35};
36
37
38} // namespace __gnu_cxx
39
40
41// The S2 class is simply a namespace for constants and static utility
42// functions related to spherical geometry, such as area calculations and edge
43// intersection tests. The name "S2" is derived from the mathematical symbol
44// for the two-dimensional unit sphere (note that the "2" refers to the
45// dimension of the surface, not the space it is embedded in).
46//
47// This class also defines a framework for decomposing the unit sphere into a
48// hierarchy of "cells". Each cell is a quadrilateral bounded by four
49// geodesics. The top level of the hierarchy is obtained by projecting the
50// six faces of a cube onto the unit sphere, and lower levels are obtained by
51// subdividing each cell into four children recursively.
52//
53// This class specifies the details of how the cube faces are projected onto
54// the unit sphere. This includes getting the face ordering and orientation
55// correct so that sequentially increasing cell ids follow a continuous
56// space-filling curve over the entire sphere, and defining the
57// transformation from cell-space to cube-space in order to make the cells
58// more uniform in size.
59//
60// This file also contains documentation of the various coordinate systems
61// and conventions used.
62//
63// This class is not thread-safe for loops and objects that use loops.
64//
65class S2 {
66 public:
67 static const bool debug;
68 // Return a unique "origin" on the sphere for operations that need a fixed
69 // reference point. In particular, this is the "point at infinity" used for
70 // point-in-polygon testing (by counting the number of edge crossings).
71 //
72 // It should *not* be a point that is commonly used in edge tests in order
73 // to avoid triggering code to handle degenerate cases. (This rules out the
74 // north and south poles.) It should also not be on the boundary of any
75 // low-level S2Cell for the same reason.
76 inline static S2Point Origin();
77
78 // Return true if the given point is approximately unit length
79 // (this is mainly useful for assertions).
80 static bool IsUnitLength(S2Point const& p);
81
82 // Return a unit-length vector that is orthogonal to "a". Satisfies
83 // Ortho(-a) = -Ortho(a) for all a.
84 static S2Point Ortho(S2Point const& a);
85
86 // Given a point "z" on the unit sphere, extend this into a right-handed
87 // coordinate frame of unit-length column vectors m = (x,y,z). Note that
88 // the vectors (x,y) are an orthonormal frame for the tangent space at "z",
89 // while "z" itself is an orthonormal frame for the normal space at "z".
90 static void GetFrame(S2Point const& z, Matrix3x3_d* m);
91
92 // Given an orthonormal basis "m" of column vectors and a point "p", return
93 // the coordinates of "p" with respect to the basis "m". The resulting
94 // point "q" satisfies the identity (m * q == p).
95 static S2Point ToFrame(Matrix3x3_d const& m, S2Point const& p);
96
97 // Given an orthonormal basis "m" of column vectors and a point "q" with
98 // respect to that basis, return the equivalent point "p" with respect to
99 // the standard axis-aligned basis. The result satisfies (p == m * q).
100 static S2Point FromFrame(Matrix3x3_d const& m, S2Point const& q);
101
102 // the coordinates of "p" with respect to the basis "m". The resulting
103 // point "r" satisfies the identity (m * r == p).
104
105 // Return true if two points are within the given distance of each other
106 // (this is mainly useful for testing).
107 static bool ApproxEquals(S2Point const& a, S2Point const& b,
108 double max_error = 1e-15);
109
110 // Return a vector "c" that is orthogonal to the given unit-length vectors
111 // "a" and "b". This function is similar to a.CrossProd(b) except that it
112 // does a better job of ensuring orthogonality when "a" is nearly parallel
113 // to "b", and it returns a non-zero result even when a == b or a == -b.
114 //
115 // It satisfies the following properties (RCP == RobustCrossProd):
116 //
117 // (1) RCP(a,b) != 0 for all a, b
118 // (2) RCP(b,a) == -RCP(a,b) unless a == b or a == -b
119 // (3) RCP(-a,b) == -RCP(a,b) unless a == b or a == -b
120 // (4) RCP(a,-b) == -RCP(a,b) unless a == b or a == -b
121 static S2Point RobustCrossProd(S2Point const& a, S2Point const& b);
122
123 // Return true if the points A, B, C are strictly counterclockwise. Return
124 // false if the points are clockwise or collinear (i.e. if they are all
125 // contained on some great circle).
126 //
127 // Due to numerical errors, situations may arise that are mathematically
128 // impossible, e.g. ABC may be considered strictly CCW while BCA is not.
129 // However, the implementation guarantees the following:
130 //
131 // If SimpleCCW(a,b,c), then !SimpleCCW(c,b,a) for all a,b,c.
132 static bool SimpleCCW(S2Point const& a, S2Point const& b, S2Point const& c);
133
134 // Returns +1 if the points A, B, C are counterclockwise, -1 if the points
135 // are clockwise, and 0 if any two points are the same. This function is
136 // essentially like taking the sign of the determinant of ABC, except that
137 // it has additional logic to make sure that the above properties hold even
138 // when the three points are coplanar, and to deal with the limitations of
139 // floating-point arithmetic.
140 //
141 // RobustCCW satisfies the following conditions:
142 //
143 // (1) RobustCCW(a,b,c) == 0 if and only if a == b, b == c, or c == a
144 // (2) RobustCCW(b,c,a) == RobustCCW(a,b,c) for all a,b,c
145 // (3) RobustCCW(c,b,a) == -RobustCCW(a,b,c) for all a,b,c
146 //
147 // In other words:
148 //
149 // (1) The result is zero if and only if two points are the same.
150 // (2) Rotating the order of the arguments does not affect the result.
151 // (3) Exchanging any two arguments inverts the result.
152 //
153 // On the other hand, note that it is not true in general that
154 // RobustCCW(-a,b,c) == -RobustCCW(a,b,c), or any similar identities
155 // involving antipodal points.
156 static int RobustCCW(S2Point const& a, S2Point const& b, S2Point const& c);
157
158 // A more efficient version of RobustCCW that allows the precomputed
159 // cross-product of A and B to be specified. (Unlike the 3 argument
160 // version this method is also inlined.)
161 inline static int RobustCCW(S2Point const& a, S2Point const& b,
162 S2Point const& c, S2Point const& a_cross_b);
163
164 // This version of RobustCCW returns +1 if the points are definitely CCW,
165 // -1 if they are definitely CW, and 0 if two points are identical or the
166 // result is uncertain. Uncertain certain cases can be resolved, if
167 // desired, by calling ExpensiveCCW.
168 //
169 // The purpose of this method is to allow additional cheap tests to be done,
170 // where possible, in order to avoid calling ExpensiveCCW unnecessarily.
171 inline static int TriageCCW(S2Point const& a, S2Point const& b,
172 S2Point const& c, S2Point const& a_cross_b);
173
174 // This function is invoked by RobustCCW() if the sign of the determinant is
175 // uncertain. It always returns a non-zero result unless two of the input
176 // points are the same. It uses a combination of multiple-precision
177 // arithmetic and symbolic perturbations to ensure that its results are
178 // always self-consistent (cf. Simulation of Simplicity, Edelsbrunner and
179 // Muecke). The basic idea is to assign an infinitesmal symbolic
180 // perturbation to every possible S2Point such that no three S2Points are
181 // collinear and no four S2Points are coplanar. These perturbations are so
182 // small that they do not affect the sign of any determinant that was
183 // non-zero before the perturbations.
184 //
185 // Unlike RobustCCW(), this method does not require the input points to be
186 // normalized.
187 static int ExpensiveCCW(S2Point const& a, S2Point const& b,
188 S2Point const& c);
189
190 // Given 4 points on the unit sphere, return true if the edges OA, OB, and
191 // OC are encountered in that order while sweeping CCW around the point O.
192 // You can think of this as testing whether A <= B <= C with respect to the
193 // CCW ordering around O that starts at A, or equivalently, whether B is
194 // contained in the range of angles (inclusive) that starts at A and extends
195 // CCW to C. Properties:
196 //
197 // (1) If OrderedCCW(a,b,c,o) && OrderedCCW(b,a,c,o), then a == b
198 // (2) If OrderedCCW(a,b,c,o) && OrderedCCW(a,c,b,o), then b == c
199 // (3) If OrderedCCW(a,b,c,o) && OrderedCCW(c,b,a,o), then a == b == c
200 // (4) If a == b or b == c, then OrderedCCW(a,b,c,o) is true
201 // (5) Otherwise if a == c, then OrderedCCW(a,b,c,o) is false
202 static bool OrderedCCW(S2Point const& a, S2Point const& b, S2Point const& c,
203 S2Point const& o);
204
205 // Return the interior angle at the vertex B in the triangle ABC. The
206 // return value is always in the range [0, Pi]. The points do not need to
207 // be normalized. Ensures that Angle(a,b,c) == Angle(c,b,a) for all a,b,c.
208 //
209 // The angle is undefined if A or C is diametrically opposite from B, and
210 // becomes numerically unstable as the length of edge AB or BC approaches
211 // 180 degrees.
212 static double Angle(S2Point const& a, S2Point const& b, S2Point const& c);
213
214 // Return the exterior angle at the vertex B in the triangle ABC. The
215 // return value is positive if ABC is counterclockwise and negative
216 // otherwise. If you imagine an ant walking from A to B to C, this is the
217 // angle that the ant turns at vertex B (positive = left, negative = right).
218 // Ensures that TurnAngle(a,b,c) == -TurnAngle(c,b,a) for all a,b,c.
219 static double TurnAngle(S2Point const& a, S2Point const& b, S2Point const& c);
220
221 // Return the area of triangle ABC. The method used is about twice as
222 // expensive as Girard's formula, but it is numerically stable for both
223 // large and very small triangles. All points should be unit length.
224 // The area is always positive.
225 //
226 // The triangle area is undefined if it contains two antipodal points, and
227 // becomes numerically unstable as the length of any edge approaches 180
228 // degrees.
229 static double Area(S2Point const& a, S2Point const& b, S2Point const& c);
230
231 // Return the area of the triangle computed using Girard's formula. All
232 // points should be unit length. This is slightly faster than the Area()
233 // method above but is not accurate for very small triangles.
234 static double GirardArea(S2Point const& a, S2Point const& b,
235 S2Point const& c);
236
237 // Like Area(), but returns a positive value for counterclockwise triangles
238 // and a negative value otherwise.
239 static double SignedArea(S2Point const& a, S2Point const& b,
240 S2Point const& c);
241
242 // About centroids:
243 // ----------------
244 //
245 // There are several notions of the "centroid" of a triangle. First, there
246 // // is the planar centroid, which is simply the centroid of the ordinary
247 // (non-spherical) triangle defined by the three vertices. Second, there is
248 // the surface centroid, which is defined as the intersection of the three
249 // medians of the spherical triangle. It is possible to show that this
250 // point is simply the planar centroid projected to the surface of the
251 // sphere. Finally, there is the true centroid (mass centroid), which is
252 // defined as the area integral over the spherical triangle of (x,y,z)
253 // divided by the triangle area. This is the point that the triangle would
254 // rotate around if it was spinning in empty space.
255 //
256 // The best centroid for most purposes is the true centroid. Unlike the
257 // planar and surface centroids, the true centroid behaves linearly as
258 // regions are added or subtracted. That is, if you split a triangle into
259 // pieces and compute the average of their centroids (weighted by triangle
260 // area), the result equals the centroid of the original triangle. This is
261 // not true of the other centroids.
262 //
263 // Also note that the surface centroid may be nowhere near the intuitive
264 // "center" of a spherical triangle. For example, consider the triangle
265 // with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere).
266 // The surface centroid of this triangle is at S=(0, 2*eps, 1), which is
267 // within a distance of 2*eps of the vertex B. Note that the median from A
268 // (the segment connecting A to the midpoint of BC) passes through S, since
269 // this is the shortest path connecting the two endpoints. On the other
270 // hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto
271 // the surface is a much more reasonable interpretation of the "center" of
272 // this triangle.
273
274 // Return the centroid of the planar triangle ABC. This can be normalized
275 // to unit length to obtain the "surface centroid" of the corresponding
276 // spherical triangle, i.e. the intersection of the three medians. However,
277 // note that for large spherical triangles the surface centroid may be
278 // nowhere near the intuitive "center" (see example above).
279 static S2Point PlanarCentroid(S2Point const& a, S2Point const& b,
280 S2Point const& c);
281
282 // Returns the true centroid of the spherical triangle ABC multiplied by the
283 // signed area of spherical triangle ABC. The reasons for multiplying by
284 // the signed area are (1) this is the quantity that needs to be summed to
285 // compute the centroid of a union or difference of triangles, and (2) it's
286 // actually easier to calculate this way.
287 static S2Point TrueCentroid(S2Point const& a, S2Point const& b,
288 S2Point const& c);
289
290 ////////////////////////// S2Cell Decomposition /////////////////////////
291 //
292 // The following methods define the cube-to-sphere projection used by
293 // the S2Cell decomposition.
294 //
295 // In the process of converting a latitude-longitude pair to a 64-bit cell
296 // id, the following coordinate systems are used:
297 //
298 // (id)
299 // An S2CellId is a 64-bit encoding of a face and a Hilbert curve position
300 // on that face. The Hilbert curve position implicitly encodes both the
301 // position of a cell and its subdivision level (see s2cellid.h).
302 //
303 // (face, i, j)
304 // Leaf-cell coordinates. "i" and "j" are integers in the range
305 // [0,(2**30)-1] that identify a particular leaf cell on the given face.
306 // The (i, j) coordinate system is right-handed on each face, and the
307 // faces are oriented such that Hilbert curves connect continuously from
308 // one face to the next.
309 //
310 // (face, s, t)
311 // Cell-space coordinates. "s" and "t" are real numbers in the range
312 // [0,1] that identify a point on the given face. For example, the point
313 // (s, t) = (0.5, 0.5) corresponds to the center of the top-level face
314 // cell. This point is also a vertex of exactly four cells at each
315 // subdivision level greater than zero.
316 //
317 // (face, si, ti)
318 // Discrete cell-space coordinates. These are obtained by multiplying
319 // "s" and "t" by 2**31 and rounding to the nearest unsigned integer.
320 // Discrete coordinates lie in the range [0,2**31]. This coordinate
321 // system can represent the edge and center positions of all cells with
322 // no loss of precision (including non-leaf cells).
323 //
324 // (face, u, v)
325 // Cube-space coordinates. To make the cells at each level more uniform
326 // in size after they are projected onto the sphere, we apply apply a
327 // nonlinear transformation of the form u=f(s), v=f(t). The (u, v)
328 // coordinates after this transformation give the actual coordinates on
329 // the cube face (modulo some 90 degree rotations) before it is projected
330 // onto the unit sphere.
331 //
332 // (x, y, z)
333 // Direction vector (S2Point). Direction vectors are not necessarily unit
334 // length, and are often chosen to be points on the biunit cube
335 // [-1,+1]x[-1,+1]x[-1,+1]. They can be be normalized to obtain the
336 // corresponding point on the unit sphere.
337 //
338 // (lat, lng)
339 // Latitude and longitude (S2LatLng). Latitudes must be between -90 and
340 // 90 degrees inclusive, and longitudes must be between -180 and 180
341 // degrees inclusive.
342 //
343 // Note that the (i, j), (s, t), (si, ti), and (u, v) coordinate systems are
344 // right-handed on all six faces.
345
346 // This is the number of levels needed to specify a leaf cell. This
347 // constant is defined here so that the S2::Metric class can be
348 // implemented without including s2cellid.h.
349 static int const kMaxCellLevel = 30;
350
351 // Convert an s or t value to the corresponding u or v value. This is
352 // a non-linear transformation from [-1,1] to [-1,1] that attempts to
353 // make the cell sizes more uniform.
354 inline static double STtoUV(double s);
355
356 // The inverse of the STtoUV transformation. Note that it is not always
357 // true that UVtoST(STtoUV(x)) == x due to numerical errors.
358 inline static double UVtoST(double u);
359
360 // Convert (face, u, v) coordinates to a direction vector (not
361 // necessarily unit length).
362 inline static S2Point FaceUVtoXYZ(int face, double u, double v);
363
364 // If the dot product of p with the given face normal is positive,
365 // set the corresponding u and v values (which may lie outside the range
366 // [-1,1]) and return true. Otherwise return false.
367 inline static bool FaceXYZtoUV(int face, S2Point const& p,
368 double* pu, double* pv);
369
370 // Convert a direction vector (not necessarily unit length) to
371 // (face, u, v) coordinates.
372 inline static int XYZtoFaceUV(S2Point const& p, double* pu, double* pv);
373
374 // Return the right-handed normal (not necessarily unit length) for an
375 // edge in the direction of the positive v-axis at the given u-value on
376 // the given face. (This vector is perpendicular to the plane through
377 // the sphere origin that contains the given edge.)
378 inline static S2Point GetUNorm(int face, double u);
379
380 // Return the right-handed normal (not necessarily unit length) for an
381 // edge in the direction of the positive u-axis at the given v-value on
382 // the given face.
383 inline static S2Point GetVNorm(int face, double v);
384
385 // Return the unit-length normal, u-axis, or v-axis for the given face.
386 inline static S2Point GetNorm(int face);
387 inline static S2Point GetUAxis(int face);
388 inline static S2Point GetVAxis(int face);
389
390 ////////////////////////////////////////////////////////////////////////
391 // The canonical Hilbert traversal order looks like an inverted 'U':
392 // the subcells are visited in the order (0,0), (0,1), (1,1), (1,0).
393 // The following tables encode the traversal order for various
394 // orientations of the Hilbert curve (axes swapped and/or directions
395 // of the axes reversed).
396
397 // Together these flags define a cell orientation. If 'kSwapMask'
398 // is true, then canonical traversal order is flipped around the
399 // diagonal (i.e. i and j are swapped with each other). If
400 // 'kInvertMask' is true, then the traversal order is rotated by 180
401 // degrees (i.e. the bits of i and j are inverted, or equivalently,
402 // the axis directions are reversed).
403 static int const kSwapMask = 0x01;
404 static int const kInvertMask = 0x02;
405
406 // kIJtoPos[orientation][ij] -> pos
407 //
408 // Given a cell orientation and the (i,j)-index of a subcell (0=(0,0),
409 // 1=(0,1), 2=(1,0), 3=(1,1)), return the order in which this subcell is
410 // visited by the Hilbert curve (a position in the range [0..3]).
411 static int const kIJtoPos[4][4];
412
413 // kPosToIJ[orientation][pos] -> ij
414 //
415 // Return the (i,j) index of the subcell at the given position 'pos' in the
416 // Hilbert curve traversal order with the given orientation. This is the
417 // inverse of the previous table:
418 //
419 // kPosToIJ[r][kIJtoPos[r][ij]] == ij
420 static int const kPosToIJ[4][4];
421
422 // kPosToOrientation[pos] -> orientation_modifier
423 //
424 // Return a modifier indicating how the orientation of the child subcell
425 // with the given traversal position [0..3] is related to the orientation
426 // of the parent cell. The modifier should be XOR-ed with the parent
427 // orientation to obtain the curve orientation in the child.
428 static int const kPosToOrientation[4];
429
430 ////////////////////////// S2Cell Metrics //////////////////////////////
431 //
432 // The following are various constants that describe the shapes and sizes of
433 // cells. They are useful for deciding which cell level to use in order to
434 // satisfy a given condition (e.g. that cell vertices must be no further
435 // than "x" apart). All of the raw constants are differential quantities;
436 // you can use the GetValue(level) method to compute the corresponding length
437 // or area on the unit sphere for cells at a given level. The minimum and
438 // maximum bounds are valid for cells at all levels, but they may be
439 // somewhat conservative for very large cells (e.g. face cells).
440
441 // Defines a cell metric of the given dimension (1 == length, 2 == area).
442 template <int dim> class Metric {
443 public:
444 explicit Metric(double deriv) : deriv_(deriv) {}
445
446 // The "deriv" value of a metric is a derivative, and must be multiplied by
447 // a length or area in (s,t)-space to get a useful value.
448 double deriv() const { return deriv_; }
449
450 // Return the value of a metric for cells at the given level. The value is
451 // either a length or an area on the unit sphere, depending on the
452 // particular metric.
453 double GetValue(int level) const { return ldexp(deriv_, - dim * level); }
454
455 // Return the level at which the metric has approximately the given
456 // value. For example, S2::kAvgEdge.GetClosestLevel(0.1) returns the
457 // level at which the average cell edge length is approximately 0.1.
458 // The return value is always a valid level.
459 int GetClosestLevel(double value) const;
460
461 // Return the minimum level such that the metric is at most the given
462 // value, or S2CellId::kMaxLevel if there is no such level. For example,
463 // S2::kMaxDiag.GetMinLevel(0.1) returns the minimum level such that all
464 // cell diagonal lengths are 0.1 or smaller. The return value is always a
465 // valid level.
466 int GetMinLevel(double value) const;
467
468 // Return the maximum level such that the metric is at least the given
469 // value, or zero if there is no such level. For example,
470 // S2::kMinWidth.GetMaxLevel(0.1) returns the maximum level such that all
471 // cells have a minimum width of 0.1 or larger. The return value is
472 // always a valid level.
473 int GetMaxLevel(double value) const;
474
475 private:
476 double const deriv_;
477 DISALLOW_EVIL_CONSTRUCTORS(Metric);
478 };
479 typedef Metric<1> LengthMetric;
480 typedef Metric<2> AreaMetric;
481
482 // Each cell is bounded by four planes passing through its four edges and
483 // the center of the sphere. These metrics relate to the angle between each
484 // pair of opposite bounding planes, or equivalently, between the planes
485 // corresponding to two different s-values or two different t-values. For
486 // example, the maximum angle between opposite bounding planes for a cell at
487 // level k is kMaxAngleSpan.GetValue(k), and the average angle span for all
488 // cells at level k is approximately kAvgAngleSpan.GetValue(k).
489 static LengthMetric const kMinAngleSpan;
490 static LengthMetric const kMaxAngleSpan;
491 static LengthMetric const kAvgAngleSpan;
492
493 // The width of geometric figure is defined as the distance between two
494 // parallel bounding lines in a given direction. For cells, the minimum
495 // width is always attained between two opposite edges, and the maximum
496 // width is attained between two opposite vertices. However, for our
497 // purposes we redefine the width of a cell as the perpendicular distance
498 // between a pair of opposite edges. A cell therefore has two widths, one
499 // in each direction. The minimum width according to this definition agrees
500 // with the classic geometric one, but the maximum width is different. (The
501 // maximum geometric width corresponds to kMaxDiag defined below.)
502 //
503 // For a cell at level k, the distance between opposite edges is at least
504 // kMinWidth.GetValue(k) and at most kMaxWidth.GetValue(k). The average
505 // width in both directions for all cells at level k is approximately
506 // kAvgWidth.GetValue(k).
507 //
508 // The width is useful for bounding the minimum or maximum distance from a
509 // point on one edge of a cell to the closest point on the opposite edge.
510 // For example, this is useful when "growing" regions by a fixed distance.
511 static LengthMetric const kMinWidth;
512 static LengthMetric const kMaxWidth;
513 static LengthMetric const kAvgWidth;
514
515 // The minimum edge length of any cell at level k is at least
516 // kMinEdge.GetValue(k), and the maximum is at most kMaxEdge.GetValue(k).
517 // The average edge length is approximately kAvgEdge.GetValue(k).
518 //
519 // The edge length metrics can also be used to bound the minimum, maximum,
520 // or average distance from the center of one cell to the center of one of
521 // its edge neighbors. In particular, it can be used to bound the distance
522 // between adjacent cell centers along the space-filling Hilbert curve for
523 // cells at any given level.
524 static LengthMetric const kMinEdge;
525 static LengthMetric const kMaxEdge;
526 static LengthMetric const kAvgEdge;
527
528 // The minimum diagonal length of any cell at level k is at least
529 // kMinDiag.GetValue(k), and the maximum is at most kMaxDiag.GetValue(k).
530 // The average diagonal length is approximately kAvgDiag.GetValue(k).
531 //
532 // The maximum diagonal also happens to be the maximum diameter of any cell,
533 // and also the maximum geometric width (see the discussion above). So for
534 // example, the distance from an arbitrary point to the closest cell center
535 // at a given level is at most half the maximum diagonal length.
536 static LengthMetric const kMinDiag;
537 static LengthMetric const kMaxDiag;
538 static LengthMetric const kAvgDiag;
539
540 // The minimum area of any cell at level k is at least kMinArea.GetValue(k),
541 // and the maximum is at most kMaxArea.GetValue(k). The average area of all
542 // cells at level k is exactly kAvgArea.GetValue(k).
543 static AreaMetric const kMinArea;
544 static AreaMetric const kMaxArea;
545 static AreaMetric const kAvgArea;
546
547 // This is the maximum edge aspect ratio over all cells at any level, where
548 // the edge aspect ratio of a cell is defined as the ratio of its longest
549 // edge length to its shortest edge length.
550 static double const kMaxEdgeAspect;
551
552 // This is the maximum diagonal aspect ratio over all cells at any level,
553 // where the diagonal aspect ratio of a cell is defined as the ratio of its
554 // longest diagonal length to its shortest diagonal length.
555 static double const kMaxDiagAspect;
556
557 private:
558 // Given a *valid* face for the given point p (meaning that dot product
559 // of p with the face normal is positive), return the corresponding
560 // u and v values (which may lie outside the range [-1,1]).
561 inline static void ValidFaceXYZtoUV(int face, S2Point const& p,
562 double* pu, double* pv);
563
564 // The value below is the maximum error in computing the determinant
565 // a.CrossProd(b).DotProd(c). To derive this, observe that computing the
566 // determinant in this way requires 14 multiplications and additions. Since
567 // all three points are normalized, none of the intermediate results in this
568 // calculation exceed 1.0 in magnitude. The maximum rounding error for an
569 // operation whose result magnitude does not exceed 1.0 (before rounding) is
570 // 2**-54 (i.e., half of the difference between 1.0 and the next
571 // representable value below 1.0). Therefore, the total error in computing
572 // the determinant does not exceed 14 * (2**-54).
573 //
574 // The C++ standard requires to initialize kMaxDetError outside of
575 // the class definition, even though GCC doesn't enforce it.
576 static double const kMaxDetError;
577
578 DISALLOW_IMPLICIT_CONSTRUCTORS(S2); // Contains only static methods.
579};
580
581// Uncomment the following line for testing purposes only. It greatly
582// increases the number of degenerate cases that need to be handled using
583// ExpensiveCCW().
584// #define S2_TEST_DEGENERACIES
585
586inline S2Point S2::Origin() {
587#ifdef S2_TEST_DEGENERACIES
588 return S2Point(0, 0, 1); // This makes polygon operations much slower.
589#else
590 return S2Point(0.00457, 1, 0.0321).Normalize();
591#endif
592}
593
594inline int S2::TriageCCW(S2Point const& a, S2Point const& b,
595 S2Point const& c, S2Point const& a_cross_b) {
596 DCHECK(IsUnitLength(a));
597 DCHECK(IsUnitLength(b));
598 DCHECK(IsUnitLength(c));
599 double det = a_cross_b.DotProd(c);
600
601 // Double-check borderline cases in debug mode.
602 DCHECK(fabs(det) < kMaxDetError ||
603 fabs(det) > 100 * kMaxDetError ||
604 det * ExpensiveCCW(a, b, c) > 0);
605
606 if (det > kMaxDetError) return 1;
607 if (det < -kMaxDetError) return -1;
608 return 0;
609}
610
611inline int S2::RobustCCW(S2Point const& a, S2Point const& b,
612 S2Point const& c, S2Point const& a_cross_b) {
613 int ccw = TriageCCW(a, b, c, a_cross_b);
614 if (ccw == 0) ccw = ExpensiveCCW(a, b, c);
615 return ccw;
616}
617
618// We have implemented three different projections from cell-space (s,t) to
619// cube-space (u,v): linear, quadratic, and tangent. They have the following
620// tradeoffs:
621//
622// Linear - This is the fastest transformation, but also produces the least
623// uniform cell sizes. Cell areas vary by a factor of about 5.2, with the
624// largest cells at the center of each face and the smallest cells in
625// the corners.
626//
627// Tangent - Transforming the coordinates via atan() makes the cell sizes
628// more uniform. The areas vary by a maximum ratio of 1.4 as opposed to a
629// maximum ratio of 5.2. However, each call to atan() is about as expensive
630// as all of the other calculations combined when converting from points to
631// cell ids, i.e. it reduces performance by a factor of 3.
632//
633// Quadratic - This is an approximation of the tangent projection that
634// is much faster and produces cells that are almost as uniform in size.
635// It is about 3 times faster than the tangent projection for converting
636// cell ids to points or vice versa. Cell areas vary by a maximum ratio of
637// about 2.1.
638//
639// Here is a table comparing the cell uniformity using each projection. "Area
640// ratio" is the maximum ratio over all subdivision levels of the largest cell
641// area to the smallest cell area at that level, "edge ratio" is the maximum
642// ratio of the longest edge of any cell to the shortest edge of any cell at
643// the same level, and "diag ratio" is the ratio of the longest diagonal of
644// any cell to the shortest diagonal of any cell at the same level. "ToPoint"
645// and "FromPoint" are the times in microseconds required to convert cell ids
646// to and from points (unit vectors) respectively. "ToPointRaw" is the time
647// to convert to a non-unit-length vector, which is all that is needed for
648// some purposes.
649//
650// Area Edge Diag ToPointRaw ToPoint FromPoint
651// Ratio Ratio Ratio (microseconds)
652// -------------------------------------------------------------------
653// Linear: 5.200 2.117 2.959 0.020 0.087 0.085
654// Tangent: 1.414 1.414 1.704 0.237 0.299 0.258
655// Quadratic: 2.082 1.802 1.932 0.033 0.096 0.108
656//
657// The worst-case cell aspect ratios are about the same with all three
658// projections. The maximum ratio of the longest edge to the shortest edge
659// within the same cell is about 1.4 and the maximum ratio of the diagonals
660// within the same cell is about 1.7.
661//
662// This data was produced using s2cell_unittest and s2cellid_unittest.
663
664#define S2_LINEAR_PROJECTION 0
665#define S2_TAN_PROJECTION 1
666#define S2_QUADRATIC_PROJECTION 2
667
668#define S2_PROJECTION S2_QUADRATIC_PROJECTION
669
670#if S2_PROJECTION == S2_LINEAR_PROJECTION
671
672inline double S2::STtoUV(double s) {
673 return 2 * s - 1;
674}
675
676inline double S2::UVtoST(double u) {
677 return 0.5 * (u + 1);
678}
679
680#elif S2_PROJECTION == S2_TAN_PROJECTION
681
682inline double S2::STtoUV(double s) {
683 // Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn't due to
684 // a flaw in the implementation of tan(), it's because the derivative of
685 // tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating
686 // point numbers on either side of the infinite-precision value of pi/4 have
687 // tangents that are slightly below and slightly above 1.0 when rounded to
688 // the nearest double-precision result.
689
690 s = tan(M_PI_2 * s - M_PI_4);
691 return s + (1.0 / (GG_LONGLONG(1) << 53)) * s;
692}
693
694inline double S2::UVtoST(double u) {
695 volatile double a = atan(u);
696 return (2 * M_1_PI) * (a + M_PI_4);
697}
698
699#elif S2_PROJECTION == S2_QUADRATIC_PROJECTION
700
701inline double S2::STtoUV(double s) {
702 if (s >= 0.5) return (1/3.) * (4*s*s - 1);
703 else return (1/3.) * (1 - 4*(1-s)*(1-s));
704}
705
706inline double S2::UVtoST(double u) {
707 if (u >= 0) return 0.5 * sqrt(1 + 3*u);
708 else return 1 - 0.5 * sqrt(1 - 3*u);
709}
710
711#else
712
713#error Unknown value for S2_PROJECTION
714
715#endif
716
717inline S2Point S2::FaceUVtoXYZ(int face, double u, double v) {
718 switch (face) {
719 case 0: return S2Point( 1, u, v);
720 case 1: return S2Point(-u, 1, v);
721 case 2: return S2Point(-u, -v, 1);
722 case 3: return S2Point(-1, -v, -u);
723 case 4: return S2Point( v, -1, -u);
724 default: return S2Point( v, u, -1);
725 }
726}
727
728inline void S2::ValidFaceXYZtoUV(int face, S2Point const& p,
729 double* pu, double* pv) {
730 DCHECK_GT(p.DotProd(FaceUVtoXYZ(face, 0, 0)), 0);
731 switch (face) {
732 case 0: *pu = p[1] / p[0]; *pv = p[2] / p[0]; break;
733 case 1: *pu = -p[0] / p[1]; *pv = p[2] / p[1]; break;
734 case 2: *pu = -p[0] / p[2]; *pv = -p[1] / p[2]; break;
735 case 3: *pu = p[2] / p[0]; *pv = p[1] / p[0]; break;
736 case 4: *pu = p[2] / p[1]; *pv = -p[0] / p[1]; break;
737 default: *pu = -p[1] / p[2]; *pv = -p[0] / p[2]; break;
738 }
739}
740
741inline int S2::XYZtoFaceUV(S2Point const& p, double* pu, double* pv) {
742 int face = p.LargestAbsComponent();
743 if (p[face] < 0) face += 3;
744 ValidFaceXYZtoUV(face, p, pu, pv);
745 return face;
746}
747
748inline bool S2::FaceXYZtoUV(int face, S2Point const& p,
749 double* pu, double* pv) {
750 if (face < 3) {
751 if (p[face] <= 0) return false;
752 } else {
753 if (p[face-3] >= 0) return false;
754 }
755 ValidFaceXYZtoUV(face, p, pu, pv);
756 return true;
757}
758
759inline S2Point S2::GetUNorm(int face, double u) {
760 switch (face) {
761 case 0: return S2Point( u, -1, 0);
762 case 1: return S2Point( 1, u, 0);
763 case 2: return S2Point( 1, 0, u);
764 case 3: return S2Point(-u, 0, 1);
765 case 4: return S2Point( 0, -u, 1);
766 default: return S2Point( 0, -1, -u);
767 }
768}
769
770inline S2Point S2::GetVNorm(int face, double v) {
771 switch (face) {
772 case 0: return S2Point(-v, 0, 1);
773 case 1: return S2Point( 0, -v, 1);
774 case 2: return S2Point( 0, -1, -v);
775 case 3: return S2Point( v, -1, 0);
776 case 4: return S2Point( 1, v, 0);
777 default: return S2Point( 1, 0, v);
778 }
779}
780
781inline S2Point S2::GetNorm(int face) {
782 return S2::FaceUVtoXYZ(face, 0, 0);
783}
784
785inline S2Point S2::GetUAxis(int face) {
786 switch (face) {
787 case 0: return S2Point( 0, 1, 0);
788 case 1: return S2Point(-1, 0, 0);
789 case 2: return S2Point(-1, 0, 0);
790 case 3: return S2Point( 0, 0, -1);
791 case 4: return S2Point( 0, 0, -1);
792 default: return S2Point( 0, 1, 0);
793 }
794}
795
796inline S2Point S2::GetVAxis(int face) {
797 switch (face) {
798 case 0: return S2Point( 0, 0, 1);
799 case 1: return S2Point( 0, 0, 1);
800 case 2: return S2Point( 0, -1, 0);
801 case 3: return S2Point( 0, -1, 0);
802 case 4: return S2Point( 1, 0, 0);
803 default: return S2Point( 1, 0, 0);
804 }
805}
806
807template <int dim>
808int S2::Metric<dim>::GetMinLevel(double value) const {
809 if (value <= 0) return S2::kMaxCellLevel;
810
811 // This code is equivalent to computing a floating-point "level"
812 // value and rounding up. frexp() returns a fraction in the
813 // range [0.5,1) and the corresponding exponent.
814 int level;
815 frexp(value / deriv_, &level);
816 level = max(0, min(S2::kMaxCellLevel, -((level - 1) >> (dim - 1))));
817 DCHECK(level == S2::kMaxCellLevel || GetValue(level) <= value);
818 DCHECK(level == 0 || GetValue(level - 1) > value);
819 return level;
820}
821
822template <int dim>
823int S2::Metric<dim>::GetMaxLevel(double value) const {
824 if (value <= 0) return S2::kMaxCellLevel;
825
826 // This code is equivalent to computing a floating-point "level"
827 // value and rounding down.
828 int level;
829 frexp(deriv_ / value, &level);
830 level = max(0, min(S2::kMaxCellLevel, (level - 1) >> (dim - 1)));
831 DCHECK(level == 0 || GetValue(level) >= value);
832 DCHECK(level == S2::kMaxCellLevel || GetValue(level + 1) < value);
833 return level;
834}
835
836template <int dim>
837int S2::Metric<dim>::GetClosestLevel(double value) const {
838 return GetMinLevel((dim == 1 ? M_SQRT2 : 2) * value);
839}
840
841#endif // UTIL_GEOMETRY_S2_H_
842