演算法一 knn 擴充套件 BBF演算法，在KD-tree上找KNN ( K-nearest neighbor)
Step1: BBF演算法，在KD-tree上找KNN。第一步做匹配咯~
1. 什麼是KD-tree（from wiki）
K-Dimension tree，實際上是一棵平衡二叉樹。
一般的KD-tree構造過程：
function kdtree (list of points pointList, int depth)
{
if pointList is empty
return nil;
else {
// Select axis based on depth so that axis cycles through all valid values
var int axis := depth mod k;
// Sort point list and choose median as pivot element
select median by axis from pointList;
// Create node and construct subtrees
var tree_node node;
node.location := median;
node.leftChild := kdtree(points in pointList before median, depth+1);
node.rightChild := kdtree(points in pointList after median, depth+1);
return node;
}
}
【例】pointList = [(2,3), (5,4), (9,6), (4,7), (8,1), (7,2)] tree = kdtree(pointList)
2. BBF演算法，在KD-tree上找KNN ( K-nearest neighbor)
BBF(Best Bin First)演算法，藉助優先佇列（這裡用最小堆）實現。從根開始，在KD-tree上找路子的時候，錯過的點先塞到優先佇列裡，自己先一個勁兒掃到leaf；然後再從佇列裡取出目前key值最小的（這裡是是ki維上的距離最小者），重複上述過程，一個勁兒掃到leaf；直到佇列找空了，或者已經重複了200遍了停止。
Step1: 將img2的features建KD-tree; kd_root = kdtree_build( feat2, n2 );。在這裡，ki是選取均方差最大的那個維度，kv是各特徵點在那個維度上的median值，features是你率領的整個兒子孫子特徵大軍，n是你兒子孫子個數。
|
/** a node in a k-d tree */ struct kd_node{ int ki; /**< partition key index */ double kv; /**< partition key value */ int leaf; /**< 1 if node is a leaf, 0 otherwise */ struct feature* features; /**< features at this node */ int n; /**< number of features */ struct kd_node* kd_left; /**< left child */ struct kd_node* kd_right; /**< right child */ }; |
Step2: 將img1的每個feat到KD-tree裡找k個最近鄰，這裡k=2。
k = kdtree_bbf_knn( kd_root, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS );
|
min_pq = minpq_init(); minpq_insert( min_pq, kd_root, 0 ); while( min_pq->n > 0 && t < max_nn_chks ) //佇列裡有東西就繼續搜，同時控制在t<200（即200步內） { expl = (struct kd_node*)minpq_extract_min( min_pq ); //取出最小的，front & pop expl = explore_to_leaf( expl, feat, min_pq ); //從該點開始，explore到leaf，路過的“有意義的點”就塞到最小佇列min_pq中。 for( i = 0; i < expl->n; i++ ) // { tree_feat = &expl->features[i]; bbf_data->old_data = tree_feat->feature_data; bbf_data->d = descr_dist_sq(feat, tree_feat); //兩feat均方差 tree_feat->feature_data = bbf_data; n += insert_into_nbr_array( tree_feat, _nbrs, n, k ); //按從小到大塞到neighbor數組裡，到時候取前k個就是 KNN 咯~ n 每次加1或0，表示目前已有的元素個數 } t++; } |
對“有意義的點”的解釋：
|
struct kd_node* explore_to_leaf( struct kd_node* kd_node, struct feature* feat, struct min_pq* min_pq )//expl, feat, min_pq { struct kd_node* unexpl, * expl = kd_node; double kv; int ki; while( expl && ! expl->leaf ) { ki = expl->ki; kv = expl->kv; if( feat->descr[ki] <= kv ) { unexpl = expl->kd_right; expl = expl->kd_left; //走左邊，右邊點將被記下來 } else { unexpl = expl->kd_left; expl = expl->kd_right; //走右邊，左邊點將被記下來 } minpq_insert( min_pq, unexpl, ABS( kv - feat->descr[ki] ) ) ;//將這些點插入進來，key鍵值為|kv - feat->descr[ki]| 即第ki維上的差值 } return expl; } |
Step3: 如果k近鄰找到了(k=2)，那麼判斷是否能作為有效特徵，d0/d1<0.49就算是咯~
|
d0 = descr_dist_sq( feat, nbrs[0] );//計算兩特徵間squared Euclidian distance d1 = descr_dist_sq( feat, nbrs[1] ); if( d0 < d1 * NN_SQ_DIST_RATIO_THR )//如果d0/d1小於閾值0.49 { pt1 = cvPoint( cvRound( feat->x ), cvRound( feat->y ) ); pt2 = cvPoint( cvRound( nbrs[0]->x ), cvRound( nbrs[0]->y ) ); pt2.y += img1->height; cvLine( stacked, pt1, pt2, CV_RGB(255,0,255), 1, 8, 0 );//畫線 m++;//matches個數 feat1[i].fwd_match = nbrs[0]; } |
Step2: 通過RANSAC演算法來消除錯配，什麼是RANSAC先？
1. RANSAC (Random Sample Consensus, 隨機抽樣一致) (from wiki)
該演算法做什麼呢？呵呵，用一堆資料去搞定一個待定模型，這裡所謂的搞定就是一反覆測試、迭代的過程，找出一個error最小的模型及其對應的同盟軍（consensus set）。用在我們的SIFT特徵匹配裡，就是說找一個變換矩陣出來，使得儘量多的特徵點間都符合這個變換關係。
演算法思想：
input:
data - a set of observations
model - a model that can be fitted to data
n - the minimum number of data required to fit the model
k - the maximum number of iterations allowed in the algorithm
t - a threshold value for determining when a datum fits a model
d - the number of close data values required to assert that a model fits well to data
output:
best_model - model parameters which best fit the data (or nil if no good model is found)
best_consensus_set - data point from which this model has been estimated
best_error - the error of this model relative to the data
iterations := 0
best_model := nil
best_consensus_set := nil
best_error := infinity
while iterations < k //進行K次迭代
maybe_inliers := n randomly selected values from data
maybe_model := model parameters fitted to maybe_inliers
consensus_set := maybe_inliers
for every point in data not in maybe_inliers
if point fits maybe_model with an error smaller than t //錯誤小於閾值t
add point to consensus_set //成為同盟，加入consensus set
if the number of elements in consensus_set is > d //同盟軍已經大於d個人，夠了
(this implies that we may have found a good model,
now test how good it is)
better_model := model parameters fitted to all points in consensus_set
this_error := a measure of how well better_model fits these points
if this_error < best_error
(we have found a model which is better than any of the previous ones,
keep it until a better one is found)
best_model := better_model
best_consensus_set := consensus_set
best_error := this_error
increment iterations
return best_model, best_consensus_set, best_error
2. RANSAC去除錯配：
H = ransac_xform( feat1, n1, FEATURE_FWD_MATCH, lsq_homog, 4, 0.01,homog_xfer_err, 3.0, NULL, NULL );
|
nm = get_matched_features( features, n, mtype, &matched ); /* initialize random number generator */ rng = gsl_rng_alloc( gsl_rng_mt19937 ); gsl_rng_set( rng, time(NULL) ); in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform ); //符合這一要求的內點至少得有多少個 p = pow( 1.0 - pow( in_frac, m ), k ); i = 0; while( p > p_badxform )//p>0.01 { sample = draw_ransac_sample( matched, nm, m, rng ); extract_corresp_pts( sample, m, mtype, &pts, &mpts ); M = xform_fn( pts, mpts, m ); if( ! M ) goto iteration_end; in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus); if( in > in_max ) { if( consensus_max ) free( consensus_max ); consensus_max = consensus; in_max = in; in_frac = (double)in_max / nm; } else free( consensus ); cvReleaseMat( &M ); iteration_end: release_mem( pts, mpts, sample ); p = pow( 1.0 - pow( in_frac, m ), ++k ); } /* calculate final transform based on best consensus set */ if( in_max >= in_min ) { extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts ); M = xform_fn( pts, mpts, in_max ); in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus); cvReleaseMat( &M ); release_mem( pts, mpts, consensus_max ); extract_corresp_pts( consensus, in, mtype, &pts, &mpts ); M = xform_fn( pts, mpts, in ); |
思考中的一些問題：
features間的對應關係，記錄在features->fwd_match裡（matching feature from forward
imge）。
1. 資料是nm個特徵點間的對應關係，由它們產生一個3*3變換矩陣（xform_fn = hsq_homog函式，此要>=4對的對應才可能計算出來咯~），此乃模型model。
2. 然後開始找同盟軍（find_consensus函式），判斷除了sample的其它對應關係是否滿足這個模型（err_fn = homog_xfer_err函式，<=err_tol就OK~），滿足則留下。
3. 一旦大於當前的in_max，那麼該模型就升級為目前最牛的模型。（最最原始的RANSAC是按錯誤率最小走的，我們這會兒已經保證了錯誤率在err_tol範圍內，按符合要求的對應數最大走，儘量多的特徵能匹配地上）
4. 重複以上3步，直到(1-wm)k <=p_badxform (即0.01)，模型就算找定~
5. 最後再把模型和同盟軍定一下，齊活兒~
宣告：以上程式碼參考Rob Hess的SIFT實現。