平面最近点对

最近点对问题：给定平面上n个点的集合S，找其中的一对点，使得在n个点组成的所有点对中，该点对间的距离最小。 

思路：

（1）先考虑一维情况下，平面上的点退化成数轴上的n个实数点，最近点实际上就是实数中相差的最小点。我们用坐标上求出S中第n/2（向上取整） 小、第(n/2（向上取整） +1)小的p，q的坐标作为分割点（确保了两个子部分的大小接近），将数轴分为两部分递归的求解，再将p,q距离与分隔的两部分算出的最少距离比较，得出实际该段的最小距离。

复杂度为 O(n)=nlgn

（2）在二维情况下，选取x=l的垂线分割平面，分割位置在所有点x坐标轴映射的第n/2（向上取整） 小、第(n/2（向上取整） +1)小的p，q两点之间，递归对其左右两个二维平面求解最近点对。对于得到当前的最近距离d，将与l距离的绝对值小于等于d的区域的点按照y值从小到大，遍历l左端的点，计算他与右端且满足y在给点+-d区域内的点的距离与当前距离最小值d比较，得该平面最近点对。

===============================================================

问题（北邮算法第二章课后）：采用平面最近点对算法，根据基站经纬度，挑选出

距离最近的2个基站

距离次最近的2个基站

要求：返回

1) 最近/次最近的2个基站间距离

2) 最近/次最近的2个基站点对（用基站ENodeBID表示）

===============================================================
#include<iostream>
#include<string>
#include<fstream>
#include<cmath>
#include<algorithm>
#define max 2000
#define EARTH_RADIUS 6378.13700

using namespace std;

typedef struct Basestation
{
int ENODEBID;
double LONGITUDE;//经度
double LATITUDE;//纬度
double K_DIST ;
}Basestation;

typedef struct Pair
{
Basestation a;
Basestation b;
double d;
Pair(){}
Pair(Basestation a, Basestation b, double d):a(a), b(b), d(d){};
}Pair;

Basestation a[max], b[max], c[max];
Pair pair1, pair2;

double Rad(double Latitude_Longitude);//把经纬度转化成弧度
double GetDistance(Basestation a, Basestation b);//获取两个基站之间的距离
Pair closestPair( int l, int r);
double min(double a, double b)
{
if(a<b)
return a;
else
return b;
}

bool SortByX(Basestation a, Basestation b)
{
if(a.LATITUDE == b.LATITUDE)
return a.LONGITUDE<b.LONGITUDE;
return a.LATITUDE<b.LATITUDE;
}
bool SortByY(const Basestation &a, const Basestation &b)
{
return a.LONGITUDE<b.LONGITUDE;
}

int main()
{
int i = 1;
int num = 0;
ifstream file;
file.open("Data_Of_Basestation.txt", ios::in);
if(file.bad())
{
cout<<"打开文件时发生错误"<<endl;
return 0;
}
while(!file.eof())
{
file>>a[i].ENODEBID>>a[i].LONGITUDE>>a[i].LATITUDE>>a[i].K_DIST;
i++;
num++;
//cout<<a[i].ENODEBID<<" "<<a[i].LONGITUDE<<" "<<a[i].LATITUDE<<" "<<a[i].K_DIST<<endl;
file.get();
if(file.peek()==EOF)
break;
}

sort(a+1, a+num, SortByX);
for(int i = 1; i<=num; i++)
{
b[i] = a[i];
}
sort(b+1, b+num, SortByY);
pair1.d = pair2.d = 100000.0;
closestPair(1, num);

cout<<"==============================================================================="<<endl;
cout<<"距离最近的两个基站：\t";
cout<<pair1.a.ENODEBID<<"\t\t"<<pair1.b.ENODEBID<<"\t\t"<<endl;
cout<<"最近距离为："<<pair1.d<<endl<<endl;
cout<<"距离次最近的两个基站：\t";

cout<<pair2.a.ENODEBID<<"\t\t"<<pair2.b.ENODEBID<<endl;
cout<<"最近距离为："<<pair2.d<<endl;
cout<<"==============================================================================="<<endl;
file.close();
return 0;
}

double Rad(double Latitude_Longitude)
{
return Latitude_Longitude*M_PI/180.0;
}

double GetDistance(Basestation a, Basestation b)
{
double RadLat1, RadLat2, RadLon1, RadLon2, dis;
RadLat1 = Rad(a.LATITUDE);
RadLon1 = Rad(a.LONGITUDE);
RadLat2 = Rad(b.LATITUDE);
RadLon2 = Rad(b.LONGITUDE);
dis = acos(cos(RadLat1)*cos(RadLat2)*cos(RadLon1-RadLon2)+sin(RadLat1)*sin(RadLat2))*EARTH_RADIUS;
dis*=1000.0;
cout<<dis<<endl;
return dis;
}

Pair closestPair(int l, int r)
{
if( r== l)//1点的情形
{
return Pair(a[l], a[r], 100000);
}
if(1 == r-l)//2点的情形
{
double dis1 = GetDistance(a[l], a[r]);
if(dis1 < pair1.d)
{
pair2 = pair1;
pair1 = Pair(a[l], a[r], dis1);
}
else if(dis1>pair1.d && dis1<pair2.d)
{
pair2 = Pair(a[l], a[r], dis1);
}
return Pair(a[l], a[r], dis1);
}

int mid = (l+r)/2;
int f = l,
g = mid+1;
Pair best = closestPair(l, mid);
Pair right = closestPair(mid+1, r);
Pair beitai;//无处不在的备胎
if(best.d > right.d)
{
beitai = best;
best = right;
right = best;
}
if(best.d < pair1.d)
{
pair2 = pair1;
pair1 = best;
}
else if(best.d > pair1.d && best.d <pair2.d)
{
pair2 = best;
}
else if(right.d >pair1.d && right.d<pair2.d)
{
pair2 = right;
}

int k = 1;
double temp = a[mid].LATITUDE;
for(int i = l; i<=r; i++)
{
if(fabs(a[i].LATITUDE - temp) <= best.d)
{
c[k++] = a[i];
}
}
k--;
sort(c+1, c+k, SortByY);
for(int i=1; i<=k; i++)
{
for(int j = i+1; j<=k && c[j].LONGITUDE - c[i].LONGITUDE<best.d; j++)
{
double dis2 = GetDistance(c[i], c[j]);
if(dis2 < pair1.d)
{
pair2 = pair1;
pair1 = Pair(c[i], c[j], dis2);
}
else if(dis2 > pair1.d && dis2<pair2.d)
{
pair2 = Pair(c[i], c[j], dis2);
}
if(dis2 < best.d)
best = Pair(c[i], c[j], dis2);
}
}
return best;
}