1 // Copyright 2015 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 "container/heap" 19 "sort" 20 ) 21 22 // RegionCoverer allows arbitrary regions to be approximated as unions of cells (CellUnion). 23 // This is useful for implementing various sorts of search and precomputation operations. 24 // 25 // Typical usage: 26 // 27 // rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 5} 28 // r := s2.Region(CapFromCenterArea(center, area)) 29 // covering := rc.Covering(r) 30 // 31 // This yields a CellUnion of at most 5 cells that is guaranteed to cover the 32 // given region (a disc-shaped region on the sphere). 33 // 34 // For covering, only cells where (level - MinLevel) is a multiple of LevelMod will be used. 35 // This effectively allows the branching factor of the S2 CellID hierarchy to be increased. 36 // Currently the only parameter values allowed are 1, 2, or 3, corresponding to 37 // branching factors of 4, 16, and 64 respectively. 38 // 39 // Note the following: 40 // 41 // - MinLevel takes priority over MaxCells, i.e. cells below the given level will 42 // never be used even if this causes a large number of cells to be returned. 43 // 44 // - For any setting of MaxCells, up to 6 cells may be returned if that 45 // is the minimum number of cells required (e.g. if the region intersects 46 // all six face cells). Up to 3 cells may be returned even for very tiny 47 // convex regions if they happen to be located at the intersection of 48 // three cube faces. 49 // 50 // - For any setting of MaxCells, an arbitrary number of cells may be 51 // returned if MinLevel is too high for the region being approximated. 52 // 53 // - If MaxCells is less than 4, the area of the covering may be 54 // arbitrarily large compared to the area of the original region even if 55 // the region is convex (e.g. a Cap or Rect). 56 // 57 // The approximation algorithm is not optimal but does a pretty good job in 58 // practice. The output does not always use the maximum number of cells 59 // allowed, both because this would not always yield a better approximation, 60 // and because MaxCells is a limit on how much work is done exploring the 61 // possible covering as well as a limit on the final output size. 62 // 63 // Because it is an approximation algorithm, one should not rely on the 64 // stability of the output. In particular, the output of the covering algorithm 65 // may change across different versions of the library. 66 // 67 // One can also generate interior coverings, which are sets of cells which 68 // are entirely contained within a region. Interior coverings can be 69 // empty, even for non-empty regions, if there are no cells that satisfy 70 // the provided constraints and are contained by the region. Note that for 71 // performance reasons, it is wise to specify a MaxLevel when computing 72 // interior coverings - otherwise for regions with small or zero area, the 73 // algorithm may spend a lot of time subdividing cells all the way to leaf 74 // level to try to find contained cells. 75 type RegionCoverer struct { 76 MinLevel int // the minimum cell level to be used. 77 MaxLevel int // the maximum cell level to be used. 78 LevelMod int // the LevelMod to be used. 79 MaxCells int // the maximum desired number of cells in the approximation. 80 } 81 82 // NewRegionCoverer returns a region coverer with the appropriate defaults. 83 func NewRegionCoverer() *RegionCoverer { 84 return &RegionCoverer{ 85 MinLevel: 0, 86 MaxLevel: MaxLevel, 87 LevelMod: 1, 88 MaxCells: 8, 89 } 90 } 91 92 type coverer struct { 93 minLevel int // the minimum cell level to be used. 94 MaxLevel int // the maximum cell level to be used. 95 levelMod int // the LevelMod to be used. 96 maxCells int // the maximum desired number of cells in the approximation. 97 region Region 98 result CellUnion 99 pq priorityQueue 100 interiorCovering bool 101 } 102 103 type candidate struct { 104 cell Cell 105 terminal bool // Cell should not be expanded further. 106 numChildren int // Number of children that intersect the region. 107 children []*candidate // Actual size may be 0, 4, 16, or 64 elements. 108 priority int // Priority of the candidate. 109 } 110 111 type priorityQueue []*candidate 112 113 func (pq priorityQueue) Len() int { 114 return len(pq) 115 } 116 117 func (pq priorityQueue) Less(i, j int) bool { 118 // We want Pop to give us the highest, not lowest, priority so we use greater than here. 119 return pq[i].priority > pq[j].priority 120 } 121 122 func (pq priorityQueue) Swap(i, j int) { 123 pq[i], pq[j] = pq[j], pq[i] 124 } 125 126 func (pq *priorityQueue) Push(x any) { 127 item := x.(*candidate) 128 *pq = append(*pq, item) 129 } 130 131 func (pq *priorityQueue) Pop() any { 132 item := (*pq)[len(*pq)-1] 133 *pq = (*pq)[:len(*pq)-1] 134 return item 135 } 136 137 func (pq *priorityQueue) Reset() { 138 *pq = (*pq)[:0] 139 } 140 141 // newCandidate returns a new candidate with no children if the cell intersects the given region. 142 // The candidate is marked as terminal if it should not be expanded further. 143 func (c *coverer) newCandidate(cell Cell) *candidate { 144 if !c.region.IntersectsCell(cell) { 145 return nil 146 } 147 cand := &candidate{cell: cell} 148 level := int(cell.level) 149 if level >= c.minLevel { 150 if c.interiorCovering { 151 if c.region.ContainsCell(cell) { 152 cand.terminal = true 153 } else if level+c.levelMod > c.MaxLevel { 154 return nil 155 } 156 } else if level+c.levelMod > c.MaxLevel || c.region.ContainsCell(cell) { 157 cand.terminal = true 158 } 159 } 160 return cand 161 } 162 163 // expandChildren populates the children of the candidate by expanding the given number of 164 // levels from the given cell. Returns the number of children that were marked "terminal". 165 func (c *coverer) expandChildren(cand *candidate, cell Cell, numLevels int) int { 166 numLevels-- 167 var numTerminals int 168 last := cell.id.ChildEnd() 169 for ci := cell.id.ChildBegin(); ci != last; ci = ci.Next() { 170 childCell := CellFromCellID(ci) 171 if numLevels > 0 { 172 if c.region.IntersectsCell(childCell) { 173 numTerminals += c.expandChildren(cand, childCell, numLevels) 174 } 175 continue 176 } 177 if child := c.newCandidate(childCell); child != nil { 178 cand.children = append(cand.children, child) 179 cand.numChildren++ 180 if child.terminal { 181 numTerminals++ 182 } 183 } 184 } 185 return numTerminals 186 } 187 188 // addCandidate adds the given candidate to the result if it is marked as "terminal", 189 // otherwise expands its children and inserts it into the priority queue. 190 // Passing an argument of nil does nothing. 191 func (c *coverer) addCandidate(cand *candidate) { 192 if cand == nil { 193 return 194 } 195 196 if cand.terminal { 197 c.result = append(c.result, cand.cell.id) 198 return 199 } 200 201 // Expand one level at a time until we hit minLevel to ensure that we don't skip over it. 202 numLevels := c.levelMod 203 level := int(cand.cell.level) 204 if level < c.minLevel { 205 numLevels = 1 206 } 207 208 numTerminals := c.expandChildren(cand, cand.cell, numLevels) 209 maxChildrenShift := uint(2 * c.levelMod) 210 if cand.numChildren == 0 { 211 return 212 } else if !c.interiorCovering && numTerminals == 1<<maxChildrenShift && level >= c.minLevel { 213 // Optimization: add the parent cell rather than all of its children. 214 // We can't do this for interior coverings, since the children just 215 // intersect the region, but may not be contained by it - we need to 216 // subdivide them further. 217 cand.terminal = true 218 c.addCandidate(cand) 219 } else { 220 // We negate the priority so that smaller absolute priorities are returned 221 // first. The heuristic is designed to refine the largest cells first, 222 // since those are where we have the largest potential gain. Among cells 223 // of the same size, we prefer the cells with the fewest children. 224 // Finally, among cells with equal numbers of children we prefer those 225 // with the smallest number of children that cannot be refined further. 226 cand.priority = -(((level<<maxChildrenShift)+cand.numChildren)<<maxChildrenShift + numTerminals) 227 heap.Push(&c.pq, cand) 228 } 229 } 230 231 // adjustLevel returns the reduced "level" so that it satisfies levelMod. Levels smaller than minLevel 232 // are not affected (since cells at these levels are eventually expanded). 233 func (c *coverer) adjustLevel(level int) int { 234 if c.levelMod > 1 && level > c.minLevel { 235 level -= (level - c.minLevel) % c.levelMod 236 } 237 return level 238 } 239 240 // adjustCellLevels ensures that all cells with level > minLevel also satisfy levelMod, 241 // by replacing them with an ancestor if necessary. Cell levels smaller 242 // than minLevel are not modified (see AdjustLevel). The output is 243 // then normalized to ensure that no redundant cells are present. 244 func (c *coverer) adjustCellLevels(cells *CellUnion) { 245 if c.levelMod == 1 { 246 return 247 } 248 249 var out int 250 for _, ci := range *cells { 251 level := ci.Level() 252 newLevel := c.adjustLevel(level) 253 if newLevel != level { 254 ci = ci.Parent(newLevel) 255 } 256 if out > 0 && (*cells)[out-1].Contains(ci) { 257 continue 258 } 259 for out > 0 && ci.Contains((*cells)[out-1]) { 260 out-- 261 } 262 (*cells)[out] = ci 263 out++ 264 } 265 *cells = (*cells)[:out] 266 } 267 268 // initialCandidates computes a set of initial candidates that cover the given region. 269 func (c *coverer) initialCandidates() { 270 // Optimization: start with a small (usually 4 cell) covering of the region's bounding cap. 271 temp := &RegionCoverer{MaxLevel: c.MaxLevel, LevelMod: 1, MaxCells: minInt(4, c.maxCells)} 272 273 cells := temp.FastCovering(c.region) 274 c.adjustCellLevels(&cells) 275 for _, ci := range cells { 276 c.addCandidate(c.newCandidate(CellFromCellID(ci))) 277 } 278 } 279 280 // coveringInternal generates a covering and stores it in result. 281 // Strategy: Start with the 6 faces of the cube. Discard any 282 // that do not intersect the shape. Then repeatedly choose the 283 // largest cell that intersects the shape and subdivide it. 284 // 285 // result contains the cells that will be part of the output, while pq 286 // contains cells that we may still subdivide further. Cells that are 287 // entirely contained within the region are immediately added to the output, 288 // while cells that do not intersect the region are immediately discarded. 289 // Therefore pq only contains cells that partially intersect the region. 290 // Candidates are prioritized first according to cell size (larger cells 291 // first), then by the number of intersecting children they have (fewest 292 // children first), and then by the number of fully contained children 293 // (fewest children first). 294 func (c *coverer) coveringInternal(region Region) { 295 c.region = region 296 297 c.initialCandidates() 298 for c.pq.Len() > 0 && (!c.interiorCovering || len(c.result) < c.maxCells) { 299 cand := heap.Pop(&c.pq).(*candidate) 300 301 // For interior covering we keep subdividing no matter how many children 302 // candidate has. If we reach MaxCells before expanding all children, 303 // we will just use some of them. 304 // For exterior covering we cannot do this, because result has to cover the 305 // whole region, so all children have to be used. 306 // candidate.numChildren == 1 case takes care of the situation when we 307 // already have more than MaxCells in result (minLevel is too high). 308 // Subdividing of the candidate with one child does no harm in this case. 309 if c.interiorCovering || int(cand.cell.level) < c.minLevel || cand.numChildren == 1 || len(c.result)+c.pq.Len()+cand.numChildren <= c.maxCells { 310 for _, child := range cand.children { 311 if !c.interiorCovering || len(c.result) < c.maxCells { 312 c.addCandidate(child) 313 } 314 } 315 } else { 316 cand.terminal = true 317 c.addCandidate(cand) 318 } 319 } 320 321 c.pq.Reset() 322 c.region = nil 323 324 // Rather than just returning the raw list of cell ids, we construct a cell 325 // union and then denormalize it. This has the effect of replacing four 326 // child cells with their parent whenever this does not violate the covering 327 // parameters specified (min_level, level_mod, etc). This significantly 328 // reduces the number of cells returned in many cases, and it is cheap 329 // compared to computing the covering in the first place. 330 c.result.Normalize() 331 if c.minLevel > 0 || c.levelMod > 1 { 332 c.result.Denormalize(c.minLevel, c.levelMod) 333 } 334 } 335 336 // newCoverer returns an instance of coverer. 337 func (rc *RegionCoverer) newCoverer() *coverer { 338 return &coverer{ 339 minLevel: maxInt(0, minInt(MaxLevel, rc.MinLevel)), 340 MaxLevel: maxInt(0, minInt(MaxLevel, rc.MaxLevel)), 341 levelMod: maxInt(1, minInt(3, rc.LevelMod)), 342 maxCells: rc.MaxCells, 343 } 344 } 345 346 // Covering returns a CellUnion that covers the given region and satisfies the various restrictions. 347 func (rc *RegionCoverer) Covering(region Region) CellUnion { 348 covering := rc.CellUnion(region) 349 covering.Denormalize(maxInt(0, minInt(MaxLevel, rc.MinLevel)), maxInt(1, minInt(3, rc.LevelMod))) 350 return covering 351 } 352 353 // InteriorCovering returns a CellUnion that is contained within the given region and satisfies the various restrictions. 354 func (rc *RegionCoverer) InteriorCovering(region Region) CellUnion { 355 intCovering := rc.InteriorCellUnion(region) 356 intCovering.Denormalize(maxInt(0, minInt(MaxLevel, rc.MinLevel)), maxInt(1, minInt(3, rc.LevelMod))) 357 return intCovering 358 } 359 360 // CellUnion returns a normalized CellUnion that covers the given region and 361 // satisfies the restrictions except for minLevel and levelMod. These criteria 362 // cannot be satisfied using a cell union because cell unions are 363 // automatically normalized by replacing four child cells with their parent 364 // whenever possible. (Note that the list of cell ids passed to the CellUnion 365 // constructor does in fact satisfy all the given restrictions.) 366 func (rc *RegionCoverer) CellUnion(region Region) CellUnion { 367 c := rc.newCoverer() 368 c.coveringInternal(region) 369 cu := c.result 370 cu.Normalize() 371 return cu 372 } 373 374 // InteriorCellUnion returns a normalized CellUnion that is contained within the given region and 375 // satisfies the restrictions except for minLevel and levelMod. These criteria 376 // cannot be satisfied using a cell union because cell unions are 377 // automatically normalized by replacing four child cells with their parent 378 // whenever possible. (Note that the list of cell ids passed to the CellUnion 379 // constructor does in fact satisfy all the given restrictions.) 380 func (rc *RegionCoverer) InteriorCellUnion(region Region) CellUnion { 381 c := rc.newCoverer() 382 c.interiorCovering = true 383 c.coveringInternal(region) 384 cu := c.result 385 cu.Normalize() 386 return cu 387 } 388 389 // FastCovering returns a CellUnion that covers the given region similar to Covering, 390 // except that this method is much faster and the coverings are not as tight. 391 // All of the usual parameters are respected (MaxCells, MinLevel, MaxLevel, and LevelMod), 392 // except that the implementation makes no attempt to take advantage of large values of 393 // MaxCells. (A small number of cells will always be returned.) 394 // 395 // This function is useful as a starting point for algorithms that 396 // recursively subdivide cells. 397 func (rc *RegionCoverer) FastCovering(region Region) CellUnion { 398 c := rc.newCoverer() 399 cu := CellUnion(region.CellUnionBound()) 400 c.normalizeCovering(&cu) 401 return cu 402 } 403 404 // IsCanonical reports whether the given CellUnion represents a valid covering 405 // that conforms to the current covering parameters. In particular: 406 // 407 // - All CellIDs must be valid. 408 // 409 // - CellIDs must be sorted and non-overlapping. 410 // 411 // - CellID levels must satisfy MinLevel, MaxLevel, and LevelMod. 412 // 413 // - If the covering has more than MaxCells, there must be no two cells with 414 // a common ancestor at MinLevel or higher. 415 // 416 // - There must be no sequence of cells that could be replaced by an 417 // ancestor (i.e. with LevelMod == 1, the 4 child cells of a parent). 418 func (rc *RegionCoverer) IsCanonical(covering CellUnion) bool { 419 return rc.newCoverer().isCanonical(covering) 420 } 421 422 // normalizeCovering normalizes the "covering" so that it conforms to the 423 // current covering parameters (maxCells, minLevel, MaxLevel, and levelMod). 424 // This method makes no attempt to be optimal. In particular, if 425 // minLevel > 0 or levelMod > 1 then it may return more than the 426 // desired number of cells even when this isn't necessary. 427 // 428 // Note that when the covering parameters have their default values, almost 429 // all of the code in this function is skipped. 430 func (c *coverer) normalizeCovering(covering *CellUnion) { 431 // If any cells are too small, or don't satisfy levelMod, then replace them with ancestors. 432 if c.MaxLevel < MaxLevel || c.levelMod > 1 { 433 for i, ci := range *covering { 434 level := ci.Level() 435 newLevel := c.adjustLevel(minInt(level, c.MaxLevel)) 436 if newLevel != level { 437 (*covering)[i] = ci.Parent(newLevel) 438 } 439 } 440 } 441 // Sort the cells and simplify them. 442 covering.Normalize() 443 444 // Make sure that the covering satisfies minLevel and levelMod, 445 // possibly at the expense of satisfying MaxCells. 446 if c.minLevel > 0 || c.levelMod > 1 { 447 covering.Denormalize(c.minLevel, c.levelMod) 448 } 449 450 // If there are too many cells and the covering is very large, use the 451 // RegionCoverer to compute a new covering. (This avoids possible O(n^2) 452 // behavior of the simpler algorithm below.) 453 excess := len(*covering) - c.maxCells 454 if excess <= 0 || c.isCanonical(*covering) { 455 return 456 } 457 if excess*len(*covering) > 10000 { 458 rc := NewRegionCoverer() 459 (*covering) = rc.Covering(covering) 460 return 461 } 462 463 // If there are still too many cells, then repeatedly replace two adjacent 464 // cells in CellID order by their lowest common ancestor. 465 for len(*covering) > c.maxCells { 466 bestIndex := -1 467 bestLevel := -1 468 for i := 0; i+1 < len(*covering); i++ { 469 level, ok := (*covering)[i].CommonAncestorLevel((*covering)[i+1]) 470 if !ok { 471 continue 472 } 473 level = c.adjustLevel(level) 474 if level > bestLevel { 475 bestLevel = level 476 bestIndex = i 477 } 478 } 479 480 if bestLevel < c.minLevel { 481 break 482 } 483 484 // Replace all cells contained by the new ancestor cell. 485 id := (*covering)[bestIndex].Parent(bestLevel) 486 (*covering) = c.replaceCellsWithAncestor(*covering, id) 487 488 // Now repeatedly check whether all children of the parent cell are 489 // present, in which case we can replace those cells with their parent. 490 for bestLevel > c.minLevel { 491 bestLevel -= c.levelMod 492 id = id.Parent(bestLevel) 493 if !c.containsAllChildren(*covering, id) { 494 break 495 } 496 (*covering) = c.replaceCellsWithAncestor(*covering, id) 497 } 498 } 499 } 500 501 // isCanonical reports whether the covering is canonical. 502 func (c *coverer) isCanonical(covering CellUnion) bool { 503 trueMax := c.MaxLevel 504 if c.levelMod != 1 { 505 trueMax = c.MaxLevel - (c.MaxLevel-c.minLevel)%c.levelMod 506 } 507 tooManyCells := len(covering) > c.maxCells 508 sameParentCount := 1 509 510 prevID := CellID(0) 511 for _, id := range covering { 512 if !id.IsValid() { 513 return false 514 } 515 516 // Check that the CellID level is acceptable. 517 level := id.Level() 518 if level < c.minLevel || level > trueMax { 519 return false 520 } 521 if c.levelMod > 1 && (level-c.minLevel)%c.levelMod != 0 { 522 return false 523 } 524 525 if prevID != 0 { 526 // Check that cells are sorted and non-overlapping. 527 if prevID.RangeMax() >= id.RangeMin() { 528 return false 529 } 530 531 lev, ok := id.CommonAncestorLevel(prevID) 532 // If there are too many cells, check that no pair of adjacent cells 533 // could be replaced by an ancestor. 534 if tooManyCells && (ok && lev >= c.minLevel) { 535 return false 536 } 537 538 // Check that there are no sequences of (4 ** level_mod) cells that all 539 // have the same parent (considering only multiples of "level_mod"). 540 pLevel := level - c.levelMod 541 if pLevel < c.minLevel || level != prevID.Level() || 542 id.Parent(pLevel) != prevID.Parent(pLevel) { 543 sameParentCount = 1 544 } else { 545 sameParentCount++ 546 if sameParentCount == 1<<uint(2*c.levelMod) { 547 return false 548 } 549 } 550 } 551 prevID = id 552 } 553 554 return true 555 } 556 557 func (c *coverer) containsAllChildren(covering []CellID, id CellID) bool { 558 pos := sort.Search(len(covering), func(i int) bool { return (covering)[i] >= id.RangeMin() }) 559 level := id.Level() + c.levelMod 560 for child := id.ChildBeginAtLevel(level); child != id.ChildEndAtLevel(level); child = child.Next() { 561 if pos == len(covering) || covering[pos] != child { 562 return false 563 } 564 pos++ 565 } 566 return true 567 } 568 569 // replaceCellsWithAncestor replaces all descendants of the given id in covering 570 // with id. This requires the covering contains at least one descendant of id. 571 func (c *coverer) replaceCellsWithAncestor(covering []CellID, id CellID) []CellID { 572 begin := sort.Search(len(covering), func(i int) bool { return covering[i] > id.RangeMin() }) 573 end := sort.Search(len(covering), func(i int) bool { return covering[i] > id.RangeMax() }) 574 575 return append(append(covering[:begin], id), covering[end:]...) 576 } 577 578 // SimpleRegionCovering returns a set of cells at the given level that cover 579 // the connected region and a starting point on the boundary or inside the 580 // region. The cells are returned in arbitrary order. 581 // 582 // Note that this method is not faster than the regular Covering 583 // method for most region types, such as Cap or Polygon, and in fact it 584 // can be much slower when the output consists of a large number of cells. 585 // Currently it can be faster at generating coverings of long narrow regions 586 // such as polylines, but this may change in the future. 587 func SimpleRegionCovering(region Region, start Point, level int) []CellID { 588 return FloodFillRegionCovering(region, cellIDFromPoint(start).Parent(level)) 589 } 590 591 // FloodFillRegionCovering returns all edge-connected cells at the same level as 592 // the given CellID that intersect the given region, in arbitrary order. 593 func FloodFillRegionCovering(region Region, start CellID) []CellID { 594 var output []CellID 595 all := map[CellID]bool{ 596 start: true, 597 } 598 frontier := []CellID{start} 599 for len(frontier) > 0 { 600 id := frontier[len(frontier)-1] 601 frontier = frontier[:len(frontier)-1] 602 if !region.IntersectsCell(CellFromCellID(id)) { 603 continue 604 } 605 output = append(output, id) 606 for _, nbr := range id.EdgeNeighbors() { 607 if !all[nbr] { 608 all[nbr] = true 609 frontier = append(frontier, nbr) 610 } 611 } 612 } 613 614 return output 615 } 616