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

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

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

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

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

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

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

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

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

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

同理，后半部分也這樣操作。如此一直遞歸，直到各個(gè)遞歸子樹(shù)的數(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）接下來(lái)對(duì){(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)行分割，最終形成的分割圖和樹(shù)狀結(jié)構(gòu)如下所示：

圖2 上例中形成的分割圖

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

1.2 k-d樹(shù)的查詢(xún)

k-d樹(shù)建立好以后，需要查詢(xún)它的最近鄰，方法如下：

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

（2）設(shè)當(dāng)前到達(dá)的葉子節(jié)點(diǎn)為目前的最近鄰（注意：可能并非真正的最近鄰），并且記錄目前的最近鄰距離。沿著來(lái)時(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）中所說(shuō)依次遍歷，直到到達(dá)根節(jié)點(diǎn)為止，查詢(xún)結(jié)束。

上面的說(shuō)法比較抽象，下面用兩個(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）的左子樹(shù)；而（5,4）的分割標(biāo)準(zhǔn)是y軸，而3.1<4，因此我們進(jìn)入（5,4）的左子樹(shù)，即找到葉子節(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為半徑畫(huà)圓，該圓與（5,4）確定的分割線(xiàn)沒(méi)有相交（即當(dāng)前最近鄰距離比查找點(diǎn)到（5,4）所確定的分割線(xiàn)距離要小），因此不需要進(jìn)入（5,4）的右子樹(shù)，繼續(xù)回溯，同理，最近鄰點(diǎn)和最近鄰距離不變，以（2.1,3.1）為原點(diǎn)，以0.1414為半徑畫(huà)圓，該圓與（7,2）所確定的分割線(xiàn)也沒(méi)有相交，因此也不需要進(jìn)入（7,2）的右子樹(shù)；回溯結(jié)束。因此（2,3）就是真正的最近鄰節(jié)點(diǎn)。

如下圖所示：

圖4 （2.1,3.1）查詢(xún)最近鄰示意圖

上面這個(gè)例子比較簡(jiǎn)單，下面我們看一個(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為半徑畫(huà)圓，可以看到該圓與（5,4）確定的分割線(xiàn)相交，因此需要遍歷（5,4）的左子樹(shù)。如下圖所示：

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

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

k-d樹(shù)的實(shí)現(xiàn)還算是比較簡(jiǎn)單的，在我的實(shí)現(xiàn)過(guò)程中遇到的問(wèn)題是開(kāi)始我沒(méi)有理解前面提到的圓與分割線(xiàn)相交的意義，所以實(shí)現(xiàn)時(shí)遇到了一些問(wèn)題，現(xiàn)在把我實(shí)現(xiàn)的kd樹(shù)的核心算法一一介紹。 （1）kd樹(shù)的結(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)算是比較簡(jiǎn)單的，只包含了數(shù)據(jù)點(diǎn)的信息（Point類(lèi)是我自己定義的），left和right是左右孩子的指針，parent是父節(jié)點(diǎn)指針，在回溯時(shí)會(huì)用到；sort_dim是記錄當(dāng)前結(jié)點(diǎn)時(shí)按照哪個(gè)維度進(jìn)行劃分的，在回溯時(shí)判斷最近鄰和查找點(diǎn)到當(dāng)前結(jié)點(diǎn)確定的分割超平面的距離哪個(gè)大時(shí)會(huì)用到。 （2）創(chuàng)建kdTree代碼

//創(chuàng)建kd樹(shù)，keypoints為點(diǎn)數(shù)據(jù)，parent表示當(dāng)前樹(shù)的雙親，默認(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)建左子樹(shù)tmp->right = createTree(rightKeyPoints, tmp);//遞歸調(diào)用，創(chuàng)建右子樹(shù)tmp->parent = parent;//記錄父節(jié)點(diǎn)return tmp;//返回當(dāng)前樹(shù)的根節(jié)點(diǎn)}

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

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

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

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

根據(jù)前面k-d樹(shù)的搜索過(guò)程我們可以知道，在搜索時(shí)首先沿著kd樹(shù)找到葉子節(jié)點(diǎn)，然后依次回溯，而回溯的路程就是前面我們查找葉子節(jié)點(diǎn)時(shí)逆序，因此進(jìn)行回溯時(shí)并沒(méi)有利用這些點(diǎn)的信息。我們接下來(lái)介紹的算法就是利用這些信息，回溯時(shí)給各個(gè)需要回溯的結(jié)點(diǎn)以?xún)?yōu)先級(jí)，這樣找到最近鄰會(huì)更快。接下來(lái)詳細(xì)介紹bbf算法的流程。 其實(shí)bbf算法的思想比較簡(jiǎn)單，通過(guò)對(duì)回溯可能需要的路過(guò)的結(jié)點(diǎn)加入隊(duì)列，并按照查找點(diǎn)到該結(jié)點(diǎn)確定的超平面的距離進(jìn)行排序，然后每次首先遍歷的是優(yōu)先級(jí)最高（即距離最短的結(jié)點(diǎn)），直到隊(duì)列為空算法結(jié)束。同時(shí)bbf算法也設(shè)立了一個(gè)時(shí)間限制，如果算法運(yùn)行時(shí)間超過(guò)該限制，不管是不是為空，一律停止運(yùn)行，返回當(dāng)前的最近鄰點(diǎn)作為結(jié)果。 bbf的算法流程如下： 輸入：kd樹(shù)，查找點(diǎn)x 輸出：kd樹(shù)種距離查找點(diǎn)最近的點(diǎn)以及最近的距離 流程：（1）若kd樹(shù)為空，則設(shè)定兩者距離為無(wú)窮大，返回；如果kd樹(shù)非空，則將kd樹(shù)的根節(jié)點(diǎn)加入到優(yōu)先級(jí)隊(duì)列中；             （2）從優(yōu)先級(jí)隊(duì)列中出隊(duì)當(dāng)前優(yōu)先級(jí)最大的結(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)先級(jí)隊(duì)列中出隊(duì)優(yōu)先級(jí)最大的結(jié)點(diǎn)；             （3）重復(fù)（1）和（2）中的操作，直到優(yōu)先級(jí)隊(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)先級(jí)隊(duì)列，查找點(diǎn)到當(dāng)前點(diǎn)確定的分割超平面距離越小優(yōu)先級(jí)越大priority_queuepri_queue;//priorityInfo類(lèi)型包含了如下信息：//(1)當(dāng)前的結(jié)點(diǎn)指針，指向kdNode類(lèi)型//(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;//這里沒(méi)有記錄時(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)確定的分割超平面的距離則不訪(fǎng)問(wè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)為根的子樹(shù)，直到葉子節(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ù)超過(guò)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è)的太低會(huì)造成匹配結(jié)果較差，需要在效率和正確率上進(jìn)行取舍。