分治法最近点对
分治法
分治法将一个难以直接解决的大问题划分成一些规模较小的子问题，分别求解各个子问题，再合并子问题的解得到原问题的解。
一般来说，分治法的求解过程由以下三个阶段组成：
- 划分：把规模为n的原问题划分为k个规模较小的子问题。
- 求解子问题：各个子问题的解法与原问题的解法通常是相同的，可以用递归的方法求解各个子问题，有时递归也可以用循环来实现。
- 合并：把各个子问题的解合并起来，合并的代价因情况不同有很大的差异，分治算法的效率很大程度上依赖于合并的实现。
最近点对问题
问题描述
设p1=(x1,y1)，p2=(x2,y2)，···，pn=(xn,yn)是平面上n个点构成的集合S，最近点对问题就是找出集合S中距离最近的点对。严格地说，最接近的点对可能多于一个，简单起见，找出一个即可。
应用实例
假设在一个金属片上钻n个大小一样的洞，如果洞的距离太近，金属就可能会断裂。可以通过任意两个洞的最小距离来估算金属断裂的概率。这种最小距离问题实际上也就是最近点对问题。
想法
最近点对的分治策略如下：
（1）划分：将集合S分成两个子集S1和S2根据平衡子问题原则，每个子集中大约有n/2个点，设集合S的最近对是pi和pj，则会出现以下三种情况：
- pi∈S1，pj∈S1，即最近对均在集合S1中；
- pi∈S2，pj∈S2，即最近对均在集合S2中；
- pi∈S1，pj∈S2，即最近对分别在集合S1和S2中。
（2）求解子问题：对于划分阶段的情况1和情况2可递归求解，如果最近对分别在两个子集中，在后边做讨论。
（3）合并：比较在划分阶段三种情况下的最近点对，取三者之中距离较小者为原问题的解。
下面讨论划分阶段的情况3：
为了将平面上的点集S分割为点的个数大致相同的两个集合，选取垂直线x=m来作分割线，其中m为S中各点x坐标的中位数。由此将S分割为S1和S2。递归地在S1和S2中求解最近对问题，分别得到S1中的最近距离d1和S2中的最近距离d2，令d=min(d1，d2。若S的最近对(p,q)之间的距离小于d，则p和q必定分属于集合S1和S2。那么，不妨设p∈S1，q∈S2，则p和q距离直线x=m的距离均小于d，所以，可以将求解限制在以x=m为中心，宽为2d的垂直带P1和P2中，垂直带之外的任何点对之间的距离都一定大于d。
假设点p(x,y)是集合P1和P2中y坐标最小的点，则点p可能在P1中也可能在P2中，现在需要找出来和点p之间小于d的点，显然这样的点的y坐标一定位于区间[y,y+d]之间，而且这样的点不会超过8个，因为P1和P2中点的相互之间的距离至少为d，如下图所示。
所以，可以将P1和P2中的点按照y坐标升序排列，顺序地处理P1和P2中的点p(x,y)，在y坐标区间[y,y+d]内最多取出8个候选点，计算它们和点p之间的距离。
算法
简单起见，假设点集S已按照x坐标升序排列，伪代码描述如下：
输入：按x坐标升序排列的n个点的集合S
输出：最近点对的距离
- 如果n等于2，则返回(x1,y1)和(x2,y2)之间的距离，算法结束；
- 划分：m=S中各点x坐标的中位数；
- d1=计算{(x1,y1)，···，(xm,ym)}的最近对距离；
- d2=计算{(xm,ym)，···，(xn,yn)}的最近对距离；
- d=min{d1,d2}；
- 依次考察集合S中的点p(x,y)，如果x<=xm并且x>=xm-d，则将点p放入集合P1中；如果x>xm并且x<=xm+d，则将点p放入集合P2中；
- 将集合P1和P2按y坐标升序排列；
- 对集合P1和P2中的每个点p(x,y)，在y坐标区间[y,y+d]内最多取出8个候选点，计算与点p的最近距离d3；
- 返回min{d，d3}；
源代码
#include<iostream>
#include<math.h>
#include<stdlib.h>
using namespace std;
int M[2], N[2], P[2], Q[2];
typedef struct point
{
double x,y;
}point;
void Qsortx(point *A, int low, int high)
{
double tempx, tempy;
int i = low, j = high;
if(low < high)
{
tempx = A[i].x;
tempy = A[i].y;
while(i != j)
{
while(j > i && A[j].x > tempx)
j--;
if(i < j)
{
A[i].x = A[j].x;
A[i].y = A[j].y;
i++;
}
while(i < j && A[i].x < tempx)
i++;
if(i<j)
{
A[j].x = A[i].x;
A[j].y = A[i].y;
j--;
}
A[i].x = tempx;
A[i].y = tempy;
Qsortx(A,low,i-1);
Qsortx(A,i+1,high);
}
}
}
void Qsorty(point *A, int low, int high)
{
double tempx, tempy;
int i = low, j = high;
if(low < high)
{
tempx = A[i].x;
tempy = A[i].y;
while(i != j)
{
while(j > i && A[j].y > tempy)
j--;
if(i<j)
{
A[i].x = A[j].x;
A[i].y = A[j].y;
i++;
}
while(i < j && A[i].y < tempy)
i++;
if(i < j)
{
A[j].x = A[i].x;
A[j].y = A[i].y;
j--;
}
A[i].x = tempx;
A[i].y = tempy;
Qsorty(A,low,i-1);
Qsorty(A,i+1,high);
}
}
}
double d(point p1, point p2)
{
double dx, dy, dm;
dx = p1.x - p2.x;
dy = p1.y - p2.y;
dm = pow(dx,2) + pow(dy,2);
return sqrt(dm);
}
void printarry(point *points, int n)
{
for(int i=0; i < n; i++)
{
cout<<"("<<points[i].x<<","<<points[i].y<<")";
}
cout<<endl;
cout<<endl;
}
double cp(point *S,point *Y,int low,int high)
{
int e, r, mid, k = 0, n=high-low+1;
double x0, min, minl, minr, minlr;
point T[n];
mid = (low + high) / 2;
x0 = S[mid].x;
if(n == 2)
{
return d(S[low], S[high]);
M[0] = low;
M[1] = high;
}
if(n == 3)
{
double d1, d2, d3;
d1 = d(S[low], S[low+1]);
d2 = d(S[low], S[low+2]);
d3 = d(S[low+1], S[low+2]);
if(d1<=d2)
{
min = d1;
M[0] = low;
M[1] = low+1;
}
else
{
min = d2;
M[0] = low;
M[1] = low+2;
}
if(min > d3)
{
min = d3;
M[0] = low + 1;
M[1] = low + 2;
}
}
else
{
minl = cp(S, Y, low, mid);
N[0] = M[0];
N[1] = M[1];
minr = cp(S, Y, mid+1, high);
P[0] = M[0];
P[1] = M[1];
if(minl <= minr)
{
min = minl;
Q[0] = N[0];
Q[1] = N[1];
}
else
{
min = minr;
Q[0] = P[0];
Q[1] = P[1];
}
for(int i = 0; i < n; i++)
{
if(fabs(Y[i].x-x0) <= min)
{
T[k].x = Y[i].x;
T[k].y = Y[i].y;
k++;
}
}
minlr = 2*min;
for(int i=0;i<k;i++)
{
int jj;
if((i+7) < k)
jj=i+7;
else jj = k;
for(int j = i+1; j < jj; j++)
{
double dd = d(T[i], T[j]);
if(dd<minlr)
{
minlr = dd;
e = i;
r = j;
}
}
}
if(minlr < min)
{
min = minlr;
for(int i = 0; i < n; i++)
{
if(T[e].x == S[i].x && T[e].y == S[i].y)
Q[0] = i;
}
for(int j = 0; j < n; j++)
{
if(T[r].x == S[j].x && T[r].y == S[j].y)
Q[1] = j;
}
}
}
return min;
}
int main()
{
int n;
int f = 1;
double mini;
cout<<"请输入点的个数:"<<endl;
while(f)
{
cin>>n;
if(n <= 3)
cout<<"非法输入！请重新输入："<<endl;
else
{
f = 0;
break;
}
}
point S[n], Y[n];
cout<<"请输入这"<<n<<"个点的坐标:"<<endl;
for(int i = 0; i < n; i++)
{
cin>>S[i].x>>S[i].y;
}
Qsortx(S, 0, n-1);
for(int i = 0; i < n; i++)
{
Y[i].x = S[i].x;
Y[i].y = S[i].y;
}
Qsorty(Y, 0, n-1);
printarry(S, n);
printarry(Y, n);
mini = cp(S, Y, 0, n-1);
cout<<"最近的兩個點是："<<endl;
cout<<"("<<S[Q[0]].x<<","<<S[Q[0]].y<<")"<<"<----->"<<"("<<S[Q[1]].x<<","<<S[Q[1]].y<<")"<<endl;
cout<<endl;
cout<<"最近点对的距离为:"<<mini<<endl;
system("pause");
return 0;
}
算法分析
应用分治算法求解含有n个点的最近对问题，由于划分阶段的情况1和情况2可递归求解，情况3的时间代价是O(n)，合并子问题解的时间代价是O(1)，则算法的时间复杂度可由下边的递推式表示：
可得T(n)=O(nlogn)。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
