k-d tree演算法
k-d樹（k-dimensional樹的簡稱），是一種分割k維資料空間的資料結構。主要應用於多維空間關鍵資料的搜尋（如：範圍搜尋和最近鄰搜尋）。
應用背景
SIFT演算法中做特徵點匹配的時候就會利用到k-d樹。而特徵點匹配實際上就是一個通過距離函式在高維向量之間進行相似性檢索的問題。針對如何快速而準確地找到查詢點的近鄰，現在提出了很多高維空間索引結構和近似查詢的演算法，k-d樹就是其中一種。
索引結構中相似性查詢有兩種基本的方式：一種是範圍查詢（range searches），另一種是K近鄰查詢（K-neighbor searches）。範圍查詢就是給定查詢點和查詢距離的閾值，從資料集中找出所有與查詢點距離小於閾值的資料；K近鄰查詢是給定查詢點及正整數K，從資料集中找到距離查詢點最近的K個數據，當K=1時，就是最近鄰查詢（nearest neighbor searches）。
特徵匹配運算元大致可以分為兩類。一類是線性掃描法，即將資料集中的點與查詢點逐一進行距離比較，也就是窮舉，缺點很明顯，就是沒有利用資料集本身蘊含的任何結構資訊，搜尋效率較低，第二類是建立資料索引，然後再進行快速匹配。因為實際資料一般都會呈現出簇狀的聚類形態，通過設計有效的索引結構可以大大加快檢索的速度。索引樹屬於第二類，其基本思想就是對搜尋空間進行層次劃分。根據劃分的空間是否有混疊可以分為Clipping和Overlapping兩種。前者劃分空間沒有重疊，其代表就是k-d樹；後者劃分空間相互有交疊，其代表為R樹。（這裡只介紹k-d樹）
例項
先以一個簡單直觀的例項來介紹k-d樹演算法。假設有6個二維資料點{（2,3），（5,4），（9,6），（4,7），（8,1），（7,2）}，資料點位於二維空間內（如圖1中黑點所示）。k-d樹演算法就是要確定圖1中這些分割空間的分割線（多維空間即為分割平面，一般為超平面）。下面就要通過一步步展示k-d樹是如何確定這些分割線的。
圖1 二維資料k-d樹空間劃分示意圖
k-d樹演算法可以分為兩大部分，一部分是有關k-d樹本身這種資料結構建立的演算法，另一部分是在建立的k-d樹上如何進行最鄰近查詢的演算法。
k-d樹構建演算法
k-d樹是一個二叉樹，每個節點表示一個空間範圍。表1給出的是k-d樹每個節點中主要包含的資料結構。
表1 k-d樹中每個節點的資料型別
| 域名 | 資料型別 | 描述 |
| Node-data | 資料向量 | 資料集中某個資料點，是n維向量（這裡也就是k維） |
| Range | 空間向量 | 該節點所代表的空間範圍 |
| split | 整數 | 垂直於分割超平面的方向軸序號 |
| Left | k-d樹 | 由位於該節點分割超平面左子空間內所有資料點所構成的k-d樹 |
| Right | k-d樹 | 由位於該節點分割超平面右子空間內所有資料點所構成的k-d樹 |
| parent | k-d樹 | 父節點 |
從上面對k-d樹節點的資料型別的描述可以看出構建k-d樹是一個逐級展開的遞迴過程。表2給出的是構建k-d樹的偽碼。
表2 構建k-d樹的偽碼
| 演算法：構建k-d樹（createKDTree） |
| 輸入：資料點集Data-set和其所在的空間Range |
| 輸出：Kd，型別為k-d tree |
| 1.If Data-set為空，則返回空的k-d tree |
| 2.呼叫節點生成程式： （1）確定split域：對於所有描述子資料（特徵向量），統計它們在每個維上的資料方差。以SURF特徵為例，描述子為64維，可計算64個方差。挑選出最大值，對應的維就是split域的值。資料方差大表明沿該座標軸方向上的資料分散得比較開，在這個方向上進行資料分割有較好的解析度； （2）確定Node-data域：資料點集Data-set按其第split域的值排序。位於正中間的那個資料點被選為Node-data。此時新的Data-set' = Data-set\Node-data（除去其中Node-data這一點）。 |
| 3.dataleft = {d屬於Data-set' && d[split] ≤ Node-data[split]} Left_Range = {Range && dataleft}dataright = {d屬於Data-set' && d[split] > Node-data[split]} Right_Range = {Range && dataright} |
| 4.left = 由（dataleft，Left_Range）建立的k-d tree，即遞迴呼叫createKDTree（dataleft，Left_ Range）。並設定left的parent域為Kd； right = 由（dataright，Right_Range）建立的k-d tree，即呼叫createKDTree（dataleft，Left_ Range）。並設定right的parent域為Kd。 |
以上述舉的例項來看，過程如下：
由於此例簡單，資料維度只有2維，所以可以簡單地給x，y兩個方向軸編號為0,1，也即split={0,1}。
（1）確定split域的首先該取的值。分別計算x，y方向上資料的方差得知x方向上的方差最大，所以split域值首先取0，也就是x軸方向；
（2）確定Node-data的域值。根據x軸方向的值2,5,9,4,8,7排序選出中值為7，所以Node-data = （7,2）。這樣，該節點的分割超平面就是通過（7,2）並垂直於split = 0（x軸）的直線x = 7；
（3）確定左子空間和右子空間。分割超平面x = 7將整個空間分為兩部分，如圖2所示。x < = 7的部分為左子空間，包含3個節點{（2,3），（5,4），（4,7）}；另一部分為右子空間，包含2個節點{（9,6），（8,1）}。
圖2 x=7將整個空間分為兩部分
如演算法所述，k-d樹的構建是一個遞迴的過程。然後對左子空間和右子空間內的資料重複根節點的過程就可以得到下一級子節點（5,4）和（9,6）（也就是左右子空間的'根'節點），同時將空間和資料集進一步細分。如此反覆直到空間中只包含一個數據點，如圖1所示。最後生成的k-d樹如圖3所示。
圖3 上述例項生成的k-d樹
注意：每一級節點旁邊的'x'和'y'表示以該節點分割左右子空間時split所取的值。
k-d樹上的最鄰近查詢演算法
在k-d樹中進行資料的查詢也是特徵匹配的重要環節，其目的是檢索在k-d樹中與查詢點距離最近的資料點。這裡先以一個簡單的例項來描述最鄰近查詢的基本思路。
星號表示要查詢的點（2.1,3.1）。通過二叉搜尋，順著搜尋路徑很快就能找到最鄰近的近似點，也就是葉子節點（2,3）。而找到的葉子節點並不一定就是最鄰近的，最鄰近肯定距離查詢點更近，應該位於以查詢點為圓心且通過葉子節點的圓域內。為了找到真正的最近鄰，還需要進行'回溯'操作：演算法沿搜尋路徑反向查詢是否有距離查詢點更近的資料點。此例中先從（7,2）點開始進行二叉查詢，然後到達（5,4），最後到達（2,3），此時搜尋路徑中的節點為<（7,2），（5,4），（2,3）>，首先以（2,3）作為當前最近鄰點，計算其到查詢點（2.1,3.1）的距離為0.1414，然後回溯到其父節點（5,4），並判斷在該父節點的其他子節點空間中是否有距離查詢點更近的資料點。以（2.1,3.1）為圓心，以0.1414為半徑畫圓，如圖4所示。發現該圓並不和超平面y = 4交割，因此不用進入（5,4）節點右子空間中去搜索。
圖4 查詢（2.1，3.1）點的兩次回溯判斷
再回溯到（7,2），以（2.1,3.1）為圓心，以0.1414為半徑的圓更不會與x = 7超平面交割，因此不用進入（7,2）右子空間進行查詢。至此，搜尋路徑中的節點已經全部回溯完，結束整個搜尋，返回最近鄰點（2,3），最近距離為0.1414。
一個複雜點了例子如查詢點為（2，4.5）。同樣先進行二叉查詢，先從（7,2）查詢到（5,4）節點，在進行查詢時是由y = 4為分割超平面的，由於查詢點為y值為4.5，因此進入右子空間查詢到（4,7），形成搜尋路徑<（7,2），（5,4），（4,7）>，取（4,7）為當前最近鄰點，計算其與目標查詢點的距離為3.202。然後回溯到（5,4），計算其與查詢點之間的距離為3.041。以（2，4.5）為圓心，以3.041為半徑作圓，如圖5所示。可見該圓和y = 4超平面交割，所以需要進入（5,4）左子空間進行查詢。此時需將（2,3）節點加入搜尋路徑中得<（7,2），（2,3）>。回溯至（2,3）葉子節點，（2,3）距離（2,4.5）比（5,4）要近，所以最近鄰點更新為（2，3），最近距離更新為1.5。回溯至（7,2），以（2,4.5）為圓心1.5為半徑作圓，並不和x = 7分割超平面交割，如圖6所示。至此，搜尋路徑回溯完。返回最近鄰點（2,3），最近距離1.5。k-d樹查詢演算法的虛擬碼如表3所示。
圖5 查詢（2，4.5）點的第一次回溯判斷
圖6 查詢（2，4.5）點的第二次回溯判斷
表3 標準k-d樹查詢演算法
| 演算法：k-d樹最鄰近查詢 |
| 輸入：Kd， //k-d tree型別 target //查詢資料點 |
| 輸出：nearest， //最鄰近資料點 dist //最鄰近資料點和查詢點間的距離 |
| 1. If Kd為NULL，則設dist為infinite並返回 |
| 2. //進行二叉查詢，生成搜尋路徑 Kd_point = &Kd； //Kd-point中儲存k-d tree根節點地址 nearest = Kd_point -> Node-data； //初始化最近鄰點 while（Kd_point） push（Kd_point）到search_path中； //search_path是一個堆疊結構，儲存著搜尋路徑節點指標 /*** If Dist（nearest，target） > Dist（Kd_point -> Node-data，target） nearest = Kd_point -> Node-data； //更新最近鄰點 Max_dist = Dist(Kd_point，target）； //更新最近鄰點與查詢點間的距離 ***/ s = Kd_point -> split； //確定待分割的方向 If target[s] <= Kd_point -> Node-data[s] //進行二叉查詢 Kd_point = Kd_point -> left； else Kd_point = Kd_point ->right； nearest = search_path中最後一個葉子節點； //注意：二叉搜尋時不比計算選擇搜尋路徑中的最鄰近點，這部分已被註釋 Max_dist = Dist（nearest，target）； //直接取最後葉子節點作為回溯前的初始最近鄰點 |
| 3. //回溯查詢 while（search_path != NULL） back_point = 從search_path取出一個節點指標； //從search_path堆疊彈棧 s = back_point -> split； //確定分割方向 If Dist（target[s]，back_point -> Node-data[s]） < Max_dist //判斷還需進入的子空間 If target[s] <= back_point -> Node-data[s] Kd_point = back_point -> right； //如果target位於左子空間，就應進入右子空間 else Kd_point = back_point -> left; //如果target位於右子空間，就應進入左子空間 將Kd_point壓入search_path堆疊； If Dist（nearest，target） > Dist（Kd_Point -> Node-data，target） nearest = Kd_point -> Node-data； //更新最近鄰點 Min_dist = Dist（Kd_point -> Node-data,target）； //更新最近鄰點與查詢點間的距離 |
上述兩次例項表明，當查詢點的鄰域與分割超平面兩側空間交割時，需要查詢另一側子空間，導致檢索過程複雜，效率下降。研究表明N個節點的K維k-d樹搜尋過程時間複雜度為：tworst=O（kN1-1/k）。
後記
以上為了介紹方便，討論的是二維情形。像實際的應用中，如SIFT特徵向量128維，SURF特徵向量64維，維度都比較大，直接利用k-d樹快速檢索（維數不超過20）的效能急劇下降。假設資料集的維數為D，一般來說要求資料的規模N滿足N»2D，才能達到高效的搜尋。所以這就引出了一系列對k-d樹演算法的改進。有待進一步研究學習。
參考
1.《影象區域性不變特性特徵與描述》王永明 王貴錦 編著 國防工業出版社
2.http://underthehood.blog.51cto.com/2531780/687160
轉載請註明：http://www.cnblogs.com/eyeszjwang/articles/2429382.html