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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 分割示意圖

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

（3）依次遞歸進(jìn)行分割，最終形成的分割圖和樹狀結(jié)構(gòu)如下所示：

圖2 上例中形成的分割圖

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

1.2 k-d樹的查詢

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

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

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

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

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

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

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

如下圖所示：

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

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

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

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

1.3 代碼實(shí)現(xiàn)

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

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

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

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

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

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

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

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

根據(jù)前面k-d樹的搜索過程我們可以知道，在搜索時(shí)首先沿著kd樹找到葉子節(jié)點(diǎn)，然后依次回溯，而回溯的路程就是前面我們查找葉子節(jié)點(diǎn)時(shí)逆序，因此進(jìn)行回溯時(shí)并沒有利用這些點(diǎn)的信息。我們接下來介紹的算法就是利用這些信息，回溯時(shí)給各個(gè)需要回溯的結(jié)點(diǎn)以優(yōu)先級，這樣找到最近鄰會更快。接下來詳細(xì)介紹bbf算法的流程。 其實(shí)bbf算法的思想比較簡單，通過對回溯可能需要的路過的結(jié)點(diǎn)加入隊(duì)列，并按照查找點(diǎn)到該結(jié)點(diǎn)確定的超平面的距離進(jìn)行排序，然后每次首先遍歷的是優(yōu)先級最高（即距離最短的結(jié)點(diǎn)），直到隊(duì)列為空算法結(jié)束。同時(shí)bbf算法也設(shè)立了一個(gè)時(shí)間限制，如果算法運(yùn)行時(shí)間超過該限制，不管是不是為空，一律停止運(yùn)行，返回當(dāng)前的最近鄰點(diǎn)作為結(jié)果。 bbf的算法流程如下： 輸入：kd樹，查找點(diǎn)x 輸出：kd樹種距離查找點(diǎn)最近的點(diǎn)以及最近的距離 流程：（1）若kd樹為空，則設(shè)定兩者距離為無窮大，返回；如果kd樹非空，則將kd樹的根節(jié)點(diǎn)加入到優(yōu)先級隊(duì)列中；             （2）從優(yōu)先級隊(duì)列中出隊(duì)當(dāng)前優(yōu)先級最大的結(jié)點(diǎn)，計(jì)算當(dāng)前的該點(diǎn)到查找點(diǎn)的距離是否比最近鄰距離小，如果是則更新最近鄰點(diǎn)和最近鄰距離。如果查找點(diǎn)在切分維坐標(biāo)小于當(dāng)前點(diǎn)的切分維坐標(biāo)，則把他的右孩子加入到隊(duì)列中，同時(shí)檢索它的左孩子，否則就把他的左孩子加入到隊(duì)列中，同時(shí)檢索它的右孩子。這樣一直重復(fù)檢索，并加入隊(duì)列，直到檢索到葉子節(jié)點(diǎn)。然后在從優(yōu)先級隊(duì)列中出隊(duì)優(yōu)先級最大的結(jié)點(diǎn)；             （3）重復(fù)（1）和（2）中的操作，直到優(yōu)先級隊(duì)列為空，或者超出規(guī)定的時(shí)間，返回當(dāng)前的最近鄰結(jié)點(diǎn)和距離。 實(shí)現(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;//最近鄰點(diǎn)在點(diǎn)集中的下標(biāo)//優(yōu)先級隊(duì)列，查找點(diǎn)到當(dāng)前點(diǎn)確定的分割超平面距離越小優(yōu)先級越大priority_queuepri_queue;//priorityInfo類型包含了如下信息：//(1)當(dāng)前的結(jié)點(diǎn)指針，指向kdNode類型//(2)當(dāng)前點(diǎn)到查找點(diǎn)的歐式距離//(3)以及查找點(diǎn)到當(dāng)前點(diǎn)確定的分割超平面的距離pri_queue.push(priorityInfo(p,calcDistance(point,root->data),                        fabs(point.data[root->sort_dim]-root->data.data[root->sort_dim])));int t = 0;//這里沒有記錄時(shí)間，使用t記錄嘗試更新最近鄰的次數(shù)while (!pri_queue.empty()){t++;priorityInfo tmp = pri_queue.top();pri_queue.pop();int sort_dim = tmp.ptr->sort_dim;//如果最近鄰距離小于查找點(diǎn)到當(dāng)前點(diǎn)確定的分割超平面的距離則不訪問該點(diǎn)的分支if (min_dis < fabs(point.data[sort_dim] - tmp.ptr->data.data[sort_dim]))continue;//記錄當(dāng)前點(diǎn)到查找點(diǎn)的歐式距離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)前點(diǎn)為根的子樹，直到葉子節(jié)點(diǎn)while (q->right != NULL || q->left != NULL){t++;int s_d = q->sort_dim;if (point.data[s_d] <= q-="">data.data[s_d])//查找點(diǎn)在分割維的大小小于當(dāng)前點(diǎn)分割維的大小{if (q->left != NULL)//進(jìn)入左孩子之前判斷左孩子是否為空{(diào)if (q->right != NULL)//把右孩子加入節(jié)點(diǎn)時(shí)判斷右孩子是否為空{(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時(shí)運(yùn)行情況已經(jīng)同暴力查找時(shí)效率相當(dāng)，如果想要加速，把閾值設(shè)的低一些。但是如果閾值設(shè)的太低會造成匹配結(jié)果較差，需要在效率和正確率上進(jìn)行取舍。