Hilbert曲線介紹以及代碼實現
阿新 • • 發佈：2017-07-07
gis 算法 空間填充曲線是指，一維曲線去包含整個二維甚至多維空間的一種函數曲線。而根據不同的排列規則，可以得到不同的空間填充曲線。
Hilbert曲線以及其離散近似表示方法都非常實用，因為其將多維空間轉換為一維空間的方法很好地保留了空間鄰近性。（x，y）是一個單元方格中的點，d代表該點在Hilbert曲線上的位置，而由於其空間的鄰近性，在單元格上近似的點，其對應Hilbert的d值也比較接近。
正因為這種鄰近性，空間填充曲線被廣泛用於計算機科學，且在多維數據庫索引中，經常用Hilbert曲線取代Z order曲線。接下來我們看如何用代碼實現Hilbert算法。
現在假設我們以U字形來訪問區域。在每個象限中，我們同樣以U字形來訪問子象限，但是要調整好U字形的朝向使得和相鄰的象限銜接起來。如果我們正確地組織了這些U字形的朝向，我們就能完全消除不連續性，不管我們選擇了什麽分辨率，都能連續地訪問整個區域，可以在完全地探訪了一個區域後才移動到下一個。這個方案不僅消除了不連續性，而且提高了總體的局域性。這就是希爾伯特曲線。
如圖4，在第一層，枚舉這些點很簡單：選定一個方向和一個起始點，環繞四個象限，用0到3給他們編號。當我們要確定訪問子象限的順序同時維護總體的鄰接屬性，困難就來了。通過檢查我們發現，子象限的曲線是原曲線的簡單變換，而且只有四種變換。自然地，這個結論也適用於子子象限，等等。對於一個給定的象限，我們在其中畫出的曲線是由象限所在大的方形的曲線以及該象限的位置決定的。只需要費一點力，我們就能構建出如下概況所有情況的表
假設我們想用這個表來確定某個點在第三層希爾伯特曲線上的位置。在這個例子中，假設點的坐標是(5,2)。 從第一層級開始，則由上圖的第一個方形開始，找到該點所在的象限。在該例中，其在右上方的象限。那麽點在希爾伯特曲線上的位置的第一部分是3（二進制是11）。接著到第二層級我們進入象限3裏面的方塊，在這個例子中，它是（圖5中的）第二個方塊。重復剛才的過程：我們的點應該落在哪個子象限？（這次是左下角，意味著位置的下一部分是1（二進制01），我們將進入的小方塊又是第二個。最後一次進入第三層級，重復這個過程，可以發現點落在右上角的子象限，因此位置的最後部分是3（二進制11）。把這些位置連接起來，我們得到點在曲線上的位置是二進制的110111，或者十進制的55。
更具體的來表達，寫出從x, y坐標到希爾伯特曲線位置轉換的方法。首先，我們要以計算機看得懂的形式，表達如下：
hilbert_map = {
’a’: {(0, 0): (0, ’d’), (0, 1): (1, ’a’), (1, 0): (3, ’b’), (1, 1): (2, ’a’)},
’b’: {(0, 0): (2, ’b’), (0, 1): (1, ’b’), (1, 0): (3, ’a’), (1, 1): (0, ’c’)},
’c’: {(0, 0): (2, ’c’), (0, 1): (3, ’d’), (1, 0): (1, ’c’), (1, 1): (0, ’b’)},
’d’: {(0, 0): (0, ’a’), (0, 1): (3, ’c’), (1, 0): (1, ’d’), (1, 1): (2, ’d’)},
}
上面的代碼中，每個hilbert_map的元素對應圖5四個方形中的一個。為了容易區分，我用一個字母來標識每個方塊：’a’是第一個方塊，’b’是第二個，等等。每個方塊的值是個字典，將(子)象限的x, y坐標映射到曲線上的位置（元組值的第一部分）以及下一個用到的方塊（元組值的第二部分）。下面的Python代碼展示了怎麽用這個來將x, y坐標轉換成希爾伯特曲線上的位置：
def point_to_hilbert(x, y, order=16):
current_square = ’a’
position = 0
for i in range(order - 1, -1, -1):
position <<= 2
quad_x = 1 if x & (1 << i) else 0
quad_y = 1 if y & (1 << i) else 0
quad_position, current_square = hilbert_map[current_square][(quad_x, quad_y)]
position |= quad_position
return position
函數的輸入是為整數的x, y坐標和曲線的階。一階曲線填充2×2的格子，二階曲線填充4×4的格子，等等。我們的x, y坐標應該先標準化到0到2order-1的區間。這個函數從最高位開始，逐步處理x, y坐標的每個比特位。在每個階段中，通過測試對應的比特位，可以確定坐標處於哪個（子）象限，還可以從我們之前定義的hilbert_map中取得在曲線上的位置以及下一個要用的方塊。在這階段取得的位置，加入到目前總的位置的最低兩位。在下一次循環的開頭，總的位置左移兩位以便給下一個位置騰出地方。
再運行一下之前的例子來檢驗一下函數寫對了沒有：
>>> point_to_hilbert(5,2,3)
55
如果這個代碼看懂了，再回頭來看
https://en.wikipedia.org/wiki/Hilbert_curve上的c語言代碼，應該就能看懂了，
且其定義的象限規則與之前完全相同。
//convert (x,y) to d
int xy2d (int n, int x, int y) {
int rx, ry, s, d=0;
for (s=n/2; s>0; s/=2) {
//由於s是2的某次方，即在該位為1，而其他位數全為0，則如果
x<s，x&s=0，x>=s ,x&s>0
rx = (x & s) > 0;
ry = (y & s) > 0;
//(3 * rx) ^ ry這一塊象限排序規則與之前的python代碼完全相同
d += s * s * ((3 * rx) ^ ry);
rot(s, &x, &y, rx, ry);
}
return d;
}
//convert d to (x,y)
void d2xy(int n, int d, int *x, int *y) {
int rx, ry, s, t=d;
*x = *y = 0;
for (s=1; s<n; s*=2) {
rx = 1 & (t/2);
ry = 1 & (t ^ rx);
rot(s, x, y, rx, ry);
*x += s * rx;
*y += s * ry;
t /= 4;
}
}
//旋轉象限
void rot(int n, int *x, int *y, int rx, int ry) {
if (ry == 0) {
if (rx == 1) {
*x = n-1 - *x;
*y = n-1 - *y;
}
//Swap x and y
int t = *x;
*x = *y;
*y = t;
}
Hilbert曲線以及其離散近似表示方法都非常實用，因為其將多維空間轉換為一維空間的方法很好地保留了空間鄰近性。（x，y）是一個單元方格中的點，d代表該點在Hilbert曲線上的位置，而由於其空間的鄰近性，在單元格上近似的點，其對應Hilbert的d值也比較接近。
正因為這種鄰近性，空間填充曲線被廣泛用於計算機科學，且在多維數據庫索引中，經常用Hilbert曲線取代Z order曲線。接下來我們看如何用代碼實現Hilbert算法。
現在假設我們以U字形來訪問區域。在每個象限中，我們同樣以U字形來訪問子象限，但是要調整好U字形的朝向使得和相鄰的象限銜接起來。如果我們正確地組織了這些U字形的朝向，我們就能完全消除不連續性，不管我們選擇了什麽分辨率，都能連續地訪問整個區域，可以在完全地探訪了一個區域後才移動到下一個。這個方案不僅消除了不連續性，而且提高了總體的局域性。這就是希爾伯特曲線。
如圖4，在第一層，枚舉這些點很簡單：選定一個方向和一個起始點，環繞四個象限，用0到3給他們編號。當我們要確定訪問子象限的順序同時維護總體的鄰接屬性，困難就來了。通過檢查我們發現，子象限的曲線是原曲線的簡單變換，而且只有四種變換。自然地，這個結論也適用於子子象限，等等。對於一個給定的象限，我們在其中畫出的曲線是由象限所在大的方形的曲線以及該象限的位置決定的。只需要費一點力，我們就能構建出如下概況所有情況的表
假設我們想用這個表來確定某個點在第三層希爾伯特曲線上的位置。在這個例子中，假設點的坐標是(5,2)。 從第一層級開始，則由上圖的第一個方形開始，找到該點所在的象限。在該例中，其在右上方的象限。那麽點在希爾伯特曲線上的位置的第一部分是3（二進制是11）。接著到第二層級我們進入象限3裏面的方塊，在這個例子中，它是（圖5中的）第二個方塊。重復剛才的過程：我們的點應該落在哪個子象限？（這次是左下角，意味著位置的下一部分是1（二進制01），我們將進入的小方塊又是第二個。最後一次進入第三層級，重復這個過程，可以發現點落在右上角的子象限，因此位置的最後部分是3（二進制11）。把這些位置連接起來，我們得到點在曲線上的位置是二進制的110111，或者十進制的55。
更具體的來表達，寫出從x, y坐標到希爾伯特曲線位置轉換的方法。首先，我們要以計算機看得懂的形式，表達如下：
hilbert_map = {
’a’: {(0, 0): (0, ’d’), (0, 1): (1, ’a’), (1, 0): (3, ’b’), (1, 1): (2, ’a’)},
’b’: {(0, 0): (2, ’b’), (0, 1): (1, ’b’), (1, 0): (3, ’a’), (1, 1): (0, ’c’)},
’c’: {(0, 0): (2, ’c’), (0, 1): (3, ’d’), (1, 0): (1, ’c’), (1, 1): (0, ’b’)},
’d’: {(0, 0): (0, ’a’), (0, 1): (3, ’c’), (1, 0): (1, ’d’), (1, 1): (2, ’d’)},
}
上面的代碼中，每個hilbert_map的元素對應圖5四個方形中的一個。為了容易區分，我用一個字母來標識每個方塊：’a’是第一個方塊，’b’是第二個，等等。每個方塊的值是個字典，將(子)象限的x, y坐標映射到曲線上的位置（元組值的第一部分）以及下一個用到的方塊（元組值的第二部分）。下面的Python代碼展示了怎麽用這個來將x, y坐標轉換成希爾伯特曲線上的位置：
def point_to_hilbert(x, y, order=16):
current_square = ’a’
position = 0
for i in range(order - 1, -1, -1):
position <<= 2
quad_x = 1 if x & (1 << i) else 0
quad_y = 1 if y & (1 << i) else 0
quad_position, current_square = hilbert_map[current_square][(quad_x, quad_y)]
position |= quad_position
return position
函數的輸入是為整數的x, y坐標和曲線的階。一階曲線填充2×2的格子，二階曲線填充4×4的格子，等等。我們的x, y坐標應該先標準化到0到2order-1的區間。這個函數從最高位開始，逐步處理x, y坐標的每個比特位。在每個階段中，通過測試對應的比特位，可以確定坐標處於哪個（子）象限，還可以從我們之前定義的hilbert_map中取得在曲線上的位置以及下一個要用的方塊。在這階段取得的位置，加入到目前總的位置的最低兩位。在下一次循環的開頭，總的位置左移兩位以便給下一個位置騰出地方。
再運行一下之前的例子來檢驗一下函數寫對了沒有：
>>> point_to_hilbert(5,2,3)
55
如果這個代碼看懂了，再回頭來看
https://en.wikipedia.org/wiki/Hilbert_curve上的c語言代碼，應該就能看懂了，
且其定義的象限規則與之前完全相同。
//convert (x,y) to d
int xy2d (int n, int x, int y) {
int rx, ry, s, d=0;
for (s=n/2; s>0; s/=2) {
//由於s是2的某次方，即在該位為1，而其他位數全為0，則如果
x<s，x&s=0，x>=s ,x&s>0
rx = (x & s) > 0;
ry = (y & s) > 0;
//(3 * rx) ^ ry這一塊象限排序規則與之前的python代碼完全相同
d += s * s * ((3 * rx) ^ ry);
rot(s, &x, &y, rx, ry);
}
return d;
}
//convert d to (x,y)
void d2xy(int n, int d, int *x, int *y) {
int rx, ry, s, t=d;
*x = *y = 0;
for (s=1; s<n; s*=2) {
rx = 1 & (t/2);
ry = 1 & (t ^ rx);
rot(s, x, y, rx, ry);
*x += s * rx;
*y += s * ry;
t /= 4;
}
}
//旋轉象限
void rot(int n, int *x, int *y, int rx, int ry) {
if (ry == 0) {
if (rx == 1) {
*x = n-1 - *x;
*y = n-1 - *y;
}
//Swap x and y
int t = *x;
*x = *y;
*y = t;
}
}
北京激光祛斑：www.6ysh.com
Hilbert曲線介紹以及代碼實現