编程之美2.11——寻找最近点对（POJ 3714）

问题：

给定平面上N个点的坐标，找出距离最近的两个点。

解法：

我们先对N个点的x坐标进行排序，排序我们使用最坏复杂度O(n*logn)的快速排序方法，在排序的过程中minDifferent会递归计算出左右两边的最小距离，再用其中的较小值minum得到以中位数点附近的带状区域[p[median+1].x-median, p[median].x+median]，对带状区域的点按照y坐标排序，对带状区域的每个点只需计算最多7个点，就能得到所有可能小于minum的点对。

#include <iostream>#include <vector>#include <algorithm>#include <cmath>using namespace std;// 顶点信息struct Point{double m_x, m_y;Point():m_x(0.0),m_y(0.0) {}Point(double x, double y):m_x(x),m_y(y){}bool operator==(const Point& p) const {return m_x==p.m_x && m_y==p.m_y;}};ostream& operator<<(ostream& os, const Point& p){return os << "(" << p.m_x << "," << p.m_y << ")";}// 插入排序template<class T, class Pr>void insert_sort(vector<T> &vec, int l, int r, Pr pred){int i, j;for (i=l+1; i<=r; i++){T tmp = vec[i];for (j=i-1; j>=l && pred(tmp,vec[j]); j--) vec[j+1]=vec[j];vec[j+1] = tmp;}}// 找到key所在的位置template<class T>int get_position(vector<T> &vec, int l, int r, T key){for (int i=l; i<=r; i++)if (key == vec[i])return i;return -1;}// 按第一个元素对vec进行划分template<class T, class Pr>int partition(vector<T> &vec, int l, int r, Pr pred){int i, j;for (i=l+1,j=l; i<=r; i++){if (pred(vec[i],vec[l])){++j;swap(vec[i],vec[j]);}}swap(vec[j],vec[l]);return j;}// 顺序统计得到第k个元素的值template<class T, class Pr>T select(vector<T> &vec, int l, int r, int k, Pr pred){int n = r-l+1;if (n==1){if (k!=0)printf("Out of Boundary!\n");return vec[l];}// 找中位数的中位数作为分割点int cnt = n/5;int tcnt = (n+4)/5;int rem = n%5;vector<T> group(tcnt);int i, j;for (i=0,j=l; i<cnt; i++,j+=5){insert_sort(vec, j, j+4, pred);group[i] = vec[j+2];}if (rem){insert_sort(vec, j, j+rem-1, pred);group[i] = vec[j+(rem-1)/2];}T key = select(group, 0, tcnt-1, (tcnt-1)/2, pred);// 找到分割点的位置int key_pos = get_position(vec, l, r, key);swap(vec[key_pos], vec[l]);// 用分割点对数组进行划分，小的在左边，大的在右边int pos = partition(vec, l, r, pred);int x = pos - l;if (x == k) return key;else if (x < k) return select(vec, pos+1, r, k-x-1, pred);elsereturn select(vec, l, pos-1, k, pred);}// 计算点a和b的距离double dist(const Point& a, const Point& b){double x = a.m_x-b.m_x;double y = a.m_y-b.m_y;return sqrt(x*x+y*y);}bool cmpX(const Point& a, const Point& b){return a.m_x < b.m_x;}bool cmpY(const Point& a, const Point& b){return a.m_y < b.m_y;}double minDifferent(vector<Point> p, int l, int r, vector<Point> &result){// 按中位数进行划分后的子区域的元素个数都会减小到2或3，不会再到1if ((r-l+1)==2){result[0] = p[l];result[1] = p[r];if (cmpX(p[r],p[l])) swap(p[l], p[r]);return dist(p[l], p[r]);}if ((r-l+1)==3){insert_sort(p, l, r, cmpX);double tmp1 = dist(p[l], p[l+1]);double tmp2 = dist(p[l+1], p[l+2]);double ret = min(tmp1, tmp2);if (tmp1 == ret){result[0] = p[l];result[1] = p[l+1];}else{result[0] = p[l+1];result[1] = p[l+2];}return ret;}// 大于3个点的情况int mid = (r+l)>>1;Point median = select(p, l, r, mid-l, cmpX);vector<Point> res1(2), res2(2);double min_l = minDifferent(p, l, mid, res1);double min_r = minDifferent(p, mid+1, r, res2);double minum = min(min_l, min_r);if (minum == min_l){result[0] = res1[0];result[1] = res1[1];}else{result[0] = res2[0];result[1] = res2[1];}// 对[p[mid+1]-minum, p[mid]+minum]的带状区域按y排序vector<Point> yvec;int i, j;for (i=mid+1; i<=r; i++)if (p[i].m_x - p[mid].m_x < minum)yvec.push_back(Point(p[i]));for (i=mid; i>=l; i--)if (p[mid+1].m_x - p[i].m_x < minum)yvec.push_back(Point(p[i]));sort(yvec.begin(), yvec.end(), cmpY);for (i=0; i<yvec.size(); i++){// 至多只有与其后最多7个点的距离会小于minumfor (j=i+1; j<yvec.size() && yvec[j].m_y-yvec[i].m_y<minum &&j<=i+7; j++){double delta = dist(yvec[i],yvec[j]);if (delta < minum){minum = delta;result[0] = yvec[i];result[1] = yvec[j];}}}return minum;}int main(){int n, i, j, x, y;vector<Point> result(2);vector<Point> input;cin >> n;for (i=0; i<n; i++){cin >> x;cin >> y;input.push_back(Point(x,y));}double minum = minDifferent(input, 0, input.size()-1, result);cout << "nearest point: " << result[0] << " and "  << result[1] << endl;cout << "distance: " << minum << endl;return 0;}

POJ 3714 问题：

平面上有两类点，计算属于不同类的顶点对的最小值。

解法：参考http://blog.csdn.net/smsmn/article/details/5963487

算法思想与上面基本相同，但编程方式上进行了改变，更好理解。下面的代码写法上完全按照参考代码的思路，只是将数组操作改为vector，但提交到POJ 3714上却会TLE，看来动态分配空间所占用的时间也不小。所以要想AC，请使用参考代码。

#include <iostream>#include <vector>#include <algorithm>#include <cmath>using namespace std;const double INF = 1e100;struct Point{double x, y;int flag;// 顶点的类别Point(){}Point(double xx, double yy):x(xx),y(yy){}};vector<Point> p;bool cmp1(const Point& a, const Point& b){return a.x < b.x;}bool cmp2(int a, int b){return p[a].y < p[b].y;}double dist(const Point& a, const Point& b){double xx = a.x - b.x;double yy = a.y - b.y;return sqrt(xx*xx+yy*yy);}// 输出属于不同类的顶点对的最小值double min_dist(vector<Point> p, int left, int right){int mid = (left+right)>>1, i,j;if (left>=right) return INF;for (i=mid; i>=left && p[mid].x<=p[i].x; i--);double minum = min_dist(p, left, i);for (i=mid; i<=right && p[i].x<=p[mid].x; i++);minum = min(minum, min_dist(p, i, right));vector<int> yp;for (i=mid; i>=left && p[mid].x-p[i].x<minum; i--)yp.push_back(i);for (i=mid+1; i<=right && p[i].x-p[mid].x<minum; i++)yp.push_back(i);// 这个方法非常巧妙，直接对顶点索引进行排序，减少了空间使用，// 代码上也更加简洁sort(yp.begin(), yp.end(), cmp2);for (i=0; i<yp.size(); i++)for (j=i+1; j<yp.size() && p[yp[j]].y-p[yp[i]].y<minum; j++)// 主要的不同之处，产生最小距离的点对必须属于不同类别if (p[yp[j]].flag != p[yp[i]].flag)minum = min(minum, dist(p[yp[j]], p[yp[i]]));return minum;}int main(){int i,j,test;cin >> test;while (test--){int xx,yy,n;cin >> n;for (i=0; i<n; i++){cin >> xx >> yy;p.push_back(Point(xx,yy));p[i].flag = 1;}for (; i<2*n; i++){cin >> xx >> yy;p.push_back(Point(xx,yy));p[i].flag = 2;}// 按照x坐标对点集进行排序sort(p.begin(), p.end(), cmp1);printf("%.3lf\n", min_dist(p, 0, p.size()-1));}}