As for the definition of neighbors, neighbors are neighbors. Neighbors can be divided into two types: edge neighbors and point neighbors. The edges are adjacent to each other in four directions, up, down, left and right. Points that are adjacent to each other also have four directions, that is, four vertices that are adjacent to each other.
As shown in the figure above, the green area is the range represented by a quadtree with a point on it, the yellow area of the figure. Now I want to find the Hilbert neighbors of the yellow points on the quadtree. The black line is a Hilbert curve running through a quadtree. The Hilbert curve starts at 0 on the upper-left grid and ends at 63 on the upper-right grid.
The four red cells are neighbors of the edge of the yellow cell, and the four blue cells are neighbors of the vertex of the yellow cell. Therefore, the neighbors of the yellow cell are 8 cells, representing points 8, 9, 54, 11, 53, 30, 31, 32 respectively. You can see that these neighbors are not adjacent to each other at this point.
So how do we find the Hilbert neighbors of any point on the quadtree?
One side neighbors
The most direct idea of edge neighbor is to first obtain the coordinates of the center point (I, j), and then obtain the coordinates of the Cell adjacent to its edge (I + 1, j), (i-1, j), (I, J-1), (I, j + 1) through the coordinate system relationship.
In practice, too. But here’s where you need to switch. So we have to convert the points on the Hilbert curve into coordinates before we can compute the edge neighbors in the same way.
For the generation and data structure of CellID, see the author’s article “How is CellID generated in Google S2?”
According to the above ideas, the code is as follows:
func (ci CellID) EdgeNeighbors(a) [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),
}
}Copy the codeThe edges are numbered 0,1,2,3 counterclockwise, bottom, right, top, left.
Let’s look at the implementation in detail.
func sizeIJ(level int) int {
return 1 << uint(maxLevel-level)
}Copy the codeSizeIJ stores the size of the current Level grid edge. This size is relative to Level 30. For example, if level = 29, its size will be 2, indicating that the edge length of a level 29 cell is composed of 2 level 30 cells, so it is 2^2^=4 small cells. If level = 28, then the side length is 4 and consists of 16 small cells. And so on.
func (ci CellID) faceIJOrientation(a) (f, i, j, orientation int) {
f = ci.Face()
orientation = f & swapMask
nbits := maxLevel - 7*lookupBits // first iteration
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 >> 2A) & ((1 << lookupBits) - 1)) < <uint(k*lookupBits)
orientation &= (swapMask | invertMask)
nbits = lookupBits // following iterations
}
// The following judgment is explained in detail
if ci.lsb()&0x1111111111111110! =0 {
orientation ^= swapMask
}
return
}Copy the codeThe method is to factor CellID back into the original I’s and j’s. The specific process here is described in the author’s article how is CellID generated in Google S2? The cellIDFromFaceIJ method is described in detail, but I won’t go into it here. CellIDFromFaceIJ method and faceIJOrientation method are inverse of each other. CellIDFromFaceIJ is passing in face, I, j as input, and it returns CellID, faceIJOrientation is breaking CellID into face, I, j, orientation. Orientation is one more orientation than cellidfrom orientation.
The important thing to explain here is how the orientation was calculated.
We know that CellID’s data structure is 3-bit face + 60-bit position + 1-bit flag bit. For level-n non-leaf nodes, there must be 2 n bits of face, followed by 2(maxlevel-n) + 1 bits starting with 1 and ending with 0 bits. MaxLevel = 30.
For example, level-16 must have 32 bits in the middle, followed by 2*(30-16) + 1 = 29 bits. The 29 bits consist of a 1 at the beginning and a 0 at the end. 3 + 32 + 29 = 64 bits. 64-bit CellID is thus formed.
When n = 30,3 + 60 + 1 = 64, so the 1 at the end doesn’t really do anything. When n is 29,3 plus 58 plus 3 is 64, then the end must be 100. Ten doesn’t do anything for the direction, and the extra zero at the end doesn’t do anything for the direction. The key is how many 00’s there are between 10 and 0. When n = 28,3 + 56 + 5 = 64, the last five digits are 10000, with a “00” between 10 and 0. “00” will affect the direction, the initial direction should be xOR 01 to get.
It’s actually easier to understand how “00” affects the original direction. CellID makes four points from the first direction, and each quarter will bring a change of direction. The direction does not change until the last four cells are changed, because between the four cells you can uniquely determine which Cell is selected. As can be seen above, the orientation of Level-30 and Level-29 remains the same, and the other levels need to change or 01 again to obtain the original orientation.
Finally, the conversion, the specific code implementation is as follows:
func cellIDFromFaceIJWrap(f, i, j int) CellID {
/ / 1.
i = clamp(i, - 1, maxSize)
j = clamp(j, - 1, maxSize)
/ / 2.
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)))
/ / 3.
f, u, v = xyzToFaceUV(faceUVToXYZ(f, u, v))
return cellIDFromFaceIJ(f, stToIJ(0.5*(u+1)), stToIJ(0.5*(v+1)))}Copy the codeThe conversion process is divided into three steps. The first step is to deal with the I and j boundary. Step two, convert I and j into u and v. Step 3, u, v go to XYZ, then go back to u, v, and finally go back to CellID.
The first step:
func clamp(x, min, max int) int {
if x < min {
return min
}
if x > max {
return max
}
return x
}Copy the codeClamp functions are used to limit the range of I and j. The range of I and j is always limited to [-1, maxSize].
The second step:
The simplest idea is to convert the (I, j) coordinates to (x, y, z) (the point is not on the boundary) and then call xyzToFaceUV to project it onto the corresponding face.
We know that when we generate CellID, when we generate stToUV, we use a quadratic transformation:
func stToUV(s float64) float64 {
if s >= 0.5 {
return (1 / 3.) * (4*s*s - 1)}return (1 / 3.) * (1 - 4* (1-s)*(1-s))
}Copy the codeBut here, we’re using a simpler transformation, a linear transformation.
u = 2 * s - 1
v = 2 * t - 1Copy the codeThe values of u and v are limited to [-1, 1]. Specific code implementation:
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)))Copy the codeStep 3: find the leaf node and convert u and V into CellID corresponding to the Level.
f, u, v = xyzToFaceUV(faceUVToXYZ(f, u, v))
return cellIDFromFaceIJ(f, stToIJ(0.5*(u+1)), stToIJ(0.5*(v+1)))Copy the codeSo you get a CellID.
Since an edge has four edges, there are four edge neighbors.
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),
}Copy the codeThe above array will load the lower neighbor, right neighbor, upper neighbor, and left neighbor of the current CellID.
If it were shown on a map, it would look like this.
CellID = 3958610196388904960, Level 10 The edge neighbors obtained according to the above method are:
3958603599319138304 // The lower neighbor
3958607997365649408 // Right neighbor
3958612395412160512 // Upper neighbor
3958599201272627200 // The left neighborCopy the codeTo show on a map:
Two. Common vertex neighbor
The common vertex neighbors here are a little different from the vertex neighbors we talked about at the beginning of this article. And there will be some seemingly strange examples below, but also in the actual coding of the author stepped on the pit, share.
A special case is that the Cell is located on the eight vertices of the earth’s tangent cube. So this point has three vertex neighbors instead of four. Because only three faces are connected to each of these eight vertices, there is only one Cell on each face that is their vertex neighbor. Every Cell except these eight points has four vertex neighbors!
j
|
| (0.1) (1.1)
| (0.0) (1.0)
|
---------------> iCopy the codeIn the above mentioned coordinate axes, if the I axis is 1, it falls in the right column of the four quadrants. If the j direction is zero, it’s going to be on the top row of the four quadrants.
If the Cell Level is not equal to 30, that is, the last flag bit 1 is followed by 0, then the end of I and j will be 1 after the Cell is converted into I and j.
The above conclusion can be proved, because when faceIJOrientation function splits the Cell, if all zeros are encountered, such as orientation = 11, and the Cell ends with zeros, then remove the last 8 bits and add orientation, 00000000 11, after lookupIJ transformation, 1111111111, so I = 1111, j = 1111, direction 11. The 00 at the end of the Cell continues the same process, so the I and j end in 1111.
Therefore, we only need to determine which quadrant of the Level given by the input parameter according to I and j, and then we can find all the neighbors of the common vertex.
If the input parameter is given a small Level, that is, a large area of the Cell, then it is necessary to determine which of the four vertices of the input parameter Cell is the common vertex of the current Cell (function caller). A Cell is a rectangle with four vertices. If the current Cell (function caller) is near any vertex, that vertex is selected as the common vertex. Then compute the four cells around the common vertex in turn.
Assuming the given Level of the input parameter is large, that is, the area of the Cell is small, it is also necessary to determine which of the four vertices of the current Cell (function caller) is the co-vertex of the input parameter Cell. A Cell is a rectangle with four vertices. If the Cell is close to the vertex, select that vertex as the common vertex. Then compute the four cells around the common vertex in turn.
Because we need to determine which one is located in the quartile of a Cell, we need to determine the position of its four children. That is to judge the relative position of level-1 children.
halfSize := sizeIJ(level + 1)
size := halfSize << 1
f, i, j, _ := ci.faceIJOrientation()
var isame, jsame bool
var ioffset, joffset intCopy the codeHalfSize is the size of the Cell’s child.
ifi&halfSize ! =0 {
// the offset is in the last column, so we need to add a grid
ioffset = size
isame = (i + size) < maxSize
} else {
// is in the left column, so the offset is subtracted by one cell
ioffset = -size
isame = (i - size) >= 0
}Copy the codeHere we judge which vertex is close to the four vertices of the rectangle according to whether the bit of halfSize is 1. Another thing to note here is that if I + size cannot exceed maxSize, if it does, it is not on the same face. In the same way, i-size cannot be less than 0. Pages less than 0 are not on the same face.
The judgment principle of j axis is exactly the same as that of I axis.
ifj&halfSize ! =0 {
// is on the top line, so the offset needs a grid
joffset = size
jsame = (j + size) < maxSize
} else {
// Is on the bottom row, so the offset is subtracted by one cell
joffset = -size
jsame = (j - size) >= 0
}Copy the codeIn the end, the input Cell is calculated first, and then the cells on the two axes around it are calculated.
results := []CellID{
ci.Parent(level),
cellIDFromFaceIJSame(f, i+ioffset, j, isame).Parent(level),
cellIDFromFaceIJSame(f, i, j+joffset, jsame).Parent(level),
}Copy the codeIf I and j are on the same face, then the common vertices are definitely not on the eight vertices of the tangent cube. Then we can compute the fourth Cell with common vertices.
if isame || jsame {
results = append(results, cellIDFromFaceIJSame(f, i+ioffset, j+joffset, isame && jsame).Parent(level))
}Copy the codeTo sum up, the complete code for calculating common vertex neighbors is as follows:
func (ci CellID) VertexNeighbors(level int) []CellID {
halfSize := sizeIJ(level + 1)
size := halfSize << 1
f, i, j, _ := ci.faceIJOrientation()
fmt.Printf("Halfsize original value = %v-%b\n", halfSize, halfSize)
var isame, jsame bool
var ioffset, joffset int
ifi&halfSize ! =0 {
// the offset is in the last column, so we need to add a grid
ioffset = size
isame = (i + size) < maxSize
} else {
// is in the left column, so the offset is subtracted by one cell
ioffset = -size
isame = (i - size) >= 0
}
ifj&halfSize ! =0 {
// is on the top line, so the offset needs a grid
joffset = size
jsame = (j + size) < maxSize
} else {
// Is on the bottom row, so the offset is subtracted by one cell
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
}Copy the codeHere are a few examples.
The first example is the same size Cell. The input and caller cells are the same Level-10.
VertexNeighbors := cellID.Parent(10).VertexNeighbors(10)
/ / 11011011101111110011110000000000000000000000000000000000000000
3958610196388904960 / / the top right corner
3958599201272627200 / / the top left corner
3958603599319138304 / / the bottom right hand corner
3958601400295882752 / / the bottom left cornerCopy the codeThe second example is not the size Cell. The caller Cell is level-30 by default.
VertexNeighbors := cellID.VertexNeighbors(10)
/ / 11011011101111110011110000000000000000000000000000000000000000
3958610196388904960 / / the bottom right hand corner
3958599201272627200 / / the bottom left corner
3958612395412160512 / / the top right corner
3958623390528438272 / / the top left cornerCopy the codeVertexNeighbors(10) returns level-10 cells, but with different directions and locations. The nature is different in their common vertices, so the generated four cells are generated in different directions.
In the C++ version, there is a restriction on finding vertex neighbors:
DCHECK_LT(level, this->level());Copy the codeThe Level of the entry must be strictly lower than the Level of the Cell to be found. In other words, the Cell size of the input Cell must be smaller than the Cell size. This is not required in the version of Go, however, and can be large or small.
In the following example, the input parameter is smaller than the Cell Level. (You can see that Chengdu has shrunk into a dot.)
VertexNeighbors := cellID.Parent(10).VertexNeighbors(5)
3957538172551823360 / / the bottom right hand corner
3955286372738138112 / / the bottom left corner
3959789972365508608 / / the top right corner
3962041772179193856 / / the top left cornerCopy the codeIn the following example, the input parameter is larger than the Cell Level. (You can see that Level 15 is already very small.)
VertexNeighbors := cellID.Parent(10).VertexNeighbors(15)
3958610197462646784 / / the bottom left corner
3958610195315163136 / / the bottom right hand corner
3958610929754570752 / / the top left corner
3958609463023239168 / / the top right cornerCopy the codeThree. The whole neighborhood
Finally, back to the question I asked at the beginning of the article. How do you find the neighbor of the Hilbert curve in a quadtree? After some foreshadows, the reader may already know what to do.
The area of the Level of the input parameter must be smaller or equal to that of the caller Cell. That is, the Level value of the input parameter cannot be less than the Level of the caller. Once the size of the neighbor Cell is small, the area of the neighbor Cell becomes huge. It is likely that a neighbor Cell will be filled with all the neighbors of the original Cell, so such search is meaningless.
For example, if the Level of an incoming parameter is smaller than the Level of the caller’s Cell. Then when searching for its whole neighbor, the following situation will be found:
AllNeighbors := cellID.Parent(10).AllNeighbors(5)Copy the code
In this case, all neighbors can be found, but overlapping cells may occur. Why such a phenomenon occurs will be analyzed below.
If the input parameter and the caller’s Cell are of the same Level, then the found full neighbor is the problem described at the beginning of this article. The ideal state is as follows:
The concrete implementation is as follows:
func (ci CellID) AllNeighbors(level int) []CellID {
var neighbors []CellID
face, i, j, _ := ci.faceIJOrientation()
// Find the coordinates of the lower-left leaf node. We need to gauge the coordinates of I and j. It is possible that the Level of the input parameter is larger than that of the caller Cell.
size := sizeIJ(ci.Level())
i &= -size
j &= -size
nbrSize := sizeIJ(level)
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
// Upper neighbor and lower neighbor
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 neighbor, right neighbor, and 2 diagonal vertex 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 the size of the loop is larger than the size of the loop, it is not necessary to check the size of the loop
if k >= size {
break}}return neighbors
}Copy the codeThe simple idea of the above code is written in comments. I’ll explain what needs to be explained now.
The first thing you need to understand is the relationship between nbrSize and size. Why nbrSize? The input Level can be different from the Level of the caller Cell. The Level of the input parameter can represent large, small, or equal cells. The final result is expressed in terms of the size of the nbrSize cell, so the loop needs to use nbrSize to control the size of the cell. Size is just the size of the original caller’s Cell.
The change in K in the loop, when K is equal to -nbrsize, the loop is only counting the left and right neighbors. The vertex neighbor on the diagonal is also a special case of the left neighbor and the right neighbor. And then K = 0, we’re going to start counting the neighbors up here and the neighbors down here. K keeps increasing until K >= size at last. In the last loop, the left and right neighbors will be calculated first and then break out.
The caller’s Cell is in the middle, so if you want to skip this Cell to the other side (up, down, or left, or right), you need to skip size. The specific code implementation is I + size and J + size.
First, check the cyclic scanning mode of the left and right neighbors.
The left neighbor is i-nbrsize, j + k, k is going around. This is how the left neighbor is generated. It generates a left neighbor column. It starts at the bottom left and goes all the way up to the top left.
The right neighbor is I + size, j + k, and k is going around. This is how the right neighbor is generated. It generates the right neighbor column. It starts at the bottom right and goes all the way up to the top right.
Now look at the cyclic scanning of the upper and lower neighbors.
The lower neighbor is I + k, j-nbrsize, k is circulating. This is how the lower neighbor is generated. It generates the next neighbor row. It starts at the left of the lower neighbor and goes all the way up to the right of the lower neighbor.
The upper neighbor is I + k, j + size, and k is circulating. This is how the upper neighbor is generated. It generates the upper neighbor row. It starts at the left of the upper neighbor and goes up to the right of the upper neighbor.
For example:
All neighbors of the middle Cell are cells of the same Level as Cell 8 in the figure above.
The generation order is identified by the need. 1,2,5,6,7,8 are generated by the left and right neighbors. 3 and 4 are generated by the upper and lower neighbors.
The above example is all generated by Level-10 cells. The total number of neighbors is exactly eight.
AllNeighbors := cellID.Parent(10).AllNeighbors(10)
3958601400295882752.3958605798342393856.3958603599319138304.3958612395412160512.3958599201272627200.3958607997365649408.3958623390528438272.3958614594435416064Copy the codeTake another example where the Level is greater than the Level of the caller’s Cell.
AllNeighbors := cellID.Parent(10).AllNeighbors(11)
3958600575662161920.3958606622976114688.3958603324441231360.3958611570778439680.3958600025906348032.3958607172731928576.3958603874197045248.3958613220045881344.3958599476150534144.3958608821999370240.3958623115650531328.3958613769801695232Copy the codeIts full neighbor generation sequence is as follows:
1,2,5,6,9,10,11,12 are left and right neighbors, and 3,4,7,8 are up and down neighbors. We can see that the left and right neighbors are created from the bottom up. The upper and lower neighbors are generated from left to right.
If the Level is larger, such as Level-15, more neighbors will be generated:
Now let’s explain what happens if the Level of the input parameter is smaller than the Level of the caller Cell.
For example, the Level parameter is 9.
AllNeighbors := cellID.Parent(10).AllNeighbors(9)
3958589305667977216.3958580509574955008.3958580509574955008.3958615693947043840.3958598101760999424.3958606897854021632.3958624490040066048.3958615693947043840Copy the codeThe generated full neighbor is as follows:
You can see there were eight neighbors, now there are only six. It’s still 8, but 2 of them are duplicated. See the two dark red cells in the figure.
Why do they overlap?
The level-10 Cell of the intermediate caller is drawn first.
Because it is Level-9, it is a quarter of the Cell in the middle.
Let’s draw the two upper-level neighbors of Level-10.
You can see that the upper neighbor Up and the vertex neighbor up-right are in the same Level-9 Cell. Therefore, the upper neighbor and the upper – right vertex neighbor are the same Level-9 Cell. So it overlaps. Similarly, the lower neighbor overlaps with the lower right vertex neighbor. So there will be overlap of 2 cells.
And there is no space for the caller Cell in the middle. Because after I + size, the scope is still in the same Level-9 Cell.
If the Level is smaller, the overlap will change again. For example, Level-5.
AllNeighbors := cellID.Parent(10).AllNeighbors(5)
3953034572924452864.3946279173483397120.3946279173483397120.3957538172551823360.3955286372738138112.3957538172551823360.3962041772179193856.3959789972365508608Copy the codeIt’s drawn on a map
The overlap has also changed.
At this point, search neighbor – related algorithms are introduced.
Spatial Search series:
How to understand n-dimensional space and N-dimensional space-time efficient multi-dimensional space point index algorithms – Geohash and Google S2 how to generate CellID in Google S2? How to find the neighbors of the Hilbert curve on the quadtree of the magical Debruine sequence?
Making Repo: Halfrost – Field
Follow: halfrost dead simple
Source: halfrost.com/go_s2_Hilbe…