1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "s2.h"
4
5// #include "base/commandlineflags.h"
6#include "base/integral_types.h"
7#include "base/logging.h"
8#include "util/math/matrix3x3-inl.h"
9#include "util/math/vector2-inl.h"
10
11// Define storage for header file constants (the values are not needed
12// here for integral constants).
13int const S2::kMaxCellLevel;
14int const S2::kSwapMask;
15int const S2::kInvertMask;
16double const S2::kMaxDetError = 0.8e-15; // 14 * (2**-54)
17
18COMPILE_ASSERT(S2::kSwapMask == 0x01 && S2::kInvertMask == 0x02,
19 masks_changed);
20
21// DEFINE_bool(s2debug, DEBUG_MODE, "Enable debugging checks in s2 code");
22bool const S2::debug = DEBUG_MODE;
23
24static const uint32 MIX32 = 0x12b9b0a1UL;
25#include<hash_set>
26namespace __gnu_cxx {
27
28
29
30// The hash function due to Bob Jenkins (see
31// http://burtleburtle.net/bob/hash/index.html).
32static inline void mix(uint32& a, uint32& b, uint32& c) { // 32bit version
33 a -= b; a -= c; a ^= (c>>13);
34 b -= c; b -= a; b ^= (a<<8);
35 c -= a; c -= b; c ^= (b>>13);
36 a -= b; a -= c; a ^= (c>>12);
37 b -= c; b -= a; b ^= (a<<16);
38 c -= a; c -= b; c ^= (b>>5);
39 a -= b; a -= c; a ^= (c>>3);
40 b -= c; b -= a; b ^= (a<<10);
41 c -= a; c -= b; c ^= (b>>15);
42}
43
44inline uint32 CollapseZero(uint32 bits) {
45 // IEEE 754 has two representations for zero, positive zero and negative
46 // zero. These two values compare as equal, and therefore we need them to
47 // hash to the same value.
48 //
49 // We handle this by simply clearing the top bit of every 32-bit value,
50 // which clears the sign bit on both big-endian and little-endian
51 // architectures. This creates some additional hash collisions between
52 // points that differ only in the sign of their components, but this is
53 // rarely a problem with real data.
54 //
55 // The obvious alternative is to explicitly map all occurrences of positive
56 // zero to negative zero (or vice versa), but this is more expensive and
57 // makes the average case slower.
58 //
59 // We also mask off the low-bit because we've seen differences in
60 // some floating point operations (specifically 'fcos' on i386)
61 // between different implementations of the same architecure
62 // (e.g. 'Xeon 5345' vs. 'Opteron 270'). It's unknown how many bits
63 // of mask are sufficient to cover real world cases, but the intent
64 // is to be as conservative as possible in discarding bits.
65
66 return bits & 0x7ffffffe;
67}
68
69size_t hash<S2Point>::operator()(S2Point const& p) const {
70 // This function is significantly faster than calling HashTo32().
71 uint32 const* data = reinterpret_cast<uint32 const*>(p.Data());
72 DCHECK_EQ((6 * sizeof(*data)), sizeof(p));
73
74 // We call CollapseZero() on every 32-bit chunk to avoid having endian
75 // dependencies.
76 uint32 a = CollapseZero(data[0]);
77 uint32 b = CollapseZero(data[1]);
78 uint32 c = CollapseZero(data[2]) + 0x12b9b0a1UL; // An arbitrary number
79 mix(a, b, c);
80 a += CollapseZero(data[3]);
81 b += CollapseZero(data[4]);
82 c += CollapseZero(data[5]);
83 mix(a, b, c);
84 return c;
85}
86
87
88} // namespace __gnu_cxx
89
90
91bool S2::IsUnitLength(S2Point const& p) {
92
93 return fabs(p.Norm2() - 1) <= 1e-15;
94}
95
96S2Point S2::Ortho(S2Point const& a) {
97#ifdef S2_TEST_DEGENERACIES
98 // Vector3::Ortho() always returns a point on the X-Y, Y-Z, or X-Z planes.
99 // This leads to many more degenerate cases in polygon operations.
100 return a.Ortho();
101#else
102 int k = a.LargestAbsComponent() - 1;
103 if (k < 0) k = 2;
104 S2Point temp(0.012, 0.0053, 0.00457);
105 temp[k] = 1;
106 return a.CrossProd(temp).Normalize();
107#endif
108}
109
110void S2::GetFrame(S2Point const& z, Matrix3x3_d* m) {
111 DCHECK(IsUnitLength(z));
112 m->SetCol(2, z);
113 m->SetCol(1, Ortho(z));
114 m->SetCol(0, m->Col(1).CrossProd(z)); // Already unit-length.
115}
116
117S2Point S2::ToFrame(Matrix3x3_d const& m, S2Point const& p) {
118 // The inverse of an orthonormal matrix is its transpose.
119 return m.Transpose() * p;
120}
121
122S2Point S2::FromFrame(Matrix3x3_d const& m, S2Point const& q) {
123 return m * q;
124}
125
126bool S2::ApproxEquals(S2Point const& a, S2Point const& b, double max_error) {
127 return a.Angle(b) <= max_error;
128}
129
130S2Point S2::RobustCrossProd(S2Point const& a, S2Point const& b) {
131 // The direction of a.CrossProd(b) becomes unstable as (a + b) or (a - b)
132 // approaches zero. This leads to situations where a.CrossProd(b) is not
133 // very orthogonal to "a" and/or "b". We could fix this using Gram-Schmidt,
134 // but we also want b.RobustCrossProd(a) == -a.RobustCrossProd(b).
135 //
136 // The easiest fix is to just compute the cross product of (b+a) and (b-a).
137 // Mathematically, this cross product is exactly twice the cross product of
138 // "a" and "b", but it has the numerical advantage that (b+a) and (b-a)
139 // are always perpendicular (since "a" and "b" are unit length). This
140 // yields a result that is nearly orthogonal to both "a" and "b" even if
141 // these two values differ only in the lowest bit of one component.
142
143 DCHECK(IsUnitLength(a));
144 DCHECK(IsUnitLength(b));
145 S2Point x = (b + a).CrossProd(b - a);
146 if (x != S2Point(0, 0, 0)) return x;
147
148 // The only result that makes sense mathematically is to return zero, but
149 // we find it more convenient to return an arbitrary orthogonal vector.
150 return Ortho(a);
151}
152
153bool S2::SimpleCCW(S2Point const& a, S2Point const& b, S2Point const& c) {
154 // We compute the signed volume of the parallelepiped ABC. The usual
155 // formula for this is (AxB).C, but we compute it here using (CxA).B
156 // in order to ensure that ABC and CBA are not both CCW. This follows
157 // from the following identities (which are true numerically, not just
158 // mathematically):
159 //
160 // (1) x.CrossProd(y) == -(y.CrossProd(x))
161 // (2) (-x).DotProd(y) == -(x.DotProd(y))
162
163 return c.CrossProd(a).DotProd(b) > 0;
164}
165
166int S2::RobustCCW(S2Point const& a, S2Point const& b, S2Point const& c) {
167 // We don't need RobustCrossProd() here because RobustCCW() does its own
168 // error estimation and calls ExpensiveCCW() if there is any uncertainty
169 // about the result.
170 return RobustCCW(a, b, c, a.CrossProd(b));
171}
172
173// Below we define two versions of ExpensiveCCW(). The first version uses
174// arbitrary-precision arithmetic (MPFloat) and the "simulation of simplicity"
175// technique. It is completely robust (i.e., it returns consistent results
176// for all possible inputs). The second version uses normal double-precision
177// arithmetic. It is numerically stable and handles many degeneracies well,
178// but it is not perfectly robust. It exists mainly for testing purposes, so
179// that we can verify that certain tests actually require the more advanced
180// techniques implemented by the first version.
181
182#define SIMULATION_OF_SIMPLICITY
183#ifdef SIMULATION_OF_SIMPLICITY
184
185// Below we define a floating-point type with enough precision so that it can
186// represent the exact determinant of any 3x3 matrix of floating-point
187// numbers. We support two options: MPFloat (which is based on MPFR and is
188// therefore subject to an LGPL license) and ExactFloat (which is based on the
189// OpenSSL Bignum library and therefore has a permissive BSD-style license).
190
191#ifdef S2_USE_EXACTFLOAT
192
193// ExactFloat only supports exact calculations with floating-point numbers.
194#include "util/math/exactfloat/exactfloat.h"
195
196#else // S2_USE_EXACTFLOAT
197
198// MPFloat requires a "maximum precision" to be specified.
199//
200// To figure out how much precision we need, first observe that any double
201// precision number can be represented as an integer by multiplying it by
202// 2**1074. This is because the minimum exponent is -1022, and denormalized
203// numbers have 52 bits after the leading "0". On the other hand, the largest
204// double precision value has the form 1.x * (2**1023), which is a 1024-bit
205// integer. Therefore any double precision value can be represented as a
206// (1074 + 1024) = 2098 bit integer.
207//
208// A 3x3 determinant is computed by adding together 6 values, each of which is
209// the product of 3 of the input values. When an m-bit integer is multiplied
210// by an n-bit integer, the result has at most (m+n) bits. When "k" m-bit
211// integers are added together, the result has at most m + ceil(log2(k)) bits.
212// Therefore the determinant of any 3x3 matrix can be represented exactly
213// using no more than (3*2098)+3 = 6297 bits.
214//
215// Note that MPFloat only uses as much precision as required to compute the
216// exact result, and that typically far fewer bits of precision are used. The
217// worst-case estimate above is only achieved for a matrix where every row
218// contains both the maximum and minimum possible double precision values
219// (i.e. approximately 1e308 and 1e-323). For randomly chosen unit-length
220// vectors, the average case uses only about 200 bits of precision.
221
222// The maximum precision must be at least (6297 + 1) so that we can assert
223// that the result of the determinant calculation is exact (by checking that
224// the actual precision of the result is less than the maximum precision
225// specified).
226
227#include "util/math/mpfloat/mpfloat.h"
228typedef MPFloat<6300> ExactFloat;
229
230#endif // S2_USE_EXACTFLOAT
231
232typedef Vector3<ExactFloat> Vector3_xf;
233
234// The following function returns the sign of the determinant of three points
235// A, B, C under a model where every possible S2Point is slightly perturbed by
236// a unique infinitesmal amount such that no three perturbed points are
237// collinear and no four points are coplanar. The perturbations are so small
238// that they do not change the sign of any determinant that was non-zero
239// before the perturbations, and therefore can be safely ignored unless the
240// determinant of three points is exactly zero (using multiple-precision
241// arithmetic).
242//
243// Since the symbolic perturbation of a given point is fixed (i.e., the
244// perturbation is the same for all calls to this method and does not depend
245// on the other two arguments), the results of this method are always
246// self-consistent. It will never return results that would correspond to an
247// "impossible" configuration of non-degenerate points.
248//
249// Requirements:
250// The 3x3 determinant of A, B, C must be exactly zero.
251// The points must be distinct, with A < B < C in lexicographic order.
252//
253// Returns:
254// +1 or -1 according to the sign of the determinant after the symbolic
255// perturbations are taken into account.
256//
257// Reference:
258// "Simulation of Simplicity" (Edelsbrunner and Muecke, ACM Transactions on
259// Graphics, 1990).
260//
261static int SymbolicallyPerturbedCCW(
262 Vector3_xf const& a, Vector3_xf const& b,
263 Vector3_xf const& c, Vector3_xf const& b_cross_c) {
264 // This method requires that the points are sorted in lexicographically
265 // increasing order. This is because every possible S2Point has its own
266 // symbolic perturbation such that if A < B then the symbolic perturbation
267 // for A is much larger than the perturbation for B.
268 //
269 // Alternatively, we could sort the points in this method and keep track of
270 // the sign of the permutation, but it is more efficient to do this before
271 // converting the inputs to the multi-precision representation, and this
272 // also lets us re-use the result of the cross product B x C.
273 DCHECK(a < b && b < c);
274
275 // Every input coordinate x[i] is assigned a symbolic perturbation dx[i].
276 // We then compute the sign of the determinant of the perturbed points,
277 // i.e.
278 // | a[0]+da[0] a[1]+da[1] a[2]+da[2] |
279 // | b[0]+db[0] b[1]+db[1] b[2]+db[2] |
280 // | c[0]+dc[0] c[1]+dc[1] c[2]+dc[2] |
281 //
282 // The perturbations are chosen such that
283 //
284 // da[2] > da[1] > da[0] > db[2] > db[1] > db[0] > dc[2] > dc[1] > dc[0]
285 //
286 // where each perturbation is so much smaller than the previous one that we
287 // don't even need to consider it unless the coefficients of all previous
288 // perturbations are zero. In fact, it is so small that we don't need to
289 // consider it unless the coefficient of all products of the previous
290 // perturbations are zero. For example, we don't need to consider the
291 // coefficient of db[1] unless the coefficient of db[2]*da[0] is zero.
292 //
293 // The follow code simply enumerates the coefficients of the perturbations
294 // (and products of perturbations) that appear in the determinant above, in
295 // order of decreasing perturbation magnitude. The first non-zero
296 // coefficient determines the sign of the result. The easiest way to
297 // enumerate the coefficients in the correct order is to pretend that each
298 // perturbation is some tiny value "eps" raised to a power of two:
299 //
300 // eps** 1 2 4 8 16 32 64 128 256
301 // da[2] da[1] da[0] db[2] db[1] db[0] dc[2] dc[1] dc[0]
302 //
303 // Essentially we can then just count in binary and test the corresponding
304 // subset of perturbations at each step. So for example, we must test the
305 // coefficient of db[2]*da[0] before db[1] because eps**12 > eps**16.
306 //
307 // Of course, not all products of these perturbations appear in the
308 // determinant above, since the determinant only contains the products of
309 // elements in distinct rows and columns. Thus we don't need to consider
310 // da[2]*da[1], db[1]*da[1], etc. Furthermore, sometimes different pairs of
311 // perturbations have the same coefficient in the determinant; for example,
312 // da[1]*db[0] and db[1]*da[0] have the same coefficient (c[2]). Therefore
313 // we only need to test this coefficient the first time we encounter it in
314 // the binary order above (which will be db[1]*da[0]).
315 //
316 // The sequence of tests below also appears in Table 4-ii of the paper
317 // referenced above, if you just want to look it up, with the following
318 // translations: [a,b,c] -> [i,j,k] and [0,1,2] -> [1,2,3]. Also note that
319 // some of the signs are different because the opposite cross product is
320 // used (e.g., B x C rather than C x B).
321
322 int det_sign = b_cross_c[2].sgn(); // da[2]
323 if (det_sign != 0) return det_sign;
324 det_sign = b_cross_c[1].sgn(); // da[1]
325 if (det_sign != 0) return det_sign;
326 det_sign = b_cross_c[0].sgn(); // da[0]
327 if (det_sign != 0) return det_sign;
328
329 det_sign = (c[0]*a[1] - c[1]*a[0]).sgn(); // db[2]
330 if (det_sign != 0) return det_sign;
331 det_sign = c[0].sgn(); // db[2] * da[1]
332 if (det_sign != 0) return det_sign;
333 det_sign = -(c[1].sgn()); // db[2] * da[0]
334 if (det_sign != 0) return det_sign;
335 det_sign = (c[2]*a[0] - c[0]*a[2]).sgn(); // db[1]
336 if (det_sign != 0) return det_sign;
337 det_sign = c[2].sgn(); // db[1] * da[0]
338 if (det_sign != 0) return det_sign;
339 // The following test is listed in the paper, but it is redundant because
340 // the previous tests guarantee that C == (0, 0, 0).
341 DCHECK_EQ(0, (c[1]*a[2] - c[2]*a[1]).sgn()); // db[0]
342
343 det_sign = (a[0]*b[1] - a[1]*b[0]).sgn(); // dc[2]
344 if (det_sign != 0) return det_sign;
345 det_sign = -(b[0].sgn()); // dc[2] * da[1]
346 if (det_sign != 0) return det_sign;
347 det_sign = b[1].sgn(); // dc[2] * da[0]
348 if (det_sign != 0) return det_sign;
349 det_sign = a[0].sgn(); // dc[2] * db[1]
350 if (det_sign != 0) return det_sign;
351 return 1; // dc[2] * db[1] * da[0]
352}
353
354int S2::ExpensiveCCW(S2Point const& a, S2Point const& b, S2Point const& c) {
355 // Return zero if and only if two points are the same. This ensures (1).
356 if (a == b || b == c || c == a) return 0;
357
358 // Sort the three points in lexicographic order, keeping track of the sign
359 // of the permutation. (Each exchange inverts the sign of the determinant.)
360 int perm_sign = 1;
361 S2Point pa = a, pb = b, pc = c;
362 if (pa > pb) { swap(pa, pb); perm_sign = -perm_sign; }
363 if (pb > pc) { swap(pb, pc); perm_sign = -perm_sign; }
364 if (pa > pb) { swap(pa, pb); perm_sign = -perm_sign; }
365 DCHECK(pa < pb && pb < pc);
366
367 // Construct multiple-precision versions of the sorted points and compute
368 // their exact 3x3 determinant.
369 Vector3_xf xa = Vector3_xf::Cast(pa);
370 Vector3_xf xb = Vector3_xf::Cast(pb);
371 Vector3_xf xc = Vector3_xf::Cast(pc);
372 Vector3_xf xb_cross_xc = xb.CrossProd(xc);
373 ExactFloat det = xa.DotProd(xb_cross_xc);
374
375 // The precision of ExactFloat is high enough that the result should always
376 // be exact (no rounding was performed).
377 DCHECK(!det.is_nan());
378 DCHECK_LT(det.prec(), det.max_prec());
379
380 // If the exact determinant is non-zero, we're done.
381 int det_sign = det.sgn();
382 if (det_sign == 0) {
383 // Otherwise, we need to resort to symbolic perturbations to resolve the
384 // sign of the determinant.
385 det_sign = SymbolicallyPerturbedCCW(xa, xb, xc, xb_cross_xc);
386 }
387 DCHECK(det_sign != 0);
388 return perm_sign * det_sign;
389}
390
391#else // SIMULATION_OF_SIMPLICITY
392
393static inline int PlanarCCW(Vector2_d const& a, Vector2_d const& b) {
394 // Return +1 if the edge AB is CCW around the origin, etc.
395 double sab = (a.DotProd(b) > 0) ? -1 : 1;
396 Vector2_d vab = a + sab * b;
397 double da = a.Norm2();
398 double db = b.Norm2();
399 double sign;
400 if (da < db || (da == db && a < b)) {
401 sign = a.CrossProd(vab) * sab;
402 } else {
403 sign = vab.CrossProd(b);
404 }
405 if (sign > 0) return 1;
406 if (sign < 0) return -1;
407 return 0;
408}
409
410static inline int PlanarOrderedCCW(Vector2_d const& a, Vector2_d const& b,
411 Vector2_d const& c) {
412 int sum = 0;
413 sum += PlanarCCW(a, b);
414 sum += PlanarCCW(b, c);
415 sum += PlanarCCW(c, a);
416 if (sum > 0) return 1;
417 if (sum < 0) return -1;
418 return 0;
419}
420
421int S2::ExpensiveCCW(S2Point const& a, S2Point const& b, S2Point const& c) {
422 // Return zero if and only if two points are the same. This ensures (1).
423 if (a == b || b == c || c == a) return 0;
424
425 // Now compute the determinant in a stable way. Since all three points are
426 // unit length and we know that the determinant is very close to zero, this
427 // means that points are very nearly collinear. Furthermore, the most common
428 // situation is where two points are nearly identical or nearly antipodal.
429 // To get the best accuracy in this situation, it is important to
430 // immediately reduce the magnitude of the arguments by computing either
431 // A+B or A-B for each pair of points. Note that even if A and B differ
432 // only in their low bits, A-B can be computed very accurately. On the
433 // other hand we can't accurately represent an arbitrary linear combination
434 // of two vectors as would be required for Gaussian elimination. The code
435 // below chooses the vertex opposite the longest edge as the "origin" for
436 // the calculation, and computes the different vectors to the other two
437 // vertices. This minimizes the sum of the lengths of these vectors.
438 //
439 // This implementation is very stable numerically, but it still does not
440 // return consistent results in all cases. For example, if three points are
441 // spaced far apart from each other along a great circle, the sign of the
442 // result will basically be random (although it will still satisfy the
443 // conditions documented in the header file). The only way to return
444 // consistent results in all cases is to compute the result using
445 // multiple-precision arithmetic. I considered using the Gnu MP library,
446 // but this would be very expensive (up to 2000 bits of precision may be
447 // needed to store the intermediate results) and seems like overkill for
448 // this problem. The MP library is apparently also quite particular about
449 // compilers and compilation options and would be a pain to maintain.
450
451 // We want to handle the case of nearby points and nearly antipodal points
452 // accurately, so determine whether A+B or A-B is smaller in each case.
453 double sab = (a.DotProd(b) > 0) ? -1 : 1;
454 double sbc = (b.DotProd(c) > 0) ? -1 : 1;
455 double sca = (c.DotProd(a) > 0) ? -1 : 1;
456 S2Point vab = a + sab * b;
457 S2Point vbc = b + sbc * c;
458 S2Point vca = c + sca * a;
459 double dab = vab.Norm2();
460 double dbc = vbc.Norm2();
461 double dca = vca.Norm2();
462
463 // Sort the difference vectors to find the longest edge, and use the
464 // opposite vertex as the origin. If two difference vectors are the same
465 // length, we break ties deterministically to ensure that the symmetry
466 // properties guaranteed in the header file will be true.
467 double sign;
468 if (dca < dbc || (dca == dbc && a < b)) {
469 if (dab < dbc || (dab == dbc && a < c)) {
470 // The "sab" factor converts A +/- B into B +/- A.
471 sign = vab.CrossProd(vca).DotProd(a) * sab; // BC is longest edge
472 } else {
473 sign = vca.CrossProd(vbc).DotProd(c) * sca; // AB is longest edge
474 }
475 } else {
476 if (dab < dca || (dab == dca && b < c)) {
477 sign = vbc.CrossProd(vab).DotProd(b) * sbc; // CA is longest edge
478 } else {
479 sign = vca.CrossProd(vbc).DotProd(c) * sca; // AB is longest edge
480 }
481 }
482 if (sign > 0) return 1;
483 if (sign < 0) return -1;
484
485 // The points A, B, and C are numerically indistinguishable from coplanar.
486 // This may be due to roundoff error, or the points may in fact be exactly
487 // coplanar. We handle this situation by perturbing all of the points by a
488 // vector (eps, eps**2, eps**3) where "eps" is an infinitesmally small
489 // positive number (e.g. 1 divided by a googolplex). The perturbation is
490 // done symbolically, i.e. we compute what would happen if the points were
491 // perturbed by this amount. It turns out that this is equivalent to
492 // checking whether the points are ordered CCW around the origin first in
493 // the Y-Z plane, then in the Z-X plane, and then in the X-Y plane.
494
495 int ccw = PlanarOrderedCCW(Vector2_d(a.y(), a.z()), Vector2_d(b.y(), b.z()),
496 Vector2_d(c.y(), c.z()));
497 if (ccw == 0) {
498 ccw = PlanarOrderedCCW(Vector2_d(a.z(), a.x()), Vector2_d(b.z(), b.x()),
499 Vector2_d(c.z(), c.x()));
500 if (ccw == 0) {
501 ccw = PlanarOrderedCCW(Vector2_d(a.x(), a.y()), Vector2_d(b.x(), b.y()),
502 Vector2_d(c.x(), c.y()));
503 // There are a few cases where "ccw" may still be zero despite our best
504 // efforts. For example, two input points may be exactly proportional
505 // to each other (where both still satisfy IsNormalized()).
506 }
507 }
508 return ccw;
509}
510
511#endif // SIMULATION_OF_SIMPLICITY
512
513double S2::Angle(S2Point const& a, S2Point const& b, S2Point const& c) {
514 return RobustCrossProd(a, b).Angle(RobustCrossProd(c, b));
515}
516
517double S2::TurnAngle(S2Point const& a, S2Point const& b, S2Point const& c) {
518 // This is a bit less efficient because we compute all 3 cross products, but
519 // it ensures that TurnAngle(a,b,c) == -TurnAngle(c,b,a) for all a,b,c.
520 double angle = RobustCrossProd(b, a).Angle(RobustCrossProd(c, b));
521 return (RobustCCW(a, b, c) > 0) ? angle : -angle;
522}
523
524double S2::Area(S2Point const& a, S2Point const& b, S2Point const& c) {
525 DCHECK(IsUnitLength(a));
526 DCHECK(IsUnitLength(b));
527 DCHECK(IsUnitLength(c));
528 // This method is based on l'Huilier's theorem,
529 //
530 // tan(E/4) = sqrt(tan(s/2) tan((s-a)/2) tan((s-b)/2) tan((s-c)/2))
531 //
532 // where E is the spherical excess of the triangle (i.e. its area),
533 // a, b, c, are the side lengths, and
534 // s is the semiperimeter (a + b + c) / 2 .
535 //
536 // The only significant source of error using l'Huilier's method is the
537 // cancellation error of the terms (s-a), (s-b), (s-c). This leads to a
538 // *relative* error of about 1e-16 * s / min(s-a, s-b, s-c). This compares
539 // to a relative error of about 1e-15 / E using Girard's formula, where E is
540 // the true area of the triangle. Girard's formula can be even worse than
541 // this for very small triangles, e.g. a triangle with a true area of 1e-30
542 // might evaluate to 1e-5.
543 //
544 // So, we prefer l'Huilier's formula unless dmin < s * (0.1 * E), where
545 // dmin = min(s-a, s-b, s-c). This basically includes all triangles
546 // except for extremely long and skinny ones.
547 //
548 // Since we don't know E, we would like a conservative upper bound on
549 // the triangle area in terms of s and dmin. It's possible to show that
550 // E <= k1 * s * sqrt(s * dmin), where k1 = 2*sqrt(3)/Pi (about 1).
551 // Using this, it's easy to show that we should always use l'Huilier's
552 // method if dmin >= k2 * s^5, where k2 is about 1e-2. Furthermore,
553 // if dmin < k2 * s^5, the triangle area is at most k3 * s^4, where
554 // k3 is about 0.1. Since the best case error using Girard's formula
555 // is about 1e-15, this means that we shouldn't even consider it unless
556 // s >= 3e-4 or so.
557
558 // We use volatile doubles to force the compiler to truncate all of these
559 // quantities to 64 bits. Otherwise it may compute a value of dmin > 0
560 // simply because it chose to spill one of the intermediate values to
561 // memory but not one of the others.
562 volatile double sa = b.Angle(c);
563 volatile double sb = c.Angle(a);
564 volatile double sc = a.Angle(b);
565 volatile double s = 0.5 * (sa + sb + sc);
566 if (s >= 3e-4) {
567 // Consider whether Girard's formula might be more accurate.
568 double s2 = s * s;
569 double dmin = s - max(sa, max(sb, sc));
570 if (dmin < 1e-2 * s * s2 * s2) {
571 // This triangle is skinny enough to consider Girard's formula.
572 double area = GirardArea(a, b, c);
573 if (dmin < s * (0.1 * area)) return area;
574 }
575 }
576 // Use l'Huilier's formula.
577 return 4 * atan(sqrt(max(0.0, tan(0.5 * s) * tan(0.5 * (s - sa)) *
578 tan(0.5 * (s - sb)) * tan(0.5 * (s - sc)))));
579}
580
581double S2::GirardArea(S2Point const& a, S2Point const& b, S2Point const& c) {
582 // This is equivalent to the usual Girard's formula but is slightly
583 // more accurate, faster to compute, and handles a == b == c without
584 // a special case. The use of RobustCrossProd() makes it much more
585 // accurate when two vertices are nearly identical or antipodal.
586
587 S2Point ab = RobustCrossProd(a, b);
588 S2Point bc = RobustCrossProd(b, c);
589 S2Point ac = RobustCrossProd(a, c);
590 return max(0.0, ab.Angle(ac) - ab.Angle(bc) + bc.Angle(ac));
591}
592
593double S2::SignedArea(S2Point const& a, S2Point const& b, S2Point const& c) {
594 return Area(a, b, c) * RobustCCW(a, b, c);
595}
596
597S2Point S2::PlanarCentroid(S2Point const& a, S2Point const& b,
598 S2Point const& c) {
599 return (1./3) * (a + b + c);
600}
601
602S2Point S2::TrueCentroid(S2Point const& a, S2Point const& b,
603 S2Point const& c) {
604 DCHECK(IsUnitLength(a));
605 DCHECK(IsUnitLength(b));
606 DCHECK(IsUnitLength(c));
607
608 // I couldn't find any references for computing the true centroid of a
609 // spherical triangle... I have a truly marvellous demonstration of this
610 // formula which this margin is too narrow to contain :)
611
612 // Use Angle() in order to get accurate results for small triangles.
613 double angle_a = b.Angle(c);
614 double angle_b = c.Angle(a);
615 double angle_c = a.Angle(b);
616 double ra = (angle_a == 0) ? 1 : (angle_a / sin(angle_a));
617 double rb = (angle_b == 0) ? 1 : (angle_b / sin(angle_b));
618 double rc = (angle_c == 0) ? 1 : (angle_c / sin(angle_c));
619
620 // Now compute a point M such that:
621 //
622 // [Ax Ay Az] [Mx] [ra]
623 // [Bx By Bz] [My] = 0.5 * det(A,B,C) * [rb]
624 // [Cx Cy Cz] [Mz] [rc]
625 //
626 // To improve the numerical stability we subtract the first row (A) from the
627 // other two rows; this reduces the cancellation error when A, B, and C are
628 // very close together. Then we solve it using Cramer's rule.
629 //
630 // TODO(user): This code still isn't as numerically stable as it could be.
631 // The biggest potential improvement is to compute B-A and C-A more
632 // accurately so that (B-A)x(C-A) is always inside triangle ABC.
633 S2Point x(a.x(), b.x() - a.x(), c.x() - a.x());
634 S2Point y(a.y(), b.y() - a.y(), c.y() - a.y());
635 S2Point z(a.z(), b.z() - a.z(), c.z() - a.z());
636 S2Point r(ra, rb - ra, rc - ra);
637 return 0.5 * S2Point(y.CrossProd(z).DotProd(r),
638 z.CrossProd(x).DotProd(r),
639 x.CrossProd(y).DotProd(r));
640}
641
642bool S2::OrderedCCW(S2Point const& a, S2Point const& b, S2Point const& c,
643 S2Point const& o) {
644 // The last inequality below is ">" rather than ">=" so that we return true
645 // if A == B or B == C, and otherwise false if A == C. Recall that
646 // RobustCCW(x,y,z) == -RobustCCW(z,y,x) for all x,y,z.
647
648 int sum = 0;
649 if (RobustCCW(b, o, a) >= 0) ++sum;
650 if (RobustCCW(c, o, b) >= 0) ++sum;
651 if (RobustCCW(a, o, c) > 0) ++sum;
652 return sum >= 2;
653}
654
655// kIJtoPos[orientation][ij] -> pos
656int const S2::kIJtoPos[4][4] = {
657 // (0,0) (0,1) (1,0) (1,1)
658 { 0, 1, 3, 2 }, // canonical order
659 { 0, 3, 1, 2 }, // axes swapped
660 { 2, 3, 1, 0 }, // bits inverted
661 { 2, 1, 3, 0 }, // swapped & inverted
662};
663
664// kPosToIJ[orientation][pos] -> ij
665int const S2::kPosToIJ[4][4] = {
666 // 0 1 2 3
667 { 0, 1, 3, 2 }, // canonical order: (0,0), (0,1), (1,1), (1,0)
668 { 0, 2, 3, 1 }, // axes swapped: (0,0), (1,0), (1,1), (0,1)
669 { 3, 2, 0, 1 }, // bits inverted: (1,1), (1,0), (0,0), (0,1)
670 { 3, 1, 0, 2 }, // swapped & inverted: (1,1), (0,1), (0,0), (1,0)
671};
672
673// kPosToOrientation[pos] -> orientation_modifier
674int const S2::kPosToOrientation[4] = {
675 kSwapMask,
676 0,
677 0,
678 kInvertMask + kSwapMask,
679};
680
681// All of the values below were obtained by a combination of hand analysis and
682// Mathematica. In general, S2_TAN_PROJECTION produces the most uniform
683// shapes and sizes of cells, S2_LINEAR_PROJECTION is considerably worse, and
684// S2_QUADRATIC_PROJECTION is somewhere in between (but generally closer to
685// the tangent projection than the linear one).
686
687S2::LengthMetric const S2::kMinAngleSpan(
688 S2_PROJECTION == S2_LINEAR_PROJECTION ? 1.0 : // 1.000
689 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI / 2 : // 1.571
690 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 4. / 3 : // 1.333
691 0);
692
693S2::LengthMetric const S2::kMaxAngleSpan(
694 S2_PROJECTION == S2_LINEAR_PROJECTION ? 2 : // 2.000
695 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI / 2 : // 1.571
696 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 1.704897179199218452 : // 1.705
697 0);
698
699S2::LengthMetric const S2::kAvgAngleSpan(M_PI / 2); // 1.571
700// This is true for all projections.
701
702S2::LengthMetric const S2::kMinWidth(
703 S2_PROJECTION == S2_LINEAR_PROJECTION ? sqrt(2. / 3) : // 0.816
704 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI / (2 * sqrt(2)) : // 1.111
705 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 2 * sqrt(2) / 3 : // 0.943
706 0);
707
708S2::LengthMetric const S2::kMaxWidth(S2::kMaxAngleSpan.deriv());
709// This is true for all projections.
710
711S2::LengthMetric const S2::kAvgWidth(
712 S2_PROJECTION == S2_LINEAR_PROJECTION ? 1.411459345844456965 : // 1.411
713 S2_PROJECTION == S2_TAN_PROJECTION ? 1.437318638925160885 : // 1.437
714 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 1.434523672886099389 : // 1.435
715 0);
716
717S2::LengthMetric const S2::kMinEdge(
718 S2_PROJECTION == S2_LINEAR_PROJECTION ? 2 * sqrt(2) / 3 : // 0.943
719 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI / (2 * sqrt(2)) : // 1.111
720 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 2 * sqrt(2) / 3 : // 0.943
721 0);
722
723S2::LengthMetric const S2::kMaxEdge(S2::kMaxAngleSpan.deriv());
724// This is true for all projections.
725
726S2::LengthMetric const S2::kAvgEdge(
727 S2_PROJECTION == S2_LINEAR_PROJECTION ? 1.440034192955603643 : // 1.440
728 S2_PROJECTION == S2_TAN_PROJECTION ? 1.461667032546739266 : // 1.462
729 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 1.459213746386106062 : // 1.459
730 0);
731
732S2::LengthMetric const S2::kMinDiag(
733 S2_PROJECTION == S2_LINEAR_PROJECTION ? 2 * sqrt(2) / 3 : // 0.943
734 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI * sqrt(2) / 3 : // 1.481
735 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 8 * sqrt(2) / 9 : // 1.257
736 0);
737
738S2::LengthMetric const S2::kMaxDiag(
739 S2_PROJECTION == S2_LINEAR_PROJECTION ? 2 * sqrt(2) : // 2.828
740 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI * sqrt(2. / 3) : // 2.565
741 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 2.438654594434021032 : // 2.439
742 0);
743
744S2::LengthMetric const S2::kAvgDiag(
745 S2_PROJECTION == S2_LINEAR_PROJECTION ? 2.031817866418812674 : // 2.032
746 S2_PROJECTION == S2_TAN_PROJECTION ? 2.063623197195635753 : // 2.064
747 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 2.060422738998471683 : // 2.060
748 0);
749
750S2::AreaMetric const S2::kMinArea(
751 S2_PROJECTION == S2_LINEAR_PROJECTION ? 4 / (3 * sqrt(3)) : // 0.770
752 S2_PROJECTION == S2_TAN_PROJECTION ? (M_PI*M_PI) / (4*sqrt(2)) : // 1.745
753 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 8 * sqrt(2) / 9 : // 1.257
754 0);
755
756S2::AreaMetric const S2::kMaxArea(
757 S2_PROJECTION == S2_LINEAR_PROJECTION ? 4 : // 4.000
758 S2_PROJECTION == S2_TAN_PROJECTION ? M_PI * M_PI / 4 : // 2.467
759 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 2.635799256963161491 : // 2.636
760 0);
761
762S2::AreaMetric const S2::kAvgArea(4 * M_PI / 6); // 2.094
763// This is true for all projections.
764
765double const S2::kMaxEdgeAspect = (
766 S2_PROJECTION == S2_LINEAR_PROJECTION ? sqrt(2) : // 1.414
767 S2_PROJECTION == S2_TAN_PROJECTION ? sqrt(2) : // 1.414
768 S2_PROJECTION == S2_QUADRATIC_PROJECTION ? 1.442615274452682920 : // 1.443
769 0);
770
771double const S2::kMaxDiagAspect = sqrt(3); // 1.732
772// This is true for all projections.
773