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")}
}