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，因此更新最近鄰點和最近鄰距離。回溯到（7，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;
};//建立kd樹，keypoints為點資料，parent表示當前樹的雙親，預設為NULL
kdNode* kdTree::createTree(vector<Point> &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;
vector<Point> leftKeyPoints(keypoints.begin(), keypoints.begin() + sort_num);
vector<Point> rightKeyPoints(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;//返回當前樹的根節點
}//通過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);
}
else
cout << "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;
}上面的描述就是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_queue<priorityInfo> pri_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;
}
else
break;
}
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;
}
else
break;
}
//更新最近鄰
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;
}