gtsocial-umbx

cellid.go (34379B)

---

 // Copyright 2014 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 (
 	"bytes"
 	"fmt"
 	"io"
 	"math"
 	"sort"
 	"strconv"
 	"strings"
 
 	"github.com/golang/geo/r1"
 	"github.com/golang/geo/r2"
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
 )
 
 // CellID uniquely identifies a cell in the S2 cell decomposition.
 // The most significant 3 bits encode the face number (0-5). The
 // remaining 61 bits encode the position of the center of this cell
 // along the Hilbert curve on that face. The zero value and the value
 // (1<<64)-1 are invalid cell IDs. The first compares less than any
 // valid cell ID, the second as greater than any valid cell ID.
 //
 // Sequentially increasing cell IDs follow a continuous space-filling curve
 // over the entire sphere. They have the following properties:
 //
 //  - The ID of a cell at level k consists of a 3-bit face number followed
 //    by k bit pairs that recursively select one of the four children of
 //    each cell. The next bit is always 1, and all other bits are 0.
 //    Therefore, the level of a cell is determined by the position of its
 //    lowest-numbered bit that is turned on (for a cell at level k, this
 //    position is 2 * (maxLevel - k)).
 //
 //  - The ID of a parent cell is at the midpoint of the range of IDs spanned
 //    by its children (or by its descendants at any level).
 //
 // Leaf cells are often used to represent points on the unit sphere, and
 // this type provides methods for converting directly between these two
 // representations. For cells that represent 2D regions rather than
 // discrete point, it is better to use Cells.
 type CellID uint64
 
 // SentinelCellID is an invalid cell ID guaranteed to be larger than any
 // valid cell ID. It is used primarily by ShapeIndex. The value is also used
 // by some S2 types when encoding data.
 // Note that the sentinel's RangeMin == RangeMax == itself.
 const SentinelCellID = CellID(^uint64(0))
 
 // sortCellIDs sorts the slice of CellIDs in place.
 func sortCellIDs(ci []CellID) {
 	sort.Sort(cellIDs(ci))
 }
 
 // cellIDs implements the Sort interface for slices of CellIDs.
 type cellIDs []CellID
 
 func (c cellIDs) Len() int           { return len(c) }
 func (c cellIDs) Swap(i, j int)      { c[i], c[j] = c[j], c[i] }
 func (c cellIDs) Less(i, j int) bool { return c[i] < c[j] }
 
 // TODO(dsymonds): Some of these constants should probably be exported.
 const (
 	faceBits = 3
 	numFaces = 6
 
 	// This is the number of levels needed to specify a leaf cell.
 	maxLevel = 30
 
 	// The extra position bit (61 rather than 60) lets us encode each cell as its
 	// Hilbert curve position at the cell center (which is halfway along the
 	// portion of the Hilbert curve that fills that cell).
 	posBits = 2*maxLevel + 1
 
 	// The maximum index of a valid leaf cell plus one. The range of valid leaf
 	// cell indices is [0..maxSize-1].
 	maxSize = 1 << maxLevel
 
 	wrapOffset = uint64(numFaces) << posBits
 )
 
 // CellIDFromFacePosLevel returns a cell given its face in the range
 // [0,5], the 61-bit Hilbert curve position pos within that face, and
 // the level in the range [0,maxLevel]. The position in the cell ID
 // will be truncated to correspond to the Hilbert curve position at
 // the center of the returned cell.
 func CellIDFromFacePosLevel(face int, pos uint64, level int) CellID {
 	return CellID(uint64(face)<<posBits + pos | 1).Parent(level)
 }
 
 // CellIDFromFace returns the cell corresponding to a given S2 cube face.
 func CellIDFromFace(face int) CellID {
 	return CellID((uint64(face) << posBits) + lsbForLevel(0))
 }
 
 // CellIDFromLatLng returns the leaf cell containing ll.
 func CellIDFromLatLng(ll LatLng) CellID {
 	return cellIDFromPoint(PointFromLatLng(ll))
 }
 
 // CellIDFromToken returns a cell given a hex-encoded string of its uint64 ID.
 func CellIDFromToken(s string) CellID {
 	if len(s) > 16 {
 		return CellID(0)
 	}
 	n, err := strconv.ParseUint(s, 16, 64)
 	if err != nil {
 		return CellID(0)
 	}
 	// Equivalent to right-padding string with zeros to 16 characters.
 	if len(s) < 16 {
 		n = n << (4 * uint(16-len(s)))
 	}
 	return CellID(n)
 }
 
 // ToToken returns a hex-encoded string of the uint64 cell id, with leading
 // zeros included but trailing zeros stripped.
 func (ci CellID) ToToken() string {
 	s := strings.TrimRight(fmt.Sprintf("%016x", uint64(ci)), "0")
 	if len(s) == 0 {
 		return "X"
 	}
 	return s
 }
 
 // IsValid reports whether ci represents a valid cell.
 func (ci CellID) IsValid() bool {
 	return ci.Face() < numFaces && (ci.lsb()&0x1555555555555555 != 0)
 }
 
 // Face returns the cube face for this cell ID, in the range [0,5].
 func (ci CellID) Face() int { return int(uint64(ci) >> posBits) }
 
 // Pos returns the position along the Hilbert curve of this cell ID, in the range [0,2^posBits-1].
 func (ci CellID) Pos() uint64 { return uint64(ci) & (^uint64(0) >> faceBits) }
 
 // Level returns the subdivision level of this cell ID, in the range [0, maxLevel].
 func (ci CellID) Level() int {
 	return maxLevel - findLSBSetNonZero64(uint64(ci))>>1
 }
 
 // IsLeaf returns whether this cell ID is at the deepest level;
 // that is, the level at which the cells are smallest.
 func (ci CellID) IsLeaf() bool { return uint64(ci)&1 != 0 }
 
 // ChildPosition returns the child position (0..3) of this cell's
 // ancestor at the given level, relative to its parent.  The argument
 // should be in the range 1..kMaxLevel.  For example,
 // ChildPosition(1) returns the position of this cell's level-1
 // ancestor within its top-level face cell.
 func (ci CellID) ChildPosition(level int) int {
 	return int(uint64(ci)>>uint64(2*(maxLevel-level)+1)) & 3
 }
 
 // lsbForLevel returns the lowest-numbered bit that is on for cells at the given level.
 func lsbForLevel(level int) uint64 { return 1 << uint64(2*(maxLevel-level)) }
 
 // Parent returns the cell at the given level, which must be no greater than the current level.
 func (ci CellID) Parent(level int) CellID {
 	lsb := lsbForLevel(level)
 	return CellID((uint64(ci) & -lsb) | lsb)
 }
 
 // immediateParent is cheaper than Parent, but assumes !ci.isFace().
 func (ci CellID) immediateParent() CellID {
 	nlsb := CellID(ci.lsb() << 2)
 	return (ci & -nlsb) | nlsb
 }
 
 // isFace returns whether this is a top-level (face) cell.
 func (ci CellID) isFace() bool { return uint64(ci)&(lsbForLevel(0)-1) == 0 }
 
 // lsb returns the least significant bit that is set.
 func (ci CellID) lsb() uint64 { return uint64(ci) & -uint64(ci) }
 
 // Children returns the four immediate children of this cell.
 // If ci is a leaf cell, it returns four identical cells that are not the children.
 func (ci CellID) Children() [4]CellID {
 	var ch [4]CellID
 	lsb := CellID(ci.lsb())
 	ch[0] = ci - lsb + lsb>>2
 	lsb >>= 1
 	ch[1] = ch[0] + lsb
 	ch[2] = ch[1] + lsb
 	ch[3] = ch[2] + lsb
 	return ch
 }
 
 func sizeIJ(level int) int {
 	return 1 << uint(maxLevel-level)
 }
 
 // EdgeNeighbors returns the four cells that are adjacent across the cell's four edges.
 // Edges 0, 1, 2, 3 are in the down, right, up, left directions in the face space.
 // All neighbors are guaranteed to be distinct.
 func (ci CellID) EdgeNeighbors() [4]CellID {
 	level := ci.Level()
 	size := sizeIJ(level)
 	f, i, j, _ := ci.faceIJOrientation()
 	return [4]CellID{
 		cellIDFromFaceIJWrap(f, i, j-size).Parent(level),
 		cellIDFromFaceIJWrap(f, i+size, j).Parent(level),
 		cellIDFromFaceIJWrap(f, i, j+size).Parent(level),
 		cellIDFromFaceIJWrap(f, i-size, j).Parent(level),
 	}
 }
 
 // VertexNeighbors returns the neighboring cellIDs with vertex closest to this cell at the given level.
 // (Normally there are four neighbors, but the closest vertex may only have three neighbors if it is one of
 // the 8 cube vertices.)
 func (ci CellID) VertexNeighbors(level int) []CellID {
 	halfSize := sizeIJ(level + 1)
 	size := halfSize << 1
 	f, i, j, _ := ci.faceIJOrientation()
 
 	var isame, jsame bool
 	var ioffset, joffset int
 	if i&halfSize != 0 {
 		ioffset = size
 		isame = (i + size) < maxSize
 	} else {
 		ioffset = -size
 		isame = (i - size) >= 0
 	}
 	if j&halfSize != 0 {
 		joffset = size
 		jsame = (j + size) < maxSize
 	} else {
 		joffset = -size
 		jsame = (j - size) >= 0
 	}
 
 	results := []CellID{
 		ci.Parent(level),
 		cellIDFromFaceIJSame(f, i+ioffset, j, isame).Parent(level),
 		cellIDFromFaceIJSame(f, i, j+joffset, jsame).Parent(level),
 	}
 
 	if isame || jsame {
 		results = append(results, cellIDFromFaceIJSame(f, i+ioffset, j+joffset, isame && jsame).Parent(level))
 	}
 
 	return results
 }
 
 // AllNeighbors returns all neighbors of this cell at the given level. Two
 // cells X and Y are neighbors if their boundaries intersect but their
 // interiors do not. In particular, two cells that intersect at a single
 // point are neighbors. Note that for cells adjacent to a face vertex, the
 // same neighbor may be returned more than once. There could be up to eight
 // neighbors including the diagonal ones that share the vertex.
 //
 // This requires level >= ci.Level().
 func (ci CellID) AllNeighbors(level int) []CellID {
 	var neighbors []CellID
 
 	face, i, j, _ := ci.faceIJOrientation()
 
 	// Find the coordinates of the lower left-hand leaf cell. We need to
 	// normalize (i,j) to a known position within the cell because level
 	// may be larger than this cell's level.
 	size := sizeIJ(ci.Level())
 	i &= -size
 	j &= -size
 
 	nbrSize := sizeIJ(level)
 
 	// We compute the top-bottom, left-right, and diagonal neighbors in one
 	// pass. The loop test is at the end of the loop to avoid 32-bit overflow.
 	for k := -nbrSize; ; k += nbrSize {
 		var sameFace bool
 		if k < 0 {
 			sameFace = (j+k >= 0)
 		} else if k >= size {
 			sameFace = (j+k < maxSize)
 		} else {
 			sameFace = true
 			// Top and bottom neighbors.
 			neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+k, j-nbrSize,
 				j-size >= 0).Parent(level))
 			neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+k, j+size,
 				j+size < maxSize).Parent(level))
 		}
 
 		// Left, right, and diagonal neighbors.
 		neighbors = append(neighbors, cellIDFromFaceIJSame(face, i-nbrSize, j+k,
 			sameFace && i-size >= 0).Parent(level))
 		neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+size, j+k,
 			sameFace && i+size < maxSize).Parent(level))
 
 		if k >= size {
 			break
 		}
 	}
 
 	return neighbors
 }
 
 // RangeMin returns the minimum CellID that is contained within this cell.
 func (ci CellID) RangeMin() CellID { return CellID(uint64(ci) - (ci.lsb() - 1)) }
 
 // RangeMax returns the maximum CellID that is contained within this cell.
 func (ci CellID) RangeMax() CellID { return CellID(uint64(ci) + (ci.lsb() - 1)) }
 
 // Contains returns true iff the CellID contains oci.
 func (ci CellID) Contains(oci CellID) bool {
 	return uint64(ci.RangeMin()) <= uint64(oci) && uint64(oci) <= uint64(ci.RangeMax())
 }
 
 // Intersects returns true iff the CellID intersects oci.
 func (ci CellID) Intersects(oci CellID) bool {
 	return uint64(oci.RangeMin()) <= uint64(ci.RangeMax()) && uint64(oci.RangeMax()) >= uint64(ci.RangeMin())
 }
 
 // String returns the string representation of the cell ID in the form "1/3210".
 func (ci CellID) String() string {
 	if !ci.IsValid() {
 		return "Invalid: " + strconv.FormatInt(int64(ci), 16)
 	}
 	var b bytes.Buffer
 	b.WriteByte("012345"[ci.Face()]) // values > 5 will have been picked off by !IsValid above
 	b.WriteByte('/')
 	for level := 1; level <= ci.Level(); level++ {
 		b.WriteByte("0123"[ci.ChildPosition(level)])
 	}
 	return b.String()
 }
 
 // cellIDFromString returns a CellID from a string in the form "1/3210".
 func cellIDFromString(s string) CellID {
 	level := len(s) - 2
 	if level < 0 || level > maxLevel {
 		return CellID(0)
 	}
 	face := int(s[0] - '0')
 	if face < 0 || face > 5 || s[1] != '/' {
 		return CellID(0)
 	}
 	id := CellIDFromFace(face)
 	for i := 2; i < len(s); i++ {
 		childPos := s[i] - '0'
 		if childPos < 0 || childPos > 3 {
 			return CellID(0)
 		}
 		id = id.Children()[childPos]
 	}
 	return id
 }
 
 // Point returns the center of the s2 cell on the sphere as a Point.
 // The maximum directional error in Point (compared to the exact
 // mathematical result) is 1.5 * dblEpsilon radians, and the maximum length
 // error is 2 * dblEpsilon (the same as Normalize).
 func (ci CellID) Point() Point { return Point{ci.rawPoint().Normalize()} }
 
 // LatLng returns the center of the s2 cell on the sphere as a LatLng.
 func (ci CellID) LatLng() LatLng { return LatLngFromPoint(Point{ci.rawPoint()}) }
 
 // ChildBegin returns the first child in a traversal of the children of this cell, in Hilbert curve order.
 //
 //    for ci := c.ChildBegin(); ci != c.ChildEnd(); ci = ci.Next() {
 //        ...
 //    }
 func (ci CellID) ChildBegin() CellID {
 	ol := ci.lsb()
 	return CellID(uint64(ci) - ol + ol>>2)
 }
 
 // ChildBeginAtLevel returns the first cell in a traversal of children a given level deeper than this cell, in
 // Hilbert curve order. The given level must be no smaller than the cell's level.
 // See ChildBegin for example use.
 func (ci CellID) ChildBeginAtLevel(level int) CellID {
 	return CellID(uint64(ci) - ci.lsb() + lsbForLevel(level))
 }
 
 // ChildEnd returns the first cell after a traversal of the children of this cell in Hilbert curve order.
 // The returned cell may be invalid.
 func (ci CellID) ChildEnd() CellID {
 	ol := ci.lsb()
 	return CellID(uint64(ci) + ol + ol>>2)
 }
 
 // ChildEndAtLevel returns the first cell after the last child in a traversal of children a given level deeper
 // than this cell, in Hilbert curve order.
 // The given level must be no smaller than the cell's level.
 // The returned cell may be invalid.
 func (ci CellID) ChildEndAtLevel(level int) CellID {
 	return CellID(uint64(ci) + ci.lsb() + lsbForLevel(level))
 }
 
 // Next returns the next cell along the Hilbert curve.
 // This is expected to be used with ChildBegin and ChildEnd,
 // or ChildBeginAtLevel and ChildEndAtLevel.
 func (ci CellID) Next() CellID {
 	return CellID(uint64(ci) + ci.lsb()<<1)
 }
 
 // Prev returns the previous cell along the Hilbert curve.
 func (ci CellID) Prev() CellID {
 	return CellID(uint64(ci) - ci.lsb()<<1)
 }
 
 // NextWrap returns the next cell along the Hilbert curve, wrapping from last to
 // first as necessary. This should not be used with ChildBegin and ChildEnd.
 func (ci CellID) NextWrap() CellID {
 	n := ci.Next()
 	if uint64(n) < wrapOffset {
 		return n
 	}
 	return CellID(uint64(n) - wrapOffset)
 }
 
 // PrevWrap returns the previous cell along the Hilbert curve, wrapping around from
 // first to last as necessary. This should not be used with ChildBegin and ChildEnd.
 func (ci CellID) PrevWrap() CellID {
 	p := ci.Prev()
 	if uint64(p) < wrapOffset {
 		return p
 	}
 	return CellID(uint64(p) + wrapOffset)
 }
 
 // AdvanceWrap advances or retreats the indicated number of steps along the
 // Hilbert curve at the current level and returns the new position. The
 // position wraps between the first and last faces as necessary.
 func (ci CellID) AdvanceWrap(steps int64) CellID {
 	if steps == 0 {
 		return ci
 	}
 
 	// We clamp the number of steps if necessary to ensure that we do not
 	// advance past the End() or before the Begin() of this level.
 	shift := uint(2*(maxLevel-ci.Level()) + 1)
 	if steps < 0 {
 		if min := -int64(uint64(ci) >> shift); steps < min {
 			wrap := int64(wrapOffset >> shift)
 			steps %= wrap
 			if steps < min {
 				steps += wrap
 			}
 		}
 	} else {
 		// Unlike Advance(), we don't want to return End(level).
 		if max := int64((wrapOffset - uint64(ci)) >> shift); steps > max {
 			wrap := int64(wrapOffset >> shift)
 			steps %= wrap
 			if steps > max {
 				steps -= wrap
 			}
 		}
 	}
 
 	// If steps is negative, then shifting it left has undefined behavior.
 	// Cast to uint64 for a 2's complement answer.
 	return CellID(uint64(ci) + (uint64(steps) << shift))
 }
 
 // Encode encodes the CellID.
 func (ci CellID) Encode(w io.Writer) error {
 	e := &encoder{w: w}
 	ci.encode(e)
 	return e.err
 }
 
 func (ci CellID) encode(e *encoder) {
 	e.writeUint64(uint64(ci))
 }
 
 // Decode decodes the CellID.
 func (ci *CellID) Decode(r io.Reader) error {
 	d := &decoder{r: asByteReader(r)}
 	ci.decode(d)
 	return d.err
 }
 
 func (ci *CellID) decode(d *decoder) {
 	*ci = CellID(d.readUint64())
 }
 
 // TODO: the methods below are not exported yet.  Settle on the entire API design
 // before doing this.  Do we want to mirror the C++ one as closely as possible?
 
 // distanceFromBegin returns the number of steps along the Hilbert curve that
 // this cell is from the first node in the S2 hierarchy at our level. (i.e.,
 // FromFace(0).ChildBeginAtLevel(ci.Level())). This is analogous to Pos(), but
 // for this cell's level.
 // The return value is always non-negative.
 func (ci CellID) distanceFromBegin() int64 {
 	return int64(ci >> uint64(2*(maxLevel-ci.Level())+1))
 }
 
 // rawPoint returns an unnormalized r3 vector from the origin through the center
 // of the s2 cell on the sphere.
 func (ci CellID) rawPoint() r3.Vector {
 	face, si, ti := ci.faceSiTi()
 	return faceUVToXYZ(face, stToUV((0.5/maxSize)*float64(si)), stToUV((0.5/maxSize)*float64(ti)))
 }
 
 // faceSiTi returns the Face/Si/Ti coordinates of the center of the cell.
 func (ci CellID) faceSiTi() (face int, si, ti uint32) {
 	face, i, j, _ := ci.faceIJOrientation()
 	delta := 0
 	if ci.IsLeaf() {
 		delta = 1
 	} else {
 		if (i^(int(ci)>>2))&1 != 0 {
 			delta = 2
 		}
 	}
 	return face, uint32(2*i + delta), uint32(2*j + delta)
 }
 
 // faceIJOrientation uses the global lookupIJ table to unfiddle the bits of ci.
 func (ci CellID) faceIJOrientation() (f, i, j, orientation int) {
 	f = ci.Face()
 	orientation = f & swapMask
 	nbits := maxLevel - 7*lookupBits // first iteration
 
 	// Each iteration maps 8 bits of the Hilbert curve position into
 	// 4 bits of "i" and "j". The lookup table transforms a key of the
 	// form "ppppppppoo" to a value of the form "iiiijjjjoo", where the
 	// letters [ijpo] represents bits of "i", "j", the Hilbert curve
 	// position, and the Hilbert curve orientation respectively.
 	//
 	// On the first iteration we need to be careful to clear out the bits
 	// representing the cube face.
 	for k := 7; k >= 0; k-- {
 		orientation += (int(uint64(ci)>>uint64(k*2*lookupBits+1)) & ((1 << uint(2*nbits)) - 1)) << 2
 		orientation = lookupIJ[orientation]
 		i += (orientation >> (lookupBits + 2)) << uint(k*lookupBits)
 		j += ((orientation >> 2) & ((1 << lookupBits) - 1)) << uint(k*lookupBits)
 		orientation &= (swapMask | invertMask)
 		nbits = lookupBits // following iterations
 	}
 
 	// The position of a non-leaf cell at level "n" consists of a prefix of
 	// 2*n bits that identifies the cell, followed by a suffix of
 	// 2*(maxLevel-n)+1 bits of the form 10*. If n==maxLevel, the suffix is
 	// just "1" and has no effect. Otherwise, it consists of "10", followed
 	// by (maxLevel-n-1) repetitions of "00", followed by "0". The "10" has
 	// no effect, while each occurrence of "00" has the effect of reversing
 	// the swapMask bit.
 	if ci.lsb()&0x1111111111111110 != 0 {
 		orientation ^= swapMask
 	}
 
 	return
 }
 
 // cellIDFromFaceIJ returns a leaf cell given its cube face (range 0..5) and IJ coordinates.
 func cellIDFromFaceIJ(f, i, j int) CellID {
 	// Note that this value gets shifted one bit to the left at the end
 	// of the function.
 	n := uint64(f) << (posBits - 1)
 	// Alternating faces have opposite Hilbert curve orientations; this
 	// is necessary in order for all faces to have a right-handed
 	// coordinate system.
 	bits := f & swapMask
 	// Each iteration maps 4 bits of "i" and "j" into 8 bits of the Hilbert
 	// curve position.  The lookup table transforms a 10-bit key of the form
 	// "iiiijjjjoo" to a 10-bit value of the form "ppppppppoo", where the
 	// letters [ijpo] denote bits of "i", "j", Hilbert curve position, and
 	// Hilbert curve orientation respectively.
 	for k := 7; k >= 0; k-- {
 		mask := (1 << lookupBits) - 1
 		bits += ((i >> uint(k*lookupBits)) & mask) << (lookupBits + 2)
 		bits += ((j >> uint(k*lookupBits)) & mask) << 2
 		bits = lookupPos[bits]
 		n |= uint64(bits>>2) << (uint(k) * 2 * lookupBits)
 		bits &= (swapMask | invertMask)
 	}
 	return CellID(n*2 + 1)
 }
 
 func cellIDFromFaceIJWrap(f, i, j int) CellID {
 	// Convert i and j to the coordinates of a leaf cell just beyond the
 	// boundary of this face.  This prevents 32-bit overflow in the case
 	// of finding the neighbors of a face cell.
 	i = clampInt(i, -1, maxSize)
 	j = clampInt(j, -1, maxSize)
 
 	// We want to wrap these coordinates onto the appropriate adjacent face.
 	// The easiest way to do this is to convert the (i,j) coordinates to (x,y,z)
 	// (which yields a point outside the normal face boundary), and then call
 	// xyzToFaceUV to project back onto the correct face.
 	//
 	// The code below converts (i,j) to (si,ti), and then (si,ti) to (u,v) using
 	// the linear projection (u=2*s-1 and v=2*t-1).  (The code further below
 	// converts back using the inverse projection, s=0.5*(u+1) and t=0.5*(v+1).
 	// Any projection would work here, so we use the simplest.)  We also clamp
 	// the (u,v) coordinates so that the point is barely outside the
 	// [-1,1]x[-1,1] face rectangle, since otherwise the reprojection step
 	// (which divides by the new z coordinate) might change the other
 	// coordinates enough so that we end up in the wrong leaf cell.
 	const scale = 1.0 / maxSize
 	limit := math.Nextafter(1, 2)
 	u := math.Max(-limit, math.Min(limit, scale*float64((i<<1)+1-maxSize)))
 	v := math.Max(-limit, math.Min(limit, scale*float64((j<<1)+1-maxSize)))
 
 	// Find the leaf cell coordinates on the adjacent face, and convert
 	// them to a cell id at the appropriate level.
 	f, u, v = xyzToFaceUV(faceUVToXYZ(f, u, v))
 	return cellIDFromFaceIJ(f, stToIJ(0.5*(u+1)), stToIJ(0.5*(v+1)))
 }
 
 func cellIDFromFaceIJSame(f, i, j int, sameFace bool) CellID {
 	if sameFace {
 		return cellIDFromFaceIJ(f, i, j)
 	}
 	return cellIDFromFaceIJWrap(f, i, j)
 }
 
 // ijToSTMin converts the i- or j-index of a leaf cell to the minimum corresponding
 // s- or t-value contained by that cell. The argument must be in the range
 // [0..2**30], i.e. up to one position beyond the normal range of valid leaf
 // cell indices.
 func ijToSTMin(i int) float64 {
 	return float64(i) / float64(maxSize)
 }
 
 // stToIJ converts value in ST coordinates to a value in IJ coordinates.
 func stToIJ(s float64) int {
 	return clampInt(int(math.Floor(maxSize*s)), 0, maxSize-1)
 }
 
 // cellIDFromPoint returns a leaf cell containing point p. Usually there is
 // exactly one such cell, but for points along the edge of a cell, any
 // adjacent cell may be (deterministically) chosen. This is because
 // s2.CellIDs are considered to be closed sets. The returned cell will
 // always contain the given point, i.e.
 //
 //   CellFromPoint(p).ContainsPoint(p)
 //
 // is always true.
 func cellIDFromPoint(p Point) CellID {
 	f, u, v := xyzToFaceUV(r3.Vector{p.X, p.Y, p.Z})
 	i := stToIJ(uvToST(u))
 	j := stToIJ(uvToST(v))
 	return cellIDFromFaceIJ(f, i, j)
 }
 
 // ijLevelToBoundUV returns the bounds in (u,v)-space for the cell at the given
 // level containing the leaf cell with the given (i,j)-coordinates.
 func ijLevelToBoundUV(i, j, level int) r2.Rect {
 	cellSize := sizeIJ(level)
 	xLo := i & -cellSize
 	yLo := j & -cellSize
 
 	return r2.Rect{
 		X: r1.Interval{
 			Lo: stToUV(ijToSTMin(xLo)),
 			Hi: stToUV(ijToSTMin(xLo + cellSize)),
 		},
 		Y: r1.Interval{
 			Lo: stToUV(ijToSTMin(yLo)),
 			Hi: stToUV(ijToSTMin(yLo + cellSize)),
 		},
 	}
 }
 
 // Constants related to the bit mangling in the Cell ID.
 const (
 	lookupBits = 4
 	swapMask   = 0x01
 	invertMask = 0x02
 )
 
 // The following lookup tables are used to convert efficiently between an
 // (i,j) cell index and the corresponding position along the Hilbert curve.
 //
 // lookupPos maps 4 bits of "i", 4 bits of "j", and 2 bits representing the
 // orientation of the current cell into 8 bits representing the order in which
 // that subcell is visited by the Hilbert curve, plus 2 bits indicating the
 // new orientation of the Hilbert curve within that subcell. (Cell
 // orientations are represented as combination of swapMask and invertMask.)
 //
 // lookupIJ is an inverted table used for mapping in the opposite
 // direction.
 //
 // We also experimented with looking up 16 bits at a time (14 bits of position
 // plus 2 of orientation) but found that smaller lookup tables gave better
 // performance. (2KB fits easily in the primary cache.)
 var (
 	ijToPos = [4][4]int{
 		{0, 1, 3, 2}, // canonical order
 		{0, 3, 1, 2}, // axes swapped
 		{2, 3, 1, 0}, // bits inverted
 		{2, 1, 3, 0}, // swapped & inverted
 	}
 	posToIJ = [4][4]int{
 		{0, 1, 3, 2}, // canonical order:    (0,0), (0,1), (1,1), (1,0)
 		{0, 2, 3, 1}, // axes swapped:       (0,0), (1,0), (1,1), (0,1)
 		{3, 2, 0, 1}, // bits inverted:      (1,1), (1,0), (0,0), (0,1)
 		{3, 1, 0, 2}, // swapped & inverted: (1,1), (0,1), (0,0), (1,0)
 	}
 	posToOrientation = [4]int{swapMask, 0, 0, invertMask | swapMask}
 	lookupIJ         [1 << (2*lookupBits + 2)]int
 	lookupPos        [1 << (2*lookupBits + 2)]int
 )
 
 func init() {
 	initLookupCell(0, 0, 0, 0, 0, 0)
 	initLookupCell(0, 0, 0, swapMask, 0, swapMask)
 	initLookupCell(0, 0, 0, invertMask, 0, invertMask)
 	initLookupCell(0, 0, 0, swapMask|invertMask, 0, swapMask|invertMask)
 }
 
 // initLookupCell initializes the lookupIJ table at init time.
 func initLookupCell(level, i, j, origOrientation, pos, orientation int) {
 	if level == lookupBits {
 		ij := (i << lookupBits) + j
 		lookupPos[(ij<<2)+origOrientation] = (pos << 2) + orientation
 		lookupIJ[(pos<<2)+origOrientation] = (ij << 2) + orientation
 		return
 	}
 
 	level++
 	i <<= 1
 	j <<= 1
 	pos <<= 2
 	r := posToIJ[orientation]
 	initLookupCell(level, i+(r[0]>>1), j+(r[0]&1), origOrientation, pos, orientation^posToOrientation[0])
 	initLookupCell(level, i+(r[1]>>1), j+(r[1]&1), origOrientation, pos+1, orientation^posToOrientation[1])
 	initLookupCell(level, i+(r[2]>>1), j+(r[2]&1), origOrientation, pos+2, orientation^posToOrientation[2])
 	initLookupCell(level, i+(r[3]>>1), j+(r[3]&1), origOrientation, pos+3, orientation^posToOrientation[3])
 }
 
 // CommonAncestorLevel returns the level of the common ancestor of the two S2 CellIDs.
 func (ci CellID) CommonAncestorLevel(other CellID) (level int, ok bool) {
 	bits := uint64(ci ^ other)
 	if bits < ci.lsb() {
 		bits = ci.lsb()
 	}
 	if bits < other.lsb() {
 		bits = other.lsb()
 	}
 
 	msbPos := findMSBSetNonZero64(bits)
 	if msbPos > 60 {
 		return 0, false
 	}
 	return (60 - msbPos) >> 1, true
 }
 
 // Advance advances or retreats the indicated number of steps along the
 // Hilbert curve at the current level, and returns the new position. The
 // position is never advanced past End() or before Begin().
 func (ci CellID) Advance(steps int64) CellID {
 	if steps == 0 {
 		return ci
 	}
 
 	// We clamp the number of steps if necessary to ensure that we do not
 	// advance past the End() or before the Begin() of this level. Note that
 	// minSteps and maxSteps always fit in a signed 64-bit integer.
 	stepShift := uint(2*(maxLevel-ci.Level()) + 1)
 	if steps < 0 {
 		minSteps := -int64(uint64(ci) >> stepShift)
 		if steps < minSteps {
 			steps = minSteps
 		}
 	} else {
 		maxSteps := int64((wrapOffset + ci.lsb() - uint64(ci)) >> stepShift)
 		if steps > maxSteps {
 			steps = maxSteps
 		}
 	}
 	return ci + CellID(steps)<<stepShift
 }
 
 // centerST return the center of the CellID in (s,t)-space.
 func (ci CellID) centerST() r2.Point {
 	_, si, ti := ci.faceSiTi()
 	return r2.Point{siTiToST(si), siTiToST(ti)}
 }
 
 // sizeST returns the edge length of this CellID in (s,t)-space at the given level.
 func (ci CellID) sizeST(level int) float64 {
 	return ijToSTMin(sizeIJ(level))
 }
 
 // boundST returns the bound of this CellID in (s,t)-space.
 func (ci CellID) boundST() r2.Rect {
 	s := ci.sizeST(ci.Level())
 	return r2.RectFromCenterSize(ci.centerST(), r2.Point{s, s})
 }
 
 // centerUV returns the center of this CellID in (u,v)-space. Note that
 // the center of the cell is defined as the point at which it is recursively
 // subdivided into four children; in general, it is not at the midpoint of
 // the (u,v) rectangle covered by the cell.
 func (ci CellID) centerUV() r2.Point {
 	_, si, ti := ci.faceSiTi()
 	return r2.Point{stToUV(siTiToST(si)), stToUV(siTiToST(ti))}
 }
 
 // boundUV returns the bound of this CellID in (u,v)-space.
 func (ci CellID) boundUV() r2.Rect {
 	_, i, j, _ := ci.faceIJOrientation()
 	return ijLevelToBoundUV(i, j, ci.Level())
 }
 
 // expandEndpoint returns a new u-coordinate u' such that the distance from the
 // line u=u' to the given edge (u,v0)-(u,v1) is exactly the given distance
 // (which is specified as the sine of the angle corresponding to the distance).
 func expandEndpoint(u, maxV, sinDist float64) float64 {
 	// This is based on solving a spherical right triangle, similar to the
 	// calculation in Cap.RectBound.
 	// Given an edge of the form (u,v0)-(u,v1), let maxV = max(abs(v0), abs(v1)).
 	sinUShift := sinDist * math.Sqrt((1+u*u+maxV*maxV)/(1+u*u))
 	cosUShift := math.Sqrt(1 - sinUShift*sinUShift)
 	// The following is an expansion of tan(atan(u) + asin(sinUShift)).
 	return (cosUShift*u + sinUShift) / (cosUShift - sinUShift*u)
 }
 
 // expandedByDistanceUV returns a rectangle expanded in (u,v)-space so that it
 // contains all points within the given distance of the boundary, and return the
 // smallest such rectangle. If the distance is negative, then instead shrink this
 // rectangle so that it excludes all points within the given absolute distance
 // of the boundary.
 //
 // Distances are measured *on the sphere*, not in (u,v)-space. For example,
 // you can use this method to expand the (u,v)-bound of an CellID so that
 // it contains all points within 5km of the original cell. You can then
 // test whether a point lies within the expanded bounds like this:
 //
 //   if u, v, ok := faceXYZtoUV(face, point); ok && bound.ContainsPoint(r2.Point{u,v}) { ... }
 //
 // Limitations:
 //
 //  - Because the rectangle is drawn on one of the six cube-face planes
 //    (i.e., {x,y,z} = +/-1), it can cover at most one hemisphere. This
 //    limits the maximum amount that a rectangle can be expanded. For
 //    example, CellID bounds can be expanded safely by at most 45 degrees
 //    (about 5000 km on the Earth's surface).
 //
 //  - The implementation is not exact for negative distances. The resulting
 //    rectangle will exclude all points within the given distance of the
 //    boundary but may be slightly smaller than necessary.
 func expandedByDistanceUV(uv r2.Rect, distance s1.Angle) r2.Rect {
 	// Expand each of the four sides of the rectangle just enough to include all
 	// points within the given distance of that side. (The rectangle may be
 	// expanded by a different amount in (u,v)-space on each side.)
 	maxU := math.Max(math.Abs(uv.X.Lo), math.Abs(uv.X.Hi))
 	maxV := math.Max(math.Abs(uv.Y.Lo), math.Abs(uv.Y.Hi))
 	sinDist := math.Sin(float64(distance))
 	return r2.Rect{
 		X: r1.Interval{expandEndpoint(uv.X.Lo, maxV, -sinDist),
 			expandEndpoint(uv.X.Hi, maxV, sinDist)},
 		Y: r1.Interval{expandEndpoint(uv.Y.Lo, maxU, -sinDist),
 			expandEndpoint(uv.Y.Hi, maxU, sinDist)}}
 }
 
 // MaxTile returns the largest cell with the same RangeMin such that
 // RangeMax < limit.RangeMin. It returns limit if no such cell exists.
 // This method can be used to generate a small set of CellIDs that covers
 // a given range (a tiling). This example shows how to generate a tiling
 // for a semi-open range of leaf cells [start, limit):
 //
 //   for id := start.MaxTile(limit); id != limit; id = id.Next().MaxTile(limit)) { ... }
 //
 // Note that in general the cells in the tiling will be of different sizes;
 // they gradually get larger (near the middle of the range) and then
 // gradually get smaller as limit is approached.
 func (ci CellID) MaxTile(limit CellID) CellID {
 	start := ci.RangeMin()
 	if start >= limit.RangeMin() {
 		return limit
 	}
 
 	if ci.RangeMax() >= limit {
 		// The cell is too large, shrink it. Note that when generating coverings
 		// of CellID ranges, this loop usually executes only once. Also because
 		// ci.RangeMin() < limit.RangeMin(), we will always exit the loop by the
 		// time we reach a leaf cell.
 		for {
 			ci = ci.Children()[0]
 			if ci.RangeMax() < limit {
 				break
 			}
 		}
 		return ci
 	}
 
 	// The cell may be too small. Grow it if necessary. Note that generally
 	// this loop only iterates once.
 	for !ci.isFace() {
 		parent := ci.immediateParent()
 		if parent.RangeMin() != start || parent.RangeMax() >= limit {
 			break
 		}
 		ci = parent
 	}
 	return ci
 }
 
 // centerFaceSiTi returns the (face, si, ti) coordinates of the center of the cell.
 // Note that although (si,ti) coordinates span the range [0,2**31] in general,
 // the cell center coordinates are always in the range [1,2**31-1] and
 // therefore can be represented using a signed 32-bit integer.
 func (ci CellID) centerFaceSiTi() (face, si, ti int) {
 	// First we compute the discrete (i,j) coordinates of a leaf cell contained
 	// within the given cell. Given that cells are represented by the Hilbert
 	// curve position corresponding at their center, it turns out that the cell
 	// returned by faceIJOrientation is always one of two leaf cells closest
 	// to the center of the cell (unless the given cell is a leaf cell itself,
 	// in which case there is only one possibility).
 	//
 	// Given a cell of size s >= 2 (i.e. not a leaf cell), and letting (imin,
 	// jmin) be the coordinates of its lower left-hand corner, the leaf cell
 	// returned by faceIJOrientation is either (imin + s/2, jmin + s/2)
 	// (imin + s/2 - 1, jmin + s/2 - 1). The first case is the one we want.
 	// We can distinguish these two cases by looking at the low bit of i or
 	// j. In the second case the low bit is one, unless s == 2 (i.e. the
 	// level just above leaf cells) in which case the low bit is zero.
 	//
 	// In the code below, the expression ((i ^ (int(id) >> 2)) & 1) is true
 	// if we are in the second case described above.
 	face, i, j, _ := ci.faceIJOrientation()
 	delta := 0
 	if ci.IsLeaf() {
 		delta = 1
 	} else if (int64(i)^(int64(ci)>>2))&1 == 1 {
 		delta = 2
 	}
 
 	// Note that (2 * {i,j} + delta) will never overflow a 32-bit integer.
 	return face, 2*i + delta, 2*j + delta
 }