Hilbert编码思路和scala 代码

需求：

使用Hilbert 曲线对遥感影像瓦片数据进行编码，获取某个区域的编码值即可

Hilbert 曲线编码方式

思路

大致可以对四个方向的数据进行归类

左下
左上
右上
右下

这个也对应着编码的顺序

思考在不同Hilbert深度（阶）情况下的四个区域的取值范围

左下

$0,4^{n-1})$
左上

$[4^{n-1}, \ 2\ \times \ 4^{n-1})$
右上

$[2\ \times\ 4^{n-1},3 \times4^{n-1})$
右下

$\ \times \ 4^{n-1}, 4^n)$

因此可以使用递归的方式来进行生成数据

主要问题

左上角和右上角不需要做旋转操作

 // 左上角，待判定点不需要做任何调整，直接递归新子区域即可case (true, false)  =>recursiveHilbertEncode(point, (Point(lat_half, lon_half), Point(bounds._1.lat,bounds._2.lon)), maxLevel - 1 ) +   Math.pow(4, maxLevel-2).toInt// 右上角case (false, false) =>recursiveHilbertEncode(point, (Point(lat_half, lon_half), bounds._2), maxLevel - 1) + 2 * Math.pow(4, maxLevel - 2).toInt

左下角和左上角需要考虑旋转问题

可以使用举例子的方式，来找到旋转后的关系

假设空间区域的左下角坐标为（10，20），右上角的坐标为（20，30），中心坐标为(15, 25)

同样以二阶图为例

这里点位置变更主要针对的是待判定的位置属于哪个区域！

Hilbert 编码区域左下区域(左侧为原始点位置，右侧为旋转后的点位置)

（18，22.5）=> （12.5，28）

（17.5，22） => (12，27.5)

Hilbert 编码区域右下区域

（12，27.5）=>（12.5, 28）

（12.5, 22） => （18, 27.5）

左下旋转后坐标生成方式

原始点位置-中心点位置 = （a,b）

结果 = 中心点位置 + （b,a)

例：（18，22.5） - （15，25） = （3，-2.5）

结果 = (15,25 ) + (-2.5, 3) = (12.5, 28)

右下旋转后坐标生成方式

原始点位置-中心点位置 = （a,b）

结果 = 中心点位置 -（b,a)

例：（12， 27.5） - （15，25） = （-3，2.5）

结果 = （15，25） -（2.5，-3） = （12.5，28）

思路捋清之后，可以考虑代码

// 左下角，case (true, true)   =>recursiveHilbertEncode(Point(lat_half + lon_reverse,lon_half + lat_reverse),   (bounds._1, Point(lat_half, lon_half)),maxLevel - 1)// 右下角case (false, true)  =>recursiveHilbertEncode(Point(lat_half - lon_reverse,lon_half - lat_reverse), ((Point(lat_half,bounds._1.lon),Point(bounds._2.lat,lon_half))), maxLevel - 1) + 3 * Math.pow(4, maxLevel-2).toInt

需要对待判定点重新生成，但区域的经纬坐标可以不变，我们不关系点坐标怎么变换，我们只是需要最终的编码结果，保证递归正常运行即可

Hilbert 源码生成

input：

正方形边界两个对角顶点的坐标
最大Hilbert 层级（阶、深度）
待判定的区域两点坐标（或者是一个点）

output：

该区域（区域也是计算区域的中心点位置）、点所在的Hilbert 编码值

递归函数

recursiveHilbertEncode

  def recursiveHilbertEncode(point: Point, bounds:(Point,Point), maxLevel:Int): Int = {if (maxLevel == 1)return 0// 边界中心点val lat_half = (bounds._1.lat + bounds._2.lat) / 2 //xval lon_half = (bounds._1.lon + bounds._2.lon) / 2 //y// 反转差值val lat_reverse = point.lat - lat_halfval lon_reverse = point.lon - lon_half(point.lat < lat_half, point.lon < lon_half) match {// 右上角 , p1,p2 不需要动，需要动的是边界和case (false, false) =>recursiveHilbertEncode(point, (Point(lat_half, lon_half), bounds._2), maxLevel - 1) + 2 * Math.pow(4, maxLevel - 2).toInt// 右下角 向左旋转case (false, true)  =>recursiveHilbertEncode(Point(lat_half - lon_reverse,lon_half - lat_reverse), ((Point(lat_half,bounds._1.lon),Point(bounds._2.lat,lon_half))), maxLevel - 1) + 3 * Math.pow(4, maxLevel-2).toInt//左上case (true, false)  =>recursiveHilbertEncode(point, (Point(lat_half, lon_half), Point(bounds._1.lat,bounds._2.lon)), maxLevel - 1 ) +   Math.pow(4, maxLevel-2).toInt// 左下角 向右旋转case (true, true)   =>recursiveHilbertEncode(Point(lat_half + lon_reverse,lon_half + lat_reverse),   (bounds._1, Point(lat_half, lon_half)),maxLevel - 1)}}

完整代码

测试代码没整理，有两个main函数，一个是零散的编码区域，一个是从0-15区域的编码结果

import HilbertCurve.getTileCodeobject HilbertCurve {case class Point(lat: Double, lon: Double)// 递归查看Hilbert编码/*** 瓦片中心点坐标* 新bounds边界（其实也是用来定义边界的中心点的）* maxLevel* @return*/def recursiveHilbertEncode(point: Point, bounds:(Point,Point), maxLevel:Int): Int = {if (maxLevel == 1)return 0// 边界中心点val lat_half = (bounds._1.lat + bounds._2.lat) / 2 //xval lon_half = (bounds._1.lon + bounds._2.lon) / 2 //y
//
//    println("边界中心点  :(" + lat_half + ", " + lon_half + ")")
//
//    println("瓦片中心点：(" + point.lat + ", " + point.lon + ")")// point 为瓦片的中心点// 反转差值val lat_reverse = point.lat - lat_halfval lon_reverse = point.lon - lon_half(point.lat < lat_half, point.lon < lon_half) match {// 右上角 , p1,p2 不需要动，需要动的是边界和case (false, false) =>recursiveHilbertEncode(point, (Point(lat_half, lon_half), bounds._2), maxLevel - 1) + 2 * Math.pow(4, maxLevel - 2).toInt//TODO: 问题// 右下角 向左旋转// 此处需要也对point 点进行旋转吗？，好像不需要，这里是沿着中轴线旋转180度case (false, true)  =>recursiveHilbertEncode(Point(lat_half - lon_reverse,lon_half - lat_reverse), ((Point(lat_half,bounds._1.lon),Point(bounds._2.lat,lon_half))), maxLevel - 1) + 3 * Math.pow(4, maxLevel-2).toInt//左上case (true, false)  =>recursiveHilbertEncode(point, (Point(lat_half, lon_half), Point(bounds._1.lat,bounds._2.lon)), maxLevel - 1 ) +   Math.pow(4, maxLevel-2).toInt//TODO： 左下角也有问题// 左下角 向右旋转case (true, true)   =>recursiveHilbertEncode(Point(lat_half + lon_reverse,lon_half + lat_reverse),   (bounds._1, Point(lat_half, lon_half)),maxLevel - 1)}}// 第一步函数：获取瓦片编码def getTileCode(p1: Point, p2: Point, bounds: (Point, Point), maxLevel: Int): Int = {// TODO：此处也默认 bounds 第一个点位于左下角// TODO: 此处默认p1 位于p2的左下角// 定义中心点 x，y值val x_half = (p1.lat + p2.lat)/2val y_half = (p1.lon + p2.lon)/2// 边界的中心点val encoder = recursiveHilbertEncode(Point(x_half,y_half), bounds, maxLevel)encoder}/*** 测试左下角的16个区域* @param args*/def main2(args: Array[String]): Unit = {val bounds = (Point(10.0, 10.0), Point(20.0, 20.0))val maxLevel4 = 4// TODO： 左下角//  15 Error//     最左下角的四个，但是此时如果是旋转过，再去查看右上角和左下角，就会有旋转问题
//val p1 = Point(10,10)val p2 = Point(11.25,11.25)val tileCode0 = getTileCode(p1,p2,bounds,maxLevel4)println(s"TIle Code0 is given region :$tileCode0")val p3 = Point(10,11.25)val p4 = Point(11.25,12.5)val tileCode1 = getTileCode(p3,p4,bounds,maxLevel4)println(s"TIle Code1 is given region :$tileCode1")val p5 = Point(11.25,11.25)val p6 = Point(12.5,12.5)val tileCode2 = getTileCode(p5,p6,bounds,maxLevel4)println(s"TIle Code2 is given region :$tileCode2")val p7 = Point(11.25,10)val p8 = Point(12.5,11.25)val tileCode3 = getTileCode(p7, p8, bounds, maxLevel4)println(s"TIle Code3 is given region :$tileCode3")println("==================")val a = Point(12.5,10)val b = Point(13.75,11.25)val tileCode4 = getTileCode(a, b, bounds, maxLevel4)println(s"TIle Code4 is given region :$tileCode4")val t1 = Point(13.75,10)val t2 = Point(15.0, 11.25)val tileCode5 = getTileCode(t1,t2,bounds,maxLevel4)println(s"Tile Code5 is given region:$tileCode5")val c = Point(13.75,11.25)val d = Point(15,12.5)val tileCode6 = getTileCode(c,d,bounds,maxLevel4)println(s"TIle Code6 is given region :$tileCode6")val b3 = Point(12.5,11.25)val b4 = Point(13.75,12.5)val tileCode7 = getTileCode(b3,b4,bounds,maxLevel4)println(s"TIle Code7 is given region :$tileCode7")println("=================")val e = Point(12.5 ,12.5)val f = Point(13.75, 13.75)val tileCode8 = getTileCode(e,f,bounds,maxLevel4)println(s"TIle Code8 is given region :$tileCode8")val g = Point(13.75, 12.5)val h = Point(15,13.75)val tileCode9 = getTileCode(g, h,bounds, maxLevel4)println(s"TIle Code9 is given region :$tileCode9")val c3 = Point(12.5,13.75)val c4 = Point(13.75,15)val tileCode11 = getTileCode(c3,c4, bounds, maxLevel4)println(s"TIle Code11 is given region :$tileCode11")
//     左下角的左上角的右上角有问题println("==========")val c1 = Point(11.25,13.75)val c2 = Point(12.5,15)val tileCode12 = getTileCode(c1,c2,bounds,maxLevel4)println(s"TIle Code12 is given region :$tileCode12")val d1 = Point(11.25,12.5)val d2 = Point(12.5,13.75)val tileCode13 = getTileCode(d1,d2,bounds,maxLevel4)println(s"TIle Code13 is given region :$tileCode13")val d3 = Point(10,12.5)val d4 = Point(11.25,13.75)val tileCode14 = getTileCode(d3,d4,bounds,maxLevel4)println(s"TIle Code14 is given region :$tileCode14")val d5 = Point(10,13.75)val d6 = Point(11.25,15)val tileCode15 = getTileCode(d5,d6,bounds,maxLevel4)println(s"TIle Code15 is given region :$tileCode15")// TODO； 右下角TODO Error// 58
//    val p74 = Point(15,10)
//    val p84 = Point(16.25,11.25)
//    val tileCode44 = getTileCode(p74, p84, bounds, maxLevel4)
//    println(s"Tile Code for given region: $tileCode44")}def main(args: Array[String]): Unit = {val bounds = (Point(10.0, 10.0), Point(20.0, 20.0))// 右上角  42val p14 = Point(18.75, 18.75)val p24 = Point(20, 20)val maxLevel4 = 4val tileCode14 = getTileCode(p14, p24, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode14")// 左上角 //21val p34 = Point(10, 18.75)val p44 = Point(11.25, 20)val tileCode24 = getTileCode(p34, p44, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode24")//   5val p54 = Point(13.75, 10)val p64 = Point(15, 11.25)val tileCode34 = getTileCode(p54, p64, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode34")// TODO 15 Errorval p545 = Point(10, 13.75)val p645 = Point(11.25, 15)val tileCode344 = getTileCode(p545, p645, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode344")// 右下角// TODO Error// 58val p74 = Point(15,10)val p84 = Point(16.25,11.25)val tileCode44 = getTileCode(p74, p84, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode44")// 34val p1 = Point(16.25, 16.25)val p2 = Point(17.5, 17.5)val maxLevel = 3val tileCode = getTileCode(p1, p2, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode")// 10val p3 = Point(13.75, 13.75)val p4 = Point(15.0, 15.0)val tileCode2 = getTileCode(p3, p4, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode2")// 8val p5 = Point(12.5, 12.5)val p6 = Point(13.75, 13.75)val tileCode3 = getTileCode(p5, p6, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode3")// 32val p7 = Point(15.0, 15.0)val p8 = Point(16.25, 16.25)val tileCode4 = getTileCode(p7, p8, bounds, maxLevel4)println(s"Tile Code for given region: $tileCode4")}
}