Google 的S2算法原理以及使⽤Java 版本--部分参考⾃《⾼效的多维空间点索引算法》
⽂章⽬录
相关资料
1.
2.
3.
4.
5.
6.
7.
8. S2jar：使⽤Maven⼯具
1.S2算法是什么？
S2其实是来⾃⼏何数学中的⼀个数学符号 S ²，它表⽰的是单位球。S2 这个库其实是被设计⽤来解决球⾯上各种⼏何问题的。
2.为什么要使⽤S2算法？
S2 来解决多维空间点的问题的
1. s2有30级，geohash只有12级。s2的层级变化较平缓，⽅便选择。
2. s2功能强⼤，解决了向量计算，⾯积计算，多边形覆盖，距离计算等问题，减少了开发⼯作量。 3. s2解决了多边形覆盖问题。
这是其与geohash功能上最本质的不同。给定不规则范围，s2可以计算出⼀个多边形近似覆盖这个范围。其覆盖⽤的格⼦数量根据精确度可控。geohash在这⽅⾯⼗分不友好，划定⼀个⼤⼀点的区域，其格⼦数可能达到数千，若减少格⼦数则丢失精度，查询区域过⼤。
如下，在min level和max level不变的情况下，只需设置可接受的max cells数值，即可控制覆盖精度。
⽽且其cell的region⼤⼩⾃动适配。geohash要在如此⼤范围实现⾼精度覆盖则会产⽣极为庞⼤的⽹格数。
另外需要注意的是，在minLevel,maxLevel,maxCells这3个参数中，不⼀定能完全满⾜.⼀般⽽⾔是优先满⾜maxLevel即最精细的⽹格⼤⼩，再尽可能控制cell数量在maxCells⾥⾯.⽽minLevel由于会合并⽹格，所以很难满⾜(在查询⼤区域的时候可能会出现⼀个⼤⽹格和很多个⼩⽹格，导致⽊桶效应.这个时候可能将⼤⽹格划分为指定等级的⼩⽹格，即最终效果为，严格遵循minLevel和maxLevel，为此牺牲maxCells)
eg：9个格⼦，maxcells为：10
eg：34个格⼦，maxcells为：45
3.S2的原理是什么？
1）球⾯坐标变换
球⾯上的⼀个点，在直⾓坐标系中，可以这样表⽰： <!--google 的S2包--> <dependency > <groupId >io .sgr </groupId > <artifactId >s2-geometry -library -java </artifactId > <version >1.0.0</version > </dependency > <!--附带的google common 组件包--> <dependency > <groupId >com .google .guava </groupId > <artifactId >guava </artifactId > <version >21.0</version > </dependency >
1
2
3
4
5
6
7脊髓损伤论坛
8
9
10
11
12
13//经纬度转弧度计算球⾯坐标S2LatLng ll = S2LatLng .fromDegrees (36.683, 117.1412);//弧度转⾓度乘以 π / 180°//⾓度转弧度乘以 180°/πS2Point point = ll .toPoint ();
1
2
3
4
5
6
7
再进⼀步，我们可以和球⾯上的经纬度联系起来。不过这⾥需要注意的是，纬度的⾓度 α 和直⾓坐标系下的球⾯坐标 θ 加起来等于 90°。所以三⾓函数要注意转换。 于是地球上任意的⼀个经纬度的点，就可以转换成 f(x,y,z)。 在 S2 中，地球半径被当成单位 1 了。所以半径不⽤考虑。x，y，z的值域都被限定在了[-1,1] x [-1,1] x [-1,1]这个区间之内了。
2）球⾯坐标转平⾯坐标（）
⾸先在地球外⾯套了⼀个外切的正⽅体，如下图:
[外链图⽚转存失败，源站可能有防盗链机制，建议将图⽚保存下来直接上传(img-FiBcpuxp-1620986757996)(C:\Users\余祥超\Desktop\notes\googleS2Img\球⾯外切正⽅形.png)]
从球⼼向外切正⽅体6个⾯分别投影。S2 是把球⾯上所有的点都投影到外切正⽅体的6个⾯上。 这⾥简单的画了⼀个投影图，上图左边的是投影到正⽅体⼀个⾯的⽰意图，实际上影响到的球⾯是右边那张图。
从侧⾯看，其中⼀个球⾯投影到正⽅体其中⼀个⾯上，边缘与圆⼼的连线相互之间的夹⾓为90°，但是和x，y，z轴的⾓度是45°。我们可以在球的6个⽅向上，把45°的辅
助圆画出来，见下图左边。
上图左边的图画了6个辅助线，蓝线是前后⼀对，红线是左右⼀对，绿线是上下⼀对。分别都是45°的地⽅和圆⼼连线与球⾯相交的点的轨迹。这样我们就可以把投影到外切正⽅体6个⾯上的球⾯画出来，见上图右边。
投影到正⽅体以后，我们就可以把这个正⽅体展开了。
在 Google S2 中，它是把地球展开成如下的样⼦：
这样第⼀步的球⾯坐标进⼀步的被转换成 f(x,y,z) -> g(face,u,v)，face是正⽅形的六个⾯，u，v对应的是六个⾯中的⼀个⾯上的x，y坐标。
remark:
1.默认定义 x 轴为0，y轴为1，z轴为2 。
2.最后 face 的值就是三个轴⾥⾯最长的轴，注意这⾥限定了他们三者都在 [0,5] 之间，所以如果是负轴就需要 + 3 进⾏修正。
3.主轴为 x 正半轴，face = 0；主轴为 y 正半轴，face = 1；主轴为 z 正半轴，face = 2；主轴为 x 负半轴，face = 3；主轴为 y 负半轴，face = 4；主轴为 z 负半轴，face = 5 。
4.如果直观的对应⼀个外切⽴⽅体的哪6个⾯，那就是 face = 0 对应的是前⾯，face = 1 对应的是右⾯，face = 2 对应的是上⾯，face = 3 对应的是后⾯，face = 4 对应的是左⾯，face = 5 对应的是下⾯。
上⼀步我们把球⾯上的球⾯矩形投影到正⽅形的某个⾯上，形成的形状类似于矩形，但是由于球⾯上⾓度的不同，最终会导致即使是投影到同⼀个⾯上，每个矩形的⾯积也不⼤
相同。
上图就表⽰出了球⾯上个⼀个球⾯矩形投影到正⽅形⼀个⾯上的情况。
经过实际计算发现，最⼤的⾯积和最⼩的⾯积相差5.2倍。见上图左边。相同的弧度区间，在不同的纬度上投影到正⽅形上的⾯积不同。
现在就需要修正各个投影出来形状的⾯积。如何选取合适的映射修正函数就成了关键。⽬标是能达到上图右边的样⼦，让各个矩形的⾯积尽量相同。
线性变换是最快的变换，但是变换⽐最⼩。tan() 变换可以使每个投影以后的矩形的⾯积更加⼀致，最⼤和最⼩的矩形⽐例仅仅只差0.414。可以说⾮常接近了。但是 tan() 函数的调⽤时间⾮常长。如果把所有点都按照这种⽅式计算的话，性能将会降低3倍。
最后⾕歌选择的是⼆次变换，这是⼀个近似切线的投影曲线。它的计算速度远远快于 tan() ，⼤概是 tan() 计算的3倍速度。⽣成的投影以后的矩形⼤⼩也类似。不过最⼤的矩形和最⼩的矩形相⽐依旧有2.082的⽐率。
上表中，ToPoint 和 FromPoint 分别是把单位向量转换到 Cell ID 所需要的毫秒数、把 Cell ID 转换回单位向量所需的毫秒数（Cell ID 就是投影到正⽅体六个⾯，某个⾯上矩形的 ID，矩形称为 Cell，它对应的 ID 称为 Cell ID）。ToPointRaw 是某种⽬的下，把 Cell ID 转换为⾮单位向量所需的毫秒数。
S2的矩形⾯积修正源码x = r * sin θ * cos φy = r * sin θ * sin φ z = r * cos θ
1
2
3//投影int face = S2Projections .xyzToFace (point );R2Vector vector = S2Projections .validFaceXyzToUv (face , point );
1
23//face 是不变的，修正double s = S2Projections .uvToST (vector .x ());double t = S2Projections .uvToST (vector .y ());
1
2
3
经过修正变换以后，u，v都变换成了s，t。值域也发⽣了变化。u，v的值域是[-1,1]，变换以后，使s，t的值域是[0,1]。
⾄此，⼩结⼀下，球⾯上的点S(lat,lng) -> f(x,y,z) -> g(face,u,v) -> h(face,s,t)。⽬前总共转换了4步，球⾯经纬度坐标转换成球⾯xyz坐标，再转换成外切正⽅体投影⾯上的坐标，最后变换成修正后的坐标。
4）点与坐标轴点相互转换
在 S2 算法中，默认划分 Cell 的等级是30，也就是说把⼀个正⽅形划分为 2^30 * 2^30个⼩的正⽅形。
那么上⼀步的s，t映射到这个正⽅形上⾯来，对应该如何转换呢？
s，t的值域是[0,1]，现在值域要扩⼤到[0,2-1]。
S2源码
到这⼀步，是h(face,s,t) -> H(face,i,j)。
5）坐标轴点与CellId 相互转换 //S2_LINEAR_PROJECTION 线性变换，未实现 //sttoUV return 2 * s - 1 //uvToST return 0.5 * (u + 1); //S2_TAN_PROJECTION tan()变换 //S2_QUADRATIC_PROJECTION ⼆次变换 public static strictfp double stToUV (double s ) { switch (S2_PROJECTION ) { case S2_LINEAR_PROJECTION : return s ; case S2_TAN_PROJECTION : s = Math .tan (0.7853981633974483D * s ); return s + 1.1102230246251565E-16D * s ; case S2_QUADRATIC_PROJECTION : if (s >= 0.0D ) { return 0.3333333333333333D * ((1.0D + s ) * (1.0D + s ) - 1.0D ); } return
0.3333333333333333D * (1.0D - (1.0D - s ) * (1.0D - s )); default : throw new IllegalStateException ("Invalid value for S2_PROJECTION"); } } public static strictfp double uvToST (double u ) { switch (S2_PROJECTION ) { case S2_LINEAR_PROJECTION : return u ; case S2_TAN_PROJECTION : return 1.2732395447351628D * Math .atan (u ); case S2_QUADRATIC_PROJECTION : if (u >= 0.0D ) { return Math .sqrt (1.0D + 3.0D * u ) - 1.0D ; } return 1.0D - Math .sqrt (1.0D - 3.0D * u ); default : throw new IllegalStateException ("Invalid value for S2_PROJECTION"); } }
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31ai ei ui 教学设计
32
33
34
35
36
37
38
39int i = S2CellId .stToIJ (s );int j = S2CellId .stToIJ (t );
1
230private static strictfp int stToIJ (double s ) { int m = 536870912; return (int )Math .max (0L , Math .min (1073741823L , Math .round (5.36870912E8D * s + 5.368709115E8D ))); }
1
2
3
4S2CellId cellid = S2CellId .fromFaceIJ (face , i , j );S2LatLng lan = cellid .toLatLng ();
1哈卡斯人
2
3
4
1.⽣成曲线
1.解释变量。
posToIJ 代表的是⼀个矩阵，⾥⾯记录了⼀些单元希尔伯特曲线的位置信息。
注意在下⾯⼀阶曲线中：横着的是i，竖着的是j
把 posToIJ 数组⾥⾯的信息⽤图表⽰出来，如下图：
从上⾯这四张图我们可以看出：
posToIJ 的四张图其实是“ U ” 字形逆时针分别旋转90°得到的。这⾥我们只能看出四张图相互之间的联系，即兄弟之间的联系，但是看不到⽗⼦图相互之间的联系。
同理，把 ijToPos 数组⾥⾯的信息⽤图表⽰出来，如下图： static int lookupBits = 4; static int swapMask = 0x01; static int invertMask = 0x02; //在整个库中，没有使⽤ static int ijToPos [][] ={ {0, 1, 3, 2}, // 标准顺序 {0, 3, 1, 2}, // 轴旋转 {2, 3, 1, 0}, // 上下
倒置 {2, 1, 3, 0}, // 轴旋转左右倒置 }; //数组存的值是ij 组合的⼆进制的数值，表⽰⽅向 static int posToIJ [][] = { {0, 1, 3, 2}, // 标准顺序: (0,0), (0,1), (1,1), (1,0) {0, 2, 3, 1}, // 轴旋转: (0,0), (1,0), (1,1), (0,1) {3, 2, 0, 1}, // 上下倒置: (1,1), (1,0), (0,0), (0,1) {3, 1, 0, 2}, // 轴旋转左右倒置: (1,1), (0,1), (0,0), (1,0) }; //posToOrientation 数组⾥⾯装了4个数字，分别是1,0,0,3，表⽰⼦格⼦的⽅向 static int posToOrientation [] = {swapMask , 0, 0, invertMask | swapMask };//1 | 2=3 //希尔伯特曲线 ID 转换成坐标轴 IJ 的转换表 static int [] lookupIJ = new int [1 << (2 * lookupBits + 2)];//1<<10=2^10=1024 //坐标轴 IJ 转换成希尔伯特曲线 ID 的转换表 static int [] lookupPos = new int [1 << (2 * lookupBits + 2)];
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
这⾥是初始化的递归函数，在希尔伯特曲线的标准顺序中可以看到是有4个格⼦，并且格⼦都有顺序的，所以初始化要遍历满所有顺序。⼊参的第4个参数origOrientation，就是从0 - 3 。
下⾯这个函数是⽣成希尔伯特曲线的。我们可以看到有⼀处对 pos << 2的操作，这⾥是把位置变换到第⼀个4个⼩格⼦中，所以位置乘以了4。
2.r[0] >> 1和r[0] & 1
r 数组来⾃于 posToIJ 数组。posToIJ 数组上⾯说过了，它⾥⾯装的其实是4个不同⽅向的“ U ”字。相当于表⽰了当前四个⼩⽅格兄弟相互之间的⽅向。r[0]、r[1]、r[2]、r[3] 取出的其实就是 00，01，10，11 这4个数。
那么 r[0]>>1 操作就是取出⼆位⼆进制位的前⼀位，即 i 位。r[0]&1 操作就是取出⼆位⼆进制位的后⼀位，即 j 位。r[1]、r[2]、r[3] 同理。
其实这个对应的就是 图0 中4个⼩⽅块接下来再划分的⽅向。
图0 中0号的位置下⼀个图的⽅向应该是图1，即01；
图0 中1号的位置下⼀个图的⽅向应该是图0，即00；
宜昌市康庄路小学
图0 中2号的位置下⼀个图的⽅向应该是图0，即00；
图0 中3号的位置下⼀个图的⽅向应该是图3，即11 。
这就是初始化 posToOrientation 数组⾥⾯的⽞机了。
posToIJ 的四张图我们只能看出兄弟之间的关系，那么 posToOrientation 的四张图让我们知道了⽗⼦之间的关系。
eg:public static void 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 );}
1
2
3
4
5
6public static void initLookupCell (int level , int i , int j , int origOrientation , int pos , int orientation ) { if (level == lookupBits ) { int 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; int 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]); }
1
2
3
4
5
6
7
8
9
10
正定方言11
12报效祖国为国争光的资料
13
14
15
16
17
18
19
20
21
22posToOrientation [] = {swapMask , 0, 0, invertMask | swapMask };//数值代⼊：posToOrientation = [4]int {1, 0, 0, 3}//根据上⼀次的原始的⽅向推算出当前的 pos 所在的⽅向。即计算⽗⼦之间关系。orientation ^posToOrientation [0]//举个例⼦，假设 orientation = 0，即图0，那么：00 ^ 01 = 0100 ^ 00 = 0000 ^ 00 = 0000 ^ 11 = 11//图0 的四个孩⼦的⽅向就被我们算出来了，01，00，00，11，1003 。和上⾯图⽚中图0展⽰的是⼀致的。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议