算法作业1 3D最近点对问题

算法思路：

数据预处理。先把点集中的所有点，按照Y坐标大小从左往右排序，得到关于Y轴的有序点集合。具体使用STL的qsort方法，自定义比较函数。

使用分治的思想， 把点集分开出来求解。先选取点集中Y坐标最中间的点，该点所在垂直于Y轴的平面把点集划分成两个子点集leftSubSet和rightSubSet，然后分别计算两个子点集中最短的点对。以此类推，分到子点集仅包含3个点以下。当子点集小于等于3个点时，使用暴力遍历法求解这些点对的最短距离。

“分”求出leftSubSet和rightSubSet的最短距离currentMin后，“治”的方法是，把Y坐标与切平面的距离小于currentMin的点取出，分成l_middle和r_middle。因为只有这些点才有可能构成两个子点集的最短距离点对。接着对l_middle中的每一个点来说，只有与r_middle中点的X坐标和Z坐标小于currentMin的点才有可能是最短距离点。分别计算他们的欧几里得距离，将距离最小的点对记录下来，更新currentMin。

为了能找出所有最短距离相等的点对，这里使用了一个stack全局变量把最短距离点对存储起来，每次计算出比stack中的点对距离还要小的点对时，把stack清空，把新点对压进stack。如果距离一样，则继续往stack中放。

代码实现：

#include<iostream>
#include<stack>
#include<vector>
#include<time.h>
#define maxn 100
using namespace std;

struct Point{
double x;
double y;
double z;
}dataset[maxn];

stack<Point> resultpairs;

int cmp(const void *a, const void *b){
return ((Point *)a)->y - ((Point *)b)->y;
}

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(int left_m, int right_m, int n){
double minDistance = INT_MAX;

//小于等于3个点时暴力解决求出最短距离
if(n <= 3){
for (int i = left_m; i <= right_m; i++){
for (int j = i + 1; j <= right_m; j++){
if(Distance(dataset[i], dataset[j]) <= minDistance){
minDistance = Distance(dataset[i], dataset[j]);

//把距离最小的点压栈到resultpairs
if(!resultpairs.empty()){
Point t1 = resultpairs.top();
resultpairs.pop();
Point t2 = resultpairs.top();
resultpairs.pop();

double now = Distance(t1, t2);
if(minDistance < now){
while(!resultpairs.empty())resultpairs.pop();
resultpairs.push(dataset[i]);
resultpairs.push(dataset[j]);
}
else {
resultpairs.push(t2);
resultpairs.push(t1);
if(minDistance == now){
resultpairs.push(dataset[i]);
resultpairs.push(dataset[j]);
}
}
}
else{
resultpairs.push(dataset[i]);
resultpairs.push(dataset[j]);
}
}
}
}
return minDistance;
}
else{

//找到在Y轴看最中间的点，记下其Y坐标，把点集分成两堆pL和pR。
int middle = (left_m + right_m) / 2;
double middleY = dataset[middle].y;

//分而治之
//找到两边最小的点对，把距离放进currentMinDistance
double lMinDistance = findMinDistance(left_m, middle, middle - left_m + 1);
double rMinDistance = findMinDistance(middle + 1, right_m, right_m - middle);
double currentMinDistance = lMinDistance < rMinDistance ? lMinDistance: rMinDistance;

//找出所有分居切面两边，且两两之间有可能成为最小距离的点放进middleP,它们的Y坐标放进middleArrayY
vector<Point> l_middle_Points;
vector<Point> r_middle_Points;
for(int i = left_m; i <= right_m; i++){
if((i <= middle) && (middleY - dataset[i].y <= currentMinDistance)){
l_middle_Points.push_back(dataset[i]);
}
if((i > middle) && (dataset[i].y - middleY <= currentMinDistance)){
r_middle_Points.push_back(dataset[i]);
}
}

vector<Point>::iterator end = l_middle_Points.end();
vector<Point>::iterator end2 = r_middle_Points.end();

double newdis = INT_MAX;
for(vector<Point>::iterator it = l_middle_Points.begin(); it != end; ++it){
for(vector<Point>::iterator it2 = r_middle_Points.begin(); it2 != end2; ++it2){

if((fabs(it->x - it2->x) <= currentMinDistance) && (fabs(it->z - it2->z) <= currentMinDistance)){
newdis = Distance(*it, *it2);
if(newdis <= currentMinDistance){
currentMinDistance = newdis;
if(!resultpairs.empty()){
Point t1 = resultpairs.top();
resultpairs.pop();
Point t2 = resultpairs.top();
resultpairs.pop();

double now = Distance(t1, t2);
if(newdis < now){
while(!resultpairs.empty())resultpairs.pop();
resultpairs.push(*it);
resultpairs.push(*it2);
}
else {
resultpairs.push(t2);
resultpairs.push(t1);
if(newdis == now){
resultpairs.push(*it);
resultpairs.push(*it2);
}
}
}
else{
resultpairs.push(*it);
resultpairs.push(*it2);
}
}
}
}
}
minDistance = currentMinDistance;
return minDistance;
}
}

int main(){

//把标准输入重定向到文件输入
freopen("dataout.txt","r",stdin);

//记录程序运行时间
clock_t start,finish;
double totaltime;
start=clock();

//数据读取
for(int i = 0; i < maxn; i++){
cin >> dataset[i].x;
cin >> dataset[i].y;
cin >> dataset[i].z;
}

//快速排序点集
qsort(dataset, maxn, sizeof(dataset[0]), cmp);

//查找最短距离，并将所有距离相等的点对压入resultpairs
double mindis = findMinDistance(0, maxn - 1, maxn);

cout<<"shortest distance is: "<<mindis<<endl;
cout<<"\nall closest pairs are: "<<endl;

Point a;
while(!resultpairs.empty()) {
a = resultpairs.top();
cout << "(" << a.x <<", " <<a.y<<", "<<a.z<<") and ";
resultpairs.pop();
a = resultpairs.top();
cout << "(" << a.x <<", " <<a.y<<", "<<a.z<<")"<<endl;
resultpairs.pop();
}

finish=clock();
totaltime=(double)(finish-start)/CLOCKS_PER_SEC;
cout<<"\nit runs "<<totaltime<<" seconds！"<<endl;
return 0;
}

View Code