基於kd樹的knn的實現原理可以參考文末的鏈接,都是一些好文章。
這裡參考了別人的代碼。用c語言寫的包括kd樹的構建與查找k近鄰的程序。
code:
1 #include<stdio.h>
2 #include<stdlib.h>
3 #include<math.h>
4 #include<time.h>
5
6 typedef struct{//數據維度
7 double x;
8 double y;
9 }data_struct;
10
11 typedef struct kd_node{
12 data_struct split_data;//數據結點
13 int split;//分裂維
14 struct kd_node *left;//由位於該結點分割超面左子空間內所有數據點構成的kd-tree
15 struct kd_node *right;//由位於該結點分割超面右子空間內所有數據點構成的kd-tree
16 }kd_struct;
17
18 //用於排序
19 int cmp1( const void *a , const void *b )
20 {
21 return (*(data_struct *)a).x > (*(data_struct *)b).x ? 1:-1;
22 }
23 //用於排序
24 int cmp2( const void *a , const void *b )
25 {
26 return (*(data_struct *)a).y > (*(data_struct *)b).y ? 1:-1;
27 }
28 //計算分裂維和分裂結點
29 void choose_split(data_struct data_set[],int size,int dimension,int *split,data_struct *split_data)
30 {
31 int i;
32 data_struct *data_temp;
33 data_temp=(data_struct *)malloc(size*sizeof(data_struct));
34 for(i=0;i<size;i++)
35 data_temp[i]=data_set[i];
36 static int count=0;//設為靜態
37 *split=(count++)%dimension;//分裂維
38 if((*split)==0) qsort(data_temp,size,sizeof(data_temp[0]),cmp1);
39 else qsort(data_temp,size,sizeof(data_temp[0]),cmp2);
40 *split_data=data_temp[(size-1)/2];//分裂結點排在中位
41 }
42 //判斷兩個數據點是否相等
43 int equal(data_struct a,data_struct b){
44 if(a.x==b.x && a.y==b.y) return 1;
45 else return 0;
46 }
47 //建立KD樹
48 kd_struct *build_kdtree(data_struct data_set[],int size,int dimension,kd_struct *T)
49 {
50 if(size==0) return NULL;//遞歸出口
51 else{
52 int sizeleft=0,sizeright=0;
53 int i,split;
54 data_struct split_data;
55 choose_split(data_set,size,dimension,&split,&split_data);
56 data_struct data_right[size];
57 data_struct data_left[size];
58
59 if (split==0){//x維
60 for(i=0;i<size;++i){
61 if(!equal(data_set[i],split_data) && data_set[i].x <= split_data.x){//比分裂結點小
62 data_left[sizeleft].x=data_set[i].x;
63 data_left[sizeleft].y=data_set[i].y;
64 sizeleft++;//位於分裂結點的左子空間的結點數
65 }
66 else if(!equal(data_set[i],split_data) && data_set[i].x > split_data.x){//比分裂結點大
67 data_right[sizeright].x=data_set[i].x;
68 data_right[sizeright].y=data_set[i].y;
69 sizeright++;//位於分裂結點的右子空間的結點數
70 }
71 }
72 }
73 else{//y維
74 for(i=0;i<size;++i){
75 if(!equal(data_set[i],split_data) && data_set[i].y <= split_data.y){
76 data_left[sizeleft].x=data_set[i].x;
77 data_left[sizeleft].y=data_set[i].y;
78 sizeleft++;
79 }
80 else if (!equal(data_set[i],split_data) && data_set[i].y > split_data.y){
81 data_right[sizeright].x = data_set[i].x;
82 data_right[sizeright].y = data_set[i].y;
83 sizeright++;
84 }
85 }
86 }
87 T=(kd_struct *)malloc(sizeof(kd_struct));
88 T->split_data.x=split_data.x;
89 T->split_data.y=split_data.y;
90 T->split=split;
91 T->left=build_kdtree(data_left,sizeleft,dimension,T->left);//左子空間
92 T->right=build_kdtree(data_right,sizeright,dimension,T->right);//右子空間
93 return T;//返回指針
94 }
95 }
96 //計算歐氏距離
97 double compute_distance(data_struct a,data_struct b){
98 double tmp=pow(a.x-b.x,2.0)+pow(a.y-b.y,2.0);
99 return sqrt(tmp);
100 }
101 //搜索1近鄰
102 void search_nearest(kd_struct *T,int size,data_struct test,data_struct *nearest_point,double *distance)
103 {
104 int path_size;//搜索路徑內的指針數目
105 kd_struct *search_path[size];//搜索路徑保存各結點的指針
106 kd_struct* psearch=T;
107 data_struct nearest;//最近鄰的結點
108 double dist;//查詢結點與最近鄰結點的距離
109 search_path[0]=psearch;//初始化搜索路徑
110 path_size=1;
111 while(psearch->left!=NULL || psearch->right!=NULL){
112 if (psearch->split==0){
113 if(test.x <= psearch->split_data.x)//如果小於就進入左子樹
114 psearch=psearch->left;
115 else
116 psearch=psearch->right;
117 }
118 else{
119 if(test.y <= psearch->split_data.y)//如果小於就進入右子樹
120 psearch=psearch->left;
121 else
122 psearch=psearch->right;
123 }
124 search_path[path_size++]=psearch;//將經過的分裂結點保存在搜索路徑中
125 }
126 //取出search_path最後一個元素,即葉子結點賦給nearest
127 nearest.x=search_path[path_size-1]->split_data.x;
128 nearest.y=search_path[path_size-1]->split_data.y;
129 path_size--;//search_path的指針數減一
130 dist=compute_distance(nearest,test);//計算與該葉子結點的距離作為初始距離
131
132 //回溯搜索路徑
133 kd_struct* pback;
134 while(path_size!=0){
135 pback=search_path[path_size-1];//取出search_path最後一個結點賦給pback
136 path_size--;//search_path的指針數減一
137
138 if(pback->left==NULL && pback->right==NULL){//如果pback為葉子結點
139 if(dist>compute_distance(pback->split_data,test)){
140 nearest=pback->split_data;
141 dist=compute_distance(pback->split_data,test);
142 }
143 }
144 else{//如果pback為分裂結點
145 int s=pback->split;
146 if(s==0){//x維
147 if(fabs(pback->split_data.x-test.x)<dist){//若以查詢點為中心的圓(球或超球),半徑為dist的圓與分割超平面相交,那麼就要跳到另一邊的子空間去搜索
148 if(dist>compute_distance(pback->split_data,test)){
149 nearest=pback->split_data;
150 dist=compute_distance(pback->split_data, test);
151 }
152 if(test.x<=pback->split_data.x)//若查詢點位於pback的左子空間,那麼就要跳到右子空間去搜索
153 psearch=pback->right;
154 else
155 psearch=pback->left;//若以查詢點位於pback的右子空間,那麼就要跳到左子空間去搜索
156 if(psearch!=NULL)
157 search_path[path_size++]=psearch;//psearch加入到search_path中
158 }
159 }
160 else {//y維
161 if(fabs(pback->split_data.y-test.y)<dist){//若以查詢點為中心的圓(球或超球),半徑為dist的圓與分割超平面相交,那麼就要跳到另一邊的子空間去搜索
162 if(dist>compute_distance(pback->split_data,test)){
163 nearest=pback->split_data;
164 dist=compute_distance(pback->split_data,test);
165 }
166 if(test.y<=pback->split_data.y)//若查詢點位於pback的左子空間,那麼就要跳到右子空間去搜索
167 psearch=pback->right;
168 else
169 psearch=pback->left;//若查詢點位於pback的的右子空間,那麼就要跳到左子空間去搜索
170 if(psearch!=NULL)
171 search_path[path_size++]=psearch;//psearch加入到search_path中
172 }
173 }
174 }
175 }
176
177 (*nearest_point).x=nearest.x;//最近鄰
178 (*nearest_point).y=nearest.y;
179 *distance=dist;//距離
180 }
181
182 int main()
183 {
184 int n=6;//數據個數
185 data_struct nearest_point;
186 double distance;
187 kd_struct *root=NULL;
188 data_struct data_set[6]={{2,3},{5,4},{9,6},{4,7},{8,1},{7,2}};//數據集
189 data_struct test={7.1,2.1};//查詢點
190 root=build_kdtree(data_set,n,2,root);
191
192 search_nearest(root,n,test,&nearest_point,&distance);
193 printf("nearest neighbor:(%.2f,%.2f)\ndistance:%.2f \n",nearest_point.x,nearest_point.y,distance);
194 return 0;
195 }
196 /* x 5,4
197 / \
198 y 2,3 7.2
199 \ / \
200 x 4,7 8.1 9.6
201 */
參考:
https://www.joinquant.com/post/2627?f=study&m=math
https://www.joinquant.com/post/2843?f=study&m=math
http://blog.csdn.net/zhl30041839/article/details/9277807