通訊！k-d樹和bbf算法 一直遞歸子樹的數據點集算法

最近還是一直在研究SIFT算法，而SIFT特征點匹配是一個比較經典的問題，使用暴力匹配的話確實可以得到結果，但是運行速度較慢。我的計算機處理是i5的二代系列，匹配兩張各檢測有2000+個SIFT特征點的圖像，通過正反匹配（即取圖像1與圖像2的匹配結果余圖像2和圖像1的匹配結果的交集），再加上OpenMP多線程加速，使用暴力匹配，大概要花20多秒，還是比較慢的。所以這一周啥也沒做，一直在實現kd樹和對應的bbf算法。下面詳細介紹下種數據結構。

一、k-d樹的介紹與實現

(相關資料圖)

1.1 k-d樹的創建

k-d樹其實就是一種樹形的數據結構，但是在創建這棵樹時有一些固定的規則。下面來講一下kd樹的創建過程

輸入：一組數據點集，n個數據點，每個點有m維

輸出：k-d樹的根結點指針

過程：（1）分別計算這n個數據點在m維中各個維度的方差，取方差最大的維度dim作為分割維度；

（2）把數據點集按照該維度中值的大小進行排列，選擇具有中間值的點作為該樹的根結點；

（3）前半部分點進行如（1）、（2）所示的遞歸操作，選出的遞歸子樹的根節點作為（2）中得到的根節點的左孩子；

同理，后半部分也這樣操作。如此一直遞歸，直到各個遞歸子樹的數據點集為空則算法截止。

例子：以2維平面上的點集為例，設有6個二維數據點{(2,3)，(5,4)，(9,6)，(4,7)，(8,1)，(7,2)}。

（1） 首先計算這6個點的橫坐標和縱坐標的方差值，橫坐標的方差值為39，縱坐標上的方差值為28.63，因此第一次分割取橫坐標上的值作為分割標準。把這些點按照橫坐標進行排序得到{(2,3)，(4,7)，(5,4)，(7,2)，(8,1)，(9,6)}，取中間點為（7,2），因此根節點為（7,2）進行分割，如下圖所示：

圖1 分割示意圖

（2）接下來對{(2,3)，(4,7)，(5,4)}和{(8,1)，(9,6)}分別進行分割，在{(2,3)，(4,7)，(5,4)}中縱坐標的方差較大，因此按縱坐標進行排序后分割，則（5,4）為（7,2）的左孩子，{(8,1)，(9,6)}中也是縱坐標方差較大，因此選縱坐標進行排序后分割，這里算則（9,6）作為（7,2）的右孩子。

（3）依次遞歸進行分割，最終形成的分割圖和樹狀結構如下所示：

圖2 上例中形成的分割圖

圖3 上例中形成的樹狀結構

1.2 k-d樹的查詢

k-d樹建立好以后，需要查詢它的最近鄰，方法如下：

（1）查詢點與k-d樹的根節點進行比較，比較兩者在根節點劃分時的維度的值的大小，若查詢點在該維的值小，則進入根節點的左子樹，否則進入右子樹。依次類推，進行查找，直到到達樹的葉子節點。

（2）設當前到達的葉子節點為目前的最近鄰（注意：可能并非真正的最近鄰），并且記錄目前的最近鄰距離。沿著來時的路向前回溯，讓目前的最近鄰距離與查找點與當前葉子節點的父節點形成的分割超平面的距離進行比較，若當前最近鄰比較小，則不用遍歷當前葉子節點的父節點的另一邊，否則需要遍歷查找以更新最近鄰距離和最近鄰節點。

（3）按照（2）中所說依次遍歷，直到到達根節點為止，查詢結束。

上面的說法比較抽象，下面用兩個博客中廣為流傳的例子進行解讀。

假設我們需要查找點（2.1,3.1）在前面中提到的二維點集中的最近鄰點。我們首先判斷（7,2）的分割標準是x軸，而2.1<7，因此查找點進入（7,2）的左子樹；而（5,4）的分割標準是y軸，而3.1<4，因此我們進入（5,4）的左子樹，即找到葉子節點（2,3）；把（2,3）作為查找點（2.1,3.1）的臨時最近鄰點，最近鄰距離為0.1414，向前回溯。

因為查找點到（5,4）的距離大于到（2,3）的距離，因此最近鄰點和最近鄰距離保持不變，因此以（2.1,3.1）為原點，以0.1414為半徑畫圓，該圓與（5,4）確定的分割線沒有相交（即當前最近鄰距離比查找點到（5,4）所確定的分割線距離要小），因此不需要進入（5,4）的右子樹，繼續回溯，同理，最近鄰點和最近鄰距離不變，以（2.1,3.1）為原點，以0.1414為半徑畫圓，該圓與（7,2）所確定的分割線也沒有相交，因此也不需要進入（7,2）的右子樹；回溯結束。因此（2,3）就是真正的最近鄰節點。

如下圖所示：

圖4 （2.1,3.1）查詢最近鄰示意圖

上面這個例子比較簡單，下面我們看一個復雜一些的例子，假設我們要查找（2,4.5）的最近鄰。

同上，首先我們判斷（7,2）的分割標準是x周，而2<7，因此到（7,2）的左孩子進行查找，而(5,4)的分割標準為y軸，而4.5>4因此3.041，因此需要到（5,4）的右孩子進行查找，找到了葉子節點（4,7）。那么我們把（4,7）作為查找點的臨時最近鄰，最近鄰距離為3.202，向前回溯，可以看到到（5,4）的距離為3.041，因此更新（5,4）為最近鄰點，最近鄰距離為3.041。然后以（2,4.5）為圓心，以3.041為半徑畫圓，可以看到該圓與（5,4）確定的分割線相交，因此需要遍歷（5,4）的左子樹。如下圖所示：

判斷（2,4.5）到（2,3）的距離為1.5，因此更新最近鄰點和最近鄰距離?；厮莸剑?，2），可以判斷不需到（7,2）的右子樹進行查找，如下圖所示：

1.3 代碼實現

k-d樹的實現還算是比較簡單的，在我的實現過程中遇到的問題是開始我沒有理解前面提到的圓與分割線相交的意義，所以實現時遇到了一些問題，現在把我實現的kd樹的核心算法一一介紹。 （1）kd樹的結點數據結構

class kdNode{public:kdNode(Point &data);~kdNode();Point data;//數據點的信息int sort_dim;//數據點的劃分維度kdNode *left;kdNode *right;kdNode *parent;};

數據結構算是比較簡單的，只包含了數據點的信息（Point類是我自己定義的），left和right是左右孩子的指針，parent是父節點指針，在回溯時會用到；sort_dim是記錄當前結點時按照哪個維度進行劃分的，在回溯時判斷最近鄰和查找點到當前結點確定的分割超平面的距離哪個大時會用到。 （2）創建kdTree代碼

//創建kd樹，keypoints為點數據，parent表示當前樹的雙親，默認為NULLkdNode* kdTree::createTree(vector&keypoints, kdNode *parent){if (keypoints.size() == 0)//若數據點集為空，則停止創建return NULL;int sort_dim = findSortDim(keypoints, parent);//確定分割的維度kdNode *tmp = findMidNode(keypoints);//找到分割結點int sort_num = keypoints.size() / 2;vectorleftKeyPoints(keypoints.begin(), keypoints.begin() + sort_num);vectorrightKeyPoints(keypoints.begin() + sort_num + 1, keypoints.end());tmp->sort_dim = sort_dim;//記錄當前結點的分割維度tmp->left = createTree(leftKeyPoints, tmp);//遞歸調用，創建左子樹tmp->right = createTree(rightKeyPoints, tmp);//遞歸調用，創建右子樹tmp->parent = parent;//記錄父節點return tmp;//返回當前樹的根節點}

這里面findMidNode函數是找到當前數據點的分割結點，在這里面對keypoints按照各點在分割維度上的大小進行了排序，因此后面直接把數據點集分成了兩部分。 （3）查找最近鄰結點

//通過kd樹查找距離指定點node最近的點//root是查找的kd樹的根節點//point是查找點nearestNodeInfo& kdTree::findNearestNode(kdNode* root, const Point& point){if (root == NULL){return nearestNodeInfo();}kdNode *p = root;//通過kd樹的二叉搜索，順著搜索路徑很快就能找到最鄰近的近似點while ((p->left != NULL) || (p->right != NULL))//只要p不是指向葉節點{int sort_dim = p->sort_dim;if (point.data[sort_dim] <= p-="">data.data[sort_dim]){if (p->left == NULL)break;p = p->left;}else{if (p->right == NULL)break;p = p->right;}}float min_dis = FLT_MAX;//距離查找點最近的距離float secmin_dis = FLT_MAX;//距離查找點的次近鄰距離int min_subscript = 0;min_dis = calcDistance(point, p->data);//計算查找點與近似鄰近葉子節點的距離min_subscript = p->data.subscript;//記錄最近鄰結點在數據點集中的下標，以便以后找到它kdNode* q = p;kdNode* tmp = q;//開始回溯while (q != root){q = tmp->parent;//當前結點距離查找點的距離float tmp_dis = calcDistance(point, q->data);//當tmp_dis小于最近鄰距離時，更新最近鄰和次近鄰if (tmp_dis < min_dis){secmin_dis = min_dis;min_dis = tmp_dis;min_subscript = q->data.subscript;}//當tmp_dis大于等于最近鄰且小于次近鄰時，更新次近鄰else if (tmp_dis == min_dis || tmp_dis < secmin_dis){secmin_dis = tmp_dis;}//查找點距離當前結點構成的區域分割線的垂直距離float sortdim_dis = std::fabs(point.data[q->sort_dim] - q->data.data[q->sort_dim]);//若垂直距離小于距離當前結點的距離//則證明以查找點為中心，以到當前結點距離為半徑畫圓，會與該結點構成的區域分割線相交if (sortdim_dis < min_dis){nearestNodeInfo tmpResult;if (tmp == q->left){tmpResult = findNearestNode(q->right, point);}else if (tmp == q->right){tmpResult = findNearestNode(q->left, point);}elsecout << "q is not parent of tmp" << endl;//tmpDis為查找點距離當前結點的另一邊的子樹的最小距離float tmp_nearest_dis = tmpResult.nearest_dis;float tmp_sec_nearest_dis = tmpResult.sec_nearest_dis;//當子樹中距離查找點的最小距離小于當前記錄的最鄰近距離時，更新最近鄰和次近鄰距離if (tmp_nearest_dis < min_dis){secmin_dis = min_dis;min_dis = tmp_nearest_dis;min_subscript = tmpResult.point_subscript;}//當子樹中距離查找點的最小距離在最近鄰和次近鄰距離之間時，更新次近鄰距離else if (tmp_nearest_dis == min_dis || tmp_nearest_dis < secmin_dis)secmin_dis = tmp_nearest_dis;//當子樹中距離查找點的次近鄰距離小于更新后的次近鄰距離時，再次更新if (tmp_sec_nearest_dis < secmin_dis)secmin_dis = tmp_sec_nearest_dis;}tmp = q;}nearestNodeInfo result(min_dis, secmin_dis, min_subscript);return result;}

這里的nearestNodeInfo表示的是最近鄰距離，次近鄰距離和最近鄰點在數據點集中的下標，為了后面的SIFT算法會用到。 上面的描述就是k-d樹的建立和利用k-d樹找最近鄰的方法了。在實際應用中k-d樹更加適合于低維的數據中，或者說如果數據量遠大于數據的維度的時候，使用k-d樹的效率與線性查找的方法相比還是有很大的提升的。但是我在實際應用時，一張圖像中通常有2000+個特征點，而SIFT特征為128維的，所以加速效果也不是很好。實際上，在我的實驗中，甚至不如暴力匹配的效率高（當然，這可能跟我的代碼質量有關）。因此也就引出了我們接下來要介紹的bbf算法。

前面講到了用k-d樹對于高維的數據進行最鄰近查詢時實際上效率并不高，這里介紹一個算法用以加速k-d樹對于高維數據的處理。

二、bbf（Best Bin First）算法介紹與實現

根據前面k-d樹的搜索過程我們可以知道，在搜索時首先沿著kd樹找到葉子節點，然后依次回溯，而回溯的路程就是前面我們查找葉子節點時逆序，因此進行回溯時并沒有利用這些點的信息。我們接下來介紹的算法就是利用這些信息，回溯時給各個需要回溯的結點以優先級，這樣找到最近鄰會更快。接下來詳細介紹bbf算法的流程。 其實bbf算法的思想比較簡單，通過對回溯可能需要的路過的結點加入隊列，并按照查找點到該結點確定的超平面的距離進行排序，然后每次首先遍歷的是優先級最高（即距離最短的結點），直到隊列為空算法結束。同時bbf算法也設立了一個時間限制，如果算法運行時間超過該限制，不管是不是為空，一律停止運行，返回當前的最近鄰點作為結果。 bbf的算法流程如下： 輸入：kd樹，查找點x 輸出：kd樹種距離查找點最近的點以及最近的距離 流程：（1）若kd樹為空，則設定兩者距離為無窮大，返回；如果kd樹非空，則將kd樹的根節點加入到優先級隊列中；             （2）從優先級隊列中出隊當前優先級最大的結點，計算當前的該點到查找點的距離是否比最近鄰距離小，如果是則更新最近鄰點和最近鄰距離。如果查找點在切分維坐標小于當前點的切分維坐標，則把他的右孩子加入到隊列中，同時檢索它的左孩子，否則就把他的左孩子加入到隊列中，同時檢索它的右孩子。這樣一直重復檢索，并加入隊列，直到檢索到葉子節點。然后在從優先級隊列中出隊優先級最大的結點；             （3）重復（1）和（2）中的操作，直到優先級隊列為空，或者超出規定的時間，返回當前的最近鄰結點和距離。 實現代碼如下：

nearestNodeInfo& kdTree::findNearestNode_bbf(kdNode* root, const Point& point){if (root == NULL)return nearestNodeInfo();kdNode *p = root;float min_dis = FLT_MAX;//最近鄰距離float sec_min_dis = FLT_MAX;//次近鄰距離int min_subscript = 0;//最近鄰點在點集中的下標//優先級隊列，查找點到當前點確定的分割超平面距離越小優先級越大priority_queuepri_queue;//priorityInfo類型包含了如下信息：//(1)當前的結點指針，指向kdNode類型//(2)當前點到查找點的歐式距離//(3)以及查找點到當前點確定的分割超平面的距離pri_queue.push(priorityInfo(p,calcDistance(point,root->data),                        fabs(point.data[root->sort_dim]-root->data.data[root->sort_dim])));int t = 0;//這里沒有記錄時間，使用t記錄嘗試更新最近鄰的次數while (!pri_queue.empty()){t++;priorityInfo tmp = pri_queue.top();pri_queue.pop();int sort_dim = tmp.ptr->sort_dim;//如果最近鄰距離小于查找點到當前點確定的分割超平面的距離則不訪問該點的分支if (min_dis < fabs(point.data[sort_dim] - tmp.ptr->data.data[sort_dim]))continue;//記錄當前點到查找點的歐式距離float tmp_dis = calcDistance(point, tmp.ptr->data);//判斷是否更新最近鄰、次近鄰距離if (tmp_dis < min_dis){sec_min_dis = min_dis;min_dis = tmp_dis;min_subscript = tmp.ptr->data.subscript;}else if (tmp_dis == min_dis || tmp_dis < sec_min_dis){sec_min_dis = tmp_dis;}kdNode* q = tmp.ptr;//遍歷以當前點為根的子樹，直到葉子節點while (q->right != NULL || q->left != NULL){t++;int s_d = q->sort_dim;if (point.data[s_d] <= q-="">data.data[s_d])//查找點在分割維的大小小于當前點分割維的大小{if (q->left != NULL)//進入左孩子之前判斷左孩子是否為空{if (q->right != NULL)//把右孩子加入節點時判斷右孩子是否為空{float distance = calcDistance(point, q->right->data);int s_t = q->right->sort_dim;pri_queue.push(priorityInfo(q->right, distance,fabs(point.data[s_t]-q->right->data.data[s_t])));}q = q->left;}elsebreak;}else{if (q->right != NULL){if (q->left != NULL){float distance = calcDistance(point, q->left->data);int s_t = q->left->sort_dim;pri_queue.push(priorityInfo(q->left, distance,fabs(point.data[s_t]-q->left->data.data[s_t])));}q = q->right;}elsebreak;}//更新最近鄰float dis = calcDistance(point, q->data);if (dis < min_dis){sec_min_dis = min_dis;min_dis = dis;min_subscript = q->data.subscript;}else if (dis == min_dis || dis < sec_min_dis)sec_min_dis = dis;}if (t > 600)//如果更新次數超過600次則直接退出循環，返回當前最近鄰結果break;}nearestNodeInfo result(min_dis, sec_min_dis, min_subscript);return result;}

這里t取600時運行情況已經同暴力查找時效率相當，如果想要加速，把閾值設的低一些。但是如果閾值設的太低會造成匹配結果較差，需要在效率和正確率上進行取舍。