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 "fmt" 19 "io" 20 "sort" 21 22 "github.com/golang/geo/s1" 23 ) 24 25 // A CellUnion is a collection of CellIDs. 26 // 27 // It is normalized if it is sorted, and does not contain redundancy. 28 // Specifically, it may not contain the same CellID twice, nor a CellID that 29 // is contained by another, nor the four sibling CellIDs that are children of 30 // a single higher level CellID. 31 // 32 // CellUnions are not required to be normalized, but certain operations will 33 // return different results if they are not (e.g. Contains). 34 type CellUnion []CellID 35 36 // CellUnionFromRange creates a CellUnion that covers the half-open range 37 // of leaf cells [begin, end). If begin == end the resulting union is empty. 38 // This requires that begin and end are both leaves, and begin <= end. 39 // To create a closed-ended range, pass in end.Next(). 40 func CellUnionFromRange(begin, end CellID) CellUnion { 41 // We repeatedly add the largest cell we can. 42 var cu CellUnion 43 for id := begin.MaxTile(end); id != end; id = id.Next().MaxTile(end) { 44 cu = append(cu, id) 45 } 46 // The output is normalized because the cells are added in order by the iteration. 47 return cu 48 } 49 50 // CellUnionFromUnion creates a CellUnion from the union of the given CellUnions. 51 func CellUnionFromUnion(cellUnions ...CellUnion) CellUnion { 52 var cu CellUnion 53 for _, cellUnion := range cellUnions { 54 cu = append(cu, cellUnion...) 55 } 56 cu.Normalize() 57 return cu 58 } 59 60 // CellUnionFromIntersection creates a CellUnion from the intersection of the given CellUnions. 61 func CellUnionFromIntersection(x, y CellUnion) CellUnion { 62 var cu CellUnion 63 64 // This is a fairly efficient calculation that uses binary search to skip 65 // over sections of both input vectors. It takes constant time if all the 66 // cells of x come before or after all the cells of y in CellID order. 67 var i, j int 68 for i < len(x) && j < len(y) { 69 iMin := x[i].RangeMin() 70 jMin := y[j].RangeMin() 71 if iMin > jMin { 72 // Either j.Contains(i) or the two cells are disjoint. 73 if x[i] <= y[j].RangeMax() { 74 cu = append(cu, x[i]) 75 i++ 76 } else { 77 // Advance j to the first cell possibly contained by x[i]. 78 j = y.lowerBound(j+1, len(y), iMin) 79 // The previous cell y[j-1] may now contain x[i]. 80 if x[i] <= y[j-1].RangeMax() { 81 j-- 82 } 83 } 84 } else if jMin > iMin { 85 // Identical to the code above with i and j reversed. 86 if y[j] <= x[i].RangeMax() { 87 cu = append(cu, y[j]) 88 j++ 89 } else { 90 i = x.lowerBound(i+1, len(x), jMin) 91 if y[j] <= x[i-1].RangeMax() { 92 i-- 93 } 94 } 95 } else { 96 // i and j have the same RangeMin(), so one contains the other. 97 if x[i] < y[j] { 98 cu = append(cu, x[i]) 99 i++ 100 } else { 101 cu = append(cu, y[j]) 102 j++ 103 } 104 } 105 } 106 107 // The output is generated in sorted order. 108 cu.Normalize() 109 return cu 110 } 111 112 // CellUnionFromIntersectionWithCellID creates a CellUnion from the intersection 113 // of a CellUnion with the given CellID. This can be useful for splitting a 114 // CellUnion into chunks. 115 func CellUnionFromIntersectionWithCellID(x CellUnion, id CellID) CellUnion { 116 var cu CellUnion 117 if x.ContainsCellID(id) { 118 cu = append(cu, id) 119 cu.Normalize() 120 return cu 121 } 122 123 idmax := id.RangeMax() 124 for i := x.lowerBound(0, len(x), id.RangeMin()); i < len(x) && x[i] <= idmax; i++ { 125 cu = append(cu, x[i]) 126 } 127 128 cu.Normalize() 129 return cu 130 } 131 132 // CellUnionFromDifference creates a CellUnion from the difference (x - y) 133 // of the given CellUnions. 134 func CellUnionFromDifference(x, y CellUnion) CellUnion { 135 // TODO(roberts): This is approximately O(N*log(N)), but could probably 136 // use similar techniques as CellUnionFromIntersectionWithCellID to be more efficient. 137 138 var cu CellUnion 139 for _, xid := range x { 140 cu.cellUnionDifferenceInternal(xid, &y) 141 } 142 143 // The output is generated in sorted order, and there should not be any 144 // cells that can be merged (provided that both inputs were normalized). 145 return cu 146 } 147 148 // The C++ constructor methods FromNormalized and FromVerbatim are not necessary 149 // since they don't call Normalize, and just set the CellIDs directly on the object, 150 // so straight casting is sufficient in Go to replicate this behavior. 151 152 // IsValid reports whether the cell union is valid, meaning that the CellIDs are 153 // valid, non-overlapping, and sorted in increasing order. 154 func (cu *CellUnion) IsValid() bool { 155 for i, cid := range *cu { 156 if !cid.IsValid() { 157 return false 158 } 159 if i == 0 { 160 continue 161 } 162 if (*cu)[i-1].RangeMax() >= cid.RangeMin() { 163 return false 164 } 165 } 166 return true 167 } 168 169 // IsNormalized reports whether the cell union is normalized, meaning that it is 170 // satisfies IsValid and that no four cells have a common parent. 171 // Certain operations such as Contains will return a different 172 // result if the cell union is not normalized. 173 func (cu *CellUnion) IsNormalized() bool { 174 for i, cid := range *cu { 175 if !cid.IsValid() { 176 return false 177 } 178 if i == 0 { 179 continue 180 } 181 if (*cu)[i-1].RangeMax() >= cid.RangeMin() { 182 return false 183 } 184 if i < 3 { 185 continue 186 } 187 if areSiblings((*cu)[i-3], (*cu)[i-2], (*cu)[i-1], cid) { 188 return false 189 } 190 } 191 return true 192 } 193 194 // Normalize normalizes the CellUnion. 195 func (cu *CellUnion) Normalize() { 196 sortCellIDs(*cu) 197 198 output := make([]CellID, 0, len(*cu)) // the list of accepted cells 199 // Loop invariant: output is a sorted list of cells with no redundancy. 200 for _, ci := range *cu { 201 // The first two passes here either ignore this new candidate, 202 // or remove previously accepted cells that are covered by this candidate. 203 204 // Ignore this cell if it is contained by the previous one. 205 // We only need to check the last accepted cell. The ordering of the 206 // cells implies containment (but not the converse), and output has no redundancy, 207 // so if this candidate is not contained by the last accepted cell 208 // then it cannot be contained by any previously accepted cell. 209 if len(output) > 0 && output[len(output)-1].Contains(ci) { 210 continue 211 } 212 213 // Discard any previously accepted cells contained by this one. 214 // This could be any contiguous trailing subsequence, but it can't be 215 // a discontiguous subsequence because of the containment property of 216 // sorted S2 cells mentioned above. 217 j := len(output) - 1 // last index to keep 218 for j >= 0 { 219 if !ci.Contains(output[j]) { 220 break 221 } 222 j-- 223 } 224 output = output[:j+1] 225 226 // See if the last three cells plus this one can be collapsed. 227 // We loop because collapsing three accepted cells and adding a higher level cell 228 // could cascade into previously accepted cells. 229 for len(output) >= 3 && areSiblings(output[len(output)-3], output[len(output)-2], output[len(output)-1], ci) { 230 // Replace four children by their parent cell. 231 output = output[:len(output)-3] 232 ci = ci.immediateParent() // checked !ci.isFace above 233 } 234 output = append(output, ci) 235 } 236 *cu = output 237 } 238 239 // IntersectsCellID reports whether this CellUnion intersects the given cell ID. 240 func (cu *CellUnion) IntersectsCellID(id CellID) bool { 241 // Find index of array item that occurs directly after our probe cell: 242 i := sort.Search(len(*cu), func(i int) bool { return id < (*cu)[i] }) 243 244 if i != len(*cu) && (*cu)[i].RangeMin() <= id.RangeMax() { 245 return true 246 } 247 return i != 0 && (*cu)[i-1].RangeMax() >= id.RangeMin() 248 } 249 250 // ContainsCellID reports whether the CellUnion contains the given cell ID. 251 // Containment is defined with respect to regions, e.g. a cell contains its 4 children. 252 // 253 // CAVEAT: If you have constructed a non-normalized CellUnion, note that groups 254 // of 4 child cells are *not* considered to contain their parent cell. To get 255 // this behavior you must use one of the call Normalize() explicitly. 256 func (cu *CellUnion) ContainsCellID(id CellID) bool { 257 // Find index of array item that occurs directly after our probe cell: 258 i := sort.Search(len(*cu), func(i int) bool { return id < (*cu)[i] }) 259 260 if i != len(*cu) && (*cu)[i].RangeMin() <= id { 261 return true 262 } 263 return i != 0 && (*cu)[i-1].RangeMax() >= id 264 } 265 266 // Denormalize replaces this CellUnion with an expanded version of the 267 // CellUnion where any cell whose level is less than minLevel or where 268 // (level - minLevel) is not a multiple of levelMod is replaced by its 269 // children, until either both of these conditions are satisfied or the 270 // maximum level is reached. 271 func (cu *CellUnion) Denormalize(minLevel, levelMod int) { 272 var denorm CellUnion 273 for _, id := range *cu { 274 level := id.Level() 275 newLevel := level 276 if newLevel < minLevel { 277 newLevel = minLevel 278 } 279 if levelMod > 1 { 280 newLevel += (MaxLevel - (newLevel - minLevel)) % levelMod 281 if newLevel > MaxLevel { 282 newLevel = MaxLevel 283 } 284 } 285 if newLevel == level { 286 denorm = append(denorm, id) 287 } else { 288 end := id.ChildEndAtLevel(newLevel) 289 for ci := id.ChildBeginAtLevel(newLevel); ci != end; ci = ci.Next() { 290 denorm = append(denorm, ci) 291 } 292 } 293 } 294 *cu = denorm 295 } 296 297 // RectBound returns a Rect that bounds this entity. 298 func (cu *CellUnion) RectBound() Rect { 299 bound := EmptyRect() 300 for _, c := range *cu { 301 bound = bound.Union(CellFromCellID(c).RectBound()) 302 } 303 return bound 304 } 305 306 // CapBound returns a Cap that bounds this entity. 307 func (cu *CellUnion) CapBound() Cap { 308 if len(*cu) == 0 { 309 return EmptyCap() 310 } 311 312 // Compute the approximate centroid of the region. This won't produce the 313 // bounding cap of minimal area, but it should be close enough. 314 var centroid Point 315 316 for _, ci := range *cu { 317 area := AvgAreaMetric.Value(ci.Level()) 318 centroid = Point{centroid.Add(ci.Point().Mul(area))} 319 } 320 321 if zero := (Point{}); centroid == zero { 322 centroid = PointFromCoords(1, 0, 0) 323 } else { 324 centroid = Point{centroid.Normalize()} 325 } 326 327 // Use the centroid as the cap axis, and expand the cap angle so that it 328 // contains the bounding caps of all the individual cells. Note that it is 329 // *not* sufficient to just bound all the cell vertices because the bounding 330 // cap may be concave (i.e. cover more than one hemisphere). 331 c := CapFromPoint(centroid) 332 for _, ci := range *cu { 333 c = c.AddCap(CellFromCellID(ci).CapBound()) 334 } 335 336 return c 337 } 338 339 // ContainsCell reports whether this cell union contains the given cell. 340 func (cu *CellUnion) ContainsCell(c Cell) bool { 341 return cu.ContainsCellID(c.id) 342 } 343 344 // IntersectsCell reports whether this cell union intersects the given cell. 345 func (cu *CellUnion) IntersectsCell(c Cell) bool { 346 return cu.IntersectsCellID(c.id) 347 } 348 349 // ContainsPoint reports whether this cell union contains the given point. 350 func (cu *CellUnion) ContainsPoint(p Point) bool { 351 return cu.ContainsCell(CellFromPoint(p)) 352 } 353 354 // CellUnionBound computes a covering of the CellUnion. 355 func (cu *CellUnion) CellUnionBound() []CellID { 356 return cu.CapBound().CellUnionBound() 357 } 358 359 // LeafCellsCovered reports the number of leaf cells covered by this cell union. 360 // This will be no more than 6*2^60 for the whole sphere. 361 func (cu *CellUnion) LeafCellsCovered() int64 { 362 var numLeaves int64 363 for _, c := range *cu { 364 numLeaves += 1 << uint64((MaxLevel-int64(c.Level()))<<1) 365 } 366 return numLeaves 367 } 368 369 // Returns true if the given four cells have a common parent. 370 // This requires that the four CellIDs are distinct. 371 func areSiblings(a, b, c, d CellID) bool { 372 // A necessary (but not sufficient) condition is that the XOR of the 373 // four cell IDs must be zero. This is also very fast to test. 374 if (a ^ b ^ c) != d { 375 return false 376 } 377 378 // Now we do a slightly more expensive but exact test. First, compute a 379 // mask that blocks out the two bits that encode the child position of 380 // "id" with respect to its parent, then check that the other three 381 // children all agree with "mask". 382 mask := d.lsb() << 1 383 mask = ^(mask + (mask << 1)) 384 idMasked := (uint64(d) & mask) 385 return ((uint64(a)&mask) == idMasked && 386 (uint64(b)&mask) == idMasked && 387 (uint64(c)&mask) == idMasked && 388 !d.isFace()) 389 } 390 391 // Contains reports whether this CellUnion contains all of the CellIDs of the given CellUnion. 392 func (cu *CellUnion) Contains(o CellUnion) bool { 393 // TODO(roberts): Investigate alternatives such as divide-and-conquer 394 // or alternating-skip-search that may be significantly faster in both 395 // the average and worst case. This applies to Intersects as well. 396 for _, id := range o { 397 if !cu.ContainsCellID(id) { 398 return false 399 } 400 } 401 402 return true 403 } 404 405 // Intersects reports whether this CellUnion intersects any of the CellIDs of the given CellUnion. 406 func (cu *CellUnion) Intersects(o CellUnion) bool { 407 for _, c := range *cu { 408 if o.IntersectsCellID(c) { 409 return true 410 } 411 } 412 413 return false 414 } 415 416 // lowerBound returns the index in this CellUnion to the first element whose value 417 // is not considered to go before the given cell id. (i.e., either it is equivalent 418 // or comes after the given id.) If there is no match, then end is returned. 419 func (cu *CellUnion) lowerBound(begin, end int, id CellID) int { 420 for i := begin; i < end; i++ { 421 if (*cu)[i] >= id { 422 return i 423 } 424 } 425 426 return end 427 } 428 429 // cellUnionDifferenceInternal adds the difference between the CellID and the union to 430 // the result CellUnion. If they intersect but the difference is non-empty, it divides 431 // and conquers. 432 func (cu *CellUnion) cellUnionDifferenceInternal(id CellID, other *CellUnion) { 433 if !other.IntersectsCellID(id) { 434 (*cu) = append((*cu), id) 435 return 436 } 437 438 if !other.ContainsCellID(id) { 439 for _, child := range id.Children() { 440 cu.cellUnionDifferenceInternal(child, other) 441 } 442 } 443 } 444 445 // ExpandAtLevel expands this CellUnion by adding a rim of cells at expandLevel 446 // around the unions boundary. 447 // 448 // For each cell c in the union, we add all cells at level 449 // expandLevel that abut c. There are typically eight of those 450 // (four edge-abutting and four sharing a vertex). However, if c is 451 // finer than expandLevel, we add all cells abutting 452 // c.Parent(expandLevel) as well as c.Parent(expandLevel) itself, 453 // as an expandLevel cell rarely abuts a smaller cell. 454 // 455 // Note that the size of the output is exponential in 456 // expandLevel. For example, if expandLevel == 20 and the input 457 // has a cell at level 10, there will be on the order of 4000 458 // adjacent cells in the output. For most applications the 459 // ExpandByRadius method below is easier to use. 460 func (cu *CellUnion) ExpandAtLevel(level int) { 461 var output CellUnion 462 levelLsb := lsbForLevel(level) 463 for i := len(*cu) - 1; i >= 0; i-- { 464 id := (*cu)[i] 465 if id.lsb() < levelLsb { 466 id = id.Parent(level) 467 // Optimization: skip over any cells contained by this one. This is 468 // especially important when very small regions are being expanded. 469 for i > 0 && id.Contains((*cu)[i-1]) { 470 i-- 471 } 472 } 473 output = append(output, id) 474 output = append(output, id.AllNeighbors(level)...) 475 } 476 sortCellIDs(output) 477 478 *cu = output 479 cu.Normalize() 480 } 481 482 // ExpandByRadius expands this CellUnion such that it contains all points whose 483 // distance to the CellUnion is at most minRadius, but do not use cells that 484 // are more than maxLevelDiff levels higher than the largest cell in the input. 485 // The second parameter controls the tradeoff between accuracy and output size 486 // when a large region is being expanded by a small amount (e.g. expanding Canada 487 // by 1km). For example, if maxLevelDiff == 4 the region will always be expanded 488 // by approximately 1/16 the width of its largest cell. Note that in the worst case, 489 // the number of cells in the output can be up to 4 * (1 + 2 ** maxLevelDiff) times 490 // larger than the number of cells in the input. 491 func (cu *CellUnion) ExpandByRadius(minRadius s1.Angle, maxLevelDiff int) { 492 minLevel := MaxLevel 493 for _, cid := range *cu { 494 minLevel = minInt(minLevel, cid.Level()) 495 } 496 497 // Find the maximum level such that all cells are at least "minRadius" wide. 498 radiusLevel := MinWidthMetric.MaxLevel(minRadius.Radians()) 499 if radiusLevel == 0 && minRadius.Radians() > MinWidthMetric.Value(0) { 500 // The requested expansion is greater than the width of a face cell. 501 // The easiest way to handle this is to expand twice. 502 cu.ExpandAtLevel(0) 503 } 504 cu.ExpandAtLevel(minInt(minLevel+maxLevelDiff, radiusLevel)) 505 } 506 507 // Equal reports whether the two CellUnions are equal. 508 func (cu CellUnion) Equal(o CellUnion) bool { 509 if len(cu) != len(o) { 510 return false 511 } 512 for i := 0; i < len(cu); i++ { 513 if cu[i] != o[i] { 514 return false 515 } 516 } 517 return true 518 } 519 520 // AverageArea returns the average area of this CellUnion. 521 // This is accurate to within a factor of 1.7. 522 func (cu *CellUnion) AverageArea() float64 { 523 return AvgAreaMetric.Value(MaxLevel) * float64(cu.LeafCellsCovered()) 524 } 525 526 // ApproxArea returns the approximate area of this CellUnion. This method is accurate 527 // to within 3% percent for all cell sizes and accurate to within 0.1% for cells 528 // at level 5 or higher within the union. 529 func (cu *CellUnion) ApproxArea() float64 { 530 var area float64 531 for _, id := range *cu { 532 area += CellFromCellID(id).ApproxArea() 533 } 534 return area 535 } 536 537 // ExactArea returns the area of this CellUnion as accurately as possible. 538 func (cu *CellUnion) ExactArea() float64 { 539 var area float64 540 for _, id := range *cu { 541 area += CellFromCellID(id).ExactArea() 542 } 543 return area 544 } 545 546 // Encode encodes the CellUnion. 547 func (cu *CellUnion) Encode(w io.Writer) error { 548 e := &encoder{w: w} 549 cu.encode(e) 550 return e.err 551 } 552 553 func (cu *CellUnion) encode(e *encoder) { 554 e.writeInt8(encodingVersion) 555 e.writeInt64(int64(len(*cu))) 556 for _, ci := range *cu { 557 ci.encode(e) 558 } 559 } 560 561 // Decode decodes the CellUnion. 562 func (cu *CellUnion) Decode(r io.Reader) error { 563 d := &decoder{r: asByteReader(r)} 564 cu.decode(d) 565 return d.err 566 } 567 568 func (cu *CellUnion) decode(d *decoder) { 569 version := d.readInt8() 570 if d.err != nil { 571 return 572 } 573 if version != encodingVersion { 574 d.err = fmt.Errorf("only version %d is supported", encodingVersion) 575 return 576 } 577 n := d.readInt64() 578 if d.err != nil { 579 return 580 } 581 const maxCells = 1000000 582 if n > maxCells { 583 d.err = fmt.Errorf("too many cells (%d; max is %d)", n, maxCells) 584 return 585 } 586 *cu = make([]CellID, n) 587 for i := range *cu { 588 (*cu)[i].decode(d) 589 } 590 } 591