[算法实现]基于分治的二维平面最近点对算法实现

摘要:

网上有很多关于分治方法求最近点对的讨论，但是没有完整的可运行代码，本文主要对于该问题介绍一完整的可运行代码，

供有兴趣者参考。

正文：

作为对比，我们也同时实现了最近点对的枚举解法，具体的主函数如下：

#include<stdio.h>#include<stdlib.h>#include "node.h"void initList(node* p){p[0].x= 2.0;p[0].y= 1.0;p[1].x= 1.0;p[1].y= 2.0;p[2].x= 1.2;p[2].y= 3.0;p[3].x= 3.0;p[3].y= 4.0;p[4].x= 5.0;p[4].y= 5.0;p[5].x= 1.5;p[5].y= 5.5;p[6].x= 2.5;p[6].y= 7.0;p[7].x= 3.5;p[7].y= 8.0;p[8].x= 4.0;p[8].y= 9.0;p[9].x = 3.9;p[9].y = 8.8;}//测试分治和暴力;int main(){node* p = (node*)malloc(sizeof(node)*10);initList(p);double ddf = force(p, 10); //9 is the number of elements;printf("the output of force is %lf\n", ddf);double dd = callMinDist(p, 10);printf("the output of divide-conquer is %lf\n", dd);getchar();return 0;}

上述中的force()是枚举的实现，callMinDist则是分治的实现，上述的initList主要对采用的测试案例进行初始化。

下面是node.h头文件的相关代码：

#ifndef __NODE__#define __NODE__#define SIZE 4#define MAX 100000.0typedef struct{double x;double y;}node;//排序部分;void mergeX(node a[], node b[], int s, int m, int t);void mergeY(node a[], node b[], int s, int m, int t);void mergeSortX(node a[], node b[], int s, int t);void mergeSortY(node a[], node b[], int s, int t);//utility;void show(node* a, int size);void initList(node* list);double dist(node* n1, node* n2);//枚举部分;double force(node* nodelist, int n);//分治部分;double combine(node* py, int n, double lx, double delta);void copynode(node* dt, node* sr, int b, int n);double minDist(node* px, node* py, int n);double callMinDist(node*p, int n);#endif

枚举部分的代码比较简单，具体看如下：

#include <math.h>#include "node.h"double dist(node* n1, node* n2){double temp1 = n1->x - n2->x;double temp2 = n1->y - n2->y;double sum = temp1*temp1 + temp2*temp2; //pow(temp1, 2)+pow(temp2, 2);return sqrt(sum);}double force(node* nodelist, int n) //n is number of elements{//这里d需要有一个初始的大值;double d = MAX;double t;//函数的主体就是下面这个双层循环;for(int i=0; i<n; i++){for(int j=i+1; j<n; j++){t = dist(&nodelist[i], &nodelist[j]);if(t<d)d = t;}}//这里最后返回d;return d;}

下面是对于分治算法的调用部分，调用之前需要分别将其中的点按x轴和按y轴进行排序操作，并且

将排完序的点放置在新的存储空间中；

double callMinDist(node* p, int n){node* px = (node*)malloc(n*sizeof(node)); //n主要是用于此处的空间申请;node* py = (node*)malloc(n*sizeof(node));//show(p, n);mergeSortX(p, px, 0, n-1); //按点的x轴值排序；copynode(px, p, 0, n-1);//show(px, n);//show(p, n);mergeSortY(p, py, 0, n-1); //按点的y轴值排序；copynode(py, p, 0, n-1);//show(py, n);double min = minDist(px, py, n);free(px);free(py);return min;}

下面就是分治算法的主体：

double minDist(node* px, node* py, int n)  //这里n是number of elements呢？还是max of index？根据下面的空间申请，应该是number of elements.{//printf("n is %d\n", n);if(n<=3){//show(px, n); //n is number of elements;double d = force(px, n); //n is number of elements;//printf("n=%d is %lf\n",n, d);return d;}int m=n/2;double fx = px[m].x;node* lx = (node*)malloc(m*sizeof(node));node* ly = (node*)malloc(m*sizeof(node));node* rx = (node*)malloc((n-m)*sizeof(node));node* ry = (node*)malloc((n-m)*sizeof(node));copynode(lx, px, 0, m-1); //对copy而言,这里的m应该是index;//show(lx, m);  //show中的n是number of elements;//printf("lx :%x\n", lx);copynode(rx, px, m, n-1); //copy这边是index;//show(rx, n-m);copynode(ly, py, 0, m-1);//show(ly, m);copynode(ry, py, m, n-1);//show(ry, n-m);double d1 = minDist(lx, ly, m); //m is number of elements;double dr = minDist(rx, ry, n-m);double delta = min(d1, dr);double d = combine(py, n, fx, delta); //对combine而言,这里的n是number of elements;//printf("lx :%x\n", lx);free(lx);free(ly);free(rx);free(ry);return min(delta, d);}

通过dl = minDist(lx, ly, m)来完成左半部分的计算；

dr = minDist(rx, ry, n-m)完成右半部分的计算；

最后通过combine(py, n, fx, delta)将两半部分的结果整合在一起；

这里的关键之处在于combine函数：

double combine(node* py, int n, double lx, double delta){int num; double d=MAX;double tempd;node* temp = (node*)malloc(n*sizeof(node));int j=0;for(int i=0; i<n; i++){if(fabs(py[i].x - lx)<= delta){ //求取在区间范围内的点；temp[j].x = py[i].x;temp[j].y = py[i].y;j++;}}num = j;  //temp中的元素for(int i=0; i<num; i++){for(j=i+1; j<(i+8) && (j<num); j++){tempd = dist(&temp[i], &temp[j]);if(tempd < d)d=tempd;}}free(temp);return d;}

根据书本上的分析，在区间中求取时，只需要计算当前点后(按y轴的值排序)的6到7

个点即可，因此此处的语句表现为：

for(j=i+1; j<(i+8)&&(j<num); j+)....

最后，我们来看下上述代码中用到的一些周边函数：

void copynode(node* dt, node* sr, int b, int n) //n is max of index;{int k=0;for(int i=b; i<=n; i++){dt[k].x = sr[i].x;dt[k].y = sr[i].y;k++;}}double min(double x, double y){if(x<=y)return x;return y;}

还有是通过归并排序对点集中的点进行排序的过程：

void mergeSortX(node a[], node b[], int s, int t){if(s == t){b[s].x = a[s].x;b[s].y = a[s].y;}else{int m = (s+t)/2;mergeSortX(a, b, s, m);mergeSortX(a, b, m+1, t);mergeX(a, b, s, m, t);}}void mergeSortY(node a[], node b[], int s, int t){if(s == t){b[s].x = a[s].x;b[s].y = a[s].y;}else{int m = (s+t)/2;mergeSortY(a, b, s, m);mergeSortY(a, b, m+1, t);mergeY(a, b, s, m, t);}}void mergeX(node* a, node* b, int s, int m, int t){int i, j, n;for(i=s, j=m+1, n=s; (i<=m)&&(j<=t); n++){if(b[i].x<=b[j].x){a[n].x = b[i].x;a[n].y = b[i].y;i++;}else{a[n].x = b[j].x;a[n].y = b[j].y;j++;}}while(i<=m){a[n].x = b[i].x;a[n++].y = b[i++].y;}while(j<=t){a[n].x = b[j].x;a[n++].y = b[j++].y;}//这里需要将a中的数据重新拷贝到b中;for(int i=s; i<=t; i++){b[i].x = a[i].x;b[i].y = a[i].y;}}void mergeY(node* a, node* b, int s, int m, int t){int i, j, n;for(i=s, j=m+1, n=s; (i<=m)&&(j<=t); n++){if(b[i].y<=b[j].y){a[n].x = b[i].x;a[n].y = b[i].y;i++;}else{a[n].x = b[j].x;a[n].y = b[j].y;j++;}}while(i<=m){a[n].x = b[i].x;a[n++].y = b[i++].y;}while(j<=t){a[n].x = b[j].x;a[n++].y = b[j++].y;}//这里需要将a中的数据重新拷贝到b中;for(int i=s; i<=t; i++){b[i].x = a[i].x;b[i].y = a[i].y;}}

注意：上述的归并排序函数写得不好，没必要用到a这个参数，完全可以函数内部堆上分配局部变量进行替换。另外代码中重复部分太多，应该可以写得更精简一点。

结论：

针对代码中的简单测试案例，分治案例结果正常；该算法的主要时间复杂度在于排序部分，复杂度为O(nlogn)，而枚举版本的复杂度为O(n2)。