1 // Copyright 2016 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 // This file contains various predicates that are guaranteed to produce 18 // correct, consistent results. They are also relatively efficient. This is 19 // achieved by computing conservative error bounds and falling back to high 20 // precision or even exact arithmetic when the result is uncertain. Such 21 // predicates are useful in implementing robust algorithms. 22 // 23 // See also EdgeCrosser, which implements various exact 24 // edge-crossing predicates more efficiently than can be done here. 25 26 import ( 27 "math" 28 "math/big" 29 30 "github.com/golang/geo/r3" 31 "github.com/golang/geo/s1" 32 ) 33 34 const ( 35 // If any other machine architectures need to be supported, these next three 36 // values will need to be updated. 37 38 // epsilon is a small number that represents a reasonable level of noise between two 39 // values that can be considered to be equal. 40 epsilon = 1e-15 41 // dblEpsilon is a smaller number for values that require more precision. 42 // This is the C++ DBL_EPSILON equivalent. 43 dblEpsilon = 2.220446049250313e-16 44 // dblError is the C++ value for S2 rounding_epsilon(). 45 dblError = 1.110223024625156e-16 46 47 // maxDeterminantError is the maximum error in computing (AxB).C where all vectors 48 // are unit length. Using standard inequalities, it can be shown that 49 // 50 // fl(AxB) = AxB + D where |D| <= (|AxB| + (2/sqrt(3))*|A|*|B|) * e 51 // 52 // where "fl()" denotes a calculation done in floating-point arithmetic, 53 // |x| denotes either absolute value or the L2-norm as appropriate, and 54 // e is a reasonably small value near the noise level of floating point 55 // number accuracy. Similarly, 56 // 57 // fl(B.C) = B.C + d where |d| <= (|B.C| + 2*|B|*|C|) * e . 58 // 59 // Applying these bounds to the unit-length vectors A,B,C and neglecting 60 // relative error (which does not affect the sign of the result), we get 61 // 62 // fl((AxB).C) = (AxB).C + d where |d| <= (3 + 2/sqrt(3)) * e 63 maxDeterminantError = 1.8274 * dblEpsilon 64 65 // detErrorMultiplier is the factor to scale the magnitudes by when checking 66 // for the sign of set of points with certainty. Using a similar technique to 67 // the one used for maxDeterminantError, the error is at most: 68 // 69 // |d| <= (3 + 6/sqrt(3)) * |A-C| * |B-C| * e 70 // 71 // If the determinant magnitude is larger than this value then we know 72 // its sign with certainty. 73 detErrorMultiplier = 3.2321 * dblEpsilon 74 ) 75 76 // Direction is an indication of the ordering of a set of points. 77 type Direction int 78 79 // These are the three options for the direction of a set of points. 80 const ( 81 Clockwise Direction = -1 82 Indeterminate Direction = 0 83 CounterClockwise Direction = 1 84 ) 85 86 // newBigFloat constructs a new big.Float with maximum precision. 87 func newBigFloat() *big.Float { return new(big.Float).SetPrec(big.MaxPrec) } 88 89 // Sign returns true if the points A, B, C are strictly counterclockwise, 90 // and returns false if the points are clockwise or collinear (i.e. if they are all 91 // contained on some great circle). 92 // 93 // Due to numerical errors, situations may arise that are mathematically 94 // impossible, e.g. ABC may be considered strictly CCW while BCA is not. 95 // However, the implementation guarantees the following: 96 // 97 // If Sign(a,b,c), then !Sign(c,b,a) for all a,b,c. 98 func Sign(a, b, c Point) bool { 99 // NOTE(dnadasi): In the C++ API the equivalent method here was known as "SimpleSign". 100 101 // We compute the signed volume of the parallelepiped ABC. The usual 102 // formula for this is (A ⨯ B) · C, but we compute it here using (C ⨯ A) · B 103 // in order to ensure that ABC and CBA are not both CCW. This follows 104 // from the following identities (which are true numerically, not just 105 // mathematically): 106 // 107 // (1) x ⨯ y == -(y ⨯ x) 108 // (2) -x · y == -(x · y) 109 return c.Cross(a.Vector).Dot(b.Vector) > 0 110 } 111 112 // RobustSign returns a Direction representing the ordering of the points. 113 // CounterClockwise is returned if the points are in counter-clockwise order, 114 // Clockwise for clockwise, and Indeterminate if any two points are the same (collinear), 115 // or the sign could not completely be determined. 116 // 117 // This function has additional logic to make sure that the above properties hold even 118 // when the three points are coplanar, and to deal with the limitations of 119 // floating-point arithmetic. 120 // 121 // RobustSign satisfies the following conditions: 122 // 123 // (1) RobustSign(a,b,c) == Indeterminate if and only if a == b, b == c, or c == a 124 // (2) RobustSign(b,c,a) == RobustSign(a,b,c) for all a,b,c 125 // (3) RobustSign(c,b,a) == -RobustSign(a,b,c) for all a,b,c 126 // 127 // In other words: 128 // 129 // (1) The result is Indeterminate if and only if two points are the same. 130 // (2) Rotating the order of the arguments does not affect the result. 131 // (3) Exchanging any two arguments inverts the result. 132 // 133 // On the other hand, note that it is not true in general that 134 // RobustSign(-a,b,c) == -RobustSign(a,b,c), or any similar identities 135 // involving antipodal points. 136 func RobustSign(a, b, c Point) Direction { 137 sign := triageSign(a, b, c) 138 if sign == Indeterminate { 139 sign = expensiveSign(a, b, c) 140 } 141 return sign 142 } 143 144 // stableSign reports the direction sign of the points in a numerically stable way. 145 // Unlike triageSign, this method can usually compute the correct determinant sign 146 // even when all three points are as collinear as possible. For example if three 147 // points are spaced 1km apart along a random line on the Earth's surface using 148 // the nearest representable points, there is only a 0.4% chance that this method 149 // will not be able to find the determinant sign. The probability of failure 150 // decreases as the points get closer together; if the collinear points are 1 meter 151 // apart, the failure rate drops to 0.0004%. 152 // 153 // This method could be extended to also handle nearly-antipodal points, but antipodal 154 // points are rare in practice so it seems better to simply fall back to 155 // exact arithmetic in that case. 156 func stableSign(a, b, c Point) Direction { 157 ab := b.Sub(a.Vector) 158 ab2 := ab.Norm2() 159 bc := c.Sub(b.Vector) 160 bc2 := bc.Norm2() 161 ca := a.Sub(c.Vector) 162 ca2 := ca.Norm2() 163 164 // Now compute the determinant ((A-C)x(B-C)).C, where the vertices have been 165 // cyclically permuted if necessary so that AB is the longest edge. (This 166 // minimizes the magnitude of cross product.) At the same time we also 167 // compute the maximum error in the determinant. 168 169 // The two shortest edges, pointing away from their common point. 170 var e1, e2, op r3.Vector 171 if ab2 >= bc2 && ab2 >= ca2 { 172 // AB is the longest edge. 173 e1, e2, op = ca, bc, c.Vector 174 } else if bc2 >= ca2 { 175 // BC is the longest edge. 176 e1, e2, op = ab, ca, a.Vector 177 } else { 178 // CA is the longest edge. 179 e1, e2, op = bc, ab, b.Vector 180 } 181 182 det := -e1.Cross(e2).Dot(op) 183 maxErr := detErrorMultiplier * math.Sqrt(e1.Norm2()*e2.Norm2()) 184 185 // If the determinant isn't zero, within maxErr, we know definitively the point ordering. 186 if det > maxErr { 187 return CounterClockwise 188 } 189 if det < -maxErr { 190 return Clockwise 191 } 192 return Indeterminate 193 } 194 195 // triageSign returns the direction sign of the points. It returns Indeterminate if two 196 // points are identical or the result is uncertain. Uncertain cases can be resolved, if 197 // desired, by calling expensiveSign. 198 // 199 // The purpose of this method is to allow additional cheap tests to be done without 200 // calling expensiveSign. 201 func triageSign(a, b, c Point) Direction { 202 det := a.Cross(b.Vector).Dot(c.Vector) 203 if det > maxDeterminantError { 204 return CounterClockwise 205 } 206 if det < -maxDeterminantError { 207 return Clockwise 208 } 209 return Indeterminate 210 } 211 212 // expensiveSign reports the direction sign of the points. It returns Indeterminate 213 // if two of the input points are the same. It uses multiple-precision arithmetic 214 // to ensure that its results are always self-consistent. 215 func expensiveSign(a, b, c Point) Direction { 216 // Return Indeterminate if and only if two points are the same. 217 // This ensures RobustSign(a,b,c) == Indeterminate if and only if a == b, b == c, or c == a. 218 // ie. Property 1 of RobustSign. 219 if a == b || b == c || c == a { 220 return Indeterminate 221 } 222 223 // Next we try recomputing the determinant still using floating-point 224 // arithmetic but in a more precise way. This is more expensive than the 225 // simple calculation done by triageSign, but it is still *much* cheaper 226 // than using arbitrary-precision arithmetic. This optimization is able to 227 // compute the correct determinant sign in virtually all cases except when 228 // the three points are truly collinear (e.g., three points on the equator). 229 detSign := stableSign(a, b, c) 230 if detSign != Indeterminate { 231 return detSign 232 } 233 234 // Otherwise fall back to exact arithmetic and symbolic permutations. 235 return exactSign(a, b, c, true) 236 } 237 238 // exactSign reports the direction sign of the points computed using high-precision 239 // arithmetic and/or symbolic perturbations. 240 func exactSign(a, b, c Point, perturb bool) Direction { 241 // Sort the three points in lexicographic order, keeping track of the sign 242 // of the permutation. (Each exchange inverts the sign of the determinant.) 243 permSign := CounterClockwise 244 pa := &a 245 pb := &b 246 pc := &c 247 if pa.Cmp(pb.Vector) > 0 { 248 pa, pb = pb, pa 249 permSign = -permSign 250 } 251 if pb.Cmp(pc.Vector) > 0 { 252 pb, pc = pc, pb 253 permSign = -permSign 254 } 255 if pa.Cmp(pb.Vector) > 0 { 256 pa, pb = pb, pa 257 permSign = -permSign 258 } 259 260 // Construct multiple-precision versions of the sorted points and compute 261 // their precise 3x3 determinant. 262 xa := r3.PreciseVectorFromVector(pa.Vector) 263 xb := r3.PreciseVectorFromVector(pb.Vector) 264 xc := r3.PreciseVectorFromVector(pc.Vector) 265 xbCrossXc := xb.Cross(xc) 266 det := xa.Dot(xbCrossXc) 267 268 // The precision of big.Float is high enough that the result should always 269 // be exact enough (no rounding was performed). 270 271 // If the exact determinant is non-zero, we're done. 272 detSign := Direction(det.Sign()) 273 if detSign == Indeterminate && perturb { 274 // Otherwise, we need to resort to symbolic perturbations to resolve the 275 // sign of the determinant. 276 detSign = symbolicallyPerturbedSign(xa, xb, xc, xbCrossXc) 277 } 278 return permSign * detSign 279 } 280 281 // symbolicallyPerturbedSign reports the sign of the determinant of three points 282 // A, B, C under a model where every possible Point is slightly perturbed by 283 // a unique infinitesmal amount such that no three perturbed points are 284 // collinear and no four points are coplanar. The perturbations are so small 285 // that they do not change the sign of any determinant that was non-zero 286 // before the perturbations, and therefore can be safely ignored unless the 287 // determinant of three points is exactly zero (using multiple-precision 288 // arithmetic). This returns CounterClockwise or Clockwise according to the 289 // sign of the determinant after the symbolic perturbations are taken into account. 290 // 291 // Since the symbolic perturbation of a given point is fixed (i.e., the 292 // perturbation is the same for all calls to this method and does not depend 293 // on the other two arguments), the results of this method are always 294 // self-consistent. It will never return results that would correspond to an 295 // impossible configuration of non-degenerate points. 296 // 297 // This requires that the 3x3 determinant of A, B, C must be exactly zero. 298 // And the points must be distinct, with A < B < C in lexicographic order. 299 // 300 // Reference: 301 // 302 // "Simulation of Simplicity" (Edelsbrunner and Muecke, ACM Transactions on 303 // Graphics, 1990). 304 func symbolicallyPerturbedSign(a, b, c, bCrossC r3.PreciseVector) Direction { 305 // This method requires that the points are sorted in lexicographically 306 // increasing order. This is because every possible Point has its own 307 // symbolic perturbation such that if A < B then the symbolic perturbation 308 // for A is much larger than the perturbation for B. 309 // 310 // Alternatively, we could sort the points in this method and keep track of 311 // the sign of the permutation, but it is more efficient to do this before 312 // converting the inputs to the multi-precision representation, and this 313 // also lets us re-use the result of the cross product B x C. 314 // 315 // Every input coordinate x[i] is assigned a symbolic perturbation dx[i]. 316 // We then compute the sign of the determinant of the perturbed points, 317 // i.e. 318 // | a.X+da.X a.Y+da.Y a.Z+da.Z | 319 // | b.X+db.X b.Y+db.Y b.Z+db.Z | 320 // | c.X+dc.X c.Y+dc.Y c.Z+dc.Z | 321 // 322 // The perturbations are chosen such that 323 // 324 // da.Z > da.Y > da.X > db.Z > db.Y > db.X > dc.Z > dc.Y > dc.X 325 // 326 // where each perturbation is so much smaller than the previous one that we 327 // don't even need to consider it unless the coefficients of all previous 328 // perturbations are zero. In fact, it is so small that we don't need to 329 // consider it unless the coefficient of all products of the previous 330 // perturbations are zero. For example, we don't need to consider the 331 // coefficient of db.Y unless the coefficient of db.Z *da.X is zero. 332 // 333 // The follow code simply enumerates the coefficients of the perturbations 334 // (and products of perturbations) that appear in the determinant above, in 335 // order of decreasing perturbation magnitude. The first non-zero 336 // coefficient determines the sign of the result. The easiest way to 337 // enumerate the coefficients in the correct order is to pretend that each 338 // perturbation is some tiny value "eps" raised to a power of two: 339 // 340 // eps** 1 2 4 8 16 32 64 128 256 341 // da.Z da.Y da.X db.Z db.Y db.X dc.Z dc.Y dc.X 342 // 343 // Essentially we can then just count in binary and test the corresponding 344 // subset of perturbations at each step. So for example, we must test the 345 // coefficient of db.Z*da.X before db.Y because eps**12 > eps**16. 346 // 347 // Of course, not all products of these perturbations appear in the 348 // determinant above, since the determinant only contains the products of 349 // elements in distinct rows and columns. Thus we don't need to consider 350 // da.Z*da.Y, db.Y *da.Y, etc. Furthermore, sometimes different pairs of 351 // perturbations have the same coefficient in the determinant; for example, 352 // da.Y*db.X and db.Y*da.X have the same coefficient (c.Z). Therefore 353 // we only need to test this coefficient the first time we encounter it in 354 // the binary order above (which will be db.Y*da.X). 355 // 356 // The sequence of tests below also appears in Table 4-ii of the paper 357 // referenced above, if you just want to look it up, with the following 358 // translations: [a,b,c] -> [i,j,k] and [0,1,2] -> [1,2,3]. Also note that 359 // some of the signs are different because the opposite cross product is 360 // used (e.g., B x C rather than C x B). 361 362 detSign := bCrossC.Z.Sign() // da.Z 363 if detSign != 0 { 364 return Direction(detSign) 365 } 366 detSign = bCrossC.Y.Sign() // da.Y 367 if detSign != 0 { 368 return Direction(detSign) 369 } 370 detSign = bCrossC.X.Sign() // da.X 371 if detSign != 0 { 372 return Direction(detSign) 373 } 374 375 detSign = newBigFloat().Sub(newBigFloat().Mul(c.X, a.Y), newBigFloat().Mul(c.Y, a.X)).Sign() // db.Z 376 if detSign != 0 { 377 return Direction(detSign) 378 } 379 detSign = c.X.Sign() // db.Z * da.Y 380 if detSign != 0 { 381 return Direction(detSign) 382 } 383 detSign = -(c.Y.Sign()) // db.Z * da.X 384 if detSign != 0 { 385 return Direction(detSign) 386 } 387 388 detSign = newBigFloat().Sub(newBigFloat().Mul(c.Z, a.X), newBigFloat().Mul(c.X, a.Z)).Sign() // db.Y 389 if detSign != 0 { 390 return Direction(detSign) 391 } 392 detSign = c.Z.Sign() // db.Y * da.X 393 if detSign != 0 { 394 return Direction(detSign) 395 } 396 397 // The following test is listed in the paper, but it is redundant because 398 // the previous tests guarantee that C == (0, 0, 0). 399 // (c.Y*a.Z - c.Z*a.Y).Sign() // db.X 400 401 detSign = newBigFloat().Sub(newBigFloat().Mul(a.X, b.Y), newBigFloat().Mul(a.Y, b.X)).Sign() // dc.Z 402 if detSign != 0 { 403 return Direction(detSign) 404 } 405 detSign = -(b.X.Sign()) // dc.Z * da.Y 406 if detSign != 0 { 407 return Direction(detSign) 408 } 409 detSign = b.Y.Sign() // dc.Z * da.X 410 if detSign != 0 { 411 return Direction(detSign) 412 } 413 detSign = a.X.Sign() // dc.Z * db.Y 414 if detSign != 0 { 415 return Direction(detSign) 416 } 417 return CounterClockwise // dc.Z * db.Y * da.X 418 } 419 420 // CompareDistances returns -1, 0, or +1 according to whether AX < BX, A == B, 421 // or AX > BX respectively. Distances are measured with respect to the positions 422 // of X, A, and B as though they were reprojected to lie exactly on the surface of 423 // the unit sphere. Furthermore, this method uses symbolic perturbations to 424 // ensure that the result is non-zero whenever A != B, even when AX == BX 425 // exactly, or even when A and B project to the same point on the sphere. 426 // Such results are guaranteed to be self-consistent, i.e. if AB < BC and 427 // BC < AC, then AB < AC. 428 func CompareDistances(x, a, b Point) int { 429 // We start by comparing distances using dot products (i.e., cosine of the 430 // angle), because (1) this is the cheapest technique, and (2) it is valid 431 // over the entire range of possible angles. (We can only use the sin^2 432 // technique if both angles are less than 90 degrees or both angles are 433 // greater than 90 degrees.) 434 sign := triageCompareCosDistances(x, a, b) 435 if sign != 0 { 436 return sign 437 } 438 439 // Optimization for (a == b) to avoid falling back to exact arithmetic. 440 if a == b { 441 return 0 442 } 443 444 // It is much better numerically to compare distances using cos(angle) if 445 // the distances are near 90 degrees and sin^2(angle) if the distances are 446 // near 0 or 180 degrees. We only need to check one of the two angles when 447 // making this decision because the fact that the test above failed means 448 // that angles "a" and "b" are very close together. 449 cosAX := a.Dot(x.Vector) 450 if cosAX > 1/math.Sqrt2 { 451 // Angles < 45 degrees. 452 sign = triageCompareSin2Distances(x, a, b) 453 } else if cosAX < -1/math.Sqrt2 { 454 // Angles > 135 degrees. sin^2(angle) is decreasing in this range. 455 sign = -triageCompareSin2Distances(x, a, b) 456 } 457 // C++ adds an additional check here using 80-bit floats. 458 // This is skipped in Go because we only have 32 and 64 bit floats. 459 460 if sign != 0 { 461 return sign 462 } 463 464 sign = exactCompareDistances(r3.PreciseVectorFromVector(x.Vector), r3.PreciseVectorFromVector(a.Vector), r3.PreciseVectorFromVector(b.Vector)) 465 if sign != 0 { 466 return sign 467 } 468 return symbolicCompareDistances(x, a, b) 469 } 470 471 // cosDistance returns cos(XY) where XY is the angle between X and Y, and the 472 // maximum error amount in the result. This requires X and Y be normalized. 473 func cosDistance(x, y Point) (cos, err float64) { 474 cos = x.Dot(y.Vector) 475 return cos, 9.5*dblError*math.Abs(cos) + 1.5*dblError 476 } 477 478 // sin2Distance returns sin**2(XY), where XY is the angle between X and Y, 479 // and the maximum error amount in the result. This requires X and Y be normalized. 480 func sin2Distance(x, y Point) (sin2, err float64) { 481 // The (x-y).Cross(x+y) trick eliminates almost all of error due to x 482 // and y being not quite unit length. This method is extremely accurate 483 // for small distances; the *relative* error in the result is O(dblError) for 484 // distances as small as dblError. 485 n := x.Sub(y.Vector).Cross(x.Add(y.Vector)) 486 sin2 = 0.25 * n.Norm2() 487 err = ((21+4*math.Sqrt(3))*dblError*sin2 + 488 32*math.Sqrt(3)*dblError*dblError*math.Sqrt(sin2) + 489 768*dblError*dblError*dblError*dblError) 490 return sin2, err 491 } 492 493 // triageCompareCosDistances returns -1, 0, or +1 according to whether AX < BX, 494 // A == B, or AX > BX by comparing the distances between them using cosDistance. 495 func triageCompareCosDistances(x, a, b Point) int { 496 cosAX, cosAXerror := cosDistance(a, x) 497 cosBX, cosBXerror := cosDistance(b, x) 498 diff := cosAX - cosBX 499 err := cosAXerror + cosBXerror 500 if diff > err { 501 return -1 502 } 503 if diff < -err { 504 return 1 505 } 506 return 0 507 } 508 509 // triageCompareSin2Distances returns -1, 0, or +1 according to whether AX < BX, 510 // A == B, or AX > BX by comparing the distances between them using sin2Distance. 511 func triageCompareSin2Distances(x, a, b Point) int { 512 sin2AX, sin2AXerror := sin2Distance(a, x) 513 sin2BX, sin2BXerror := sin2Distance(b, x) 514 diff := sin2AX - sin2BX 515 err := sin2AXerror + sin2BXerror 516 if diff > err { 517 return 1 518 } 519 if diff < -err { 520 return -1 521 } 522 return 0 523 } 524 525 // exactCompareDistances returns -1, 0, or 1 after comparing using the values as 526 // PreciseVectors. 527 func exactCompareDistances(x, a, b r3.PreciseVector) int { 528 // This code produces the same result as though all points were reprojected 529 // to lie exactly on the surface of the unit sphere. It is based on testing 530 // whether x.Dot(a.Normalize()) < x.Dot(b.Normalize()), reformulated 531 // so that it can be evaluated using exact arithmetic. 532 cosAX := x.Dot(a) 533 cosBX := x.Dot(b) 534 535 // If the two values have different signs, we need to handle that case now 536 // before squaring them below. 537 aSign := cosAX.Sign() 538 bSign := cosBX.Sign() 539 if aSign != bSign { 540 // If cos(AX) > cos(BX), then AX < BX. 541 if aSign > bSign { 542 return -1 543 } 544 return 1 545 } 546 cosAX2 := newBigFloat().Mul(cosAX, cosAX) 547 cosBX2 := newBigFloat().Mul(cosBX, cosBX) 548 cmp := newBigFloat().Sub(cosBX2.Mul(cosBX2, a.Norm2()), cosAX2.Mul(cosAX2, b.Norm2())) 549 return aSign * cmp.Sign() 550 } 551 552 // symbolicCompareDistances returns -1, 0, or +1 given three points such that AX == BX 553 // (exactly) according to whether AX < BX, AX == BX, or AX > BX after symbolic 554 // perturbations are taken into account. 555 func symbolicCompareDistances(x, a, b Point) int { 556 // Our symbolic perturbation strategy is based on the following model. 557 // Similar to "simulation of simplicity", we assign a perturbation to every 558 // point such that if A < B, then the symbolic perturbation for A is much, 559 // much larger than the symbolic perturbation for B. We imagine that 560 // rather than projecting every point to lie exactly on the unit sphere, 561 // instead each point is positioned on its own tiny pedestal that raises it 562 // just off the surface of the unit sphere. This means that the distance AX 563 // is actually the true distance AX plus the (symbolic) heights of the 564 // pedestals for A and X. The pedestals are infinitesmally thin, so they do 565 // not affect distance measurements except at the two endpoints. If several 566 // points project to exactly the same point on the unit sphere, we imagine 567 // that they are placed on separate pedestals placed close together, where 568 // the distance between pedestals is much, much less than the height of any 569 // pedestal. (There are a finite number of Points, and therefore a finite 570 // number of pedestals, so this is possible.) 571 // 572 // If A < B, then A is on a higher pedestal than B, and therefore AX > BX. 573 switch a.Cmp(b.Vector) { 574 case -1: 575 return 1 576 case 1: 577 return -1 578 default: 579 return 0 580 } 581 } 582 583 var ( 584 // ca45Degrees is a predefined ChordAngle representing (approximately) 45 degrees. 585 ca45Degrees = s1.ChordAngleFromSquaredLength(2 - math.Sqrt2) 586 ) 587 588 // CompareDistance returns -1, 0, or +1 according to whether the distance XY is 589 // respectively less than, equal to, or greater than the provided chord angle. Distances are measured 590 // with respect to the positions of all points as though they are projected to lie 591 // exactly on the surface of the unit sphere. 592 func CompareDistance(x, y Point, r s1.ChordAngle) int { 593 // As with CompareDistances, we start by comparing dot products because 594 // the sin^2 method is only valid when the distance XY and the limit "r" are 595 // both less than 90 degrees. 596 sign := triageCompareCosDistance(x, y, float64(r)) 597 if sign != 0 { 598 return sign 599 } 600 601 // Unlike with CompareDistances, it's not worth using the sin^2 method 602 // when the distance limit is near 180 degrees because the ChordAngle 603 // representation itself has has a rounding error of up to 2e-8 radians for 604 // distances near 180 degrees. 605 if r < ca45Degrees { 606 sign = triageCompareSin2Distance(x, y, float64(r)) 607 if sign != 0 { 608 return sign 609 } 610 } 611 return exactCompareDistance(r3.PreciseVectorFromVector(x.Vector), r3.PreciseVectorFromVector(y.Vector), big.NewFloat(float64(r)).SetPrec(big.MaxPrec)) 612 } 613 614 // triageCompareCosDistance returns -1, 0, or +1 according to whether the distance XY is 615 // less than, equal to, or greater than r2 respectively using cos distance. 616 func triageCompareCosDistance(x, y Point, r2 float64) int { 617 cosXY, cosXYError := cosDistance(x, y) 618 cosR := 1.0 - 0.5*r2 619 cosRError := 2.0 * dblError * cosR 620 diff := cosXY - cosR 621 err := cosXYError + cosRError 622 if diff > err { 623 return -1 624 } 625 if diff < -err { 626 return 1 627 } 628 return 0 629 } 630 631 // triageCompareSin2Distance returns -1, 0, or +1 according to whether the distance XY is 632 // less than, equal to, or greater than r2 respectively using sin^2 distance. 633 func triageCompareSin2Distance(x, y Point, r2 float64) int { 634 // Only valid for distance limits < 90 degrees. 635 sin2XY, sin2XYError := sin2Distance(x, y) 636 sin2R := r2 * (1.0 - 0.25*r2) 637 sin2RError := 3.0 * dblError * sin2R 638 diff := sin2XY - sin2R 639 err := sin2XYError + sin2RError 640 if diff > err { 641 return 1 642 } 643 if diff < -err { 644 return -1 645 } 646 return 0 647 } 648 649 var ( 650 bigOne = big.NewFloat(1.0).SetPrec(big.MaxPrec) 651 bigHalf = big.NewFloat(0.5).SetPrec(big.MaxPrec) 652 ) 653 654 // exactCompareDistance returns -1, 0, or +1 after comparing using PreciseVectors. 655 func exactCompareDistance(x, y r3.PreciseVector, r2 *big.Float) int { 656 // This code produces the same result as though all points were reprojected 657 // to lie exactly on the surface of the unit sphere. It is based on 658 // comparing the cosine of the angle XY (when both points are projected to 659 // lie exactly on the sphere) to the given threshold. 660 cosXY := x.Dot(y) 661 cosR := newBigFloat().Sub(bigOne, newBigFloat().Mul(bigHalf, r2)) 662 663 // If the two values have different signs, we need to handle that case now 664 // before squaring them below. 665 xySign := cosXY.Sign() 666 rSign := cosR.Sign() 667 if xySign != rSign { 668 if xySign > rSign { 669 return -1 670 } 671 return 1 // If cos(XY) > cos(r), then XY < r. 672 } 673 cmp := newBigFloat().Sub( 674 newBigFloat().Mul( 675 newBigFloat().Mul(cosR, cosR), newBigFloat().Mul(x.Norm2(), y.Norm2())), 676 newBigFloat().Mul(cosXY, cosXY)) 677 return xySign * cmp.Sign() 678 } 679 680 // TODO(roberts): Differences from C++ 681 // CompareEdgeDistance 682 // CompareEdgeDirections 683 // EdgeCircumcenterSign 684 // GetVoronoiSiteExclusion 685 // GetClosestVertex 686 // TriageCompareLineSin2Distance 687 // TriageCompareLineCos2Distance 688 // TriageCompareLineDistance 689 // TriageCompareEdgeDistance 690 // ExactCompareLineDistance 691 // ExactCompareEdgeDistance 692 // TriageCompareEdgeDirections 693 // ExactCompareEdgeDirections 694 // ArePointsAntipodal 695 // ArePointsLinearlyDependent 696 // GetCircumcenter 697 // TriageEdgeCircumcenterSign 698 // ExactEdgeCircumcenterSign 699 // UnperturbedSign 700 // SymbolicEdgeCircumcenterSign 701 // ExactVoronoiSiteExclusion 702