本文转载自:http://blog.csdn.NET/luoshixian099/article/details/47606159
继上一篇中已经介绍了SIFT原理与C源码剖析,最后得到了一系列特征点,每个特征点对应一个128维向量。假如现在有两副图片都已经提取到特征点,现在要做的就是匹配上相似的特征点。
相似性查询有两种基本方式:1.范围查询:即给点查询点和查询阈值,从数据集中找出所有与查询点距离小于阈值的点。
2.K近邻查询:给点查询点及正整数K,从数据集中找到与查询点最近的K个数据,当K=1时,就是最近邻查询。
特征匹配算子可以分为两类:1.穷举法:即将数据集中的点与查询点逐一计算距离,如果图1提取到N1个特征点,图2中提取到N2个特征点,用穷举法匹配,要做N1×N2运算,这种方法显然效率低下。
2.建立数据索引:对数据进行分析,对搜索数据空间进行划分,按划分时是否有重叠,分为KD树和R树。KD树是对空间划分时没有重叠的一种。
一个三维k-d树。第一次划分(红色)把根节点(白色)划分成两个节点,然后它们分别再次被划分(绿色)为两个子节点。最后这四个子节点的每一个都被划分(蓝色)为两个子节点。因为没有更进一步的划分,最后得到的八个节点称为叶子节点。
KD树的构建:KD树是一个二叉树,对数据空间空间进行划分,每一个结点对应一个空间范围。如上图所示,三维空间的划分方式。首先确定在数据集上对应方差最大的维度ki,并找到在ki维度上的数据集的中值kv(并作为根节点),即第一步把空间划分成两部分,在第ki维上小于kv的为一部分称为左子节点,大于kv的为另外一部分对应右子节点,,然后再利用同样的方法,对左子结点和右子节点继续构建二叉树,直所剩数据集为空。
举个例子:有5个数据,每个数据都是5维,建立KD树,A<7,5,7,3,8>;B<3,4,1,2,7>;C<5,2,6,6,9>;D<9,3,2,4,1>,E<2,1,5,1,4>,首先在计算在5个维度上的方差为6.56;2;5.36;2.96;8.56;可见在第5维度上方差最大,继续在第5个维度上找到中值为7,即B点,在第5维度上值小于7的作为左子树数据(A,C),大于7的作为右子树(D,E),然后继续在A,C,两点上计算方差最大的维度,继续划分。D,E也是如此。如下图,ki表示维度,kv表示该维度上的值。
KD树的查询:从根节点开始沿二叉树搜索,直到叶子结点为止,此时该叶节点并不一定是最近的点,但是一定是叶子结点附近的点。所以一定要有回溯的过程,回溯到根节点为止。例如:查询与M<5,4,1,3,6>点的最近邻点,查询路径为B,A,C,计算完MC的距离后,逆序向上,查询A及A的右子树,再次回溯B及B左子树,最后得到最近的距离,MB点最近。
假如数据集是维数是D,一般来说要求数据的规模N需要满足N>>2^D条件,才能达到高效的搜索,一般来说用标准的KD树时数据集的维数不超过20,但是像SIFT特征描述子128为,SURF描述子为64维,所以要对现有的KD树进行改进。
BBF:上述回溯的过程,完全是按照查询时路径决定的,没有考虑查询路径上的数据性质,BBF(Best-Bin-First)查询机制能确保优先包含最近邻点的空间,即BBF维护了一个优先队列,每一次查询到左子树或右子树的过程中,同时计算查询点在该维度的中值的距离差保存在优先队列里,同时另一个孩子节点地址也存入队列里,回溯的过程即从优先队列按(差值)从小到大的顺序依次回溯。如上一个例子,首先把B保存在优先队列里,然后开始从优先队列里取数据,取出B,发现要到左孩子A节点里继续查询,这时,要把右孩子节点D保存在优先队列里,同时加上距离属性ki=5,kv=7,所以d=7-6=1,这时优先队列里简记为D(1);同理,如果A有右孩子,也要存入优先队列,加上属性ki=2,kv=5,d=5-4=1;(例子不太恰当,o(╯□╰)o),回溯的过程是按照优先队列的距离逐个回溯,直到优先队列为空,或者超时,停止;BBF设置了超时机制,为了在高维数据上,满足检索速度的需要以精度换取时间,获得快速查询。这样可知,BBF机制找到的最近邻是近似的,并非是最近的,只能说是离最近点比较近而已。超时机制在算法的实现上,限定了从优先队列中提取数据的次数。
下面从算法上解析:
构建KD树:
- struct kd_node* kdtree_build( struct feature* features, int n )
- {
- struct kd_node* kd_root;
-
- if( ! features || n <= 0 )
- {
- fprintf( stderr, "Warning: kdtree_build(): no features, %s, line %d\n",
- __FILE__, __LINE__ );
- return NULL;
- }
-
- kd_root = kd_node_init( features, n );
- expand_kd_node_subtree( kd_root );
-
- return kd_root;
- }
- static struct kd_node* kd_node_init( struct feature* features, int n )
- {
- struct kd_node* kd_node;
-
- kd_node = malloc( sizeof( struct kd_node ) );
- memset( kd_node, 0, sizeof( struct kd_node ) );
- kd_node->ki = -1;
- kd_node->features = features;
- kd_node->n = n;
-
- return kd_node;
- }
扩展KD树:以当前结点的最大方差的维数为对应的中值为基准,把所有数据分成左右子树的结点数据,并以此递归下去,直到叶子结点的创建即返回。
- static void expand_kd_node_subtree( struct kd_node* kd_node )
- {
-
- if( kd_node->n == 1 || kd_node->n == 0 )
- {
- kd_node->leaf = 1;
- return;
- }
-
- assign_part_key( kd_node );
- partition_features( kd_node );
-
- if( kd_node->kd_left )
- expand_kd_node_subtree( kd_node->kd_left );
- if( kd_node->kd_right )
- expand_kd_node_subtree( kd_node->kd_right );
- }
计算最大方差对应的维数ki,与中值kv,取中值时,采用了最坏情况也是线性时间的选择算法,我的博客之前写过,这里不再分析<中位数排序>
点击打开链接
- static void assign_part_key( struct kd_node* kd_node )
- {
- struct feature* features;
- double kv, x, mean, var, var_max = 0;
- double* tmp;
- int d, n, i, j, ki = 0;
-
- features = kd_node->features;
- n = kd_node->n;
- d = features[0].d;
-
-
- for( j = 0; j < d; j++ )
- {
- mean = var = 0;
- for( i = 0; i < n; i++ )
- mean += features[i].descr[j];
- mean /= n;
- for( i = 0; i < n; i++ )
- {
- x = features[i].descr[j] - mean;
- var += x * x;
- }
- var /= n;
-
- if( var > var_max )
- {
- ki = j;
- var_max = var;
- }
- }
-
-
- tmp = calloc( n, sizeof( double ) );
- for( i = 0; i < n; i++ )
- tmp[i] = features[i].descr[ki];
- kv = median_select( tmp, n );
- free( tmp );
-
- kd_node->ki = ki;
- kd_node->kv = kv;
- }
按ki维上kv值,把特征点排序,小于等于kv为作为左子树数据,大于kv作为右子树数据
- static void partition_features( struct kd_node* kd_node )
- {
- struct feature* features, tmp;
- double kv;
- int n, ki, p, i, j = -1;
-
- features = kd_node->features;
- n = kd_node->n;
- ki = kd_node->ki;
- kv = kd_node->kv;
- for( i = 0; i < n; i++ )
- if( features[i].descr[ki] <= kv )
- {
- tmp = features[++j];
- features[j] = features[i];
- features[i] = tmp;
- if( features[j].descr[ki] == kv )
- p = j;
- }
- tmp = features[p];
- features[p] = features[j];
- features[j] = tmp;
-
-
- if( j == n - 1 )
- {
- kd_node->leaf = 1;
- return;
- }
-
- kd_node->kd_left = kd_node_init( features, j + 1 );
- kd_node->kd_right = kd_node_init( features + ( j + 1 ), ( n - j - 1 ) );
- }
☆KD树已经创建完毕,现在要做的是查询,查询与特征点最近邻的K个特征点,首先把根节点插入到优先队列,然后开始从有优先队列中取元素,遍历到叶节点,同时路径过程中,未查询的另一个结点的加入优先队列(按ki维上的数值与kv的差值的绝对值大小),然后再次从优先队列中取结点,再次遍历到叶节点,如此反复...直到遇到超时限制,或者遍历完所有节点为止。
-
-
-
-
-
-
-
- int kdtree_bbf_knn( struct kd_node* kd_root, struct feature* feat, int k,
- struct feature*** nbrs, int max_nn_chks )
- {
- struct kd_node* expl;
- struct min_pq* min_pq;
- struct feature* tree_feat, ** _nbrs;
- struct bbf_data* bbf_data;
- int i, t = 0, n = 0;
-
- if( ! nbrs || ! feat || ! kd_root )
- {
- fprintf( stderr, "Warning: NULL pointer error, %s, line %d\n",
- __FILE__, __LINE__ );
- return -1;
- }
-
- _nbrs = calloc( k, sizeof( struct feature* ) );
- min_pq = minpq_init();
- minpq_insert( min_pq, kd_root, 0 );
- while( min_pq->n > 0 && t < max_nn_chks )
- {
- expl = (struct kd_node*)minpq_extract_min( min_pq );
- if( ! expl )
- {
- fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n",
- __FILE__, __LINE__ );
- goto fail;
- }
-
- expl = explore_to_leaf( expl, feat, min_pq );
- if( ! expl )
- {
- fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n",
- __FILE__, __LINE__ );
- goto fail;
- }
-
- for( i = 0; i < expl->n; i++ )
- {
-
- tree_feat = &expl->features[i];
- bbf_data = malloc( sizeof( struct bbf_data ) );
- if( ! bbf_data )
- {
- fprintf( stderr, "Warning: unable to allocate memory,"
- " %s line %d\n", __FILE__, __LINE__ );
- goto fail;
- }
-
- bbf_data->old_data = tree_feat->feature_data;
- printf("%x",bbf_data->old_data);
- bbf_data->d = descr_dist_sq(feat, tree_feat);
- tree_feat->feature_data = bbf_data;
- n += insert_into_nbr_array( tree_feat, _nbrs, n, k );
- }
- t++;
- }
-
- minpq_release( &min_pq );
- for( i = 0; i < n; i++ )
- {
- bbf_data = _nbrs[i]->feature_data;
- _nbrs[i]->feature_data = bbf_data->old_data;
- free( bbf_data );
- }
- *nbrs = _nbrs;
- return n;
-
- fail:
- minpq_release( &min_pq );
- for( i = 0; i < n; i++ )
- {
- bbf_data = _nbrs[i]->feature_data;
- _nbrs[i]->feature_data = bbf_data->old_data;
- free( bbf_data );
- }
- free( _nbrs );
- *nbrs = NULL;
- return -1;
- }
这里创建优先队列采用了堆排序的思想,堆排序一个重要的应用就是优先队列,在一个包含n个元素的堆中,所有优先队列的操作都可以在lgn时间内完成。优先队列也有两种形式,最大优先队列和最小优先队列,这里使用的是最小优先队列,即key值越小优先级越高,关于堆排序的原理,可以看之前的堆排序的算法<堆排序>
点击打开链接
。
- struct min_pq* minpq_init()
- {
- struct min_pq* min_pq;
-
- min_pq = malloc( sizeof( struct min_pq ) );
- min_pq->pq_array = calloc( MINPQ_INIT_NALLOCD, sizeof( struct pq_node ) );
- min_pq->nallocd = MINPQ_INIT_NALLOCD;
- min_pq->n = 0;
-
- return min_pq;
- }
- int minpq_insert( struct min_pq* min_pq, void* data, int key )
- {
- int n = min_pq->n;
-
-
- if( min_pq->nallocd == n )
- {
- min_pq->nallocd = array_double( (void**)&min_pq->pq_array,
- min_pq->nallocd,
- sizeof( struct pq_node ) );
- if( ! min_pq->nallocd )
- {
- fprintf( stderr, "Warning: unable to allocate memory, %s, line %d\n",
- __FILE__, __LINE__ );
- return 1;
- }
- }
- min_pq->pq_array[n].data = data;
- min_pq->pq_array[n].key = INT_MAX;
- decrease_pq_node_key( min_pq->pq_array, min_pq->n, key );
- min_pq->n++;
- return 0;
- }
从队列中取出一个节点,沿着结点遍历到叶节点为止,同时未查询的加入优先队列。
- static struct kd_node* explore_to_leaf( struct kd_node* kd_node,
- struct feature* feat,
- struct min_pq* min_pq )
- {
- struct kd_node* unexpl, * expl = kd_node;
- double kv;
- int ki;
-
- while( expl && ! expl->leaf )
- {
- ki = expl->ki;
- kv = expl->kv;
-
- if( ki >= feat->d )
- {
- fprintf( stderr, "Warning: comparing imcompatible descriptors, %s" \
- " line %d\n", __FILE__, __LINE__ );
- return NULL;
- }
- if( feat->descr[ki] <= kv )
- {
- unexpl = expl->kd_right;
- expl = expl->kd_left;
- }
- else
- {
- unexpl = expl->kd_left;
- expl = expl->kd_right;
- }
-
- if( minpq_insert( min_pq, unexpl, ABS( kv - feat->descr[ki] ) ) )
- {
- fprintf( stderr, "Warning: unable to insert into PQ, %s, line %d\n",
- __FILE__, __LINE__ );
- return NULL;
- }
- }
-
- return expl;
- }
计算两特征点的欧式距离:
- double descr_dist_sq( struct feature* f1, struct feature* f2 )
- {
- double diff, dsq = 0;
- double* descr1, * descr2;
- int i, d;
-
- d = f1->d;
- if( f2->d != d )
- return DBL_MAX;
- descr1 = f1->descr;
- descr2 = f2->descr;
-
- for( i = 0; i < d; i++ )
- {
- diff = descr1[i] - descr2[i];
- dsq += diff*diff;
- }
- return dsq;
- }
找到的特征点距离插入nbrs队列,作为输出。nbrs也是按照从小到大的顺序存储距离,
假如要插入的新的距离为D,队列此时最后一个元素d
1.如果队列未满的情况(即要取得K近邻还没有找到K个)(1) D>=d,直接插入到队列后面;(2) D<d,则找到要插入的位置,然后比D大的元素都往后移动一位;
2.如果队列已经满了:(1)D>=d,直接丢弃当前的距离;(2)D<d,则找到要插入的位置,然后比D大的元素都往后移动一位;最后一个元素丢弃。
- static int insert_into_nbr_array( struct feature* feat, struct feature** nbrs,
- int n, int k )
- {
- struct bbf_data* fdata, * ndata;
- double dn, df;
- int i, ret = 0;
-
- if( n == 0 )
- {
- nbrs[0] = feat;
- return 1;
- }
-
-
- fdata = (struct bbf_data*)feat->feature_data;
- df = fdata->d;
- ndata = (struct bbf_data*)nbrs[n-1]->feature_data;
- dn = ndata->d;
- if( df >= dn )
- {
- if( n == k )
- {
- feat->feature_data = fdata->old_data;
- free( fdata );
- return 0;
- }
- nbrs[n] = feat;
- return 1;
- }
-
-
-
-
- if( n < k )
- {
- nbrs[n] = nbrs[n-1];
- ret = 1;
- }
- else
- {
- nbrs[n-1]->feature_data = ndata->old_data;
- free( ndata );
- }
- i = n-2;
- while( i >= 0 )
- {
- ndata = (struct bbf_data*)nbrs[i]->feature_data;
- dn = ndata->d;
- if( dn <= df )
- break;
- nbrs[i+1] = nbrs[i];
- i--;
- }
- i++;
- nbrs[i] = feat;
-
- return ret;
- }