diff options
Diffstat (limited to 'vendor/github.com/golang/geo/s2/cellid.go')
| -rw-r--r-- | vendor/github.com/golang/geo/s2/cellid.go | 944 |
1 files changed, 0 insertions, 944 deletions
diff --git a/vendor/github.com/golang/geo/s2/cellid.go b/vendor/github.com/golang/geo/s2/cellid.go deleted file mode 100644 index c6cbaf2db..000000000 --- a/vendor/github.com/golang/geo/s2/cellid.go +++ /dev/null @@ -1,944 +0,0 @@ -// 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 -} |