最邻近点问题----分治法

分治法

划分问题：把问题的实例划分成子问题

递归求解：递归解决子问题

合并求解：合并子问题的解得到原问题的解

“关键”步骤在与合并。

最邻近点问题

问题描述：

输入：空间平面上点集Q   输出：距离最近的两个点对

以下分为一维、二维、和三维使用分治法来分析最邻近点问题。

一维：如果是在一个直线上找最近的点对，则可以使用排序，之后找最近最近点。

分治思路：

Divide 将其划分成两个部分，Q1，Q2, T(n) = O(n)

Conquer 分别找最近点对<p1,p2>,<q1,q2> T(n) = 2T(n/2)

Merge 比较分开点附近的两个点距离<p2,q3>和找出的<p1,p2> <q1,q2>的距离T(n) = O(n)

时间复杂度：T(n) = O(1)   n=2

T(n) = 2T(n/2) + O(n) = O(nlogn) n >= 3

二维空间问题：

Divide阶段和Conquer阶段都和一维问题类似，使用x=m的直线划分为两个子部分，可以找出最近点对<p1,p2>,<q1,q2>

令d=min{<p1,p2>,<q1,q2>}，则merge阶段的点一定在[m-d,m+d]范围内。

似乎是优化很多，但是最糟糕的情况是在改范围内的点很多则对应起来还是会有n^2级复杂度。

如图，在L:x=m左侧区域中的点P(x0,y0)，则距离其小于d的点在以P的圆心，半径为d的圆内，

则在L:x=m右侧的部分圆一定在以(m,y0-d)的左下角顶点，(m+d,y0+d)为右上角顶点的矩形内，即图中红色矩形区域内。

在对该矩形进行分6块，长为2d/3,宽为d/2，则每个小矩形的内的点最长距离为5d/6(对角线)，即每个小矩形最多有一个点。

因为如果有两个点在同一个矩形则其最小距离不会是d。

因此对[m-d,m+d]范围内每个顶点最多有6个点需要判断。

需要注意的事项：

1.      要设定好数据结构便于处理

2.      要进行预处理，对点集分别按x，y排序，维护两个排序排序结合。

时间复杂度：T(n)=O(1)    n=2

T(n)=O(n)+2T(n/2)+O(n)=O(nlogn)  n≥3

代码如下：

#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>//定义点和两个数组struct POINT {double x,y;}point[10000],temp[10000];//两点之间距离double dis(struct POINT p1, struct POINT p2){return sqrt((p1.x - p2.x) * (p1.x - p2.x) +(p1.y - p2.y) * (p1.y - p2.y) );}//排序判断函数，按x坐标排序int cmp(const void * a, const void * b){struct POINT * c = (struct POINT *)a;struct POINT * d = (struct POINT *)b;if (c->x != d->x){return c->x > d->x;}else return c->y > d->y;}//排序判断函数，按y坐标排序int cmp1(const void * a, const void * b){struct POINT * c = (struct  POINT *)a;struct POINT * d = (struct  POINT *)b;if (c->y != d->y){return c->y > d->y;}else return c->x > d->x;}double findMin(int l, int r){if (l == r){return 10010;}if (l == r - 1){return dis(point[l], point[r]);}double tmp1 = findMin(l,(l + r) >> 1);double tmp2 = findMin(((l + r) >> 1) + 1, r);double Mindis,tmp, mid;mid = point[(l + r) >> 1].x;/*mid = (point[l].x + point[r].x) / 2.0;*/int i,j,cnt = 0;if (tmp1 < tmp2){Mindis = tmp1;}elseMindis = tmp2;for (i = l; i <= r; ++ i){if (fabs(point[i].x - mid) < Mindis){temp[cnt ++] = point[i];}}qsort(temp, cnt, sizeof(temp[0]), cmp1);for (i = 0; i < cnt - 1; ++ i){/*for (j = i + 1; j < cnt; ++ j)*/for (j = i + 1; j < i + 7 && j < cnt; ++ j){tmp = dis(temp[i], temp[j]);if (tmp < Mindis){Mindis = tmp;}}}return Mindis;}int main(){int n,i,j;double minDis;while (scanf("%d", &n)==1 && n){for (i = 0; i < n; ++ i){scanf("%lf%lf", &point[i].x, &point[i].y);}qsort(point, n, sizeof(point[0]), cmp);minDis = findMin(0, n-1);if (minDis > 10000){printf("INFINITY\n");}elseprintf("%.4lf\n", minDis);}return 0;}

三维问题：

#include <iostream>#include <cstdio>#include <cstdlib>#include <cmath>using namespace std;class Point{public:double x;double y;double z;};// [low, high]int randomInteger(int low, int high){return low+rand()%(high-low+1);}void swap(double & a, double & b){double tmp=a;a=b;b=tmp;}int quickSortPartition(Point * p, int low, int high){int random=randomInteger(low, high);swap(p[random], p[low]);int i, j;i=low, j=high;double pivot=p[low].y;for(j=low; j<=high; j++){if(pivot>p[j].y){i=i+1;swap(p[i], p[j]);}}swap(p[i], p[low]);return i;}int partition(double * p, int len, int low, int high){int i, j;int random=randomInteger(low, high);swap(p[low], p[random]);i=low, j=high;double pivot=p[low];for(j=low; j<=high; j++){if(pivot>p[j]){i=i+1;}}swap(p[random], p[low]);return i;}void quickSort(Point * p, int low, int high){if(low<high){int par=quickSortPartition(p, low, high);quickSort(p, low, par-1);quickSort(p, par+1, high);}}double findMiddle(double * p, int len, int ith){if(len==1)return p[0];int random=partition(p, len, 0, len-1);if(random == ith){return p[random];}else if(random < ith){return findMiddle(p+random+1, len-random-1, ith-random-1);}else{return findMiddle(p, random, ith);}}double Distance(const Point & p1, const Point & p2){return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) + (p1.z-p2.z)*(p1.z-p2.z) );}double findMinDistance(Point * p, double * arrayX, double * arrayY, double * arrayZ, int n, Point & des1, Point & des2){double minDistance=INT_MAX;if(n<=3){for(int i=0; i<n; i++){for(int j=0; j<n; j++){if(j != i){if(Distance(p[i], p[j])<minDistance){minDistance=Distance(p[i], p[j]);des1=p[i], des2=p[j];}}}return minDistance;}}else{double middleY=findMiddle(arrayY, n, (int)floor((double)(n/2)));// partition P set into PL and PR according to the middle value of yPoint * pL=new Point[n];Point * pR=new Point[n];int lIndex=0, rIndex=0;double * leftArrayX=new double[n];double * leftArrayY=new double[n];double * leftArrayZ=new double[n];double * rightArrayX=new double[n];double * rightArrayY=new double[n];double * rightArrayZ=new double[n];for(int i=0; i<n; i++){if(p[i].y<middleY){pL[lIndex]=p[i];leftArrayX[lIndex]=p[i].x;leftArrayY[lIndex]=p[i].y;leftArrayZ[lIndex]=p[i].z;lIndex++;}else{pR[rIndex]=p[i];rightArrayX[rIndex]=p[i].x;rightArrayY[rIndex]=p[i].y;rightArrayZ[rIndex]=p[i].z;rIndex++;}}// divide and conquerdouble lMinDistance=findMinDistance(pL, leftArrayX, leftArrayY, leftArrayZ, lIndex, des1, des2);double rMinDistance=findMinDistance(pR, leftArrayX, leftArrayY, leftArrayZ, rIndex, des1, des2);double currentMinDistance=lMinDistance<rMinDistance ? lMinDistance : rMinDistance;Point * middleP=new Point[n];double * middleArrayY=new double[n];int middleIndex=0;for(int i=0; i<n; i++){if(fabs(p[i].y-middleY)<=currentMinDistance){middleP[middleIndex]=p[i];middleArrayY[middleIndex]=p[i].y;}}for(int i=0; i<middleIndex; i++){for(int j=1; j<=15; j++){if(Distance(p[i], p[i+j])<currentMinDistance){currentMinDistance=Distance(p[i], p[i+j]);des1=p[i];des2=p[i+j];}}}minDistance=currentMinDistance;delete [] pL;delete [] pR;delete [] leftArrayX;delete [] leftArrayY;delete [] leftArrayZ;delete [] rightArrayX;delete [] rightArrayY;delete [] rightArrayZ;delete [] middleP;delete [] middleArrayY;return minDistance;}}int main(){Point * p=new Point[4];for(int i=0; i<4; i++){p[i].x=4-i;p[i].y=4-i;p[i].z=4-i;}quickSort(p, 0, 3);double * arrayX, * arrayY, * arrayZ;arrayX=new double[4];arrayY=new double[4];arrayZ=new double[4];for(int i=0; i<4; i++){arrayX[i]=p[i].x;arrayY[i]=p[i].y;arrayZ[i]=p[i].z;}for(int i=0; i<4; i++){cout<<p[i].x<<" "<<p[i].y<<" "<<p[i].z<<endl;}double minDistance;Point des1, des2;minDistance=findMinDistance(p, arrayX, arrayY, arrayZ, 4, des1, des2);cout << "The minimal distance is:" << minDistance <<endl;cout << "The two points that have the minimal distance is:" << endl;cout << "("<<des1.x<<","<<des1.y<<","<<des1.z<<")"<<endl;cout << "("<<des2.x<<","<<des2.y<<","<<des2.z<<")"<<endl;getchar();return 0;}