paddedcell.go

Documentation: github.com/golang/geo/s2

     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  import (
    18  	"github.com/golang/geo/r1"
    19  	"github.com/golang/geo/r2"
    20  )
    21  
    22  // PaddedCell represents a Cell whose (u,v)-range has been expanded on
    23  // all sides by a given amount of "padding". Unlike Cell, its methods and
    24  // representation are optimized for clipping edges against Cell boundaries
    25  // to determine which cells are intersected by a given set of edges.
    26  type PaddedCell struct {
    27  	id          CellID
    28  	padding     float64
    29  	bound       r2.Rect
    30  	middle      r2.Rect // A rect in (u, v)-space that belongs to all four children.
    31  	iLo, jLo    int     // Minimum (i,j)-coordinates of this cell before padding
    32  	orientation int     // Hilbert curve orientation of this cell.
    33  	level       int
    34  }
    35  
    36  // PaddedCellFromCellID constructs a padded cell with the given padding.
    37  func PaddedCellFromCellID(id CellID, padding float64) *PaddedCell {
    38  	p := &PaddedCell{
    39  		id:      id,
    40  		padding: padding,
    41  		middle:  r2.EmptyRect(),
    42  	}
    43  
    44  	// Fast path for constructing a top-level face (the most common case).
    45  	if id.isFace() {
    46  		limit := padding + 1
    47  		p.bound = r2.Rect{r1.Interval{-limit, limit}, r1.Interval{-limit, limit}}
    48  		p.middle = r2.Rect{r1.Interval{-padding, padding}, r1.Interval{-padding, padding}}
    49  		p.orientation = id.Face() & 1
    50  		return p
    51  	}
    52  
    53  	_, p.iLo, p.jLo, p.orientation = id.faceIJOrientation()
    54  	p.level = id.Level()
    55  	p.bound = ijLevelToBoundUV(p.iLo, p.jLo, p.level).ExpandedByMargin(padding)
    56  	ijSize := sizeIJ(p.level)
    57  	p.iLo &= -ijSize
    58  	p.jLo &= -ijSize
    59  
    60  	return p
    61  }
    62  
    63  // PaddedCellFromParentIJ constructs the child of parent with the given (i,j) index.
    64  // The four child cells have indices of (0,0), (0,1), (1,0), (1,1), where the i and j
    65  // indices correspond to increasing u- and v-values respectively.
    66  func PaddedCellFromParentIJ(parent *PaddedCell, i, j int) *PaddedCell {
    67  	// Compute the position and orientation of the child incrementally from the
    68  	// orientation of the parent.
    69  	pos := ijToPos[parent.orientation][2*i+j]
    70  
    71  	p := &PaddedCell{
    72  		id:          parent.id.Children()[pos],
    73  		padding:     parent.padding,
    74  		bound:       parent.bound,
    75  		orientation: parent.orientation ^ posToOrientation[pos],
    76  		level:       parent.level + 1,
    77  		middle:      r2.EmptyRect(),
    78  	}
    79  
    80  	ijSize := sizeIJ(p.level)
    81  	p.iLo = parent.iLo + i*ijSize
    82  	p.jLo = parent.jLo + j*ijSize
    83  
    84  	// For each child, one corner of the bound is taken directly from the parent
    85  	// while the diagonally opposite corner is taken from middle().
    86  	middle := parent.Middle()
    87  	if i == 1 {
    88  		p.bound.X.Lo = middle.X.Lo
    89  	} else {
    90  		p.bound.X.Hi = middle.X.Hi
    91  	}
    92  	if j == 1 {
    93  		p.bound.Y.Lo = middle.Y.Lo
    94  	} else {
    95  		p.bound.Y.Hi = middle.Y.Hi
    96  	}
    97  
    98  	return p
    99  }
   100  
   101  // CellID returns the CellID this padded cell represents.
   102  func (p PaddedCell) CellID() CellID {
   103  	return p.id
   104  }
   105  
   106  // Padding returns the amount of padding on this cell.
   107  func (p PaddedCell) Padding() float64 {
   108  	return p.padding
   109  }
   110  
   111  // Level returns the level this cell is at.
   112  func (p PaddedCell) Level() int {
   113  	return p.level
   114  }
   115  
   116  // Center returns the center of this cell.
   117  func (p PaddedCell) Center() Point {
   118  	ijSize := sizeIJ(p.level)
   119  	si := uint32(2*p.iLo + ijSize)
   120  	ti := uint32(2*p.jLo + ijSize)
   121  	return Point{faceSiTiToXYZ(p.id.Face(), si, ti).Normalize()}
   122  }
   123  
   124  // Middle returns the rectangle in the middle of this cell that belongs to
   125  // all four of its children in (u,v)-space.
   126  func (p *PaddedCell) Middle() r2.Rect {
   127  	// We compute this field lazily because it is not needed the majority of the
   128  	// time (i.e., for cells where the recursion terminates).
   129  	if p.middle.IsEmpty() {
   130  		ijSize := sizeIJ(p.level)
   131  		u := stToUV(siTiToST(uint32(2*p.iLo + ijSize)))
   132  		v := stToUV(siTiToST(uint32(2*p.jLo + ijSize)))
   133  		p.middle = r2.Rect{
   134  			r1.Interval{u - p.padding, u + p.padding},
   135  			r1.Interval{v - p.padding, v + p.padding},
   136  		}
   137  	}
   138  	return p.middle
   139  }
   140  
   141  // Bound returns the bounds for this cell in (u,v)-space including padding.
   142  func (p PaddedCell) Bound() r2.Rect {
   143  	return p.bound
   144  }
   145  
   146  // ChildIJ returns the (i,j) coordinates for the child cell at the given traversal
   147  // position. The traversal position corresponds to the order in which child
   148  // cells are visited by the Hilbert curve.
   149  func (p PaddedCell) ChildIJ(pos int) (i, j int) {
   150  	ij := posToIJ[p.orientation][pos]
   151  	return ij >> 1, ij & 1
   152  }
   153  
   154  // EntryVertex return the vertex where the space-filling curve enters this cell.
   155  func (p PaddedCell) EntryVertex() Point {
   156  	// The curve enters at the (0,0) vertex unless the axis directions are
   157  	// reversed, in which case it enters at the (1,1) vertex.
   158  	i := p.iLo
   159  	j := p.jLo
   160  	if p.orientation&invertMask != 0 {
   161  		ijSize := sizeIJ(p.level)
   162  		i += ijSize
   163  		j += ijSize
   164  	}
   165  	return Point{faceSiTiToXYZ(p.id.Face(), uint32(2*i), uint32(2*j)).Normalize()}
   166  }
   167  
   168  // ExitVertex returns the vertex where the space-filling curve exits this cell.
   169  func (p PaddedCell) ExitVertex() Point {
   170  	// The curve exits at the (1,0) vertex unless the axes are swapped or
   171  	// inverted but not both, in which case it exits at the (0,1) vertex.
   172  	i := p.iLo
   173  	j := p.jLo
   174  	ijSize := sizeIJ(p.level)
   175  	if p.orientation == 0 || p.orientation == swapMask+invertMask {
   176  		i += ijSize
   177  	} else {
   178  		j += ijSize
   179  	}
   180  	return Point{faceSiTiToXYZ(p.id.Face(), uint32(2*i), uint32(2*j)).Normalize()}
   181  }
   182  
   183  // ShrinkToFit returns the smallest CellID that contains all descendants of this
   184  // padded cell whose bounds intersect the given rect. For algorithms that use
   185  // recursive subdivision to find the cells that intersect a particular object, this
   186  // method can be used to skip all of the initial subdivision steps where only
   187  // one child needs to be expanded.
   188  //
   189  // Note that this method is not the same as returning the smallest cell that contains
   190  // the intersection of this cell with rect. Because of the padding, even if one child
   191  // completely contains rect it is still possible that a neighboring child may also
   192  // intersect the given rect.
   193  //
   194  // The provided Rect must intersect the bounds of this cell.
   195  func (p *PaddedCell) ShrinkToFit(rect r2.Rect) CellID {
   196  	// Quick rejection test: if rect contains the center of this cell along
   197  	// either axis, then no further shrinking is possible.
   198  	if p.level == 0 {
   199  		// Fast path (most calls to this function start with a face cell).
   200  		if rect.X.Contains(0) || rect.Y.Contains(0) {
   201  			return p.id
   202  		}
   203  	}
   204  
   205  	ijSize := sizeIJ(p.level)
   206  	if rect.X.Contains(stToUV(siTiToST(uint32(2*p.iLo+ijSize)))) ||
   207  		rect.Y.Contains(stToUV(siTiToST(uint32(2*p.jLo+ijSize)))) {
   208  		return p.id
   209  	}
   210  
   211  	// Otherwise we expand rect by the given padding on all sides and find
   212  	// the range of coordinates that it spans along the i- and j-axes. We then
   213  	// compute the highest bit position at which the min and max coordinates
   214  	// differ. This corresponds to the first cell level at which at least two
   215  	// children intersect rect.
   216  
   217  	// Increase the padding to compensate for the error in uvToST.
   218  	// (The constant below is a provable upper bound on the additional error.)
   219  	padded := rect.ExpandedByMargin(p.padding + 1.5*dblEpsilon)
   220  	iMin, jMin := p.iLo, p.jLo // Min i- or j- coordinate spanned by padded
   221  	var iXor, jXor int         // XOR of the min and max i- or j-coordinates
   222  
   223  	if iMin < stToIJ(uvToST(padded.X.Lo)) {
   224  		iMin = stToIJ(uvToST(padded.X.Lo))
   225  	}
   226  	if a, b := p.iLo+ijSize-1, stToIJ(uvToST(padded.X.Hi)); a <= b {
   227  		iXor = iMin ^ a
   228  	} else {
   229  		iXor = iMin ^ b
   230  	}
   231  
   232  	if jMin < stToIJ(uvToST(padded.Y.Lo)) {
   233  		jMin = stToIJ(uvToST(padded.Y.Lo))
   234  	}
   235  	if a, b := p.jLo+ijSize-1, stToIJ(uvToST(padded.Y.Hi)); a <= b {
   236  		jXor = jMin ^ a
   237  	} else {
   238  		jXor = jMin ^ b
   239  	}
   240  
   241  	// Compute the highest bit position where the two i- or j-endpoints differ,
   242  	// and then choose the cell level that includes both of these endpoints. So
   243  	// if both pairs of endpoints are equal we choose MaxLevel; if they differ
   244  	// only at bit 0, we choose (MaxLevel - 1), and so on.
   245  	levelMSB := uint64(((iXor | jXor) << 1) + 1)
   246  	level := MaxLevel - findMSBSetNonZero64(levelMSB)
   247  	if level <= p.level {
   248  		return p.id
   249  	}
   250  
   251  	return cellIDFromFaceIJ(p.id.Face(), iMin, jMin).Parent(level)
   252  }
   253
View as plain text