mirror of
https://github.com/superseriousbusiness/gotosocial.git
synced 2024-12-05 10:02:47 +00:00
98263a7de6
* start fixing up tests * fix up tests + automate with drone * fiddle with linting * messing about with drone.yml * some more fiddling * hmmm * add cache * add vendor directory * verbose * ci updates * update some little things * update sig
945 lines
34 KiB
Go
945 lines
34 KiB
Go
// 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
|
|
}
|