最鄰近點問題----分治法
阿新 • • 發佈：2019-01-22
分治法
劃分問題：把問題的例項劃分成子問題
遞迴求解：遞迴解決子問題
合併求解：合併子問題的解得到原問題的解
“關鍵”步驟在與合併。
最鄰近點問題
問題描述：
輸入：空間平面上點集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=2T(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;
}
else
Mindis = 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");
}
else
printf("%.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 y
Point * 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 conquer
double 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;
}