// Copyright 2015 Google Inc. All rights reserved. // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. // You may obtain a copy of the License at // // http://www.apache.org/licenses/LICENSE-2.0 // // Unless required by applicable law or agreed to in writing, software // distributed under the License is distributed on an "AS IS" BASIS, // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. package s2 import ( "container/heap" ) // RegionCoverer allows arbitrary regions to be approximated as unions of cells (CellUnion). // This is useful for implementing various sorts of search and precomputation operations. // // Typical usage: // // rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 5} // r := s2.Region(CapFromCenterArea(center, area)) // covering := rc.Covering(r) // // This yields a CellUnion of at most 5 cells that is guaranteed to cover the // given region (a disc-shaped region on the sphere). // // For covering, only cells where (level - MinLevel) is a multiple of LevelMod will be used. // This effectively allows the branching factor of the S2 CellID hierarchy to be increased. // Currently the only parameter values allowed are 1, 2, or 3, corresponding to // branching factors of 4, 16, and 64 respectively. // // Note the following: // // - MinLevel takes priority over MaxCells, i.e. cells below the given level will // never be used even if this causes a large number of cells to be returned. // // - For any setting of MaxCells, up to 6 cells may be returned if that // is the minimum number of cells required (e.g. if the region intersects // all six face cells). Up to 3 cells may be returned even for very tiny // convex regions if they happen to be located at the intersection of // three cube faces. // // - For any setting of MaxCells, an arbitrary number of cells may be // returned if MinLevel is too high for the region being approximated. // // - If MaxCells is less than 4, the area of the covering may be // arbitrarily large compared to the area of the original region even if // the region is convex (e.g. a Cap or Rect). // // The approximation algorithm is not optimal but does a pretty good job in // practice. The output does not always use the maximum number of cells // allowed, both because this would not always yield a better approximation, // and because MaxCells is a limit on how much work is done exploring the // possible covering as well as a limit on the final output size. // // Because it is an approximation algorithm, one should not rely on the // stability of the output. In particular, the output of the covering algorithm // may change across different versions of the library. // // One can also generate interior coverings, which are sets of cells which // are entirely contained within a region. Interior coverings can be // empty, even for non-empty regions, if there are no cells that satisfy // the provided constraints and are contained by the region. Note that for // performance reasons, it is wise to specify a MaxLevel when computing // interior coverings - otherwise for regions with small or zero area, the // algorithm may spend a lot of time subdividing cells all the way to leaf // level to try to find contained cells. type RegionCoverer struct { MinLevel int // the minimum cell level to be used. MaxLevel int // the maximum cell level to be used. LevelMod int // the LevelMod to be used. MaxCells int // the maximum desired number of cells in the approximation. } type coverer struct { minLevel int // the minimum cell level to be used. maxLevel int // the maximum cell level to be used. levelMod int // the LevelMod to be used. maxCells int // the maximum desired number of cells in the approximation. region Region result CellUnion pq priorityQueue interiorCovering bool } type candidate struct { cell Cell terminal bool // Cell should not be expanded further. numChildren int // Number of children that intersect the region. children []*candidate // Actual size may be 0, 4, 16, or 64 elements. priority int // Priority of the candidate. } type priorityQueue []*candidate func (pq priorityQueue) Len() int { return len(pq) } func (pq priorityQueue) Less(i, j int) bool { // We want Pop to give us the highest, not lowest, priority so we use greater than here. return pq[i].priority > pq[j].priority } func (pq priorityQueue) Swap(i, j int) { pq[i], pq[j] = pq[j], pq[i] } func (pq *priorityQueue) Push(x interface{}) { item := x.(*candidate) *pq = append(*pq, item) } func (pq *priorityQueue) Pop() interface{} { item := (*pq)[len(*pq)-1] *pq = (*pq)[:len(*pq)-1] return item } func (pq *priorityQueue) Reset() { *pq = (*pq)[:0] } // newCandidate returns a new candidate with no children if the cell intersects the given region. // The candidate is marked as terminal if it should not be expanded further. func (c *coverer) newCandidate(cell Cell) *candidate { if !c.region.IntersectsCell(cell) { return nil } cand := &candidate{cell: cell} level := int(cell.level) if level >= c.minLevel { if c.interiorCovering { if c.region.ContainsCell(cell) { cand.terminal = true } else if level+c.levelMod > c.maxLevel { return nil } } else if level+c.levelMod > c.maxLevel || c.region.ContainsCell(cell) { cand.terminal = true } } return cand } // expandChildren populates the children of the candidate by expanding the given number of // levels from the given cell. Returns the number of children that were marked "terminal". func (c *coverer) expandChildren(cand *candidate, cell Cell, numLevels int) int { numLevels-- var numTerminals int last := cell.id.ChildEnd() for ci := cell.id.ChildBegin(); ci != last; ci = ci.Next() { childCell := CellFromCellID(ci) if numLevels > 0 { if c.region.IntersectsCell(childCell) { numTerminals += c.expandChildren(cand, childCell, numLevels) } continue } if child := c.newCandidate(childCell); child != nil { cand.children = append(cand.children, child) cand.numChildren++ if child.terminal { numTerminals++ } } } return numTerminals } // addCandidate adds the given candidate to the result if it is marked as "terminal", // otherwise expands its children and inserts it into the priority queue. // Passing an argument of nil does nothing. func (c *coverer) addCandidate(cand *candidate) { if cand == nil { return } if cand.terminal { c.result = append(c.result, cand.cell.id) return } // Expand one level at a time until we hit minLevel to ensure that we don't skip over it. numLevels := c.levelMod level := int(cand.cell.level) if level < c.minLevel { numLevels = 1 } numTerminals := c.expandChildren(cand, cand.cell, numLevels) maxChildrenShift := uint(2 * c.levelMod) if cand.numChildren == 0 { return } else if !c.interiorCovering && numTerminals == 1<= c.minLevel { // Optimization: add the parent cell rather than all of its children. // We can't do this for interior coverings, since the children just // intersect the region, but may not be contained by it - we need to // subdivide them further. cand.terminal = true c.addCandidate(cand) } else { // We negate the priority so that smaller absolute priorities are returned // first. The heuristic is designed to refine the largest cells first, // since those are where we have the largest potential gain. Among cells // of the same size, we prefer the cells with the fewest children. // Finally, among cells with equal numbers of children we prefer those // with the smallest number of children that cannot be refined further. cand.priority = -(((level< 1 && level > c.minLevel { level -= (level - c.minLevel) % c.levelMod } return level } // adjustCellLevels ensures that all cells with level > minLevel also satisfy levelMod, // by replacing them with an ancestor if necessary. Cell levels smaller // than minLevel are not modified (see AdjustLevel). The output is // then normalized to ensure that no redundant cells are present. func (c *coverer) adjustCellLevels(cells *CellUnion) { if c.levelMod == 1 { return } var out int for _, ci := range *cells { level := ci.Level() newLevel := c.adjustLevel(level) if newLevel != level { ci = ci.Parent(newLevel) } if out > 0 && (*cells)[out-1].Contains(ci) { continue } for out > 0 && ci.Contains((*cells)[out-1]) { out-- } (*cells)[out] = ci out++ } *cells = (*cells)[:out] } // initialCandidates computes a set of initial candidates that cover the given region. func (c *coverer) initialCandidates() { // Optimization: start with a small (usually 4 cell) covering of the region's bounding cap. temp := &RegionCoverer{MaxLevel: c.maxLevel, LevelMod: 1, MaxCells: minInt(4, c.maxCells)} cells := temp.FastCovering(c.region) c.adjustCellLevels(&cells) for _, ci := range cells { c.addCandidate(c.newCandidate(CellFromCellID(ci))) } } // coveringInternal generates a covering and stores it in result. // Strategy: Start with the 6 faces of the cube. Discard any // that do not intersect the shape. Then repeatedly choose the // largest cell that intersects the shape and subdivide it. // // result contains the cells that will be part of the output, while pq // contains cells that we may still subdivide further. Cells that are // entirely contained within the region are immediately added to the output, // while cells that do not intersect the region are immediately discarded. // Therefore pq only contains cells that partially intersect the region. // Candidates are prioritized first according to cell size (larger cells // first), then by the number of intersecting children they have (fewest // children first), and then by the number of fully contained children // (fewest children first). func (c *coverer) coveringInternal(region Region) { c.region = region c.initialCandidates() for c.pq.Len() > 0 && (!c.interiorCovering || len(c.result) < c.maxCells) { cand := heap.Pop(&c.pq).(*candidate) // For interior covering we keep subdividing no matter how many children // candidate has. If we reach MaxCells before expanding all children, // we will just use some of them. // For exterior covering we cannot do this, because result has to cover the // whole region, so all children have to be used. // candidate.numChildren == 1 case takes care of the situation when we // already have more than MaxCells in result (minLevel is too high). // Subdividing of the candidate with one child does no harm in this case. if c.interiorCovering || int(cand.cell.level) < c.minLevel || cand.numChildren == 1 || len(c.result)+c.pq.Len()+cand.numChildren <= c.maxCells { for _, child := range cand.children { if !c.interiorCovering || len(c.result) < c.maxCells { c.addCandidate(child) } } } else { cand.terminal = true c.addCandidate(cand) } } c.pq.Reset() c.region = nil } // newCoverer returns an instance of coverer. func (rc *RegionCoverer) newCoverer() *coverer { return &coverer{ minLevel: maxInt(0, minInt(maxLevel, rc.MinLevel)), maxLevel: maxInt(0, minInt(maxLevel, rc.MaxLevel)), levelMod: maxInt(1, minInt(3, rc.LevelMod)), maxCells: rc.MaxCells, } } // Covering returns a CellUnion that covers the given region and satisfies the various restrictions. func (rc *RegionCoverer) Covering(region Region) CellUnion { covering := rc.CellUnion(region) covering.Denormalize(maxInt(0, minInt(maxLevel, rc.MinLevel)), maxInt(1, minInt(3, rc.LevelMod))) return covering } // InteriorCovering returns a CellUnion that is contained within the given region and satisfies the various restrictions. func (rc *RegionCoverer) InteriorCovering(region Region) CellUnion { intCovering := rc.InteriorCellUnion(region) intCovering.Denormalize(maxInt(0, minInt(maxLevel, rc.MinLevel)), maxInt(1, minInt(3, rc.LevelMod))) return intCovering } // CellUnion returns a normalized CellUnion that covers the given region and // satisfies the restrictions except for minLevel and levelMod. These criteria // cannot be satisfied using a cell union because cell unions are // automatically normalized by replacing four child cells with their parent // whenever possible. (Note that the list of cell ids passed to the CellUnion // constructor does in fact satisfy all the given restrictions.) func (rc *RegionCoverer) CellUnion(region Region) CellUnion { c := rc.newCoverer() c.coveringInternal(region) cu := c.result cu.Normalize() return cu } // InteriorCellUnion returns a normalized CellUnion that is contained within the given region and // satisfies the restrictions except for minLevel and levelMod. These criteria // cannot be satisfied using a cell union because cell unions are // automatically normalized by replacing four child cells with their parent // whenever possible. (Note that the list of cell ids passed to the CellUnion // constructor does in fact satisfy all the given restrictions.) func (rc *RegionCoverer) InteriorCellUnion(region Region) CellUnion { c := rc.newCoverer() c.interiorCovering = true c.coveringInternal(region) cu := c.result cu.Normalize() return cu } // FastCovering returns a CellUnion that covers the given region similar to Covering, // except that this method is much faster and the coverings are not as tight. // All of the usual parameters are respected (MaxCells, MinLevel, MaxLevel, and LevelMod), // except that the implementation makes no attempt to take advantage of large values of // MaxCells. (A small number of cells will always be returned.) // // This function is useful as a starting point for algorithms that // recursively subdivide cells. func (rc *RegionCoverer) FastCovering(region Region) CellUnion { c := rc.newCoverer() cu := CellUnion(region.CellUnionBound()) c.normalizeCovering(&cu) return cu } // normalizeCovering normalizes the "covering" so that it conforms to the current covering // parameters (MaxCells, minLevel, maxLevel, and levelMod). // This method makes no attempt to be optimal. In particular, if // minLevel > 0 or levelMod > 1 then it may return more than the // desired number of cells even when this isn't necessary. // // Note that when the covering parameters have their default values, almost // all of the code in this function is skipped. func (c *coverer) normalizeCovering(covering *CellUnion) { // If any cells are too small, or don't satisfy levelMod, then replace them with ancestors. if c.maxLevel < maxLevel || c.levelMod > 1 { for i, ci := range *covering { level := ci.Level() newLevel := c.adjustLevel(minInt(level, c.maxLevel)) if newLevel != level { (*covering)[i] = ci.Parent(newLevel) } } } // Sort the cells and simplify them. covering.Normalize() // If there are still too many cells, then repeatedly replace two adjacent // cells in CellID order by their lowest common ancestor. for len(*covering) > c.maxCells { bestIndex := -1 bestLevel := -1 for i := 0; i+1 < len(*covering); i++ { level, ok := (*covering)[i].CommonAncestorLevel((*covering)[i+1]) if !ok { continue } level = c.adjustLevel(level) if level > bestLevel { bestLevel = level bestIndex = i } } if bestLevel < c.minLevel { break } (*covering)[bestIndex] = (*covering)[bestIndex].Parent(bestLevel) covering.Normalize() } // Make sure that the covering satisfies minLevel and levelMod, // possibly at the expense of satisfying MaxCells. if c.minLevel > 0 || c.levelMod > 1 { covering.Denormalize(c.minLevel, c.levelMod) } } // SimpleRegionCovering returns a set of cells at the given level that cover // the connected region and a starting point on the boundary or inside the // region. The cells are returned in arbitrary order. // // Note that this method is not faster than the regular Covering // method for most region types, such as Cap or Polygon, and in fact it // can be much slower when the output consists of a large number of cells. // Currently it can be faster at generating coverings of long narrow regions // such as polylines, but this may change in the future. func SimpleRegionCovering(region Region, start Point, level int) []CellID { return FloodFillRegionCovering(region, cellIDFromPoint(start).Parent(level)) } // FloodFillRegionCovering returns all edge-connected cells at the same level as // the given CellID that intersect the given region, in arbitrary order. func FloodFillRegionCovering(region Region, start CellID) []CellID { var output []CellID all := map[CellID]bool{ start: true, } frontier := []CellID{start} for len(frontier) > 0 { id := frontier[len(frontier)-1] frontier = frontier[:len(frontier)-1] if !region.IntersectsCell(CellFromCellID(id)) { continue } output = append(output, id) for _, nbr := range id.EdgeNeighbors() { if !all[nbr] { all[nbr] = true frontier = append(frontier, nbr) } } } return output } // TODO(roberts): The differences from the C++ version // finish up FastCovering to match C++ // IsCanonical // CanonicalizeCovering // containsAllChildren // replaceCellsWithAncestor