VLFeat库实现KD-Tree算法
VLFeat库实现KD-Tree算法
K-D树(K-Dimensional Tree,即K维二叉树),K=1时,即是一棵普通的二叉树。常被用于高维空间中的搜索,比如范围搜索和最近邻搜索。考虑这样一种情况:
平面上有一堆散点,若想要找到与某一点最近的点,最朴素的办法就是逐一计算该点到其他点的距离,但算法会随着样本点的增大变得非常低效。使用2-D Tree能够比暴力搜索更快,并且随着散点集维度增加,KDTree的优势将更加明显。
本文档将介绍使用开源库VLFeat实现KDTree算法。
1 VLFeat库简介及配置
VLFeat,知名度略低于著名的计算机视觉/图像处理开源库OpenCV,VLFeat 开源库实现了大量的机器视觉算法,专门用于图像理解和局部特征提取和匹配。 算法包括 Fisher Vector、VLAD、SIFT、MSER、K-Means、分层 K-Means等。其使用 C 语言编写,提供 C 语言和 MATLAB 两种接口和详细的文档,可以在 windows,Mac 和 Linux 上使用。
首先,我们需要进入官网下载(VLFeat - 首页),并解压到本地。
然后,新建一个工程文件,随后开始配置:
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705173457176-1232267968.png)
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705173634663-1353715323.png)
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705173725794-81111299.png)
输入的命令行为:
copy D:\xxxxxx\vlfeat-0.9.21\bin\win64\vl.dll $(SolutionDir)x64\$(ConfigurationName)
这样就配置完成了,我们在main.cpp中输入测试代码:
#include "vl\generic.h"
int main() {
VL_PRINT("Hello world!\n");
system("pause");
return 0;
}
输出为:
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705173809412-262532049.png)
2 算法实现(Matlab)
(1)生成随机散点云、定义目标点
简单地说就是:在空间中有一堆散点,而我们要在茫茫多的散点中找到距离目标点(*)最近的几个。
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705173916408-1791018160.png)
rand('state',0) ;
Q = single(rand(2,1)) ; %目标查询点
X = single(rand(2, 100)) ;%生成随机散点
(2)建树
建立KDTree实际上是一个不断划分的过程,选择方差最大的维度进行切分,然后选择这一维度上的中点/中位数(median)或者均值(mean)作为分割点。因此,只要数据散点确定,这棵树也就被确定下来了。使用上述分裂规则,将当前空间内数据一分为二,再在新的子集中递归切分过程,直到分割到叶子节点。
图中展示了二维/三维数据点的建树结果(按红、绿、蓝的顺序进行切分),这样一来,KD树将平面/空间划分成了许多子平面/空间,随后我们便能根据待查询的点、左/右子空间、以及分割点计算距离,从而规避逐一计算两点距离的庞大计算量。
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705173948686-1623868865.png)
kdforest = vl_kdtreebuild (X, 'verbose','thresholdmethod', 'median') ;
(3)搜索&回溯
比如我们的目标点是(2,5),寻找它的最近邻点,得到的搜索路径为(7,2)\(\to\) (5,4)\(\to\)(4,7),因此我们得到的最近邻点为(4,7)。但是很明显,(2,3)这个点要更近,不过它被绿色这条分割线划到了对面(图①)。因此仅靠二叉树分割区域,然后把区域的上级分割点作为最邻近点是不准确的,需要通过回溯进行排查。需要注意,只有叶子节点才能作为最终的临近点去计算距离,也即搜索路径必须从根节点\(\to\)叶子节点,回溯本质是逆向\(\to\)正向的调校路径过程。
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705180021028-144652854.png)
%回溯/查询
[i, d] = vl_kdtreequery (kdforest, X, Q, 'numneighbors', 10, 'verbose') ;
% 输出i为样本点的index,d为优选点距目标点的距离distance
在整个过程完成之后,我们能够得到搜索结果。官方提供的教程中,对不同数据类型和分割规则作了演示:
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705174130148-1928781613.png)
3 转换至C/C++
VFLeat虽提供了开源C代码,好比给了一整箱零件,但离掌握函数的调用规则、搭建起一个完整的算法,还缺一张图纸。而这个图纸便是由上一章中由Matlab实现的部分。我们知道这套Matlab代码调用了开源库中的C零件,但是具体调用了哪些,怎么连接的,就需要用到MEX。
MEX搭建了Matlab与C之间的桥梁。它能够编译一个或多个用 MATLAB 数据 API 编写的 C++ 源文件并将其链接到当前文件夹中的二进制 MEX 文件中。通过MEX函数我们可以像调用内置函数一样,从 MATLAB 命令行调用 C 或 C++ 程序。比如,建树的函数vl_kdtreebuild:
-
Matlab "vlfeat-0.9.21\toolbox\misc\vl_kdtreebuild.m";
-
MEX "vlfeat-0.9.21\toolbox\mex\mexw64\vl_kdtreebuild.mexw64";
-
C "vlfeat-0.9.21\toolbox\misc\vl_kdtreebuild.c";
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705174203662-814578114.png)
可以简单的将MEX理解为一种C与Matlab都能使用的语言,它有其独特的语法。比如,MEX函数库里的mexFunction()函数,它相当于C语言中的main()函数;用mxArray表示数组。更多关于MEX的介绍和用法可访问https://ww2.mathworks.cn/help/matlab/call-mex-files-1.html。
原vl_kdtreebuild.c代码是用MEX书写:
void mexFunction(int nout, mxArray *out[], //输入及格式,
int nin, const mxArray *in[]) //输出及格式
{
enum {IN_DATA = 0, IN_END} ;
enum {OUT_TREE = 0} ;
int verbose = 0 ;
int opt ;
int next = IN_END ;
...
forest = vl_kdforest_new (dataType, dimension, numTrees, distance) ; //C函数,存在kdtree.c中
vl_kdforest_set_thresholding_method (forest, thresholdingMethod) ;
...
}
于是,我们可以根据MEX代码,根据其特定的数据格式与语法,转换成C/C++语言,从而仿制出一个C版的Demo。
*之后跳过MEX-->C/C++的转换、翻译过程,直接展示可使用的C/C++代码。
(1)生成随机散点云、定义目标点
在实际运用过程中,这部分通常是从外部输入的文件,此处我们直接用rand()函数生成10000组随机散点。并设置待搜寻的目标点为(30,50)。
#define NUM 10000
srand(7);
double *xpoint = (double*)calloc(NUM,sizeof(double));
double *ypoint = (double*)calloc(NUM, sizeof(double));
for(int i = 0; i<NUM; i++){
xpoint[i] = ((double)rand() / RAND_MAX) * 200.0 - 100.0; //生成随机二维散点
ypoint[i] = ((double)rand() / RAND_MAX) * 200.0 - 100.0;
}
double xdist = 30.0; //目标位置
double ydist = 50.0;
(2)建树
在mex代码中,提供了许多可调的参数。我们可以根据后续地实际需要增加判断语句去一一实现这些参数。首先,声明一个接受建树参数的结构体:
/*
vlmxOption options [] = { //带引号的是Matlab的用法,使用变量名的则是C
{"Verbose", 0, opt_verbose },
{"ThresholdMethod", 1, opt_threshold_method },
{"NumTrees", 1, opt_num_trees },
{"Distance", 1, opt_distance },
{0, 0, 0 }
} ;
*/
typedef struct BuildKdtreeParametert{
int num_trees;
}BuildKdtreeParametert;
随后借助Vl方法,创建可供C/C++调用的建树函数:
/**********************************************
函数功能:根据现有数据创建KD树(2D)
输入:
xpoint X坐标
ypoint Y坐标
nsample 样本个数
bkparameter 参数列表
输出:
forest VlKD树对象
备注:
函数中出现的所有以Vl为开头的关键字为VLFeat提供的类型或方法
**************************************************/
VlKDForest *kdtreebuild2D(double *xpoint, double *ypoint, unsigned int nsample, BuildKdtreeParametert bkparameter) {
if((xpoint==NULL)||(ypoint==NULL)||(nsample==0)){
return NULL;
}
vl_size numData = nsample; // 样本个数
vl_size Dimension = 2; // 样本维度
vl_size dataType = VL_TYPE_DOUBLE; // 数据类型
vl_size numTrees = 1; // 树棵数
VlVectorComparisonType distance = VlDistanceL2; //L2距离
int thresholdingMethod = VL_KDTREE_MEDIAN; //中位数作为分割点
double *data = NULL;
data=(double*)calloc(nsample*2,sizeof(double));
for(int i = 0; i < nsample; i++){ //两列数据合并为一列,为了满足后续接口
data[2 * i] = xpoint[i];
data[2 * i + 1] = ypoint[i];
}//end for(int i = 0; i < nsample; i++)
VlKDForest *forest;
forest = vl_kdforest_new(dataType, Dimension, numTrees, distance); //创建一颗新的KDTree
vl_kdforest_set_thresholding_method(forest, VL_KDTREE_MEDIAN); // 按照规则分割树
vl_kdforest_build(forest, numData, data); //建树
return forest;
}
(3)搜索&回溯
同样的,我们先声明参数列表,并创建一个最近邻属性结构体:
typedef struct FindKdtreeParameter{
int numneighbors;
}FindKdtreeParameter;
typedef struct NearstPointSet{ //Nearst
vector<double> xPoint;
vector<double> yPoint;
vector<double> distance;
int nSample;
}NearstPointSet;
移植搜索方法:
/**********************************************
函数功能:根据建好的KD树(2D)进行回溯,寻找最近的n个点
输入:
forest VlKD树对象,即上一步建好的树
xpoint X坐标
ypoint Y坐标
numpoint 样本个数
xdist 目标点X坐标
ydist 目标点Y坐标
findkdtreeparameter 参数列表
输出:
nps NearstPointSet对象,最近的n个点的属性
**************************************************/
NearstPointSet kdtreefind2D(VlKDForest *forest, double *xpoint, double *ypoint, int numpoint, double xdist, double ydist, FindKdtreeParameter findkdtreeparameter) {
NearstPointSet nps; //定义一些基础变量
VlKDForestNeighbor *neighbors = NULL;
double *qurey = NULL;
int numneighbors = findkdtreeparameter.numneighbors;
int indexTemp = 0;
double distanceTemp = 0.0;
qurey=(double*)calloc(2,sizeof(double));
neighbors = (VlKDForestNeighbor*)calloc(numneighbors,sizeof(VlKDForestNeighbor));
qurey[0] = xdist; //查询点(或目标点)
qurey[1] = ydist;
vl_kdforest_query(forest, neighbors, numneighbors, qurey); //执行搜索操作
for (int i = 0; i < numneighbors; i++){ //按指定需求输出最近的n个点
indexTemp = neighbors[i].index;
distanceTemp =neighbors[i].distance;
if((indexTemp>=0)&&(indexTemp<numpoint)){ //按索引存回坐标与距离
nps.xPoint.push_back(xpoint[indexTemp]);
nps.yPoint.push_back(ypoint[indexTemp]);
nps.distance.push_back(distanceTemp);
}//end if((indexTemp>=0)&&(indexTemp<numpoint)
}//end for (int i = 0; i < numneighbors; i++)
nps.nSample = nps.xPoint.size(); //搜索到的满足条件的点数
return nps;
}
(4)执行程序
BuildKdtreeParameter bkp;
bkp.num_trees = 1; //树棵树
VlKDForest *forest = kdtreebuild2D(xpoint, ypoint, NUM,bkp);
FindKdtreeParameter fkp;
fkp.numneighbors = 10; //挑选最近的10个点
NearstPointSet p1 = kdtreefind2D(forest, xpoint, ypoint, NUM, xdist, ydist,fkp);
//打印最近的那个点坐标和距离
printf("The coordinates of the nearest point are(%f,%f),and the distance=%f\n", p1.xPoint[0], p1.yPoint[0],p1.distance[0]);
将得到以下输出:
![](https://img2022.cnblogs.com/blog/1602810/202207/1602810-20220705174255283-23393948.png)
完整代码
I KD-Tree.h
#pragma once
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<algorithm>
#include<iostream>
#include<vector>
#define PI 3.141592654
#define EPS 0.000001
#define N 2
#include "vl\generic.h"
#include"vl\kdtree.h"
using namespace std;
typedef struct BuildKdtreeParameter{
int num_trees;
}BuildKdtreeParameter;
typedef struct FindKdtreeParameter{
int numneighbors;
}FindKdtreeParameter;
typedef struct NearstPointSet{ //Nearst
vector<double> xPoint;
vector<double> yPoint;
vector<double> distance;
int nSample;
}NearstPointSet;
VlKDForest *kdtreebuild2D(double *xpoint,double *ypoint,unsigned int nsample, BuildKdtreeParameter bkparameter); //build Kdtree Object
NearstPointSet kdtreefind2D(VlKDForest *forest,double *xpoint,double *ypoint,int numpoint,double xdist,double ydist, FindKdtreeParameter findkdtreeparameter);//estimate Nearest Point Set Of aim point ;
II KD-Tree.cpp
#include"KD-Tree.h"
//建树
VlKDForest *kdtreebuild2D(double *xpoint, double *ypoint, unsigned int nsample, BuildKdtreeParameter bkparameter) {
if((xpoint==NULL)||(ypoint==NULL)||(nsample==0)){
return NULL;
}
vl_size numData = nsample; // number of sample
vl_size Dimension = 2; // number of point dimension
vl_size dataType = VL_TYPE_DOUBLE; // Define size of data
vl_size numTrees = 1; // number of kd-tree
VlVectorComparisonType distance = VlDistanceL2; //L2 distance
int thresholdingMethod = VL_KDTREE_MEDIAN;
double *data = NULL;
data=(double*)calloc(nsample*2,sizeof(double));
for(int i = 0; i < nsample; i++){ //data to point
data[2 * i] = xpoint[i];
data[2 * i + 1] = ypoint[i];
}//end for(int i = 0; i < nsample; i++)
VlKDForest *forest;
forest = vl_kdforest_new(dataType, Dimension, numTrees, distance);
vl_kdforest_set_thresholding_method(forest, VL_KDTREE_MEDIAN);
vl_kdforest_build(forest, numData, data);
return forest;
}
//回溯
NearstPointSet kdtreefind2D(VlKDForest *forest, double *xpoint, double *ypoint, int numpoint, double xdist, double ydist, FindKdtreeParameter findkdtreeparameter) {
NearstPointSet nps;
VlKDForestNeighbor *neighbors = NULL;
double *qurey = NULL;
int numneighbors = findkdtreeparameter.numneighbors;
int indexTemp = 0;
double distanceTemp = 0.0;
qurey=(double*)calloc(2,sizeof(double));
neighbors = (VlKDForestNeighbor*)calloc(numneighbors,sizeof(VlKDForestNeighbor));
qurey[0] = xdist;
qurey[1] = ydist;
vl_kdforest_query(forest, neighbors, numneighbors, qurey);
for (int i = 0; i < numneighbors; i++){
indexTemp = neighbors[i].index;
distanceTemp =neighbors[i].distance;
if((indexTemp>=0)&&(indexTemp<numpoint)){
nps.xPoint.push_back(xpoint[indexTemp]);
nps.yPoint.push_back(ypoint[indexTemp]);
nps.distance.push_back(distanceTemp);
}//end if((indexTemp>=0)&&(indexTemp<numpoint)
}//end for (int i = 0; i < numneighbors; i++)
nps.nSample = nps.xPoint.size();
return nps;
}
III main.cpp
#include"KD-Tree.h"
#define NUM 10000
int main() {
srand(7);
double *xpoint = (double*)calloc(NUM,sizeof(double));
double *ypoint = (double*)calloc(NUM, sizeof(double));
for(int i = 0; i<NUM; i++){
xpoint[i] = ((double)rand() / RAND_MAX) * 200.0 - 100.0;
ypoint[i] = ((double)rand() / RAND_MAX) * 200.0 - 100.0;
}
double xdist = 30.0;
double ydist = 50.0;
BuildKdtreeParameter bkp;
bkp.num_trees = 10;
VlKDForest *forest = kdtreebuild2D(xpoint, ypoint, NUM, bkp);
FindKdtreeParameter fkp;
fkp.numneighbors = 10;
NearstPointSet p1 = kdtreefind2D(forest, xpoint, ypoint, NUM, xdist, ydist, fkp);
printf("The coordinates of the nearest point are(%f,%f),and the distance=%f\n", p1.xPoint[0], p1.yPoint[0],p1.distance[0]);
return 0;
}
参考: