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

一、k-d樹的介紹與實現(xiàn)

(相關(guān)資料圖)

1.1 k-d樹的創(chuàng)建

k-d樹其實就是一種樹形的數(shù)據(jù)結(jié)構(gòu)，但是在創(chuàng)建這棵樹時有一些固定的規(guī)則。下面來講一下kd樹的創(chuàng)建過程

輸入：一組數(shù)據(jù)點集，n個數(shù)據(jù)點，每個點有m維

輸出：k-d樹的根結(jié)點指針

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

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

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

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

例子：以2維平面上的點集為例，設(shè)有6個二維數(shù)據(jù)點{(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），因此根節(jié)點為（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）依次遞歸進行分割，最終形成的分割圖和樹狀結(jié)構(gòu)如下所示：

圖2 上例中形成的分割圖

圖3 上例中形成的樹狀結(jié)構(gòu)

1.2 k-d樹的查詢

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

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

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

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

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

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

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

如下圖所示：

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

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

同上，首先我們判斷（7,2）的分割標準是x周，而2<7，因此到（7,2）的左孩子進行查找，而(5,4)的分割標準為y軸，而4.5>4因此3.041，因此需要到（5,4）的右孩子進行查找，找到了葉子節(jié)點（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，因此更新最近鄰點和最近鄰距離。回溯到（7，2），可以判斷不需到（7,2）的右子樹進行查找，如下圖所示：

1.3 代碼實現(xiàn)

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

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

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

//創(chuàng)建kd樹，keypoints為點數(shù)據(jù)，parent表示當(dāng)前樹的雙親，默認為NULLkdNode* kdTree::createTree(vector&keypoints, kdNode *parent){if (keypoints.size() == 0)//若數(shù)據(jù)點集為空，則停止創(chuàng)建return NULL;int sort_dim = findSortDim(keypoints, parent);//確定分割的維度kdNode *tmp = findMidNode(keypoints);//找到分割結(jié)點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;//記錄當(dāng)前結(jié)點的分割維度tmp->left = createTree(leftKeyPoints, tmp);//遞歸調(diào)用，創(chuàng)建左子樹tmp->right = createTree(rightKeyPoints, tmp);//遞歸調(diào)用，創(chuàng)建右子樹tmp->parent = parent;//記錄父節(jié)點return tmp;//返回當(dāng)前樹的根節(jié)點}

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

//通過kd樹查找距離指定點node最近的點//root是查找的kd樹的根節(jié)點//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不是指向葉節(jié)點{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);//計算查找點與近似鄰近葉子節(jié)點的距離min_subscript = p->data.subscript;//記錄最近鄰結(jié)點在數(shù)據(jù)點集中的下標，以便以后找到它kdNode* q = p;kdNode* tmp = q;//開始回溯while (q != root){q = tmp->parent;//當(dāng)前結(jié)點距離查找點的距離float tmp_dis = calcDistance(point, q->data);//當(dāng)tmp_dis小于最近鄰距離時，更新最近鄰和次近鄰if (tmp_dis < min_dis){secmin_dis = min_dis;min_dis = tmp_dis;min_subscript = q->data.subscript;}//當(dāng)tmp_dis大于等于最近鄰且小于次近鄰時，更新次近鄰else if (tmp_dis == min_dis || tmp_dis < secmin_dis){secmin_dis = tmp_dis;}//查找點距離當(dāng)前結(jié)點構(gòu)成的區(qū)域分割線的垂直距離float sortdim_dis = std::fabs(point.data[q->sort_dim] - q->data.data[q->sort_dim]);//若垂直距離小于距離當(dāng)前結(jié)點的距離//則證明以查找點為中心，以到當(dāng)前結(jié)點距離為半徑畫圓，會與該結(jié)點構(gòu)成的區(qū)域分割線相交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為查找點距離當(dāng)前結(jié)點的另一邊的子樹的最小距離float tmp_nearest_dis = tmpResult.nearest_dis;float tmp_sec_nearest_dis = tmpResult.sec_nearest_dis;//當(dāng)子樹中距離查找點的最小距離小于當(dāng)前記錄的最鄰近距離時，更新最近鄰和次近鄰距離if (tmp_nearest_dis < min_dis){secmin_dis = min_dis;min_dis = tmp_nearest_dis;min_subscript = tmpResult.point_subscript;}//當(dāng)子樹中距離查找點的最小距離在最近鄰和次近鄰距離之間時，更新次近鄰距離else if (tmp_nearest_dis == min_dis || tmp_nearest_dis < secmin_dis)secmin_dis = tmp_nearest_dis;//當(dāng)子樹中距離查找點的次近鄰距離小于更新后的次近鄰距離時，再次更新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表示的是最近鄰距離，次近鄰距離和最近鄰點在數(shù)據(jù)點集中的下標，為了后面的SIFT算法會用到。 上面的描述就是k-d樹的建立和利用k-d樹找最近鄰的方法了。在實際應(yīng)用中k-d樹更加適合于低維的數(shù)據(jù)中，或者說如果數(shù)據(jù)量遠大于數(shù)據(jù)的維度的時候，使用k-d樹的效率與線性查找的方法相比還是有很大的提升的。但是我在實際應(yīng)用時，一張圖像中通常有2000+個特征點，而SIFT特征為128維的，所以加速效果也不是很好。實際上，在我的實驗中，甚至不如暴力匹配的效率高（當(dāng)然，這可能跟我的代碼質(zhì)量有關(guān)）。因此也就引出了我們接下來要介紹的bbf算法。

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

二、bbf（Best Bin First）算法介紹與實現(xiàn)

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

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;//最近鄰點在點集中的下標//優(yōu)先級隊列，查找點到當(dāng)前點確定的分割超平面距離越小優(yōu)先級越大priority_queuepri_queue;//priorityInfo類型包含了如下信息：//(1)當(dāng)前的結(jié)點指針，指向kdNode類型//(2)當(dāng)前點到查找點的歐式距離//(3)以及查找點到當(dāng)前點確定的分割超平面的距離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記錄嘗試更新最近鄰的次數(shù)while (!pri_queue.empty()){t++;priorityInfo tmp = pri_queue.top();pri_queue.pop();int sort_dim = tmp.ptr->sort_dim;//如果最近鄰距離小于查找點到當(dāng)前點確定的分割超平面的距離則不訪問該點的分支if (min_dis < fabs(point.data[sort_dim] - tmp.ptr->data.data[sort_dim]))continue;//記錄當(dāng)前點到查找點的歐式距離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;//遍歷以當(dāng)前點為根的子樹，直到葉子節(jié)點while (q->right != NULL || q->left != NULL){t++;int s_d = q->sort_dim;if (point.data[s_d] <= q-="">data.data[s_d])//查找點在分割維的大小小于當(dāng)前點分割維的大小{if (q->left != NULL)//進入左孩子之前判斷左孩子是否為空{(diào)if (q->right != NULL)//把右孩子加入節(jié)點時判斷右孩子是否為空{(diào)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)//如果更新次數(shù)超過600次則直接退出循環(huán)，返回當(dāng)前最近鄰結(jié)果break;}nearestNodeInfo result(min_dis, sec_min_dis, min_subscript);return result;}

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