1 // Copyright 2014 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 import ( 18 "flag" 19 "io" 20 "math" 21 "math/rand" 22 "os" 23 24 "github.com/golang/geo/r1" 25 "github.com/golang/geo/r2" 26 "github.com/golang/geo/s1" 27 ) 28 29 var ( 30 // To set in testing add "--benchmark_brute_force=true" to your test command. 31 benchmarkBruteForce = flag.Bool("benchmark_brute_force", false, 32 "When set, use brute force algorithms in benchmarking.") 33 ) 34 35 // float64Eq reports whether the two values are within the default epsilon. 36 func float64Eq(x, y float64) bool { return float64Near(x, y, epsilon) } 37 38 // float64Near reports whether the two values are within the given epsilon. 39 func float64Near(x, y, ε float64) bool { 40 return math.Abs(x-y) <= ε 41 } 42 43 // TODO(roberts): Add in flag to allow specifying the random seed for repeatable tests. 44 45 // The Earth's mean radius in kilometers (according to NASA). 46 const earthRadiusKm = 6371.01 47 48 // kmToAngle converts a distance on the Earth's surface to an angle. 49 func kmToAngle(km float64) s1.Angle { 50 return s1.Angle(km / earthRadiusKm) 51 } 52 53 // randomBits returns a 64-bit random unsigned integer whose lowest "num" are random, and 54 // whose other bits are zero. 55 func randomBits(num uint32) uint64 { 56 // Make sure the request is for not more than 63 bits. 57 if num > 63 { 58 num = 63 59 } 60 return uint64(rand.Int63()) & ((1 << num) - 1) 61 } 62 63 // Return a uniformly distributed 64-bit unsigned integer. 64 func randomUint64() uint64 { 65 return uint64(rand.Int63() | (rand.Int63() << 63)) 66 } 67 68 // Return a uniformly distributed 32-bit unsigned integer. 69 func randomUint32() uint32 { 70 return uint32(randomBits(32)) 71 } 72 73 // randomFloat64 returns a uniformly distributed value in the range [0,1). 74 // Note that the values returned are all multiples of 2**-53, which means that 75 // not all possible values in this range are returned. 76 func randomFloat64() float64 { 77 const randomFloatBits = 53 78 return math.Ldexp(float64(randomBits(randomFloatBits)), -randomFloatBits) 79 } 80 81 // randomUniformInt returns a uniformly distributed integer in the range [0,n). 82 // NOTE: This is replicated here to stay in sync with how the C++ code generates 83 // uniform randoms. (instead of using Go's math/rand package directly). 84 func randomUniformInt(n int) int { 85 return int(randomFloat64() * float64(n)) 86 } 87 88 // randomUniformFloat64 returns a uniformly distributed value in the range [min, max). 89 func randomUniformFloat64(min, max float64) float64 { 90 return min + randomFloat64()*(max-min) 91 } 92 93 // oneIn returns true with a probability of 1/n. 94 func oneIn(n int) bool { 95 return randomUniformInt(n) == 0 96 } 97 98 // randomPoint returns a random unit-length vector. 99 func randomPoint() Point { 100 return PointFromCoords(randomUniformFloat64(-1, 1), 101 randomUniformFloat64(-1, 1), randomUniformFloat64(-1, 1)) 102 } 103 104 // randomFrame returns a right-handed coordinate frame (three orthonormal vectors) for 105 // a randomly generated point. 106 func randomFrame() *matrix3x3 { 107 return randomFrameAtPoint(randomPoint()) 108 } 109 110 // randomFrameAtPoint returns a right-handed coordinate frame using the given 111 // point as the z-axis. The x- and y-axes are computed such that (x,y,z) is a 112 // right-handed coordinate frame (three orthonormal vectors). 113 func randomFrameAtPoint(z Point) *matrix3x3 { 114 x := Point{z.Cross(randomPoint().Vector).Normalize()} 115 y := Point{z.Cross(x.Vector).Normalize()} 116 117 m := &matrix3x3{} 118 m.setCol(0, x) 119 m.setCol(1, y) 120 m.setCol(2, z) 121 return m 122 } 123 124 // randomCellIDForLevel returns a random CellID at the given level. 125 // The distribution is uniform over the space of cell ids, but only 126 // approximately uniform over the surface of the sphere. 127 func randomCellIDForLevel(level int) CellID { 128 face := randomUniformInt(NumFaces) 129 pos := randomUint64() & uint64((1<<PosBits)-1) 130 return CellIDFromFacePosLevel(face, pos, level) 131 } 132 133 // randomCellID returns a random CellID at a randomly chosen 134 // level. The distribution is uniform over the space of cell ids, 135 // but only approximately uniform over the surface of the sphere. 136 func randomCellID() CellID { 137 return randomCellIDForLevel(randomUniformInt(MaxLevel + 1)) 138 } 139 140 // randomCellUnion returns a CellUnion of the given size of randomly selected cells. 141 func randomCellUnion(n int) CellUnion { 142 var cu CellUnion 143 for i := 0; i < n; i++ { 144 cu = append(cu, randomCellID()) 145 } 146 return cu 147 } 148 149 // concentricLoopsPolygon constructs a polygon with the specified center as a 150 // number of concentric loops and vertices per loop. 151 func concentricLoopsPolygon(center Point, numLoops, verticesPerLoop int) *Polygon { 152 var loops []*Loop 153 for li := 0; li < numLoops; li++ { 154 radius := s1.Angle(0.005 * float64(li+1) / float64(numLoops)) 155 loops = append(loops, RegularLoop(center, radius, verticesPerLoop)) 156 } 157 return PolygonFromLoops(loops) 158 } 159 160 // skewedInt returns a number in the range [0,2^max_log-1] with bias towards smaller numbers. 161 func skewedInt(maxLog int) int { 162 base := uint32(rand.Int31n(int32(maxLog + 1))) 163 return int(randomBits(31) & ((1 << base) - 1)) 164 } 165 166 // randomCap returns a cap with a random axis such that the log of its area is 167 // uniformly distributed between the logs of the two given values. The log of 168 // the cap angle is also approximately uniformly distributed. 169 func randomCap(minArea, maxArea float64) Cap { 170 capArea := maxArea * math.Pow(minArea/maxArea, randomFloat64()) 171 return CapFromCenterArea(randomPoint(), capArea) 172 } 173 174 // latLngsApproxEqual reports of the two LatLngs are within the given epsilon. 175 func latLngsApproxEqual(a, b LatLng, epsilon float64) bool { 176 return float64Near(float64(a.Lat), float64(b.Lat), epsilon) && 177 float64Near(float64(a.Lng), float64(b.Lng), epsilon) 178 } 179 180 // pointsApproxEqual reports whether the two points are within the given distance 181 // of each other. This is the same as Point.ApproxEqual but permits specifying 182 // the epsilon. 183 func pointsApproxEqual(a, b Point, epsilon float64) bool { 184 return float64(a.Vector.Angle(b.Vector)) <= epsilon 185 } 186 187 // pointSlicesApproxEqual reports whether corresponding elements of each slice are approximately equal. 188 func pointSlicesApproxEqual(a, b []Point, epsilon float64) bool { 189 if len(a) != len(b) { 190 return false 191 } 192 193 for i := range a { 194 if !pointsApproxEqual(a[i], b[i], epsilon) { 195 return false 196 } 197 } 198 return true 199 } 200 201 // r1IntervalsApproxEqual reports whether the two r1.Intervals are within the given 202 // epsilon of each other. This adds a test changeable value for epsilon 203 func r1IntervalsApproxEqual(a, b r1.Interval, epsilon float64) bool { 204 if a.IsEmpty() { 205 return b.Length() <= 2*epsilon 206 } 207 if b.IsEmpty() { 208 return a.Length() <= 2*epsilon 209 } 210 return math.Abs(b.Lo-a.Lo) <= epsilon && 211 math.Abs(b.Hi-a.Hi) <= epsilon 212 } 213 214 var ( 215 rectErrorLat = 10 * dblEpsilon 216 rectErrorLng = dblEpsilon 217 ) 218 219 // r2PointsApproxEqual reports whether the two points are within the given epsilon. 220 func r2PointsApproxEqual(a, b r2.Point, epsilon float64) bool { 221 return float64Near(a.X, b.X, epsilon) && float64Near(a.Y, b.Y, epsilon) 222 } 223 224 // r2PointSlicesApproxEqual reports whether corresponding elements of the slices are approximately equal. 225 func r2PointSlicesApproxEqual(a, b []r2.Point, epsilon float64) bool { 226 if len(a) != len(b) { 227 return false 228 } 229 230 for i := range a { 231 if !r2PointsApproxEqual(a[i], b[i], epsilon) { 232 return false 233 } 234 } 235 return true 236 } 237 238 // rectsApproxEqual reports whether the two rect are within the given tolerances 239 // at each corner from each other. The tolerances are specific to each axis. 240 func rectsApproxEqual(a, b Rect, tolLat, tolLng float64) bool { 241 return math.Abs(a.Lat.Lo-b.Lat.Lo) < tolLat && 242 math.Abs(a.Lat.Hi-b.Lat.Hi) < tolLat && 243 math.Abs(a.Lng.Lo-b.Lng.Lo) < tolLng && 244 math.Abs(a.Lng.Hi-b.Lng.Hi) < tolLng 245 } 246 247 // matricesApproxEqual reports whether all cells in both matrices are equal within 248 // the default floating point epsilon. 249 func matricesApproxEqual(m1, m2 *matrix3x3) bool { 250 return float64Eq(m1[0][0], m2[0][0]) && 251 float64Eq(m1[0][1], m2[0][1]) && 252 float64Eq(m1[0][2], m2[0][2]) && 253 254 float64Eq(m1[1][0], m2[1][0]) && 255 float64Eq(m1[1][1], m2[1][1]) && 256 float64Eq(m1[1][2], m2[1][2]) && 257 258 float64Eq(m1[2][0], m2[2][0]) && 259 float64Eq(m1[2][1], m2[2][1]) && 260 float64Eq(m1[2][2], m2[2][2]) 261 } 262 263 // samplePointFromRect returns a point chosen uniformly at random (with respect 264 // to area on the sphere) from the given rectangle. 265 func samplePointFromRect(rect Rect) Point { 266 // First choose a latitude uniformly with respect to area on the sphere. 267 sinLo := math.Sin(rect.Lat.Lo) 268 sinHi := math.Sin(rect.Lat.Hi) 269 lat := math.Asin(randomUniformFloat64(sinLo, sinHi)) 270 271 // Now choose longitude uniformly within the given range. 272 lng := rect.Lng.Lo + randomFloat64()*rect.Lng.Length() 273 274 return PointFromLatLng(LatLng{s1.Angle(lat), s1.Angle(lng)}.Normalized()) 275 } 276 277 // samplePointFromCap returns a point chosen uniformly at random (with respect 278 // to area) from the given cap. 279 func samplePointFromCap(c Cap) Point { 280 // We consider the cap axis to be the "z" axis. We choose two other axes to 281 // complete the coordinate frame. 282 m := getFrame(c.Center()) 283 284 // The surface area of a spherical cap is directly proportional to its 285 // height. First we choose a random height, and then we choose a random 286 // point along the circle at that height. 287 h := randomFloat64() * c.Height() 288 theta := 2 * math.Pi * randomFloat64() 289 r := math.Sqrt(h * (2 - h)) 290 291 // The result should already be very close to unit-length, but we might as 292 // well make it accurate as possible. 293 return Point{fromFrame(m, PointFromCoords(math.Cos(theta)*r, math.Sin(theta)*r, 1-h)).Normalize()} 294 } 295 296 // perturbATowardsB returns a point that has been shifted some distance towards the 297 // second point based on a random number. 298 func perturbATowardsB(a, b Point) Point { 299 choice := randomFloat64() 300 if choice < 0.1 { 301 return a 302 } 303 if choice < 0.3 { 304 // Return a point that is exactly proportional to A and that still 305 // satisfies IsUnitLength(). 306 for { 307 b := Point{a.Mul(2 - a.Norm() + 5*(randomFloat64()-0.5)*dblEpsilon)} 308 if !b.ApproxEqual(a) && b.IsUnit() { 309 return b 310 } 311 } 312 } 313 if choice < 0.5 { 314 // Return a point such that the distance squared to A will underflow. 315 return InterpolateAtDistance(1e-300, a, b) 316 } 317 // Otherwise return a point whose distance from A is near dblEpsilon such 318 // that the log of the pdf is uniformly distributed. 319 distance := dblEpsilon * 1e-5 * math.Pow(1e6, randomFloat64()) 320 return InterpolateAtDistance(s1.Angle(distance), a, b) 321 } 322 323 // perturbedCornerOrMidpoint returns a Point from a line segment whose endpoints are 324 // difficult to handle correctly. Given two adjacent cube vertices P and Q, 325 // it returns either an edge midpoint, face midpoint, or corner vertex that is 326 // in the plane of PQ and that has been perturbed slightly. It also sometimes 327 // returns a random point from anywhere on the sphere. 328 func perturbedCornerOrMidpoint(p, q Point) Point { 329 a := p.Mul(float64(randomUniformInt(3) - 1)).Add(q.Mul(float64(randomUniformInt(3) - 1))) 330 if oneIn(10) { 331 // This perturbation often has no effect except on coordinates that are 332 // zero, in which case the perturbed value is so small that operations on 333 // it often result in underflow. 334 a = a.Add(randomPoint().Mul(math.Pow(1e-300, randomFloat64()))) 335 } else if oneIn(2) { 336 // For coordinates near 1 (say > 0.5), this perturbation yields values 337 // that are only a few representable values away from the initial value. 338 a = a.Add(randomPoint().Mul(4 * dblEpsilon)) 339 } else { 340 // A perturbation whose magnitude is in the range [1e-25, 1e-10]. 341 a = a.Add(randomPoint().Mul(1e-10 * math.Pow(1e-15, randomFloat64()))) 342 } 343 344 if a.Norm2() < math.SmallestNonzeroFloat64 { 345 // If a.Norm2() is denormalized, Normalize() loses too much precision. 346 return perturbedCornerOrMidpoint(p, q) 347 } 348 return Point{a} 349 } 350 351 // readLoops returns a slice of Loops read from a file encoded using Loops Encode. 352 func readLoops(filename string) ([]*Loop, error) { 353 f, err := os.Open(filename) 354 if err != nil { 355 return nil, err 356 } 357 defer f.Close() 358 359 var loops []*Loop 360 361 // Test loop files are expected to be a direct concatenation of encoded loops with 362 // no separator tokens. Because there is no way of knowing a priori how many items 363 // or how many bytes ahead of time, the only way to end the loop is when we hit EOF. 364 for { 365 l := &Loop{} 366 err := l.Decode(f) 367 if err == io.EOF { 368 break 369 } 370 if err != nil { 371 return nil, err 372 } 373 loops = append(loops, l) 374 } 375 376 return loops, nil 377 } 378 379 // fractal is a simple type that generates "Koch snowflake" fractals (see Wikipedia 380 // for an introduction). There is an option to control the fractal dimension 381 // (between 1.0 and 2.0); values between 1.02 and 1.50 are reasonable simulations 382 // of various coastlines. The default dimension (about 1.26) corresponds to the 383 // standard Koch snowflake. (The west coast of Britain has a fractal dimension 384 // of approximately 1.25.) 385 // 386 // The fractal is obtained by starting with an equilateral triangle and 387 // recursively subdividing each edge into four segments of equal length. 388 // Therefore the shape at level "n" consists of 3*(4**n) edges. Multi-level 389 // fractals are also supported: if you set MinLevel to a non-negative 390 // value, then the recursive subdivision has an equal probability of 391 // stopping at any of the levels between the given min and max (inclusive). 392 // This yields a fractal where the perimeter of the original triangle is 393 // approximately equally divided between fractals at the various possible 394 // levels. If there are k distinct levels {min,..,max}, the expected number 395 // of edges at each level "i" is approximately 3*(4**i)/k. 396 type fractal struct { 397 maxLevel int 398 // minLevel defines the minimum subdivision level of the fractal 399 // A min level of 0 should be avoided since this creates a 400 // significant chance that none of the edges will end up subdivided. 401 minLevel int 402 403 // dimension of the fractal. A value of approximately 1.26 corresponds 404 // to the stardard Koch curve. The value must lie in the range [1.0, 2.0). 405 dimension float64 406 407 // The ratio of the sub-edge length to the original edge length at each 408 // subdivision step. 409 edgeFraction float64 410 411 // The distance from the original edge to the middle vertex at each 412 // subdivision step, as a fraction of the original edge length. 413 offsetFraction float64 414 } 415 416 // newFractal returns a new instance of the fractal type with appropriate defaults. 417 func newFractal() *fractal { 418 return &fractal{ 419 maxLevel: -1, 420 minLevel: -1, 421 dimension: math.Log(4) / math.Log(3), // standard Koch curve 422 edgeFraction: 0, 423 offsetFraction: 0, 424 } 425 } 426 427 // setLevelForApproxMinEdges sets the min level to produce approximately the 428 // given number of edges. The values are rounded to a nearby value of 3*(4**n). 429 func (f *fractal) setLevelForApproxMinEdges(minEdges int) { 430 // Map values in the range [3*(4**n)/2, 3*(4**n)*2) to level n. 431 f.minLevel = int(math.Round(0.5 * math.Log2(float64(minEdges)/3))) 432 } 433 434 // setLevelForApproxMaxEdges sets the max level to produce approximately the 435 // given number of edges. The values are rounded to a nearby value of 3*(4**n). 436 func (f *fractal) setLevelForApproxMaxEdges(maxEdges int) { 437 // Map values in the range [3*(4**n)/2, 3*(4**n)*2) to level n. 438 f.maxLevel = int(math.Round(0.5 * math.Log2(float64(maxEdges)/3))) 439 } 440 441 // minRadiusFactor returns a lower bound on the ratio (Rmin / R), where "R" is the 442 // radius passed to makeLoop and Rmin is the minimum distance from the 443 // fractal boundary to its center, where all distances are measured in the 444 // tangent plane at the fractal's center. This can be used to inscribe 445 // another geometric figure within the fractal without intersection. 446 func (f *fractal) minRadiusFactor() float64 { 447 // The minimum radius is attained at one of the vertices created by the 448 // first subdivision step as long as the dimension is not too small (at 449 // least minDimensionForMinRadiusAtLevel1, see below). Otherwise we fall 450 // back on the incircle radius of the original triangle, which is always a 451 // lower bound (and is attained when dimension = 1). 452 // 453 // The value below was obtained by letting AE be an original triangle edge, 454 // letting ABCDE be the corresponding polyline after one subdivision step, 455 // and then letting BC be tangent to the inscribed circle at the center of 456 // the fractal O. This gives rise to a pair of similar triangles whose edge 457 // length ratios can be used to solve for the corresponding "edge fraction". 458 // This method is slightly conservative because it is computed using planar 459 // rather than spherical geometry. The value below is equal to 460 // -math.Log(4)/math.Log((2 + math.Cbrt(2) - math.Cbrt(4))/6). 461 const minDimensionForMinRadiusAtLevel1 = 1.0852230903040407 462 if f.dimension >= minDimensionForMinRadiusAtLevel1 { 463 return math.Sqrt(1 + 3*f.edgeFraction*(f.edgeFraction-1)) 464 } 465 return 0.5 466 } 467 468 // maxRadiusFactor returns the ratio (Rmax / R), where "R" is the radius passed 469 // to makeLoop and Rmax is the maximum distance from the fractal boundary 470 // to its center, where all distances are measured in the tangent plane at 471 // the fractal's center. This can be used to inscribe the fractal within 472 // some other geometric figure without intersection. 473 func (f *fractal) maxRadiusFactor() float64 { 474 // The maximum radius is always attained at either an original triangle 475 // vertex or at a middle vertex from the first subdivision step. 476 return math.Max(1.0, f.offsetFraction*math.Sqrt(3)+0.5) 477 } 478 479 // r2VerticesHelper recursively subdivides the edge to the desired level given 480 // the two endpoints (v0,v4) of an edge, and returns all vertices of the resulting 481 // curve up to but not including the endpoint v4. 482 func (f *fractal) r2VerticesHelper(v0, v4 r2.Point, level int) []r2.Point { 483 if level >= f.minLevel && oneIn(f.maxLevel-level+1) { 484 // stop at this level 485 return []r2.Point{v0} 486 } 487 488 var vertices []r2.Point 489 490 // Otherwise compute the intermediate vertices v1, v2, and v3. 491 dir := v4.Sub(v0) 492 v1 := v0.Add(dir.Mul(f.edgeFraction)) 493 v2 := v0.Add(v4).Mul(0.5).Sub(dir.Ortho().Mul(f.offsetFraction)) 494 v3 := v4.Sub(dir.Mul(f.edgeFraction)) 495 496 // And recurse on the four sub-edges. 497 vertices = append(vertices, f.r2VerticesHelper(v0, v1, level+1)...) 498 vertices = append(vertices, f.r2VerticesHelper(v1, v2, level+1)...) 499 vertices = append(vertices, f.r2VerticesHelper(v2, v3, level+1)...) 500 vertices = append(vertices, f.r2VerticesHelper(v3, v4, level+1)...) 501 502 return vertices 503 } 504 505 // generateR2Vertices returns the set of r2 plane vertices for the fractal 506 // based on its current settings. 507 func (f *fractal) generateR2Vertices() []r2.Point { 508 var vertices []r2.Point 509 510 // The Koch "snowflake" consists of three Koch curves whose initial edges 511 // form an equilateral triangle. 512 v0 := r2.Point{1.0, 0.0} 513 v1 := r2.Point{-0.5, math.Sqrt(3) / 2} 514 v2 := r2.Point{-0.5, -math.Sqrt(3) / 2} 515 vertices = append(vertices, f.r2VerticesHelper(v0, v1, 0)...) 516 vertices = append(vertices, f.r2VerticesHelper(v1, v2, 0)...) 517 vertices = append(vertices, f.r2VerticesHelper(v2, v0, 0)...) 518 519 return vertices 520 } 521 522 // makeLoop returns a fractal loop centered around the z-axis of the given 523 // coordinate frame, with the first vertex in the direction of the 524 // positive x-axis. In order to avoid self-intersections, the fractal is 525 // generated by first drawing it in a 2D tangent plane to the unit sphere 526 // (touching at the fractal's center point) and then projecting the edges 527 // onto the sphere. This has the side effect of shrinking the fractal 528 // slightly compared to its nominal radius. 529 func (f *fractal) makeLoop(frame *matrix3x3, nominalRadius s1.Angle) *Loop { 530 // update dependent values before making the loop. 531 if f.minLevel < 0 || f.minLevel > f.maxLevel { 532 f.minLevel = f.maxLevel 533 } 534 535 f.edgeFraction = math.Pow(4.0, -1.0/f.dimension) 536 f.offsetFraction = math.Sqrt(f.edgeFraction - 0.25) 537 538 r2pts := f.generateR2Vertices() 539 verts := make([]Point, 0, len(r2pts)) 540 r := nominalRadius.Radians() 541 542 for _, pt := range r2pts { 543 verts = append(verts, fromFrame(*frame, PointFromCoords(pt.X*r, pt.Y*r, 1))) 544 } 545 return LoopFromPoints(verts) 546 } 547 548 // TODO(roberts): Items remaining to port: 549 // CheckCovering 550 // CheckResultSet 551 // CheckDistanceResults 552