算法作业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中放。

代码实现：

  1 #include<iostream>
  2 #include<stack>
  3 #include<vector>
  4 #include<time.h>
  5 #define maxn 100
  6 using namespace std;
  7 
  8 struct Point{
  9     double x;
 10     double y;
 11     double z;
 12 }dataset[maxn];
 13 
 14 stack<Point> resultpairs;
 15 
 16 int cmp(const void *a, const void *b){
 17     return ((Point *)a)->y - ((Point *)b)->y;
 18 }
 19 
 20 double Distance(const Point &p1, const Point &p2){
 21     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));
 22 }
 23 
 24 double findMinDistance(int left_m, int right_m, int n){
 25     double minDistance = INT_MAX;
 26 
 27     //小于等于3个点时暴力解决求出最短距离
 28     if(n <= 3){
 29         for (int i = left_m; i <= right_m; i++){
 30             for (int j = i + 1; j <= right_m; j++){
 31                 if(Distance(dataset[i], dataset[j]) <= minDistance){
 32                     minDistance = Distance(dataset[i], dataset[j]);
 33 
 34                     //把距离最小的点压栈到resultpairs
 35                     if(!resultpairs.empty()){
 36                         Point t1 = resultpairs.top();
 37                         resultpairs.pop();
 38                         Point t2 = resultpairs.top();
 39                         resultpairs.pop();
 40 
 41                         double now = Distance(t1, t2);
 42                         if(minDistance < now){
 43                             while(!resultpairs.empty())resultpairs.pop();
 44                             resultpairs.push(dataset[i]);
 45                             resultpairs.push(dataset[j]);
 46                         }
 47                         else {
 48                             resultpairs.push(t2);
 49                             resultpairs.push(t1);
 50                             if(minDistance == now){
 51                                 resultpairs.push(dataset[i]);
 52                                 resultpairs.push(dataset[j]);
 53                             }
 54                         }
 55                     }
 56                     else{
 57                         resultpairs.push(dataset[i]);
 58                         resultpairs.push(dataset[j]);
 59                     }
 60                 }
 61             }
 62         }
 63         return minDistance;
 64     }
 65     else{
 66 
 67         //找到在Y轴看最中间的点，记下其Y坐标，把点集分成两堆pL和pR。
 68         int middle = (left_m + right_m) / 2;
 69         double middleY = dataset[middle].y;
 70 
 71         //分而治之
 72         //找到两边最小的点对，把距离放进currentMinDistance
 73         double lMinDistance = findMinDistance(left_m, middle, middle - left_m + 1);
 74         double rMinDistance = findMinDistance(middle + 1, right_m, right_m - middle);
 75         double currentMinDistance = lMinDistance < rMinDistance ? lMinDistance: rMinDistance;
 76 
 77         //找出所有分居切面两边，且两两之间有可能成为最小距离的点放进middleP,它们的Y坐标放进middleArrayY
 78         vector<Point> l_middle_Points;
 79         vector<Point> r_middle_Points;
 80         for(int i = left_m; i <= right_m; i++){
 81             if((i <= middle) && (middleY - dataset[i].y <= currentMinDistance)){
 82                 l_middle_Points.push_back(dataset[i]);
 83             }
 84             if((i > middle) && (dataset[i].y - middleY <= currentMinDistance)){
 85                 r_middle_Points.push_back(dataset[i]);
 86             }
 87         }
 88 
 89         vector<Point>::iterator end = l_middle_Points.end();
 90         vector<Point>::iterator end2 = r_middle_Points.end();
 91 
 92         double newdis = INT_MAX;
 93         for(vector<Point>::iterator it = l_middle_Points.begin(); it != end; ++it){
 94             for(vector<Point>::iterator it2 = r_middle_Points.begin(); it2 != end2; ++it2){
 95 
 96                 if((fabs(it->x - it2->x) <= currentMinDistance) && (fabs(it->z - it2->z) <= currentMinDistance)){
 97                     newdis = Distance(*it, *it2);
 98                     if(newdis <= currentMinDistance){
 99                         currentMinDistance = newdis;
100                         if(!resultpairs.empty()){
101                             Point t1 = resultpairs.top();
102                             resultpairs.pop();
103                             Point t2 = resultpairs.top();
104                             resultpairs.pop();
105 
106                             double now = Distance(t1, t2);
107                             if(newdis < now){
108                                 while(!resultpairs.empty())resultpairs.pop();
109                                 resultpairs.push(*it);
110                                 resultpairs.push(*it2);
111                             }
112                             else {
113                                 resultpairs.push(t2);
114                                 resultpairs.push(t1);
115                                 if(newdis == now){
116                                     resultpairs.push(*it);
117                                     resultpairs.push(*it2);
118                                 }
119                             }
120                         }
121                         else{
122                             resultpairs.push(*it);
123                             resultpairs.push(*it2);
124                         }
125                     }
126                 }
127             }
128         }
129         minDistance = currentMinDistance;
130         return minDistance;
131     }
132 }
133 
134 int main(){
135 
136     //把标准输入重定向到文件输入
137     freopen("dataout.txt","r",stdin);  
138 
139     //记录程序运行时间
140     clock_t start,finish;  
141     double totaltime;  
142     start=clock();  
143 
144     //数据读取
145     for(int i = 0; i < maxn; i++){
146         cin >> dataset[i].x;
147         cin >> dataset[i].y;
148         cin >> dataset[i].z;
149     }
150 
151     //快速排序点集
152     qsort(dataset, maxn, sizeof(dataset[0]), cmp);
153 
154     //查找最短距离，并将所有距离相等的点对压入resultpairs
155     double mindis = findMinDistance(0, maxn - 1, maxn);
156 
157     cout<<"shortest distance is: "<<mindis<<endl;
158     cout<<"
all closest pairs are: "<<endl;
159     
160     Point a;
161     while(!resultpairs.empty()) {
162         a = resultpairs.top();
163         cout << "(" << a.x <<", " <<a.y<<", "<<a.z<<") and ";
164         resultpairs.pop();
165         a = resultpairs.top();
166         cout << "(" << a.x <<", " <<a.y<<", "<<a.z<<")"<<endl;
167         resultpairs.pop();
168     }
169 
170 
171     finish=clock();  
172     totaltime=(double)(finish-start)/CLOCKS_PER_SEC;  
173     cout<<"
it runs "<<totaltime<<" seconds！"<<endl; 
174     return 0;  
175 }

View Code