ACM / 数学 / 计算几何 · 2024-03-20

平面最近点对

分治唯一真神,拒绝人类智慧

考虑分治方法:与常规的分治算法一样,我们将这个有 个点的集合拆分成两个大小相同的集合 ,递归的解决最近点对问题。

这样分成了三个问题:

  1. 两个点都在中的最近点对
  2. 两个点分别在的点对
  3. 两个点都在中的最近点对

情况在集合大小为 时可以轻松解决,这样就只用考虑分治过程中合并的问题:先以 为关键字排序,进行分治。分治过程中以 为轴,将 的点加入待处理集合 。再按 为关键字,在递归过程中归并排序,算出满足 的每对点的距离更新答案。(可以证明在大小为 的矩形区域内的点不超过5个)

时间复杂度
分治方法Ⅰ:比起第二种较慢但是码量少一点

struct Point
{
    double x, y;
};
typedef vector<Point>::iterator Iter;
bool cmpx(const Point a, const Point b) { return a.x < b.x; }
bool cmpy(const Point a, const Point b) { return a.y < b.y; }
double dis(const Point a, const Point b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
void slv(const Iter l, const Iter r, double &d)
{
    if (r - l <= 1) return;
    vector<Point> Q;
    Iter t = l + (r - l) / 2;
    double w = t->x;
    slv(l, t, d), slv(t, r, d), inplace_merge(l, t, r, cmpy);
    for (Iter x = l; x != r; ++x) if (abs(w - x->x) <= d) Q.push_back(*x);
    for (Iter x = Q.begin(), y = x; x != Q.end(); ++x)
    {
        while (y != Q.end() && y->y <= x->y + d) ++y;
        for (Iter z = x + 1; z != y; ++z) d = min(d, dis(*x, *z));
    }
}
vector<Point> Poi; int n;
double preparata()
{
    double ans = 1e18;
    sort(Poi.begin(),Poi.end(), cmpx);
    slv(Poi.begin(),Poi.end(), ans);
    return ans;
}

分治方法Ⅱ:比第一种更快但是码量大一点

struct pt
{
    double x, y;
    int id;
};
bool cmpx(const pt a, const pt b) { return a.x < b.x; }
bool cmpy(const pt a, const pt b) { return a.y < b.y; }
int n, ansa, ansb; // 答案在这
vector<pt> a;     // 下标0开始
double mindist;
void upd_ans(const pt &a, const pt &b)
{
    double dist = sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
    if (dist < mindist)
        mindist = dist, ansa = a.id, ansb = b.id;
}
void rec(int l, int r)
{
    if (r - l <= 3) //点数小于3直接暴力找
    {
        for (int i = l; i <= r; ++i)
            for (int j = i + 1; j <= r; ++j)
                upd_ans(a[i], a[j]);
        sort(a.begin() + l, a.begin() + r + 1, cmpy);
        return;
    }
    int m = (l + r) >> 1;
    double midx = a[m].x;
    rec(l, m), rec(m + 1, r);
    //归并排序
    inplace_merge(a.begin() + l, a.begin() + m + 1, a.begin() + r + 1, cmpy);
    static pt t[N]; // 缓存数组
    int tsz = 0;
    for (int i = l; i <= r; ++i)
        if (abs(a[i].x - midx) < mindist) //可能是最近点对
        {
            for (int j = tsz - 1; j >= 0 && a[i].y - t[j].y < mindist; --j) upd_ans(a[i], t[j]);
            t[tsz++] = a[i];
        }
}
void preparata()
{
    sort(a.begin(), a.end(), cmpx);
    mindist = (ll)1e18;
    rec(0, n - 1);
}

非分治:更快(常数小)

int n;
double ans = 1e20;
struct Point
{
    double x, y;
    Point(double x = 0, double y = 0) : x(x), y(y) {}
};
struct cmp_y { bool operator()(const Point &a, const Point &b) const { return a.y < b.y; } };
void upd_ans(const Point &a, const Point &b)
{
    double dist = sqrt(pow((a.x - b.x), 2) + pow((a.y - b.y), 2));
    if (ans > dist) ans = dist;
}
Point a[N];//下标0开始
multiset<Point, cmp_y> s;
void preparata()
{
    sort(a, a + n, [&](Point a, Point b) -> bool { return a.x < b.x || (a.x == b.x && a.y < b.y); });
    for (int i = 0, l = 0; i < n; i++)
    {
        while (l < i && a[i].x - a[l].x >= ans) s.erase(s.find(a[l++]));
        for (auto it=s.lower_bound({a[i].x, a[i].y - ans}); it != s.end() && it->y - a[i].y < ans; it++)
            upd_ans(*it, a[i]);
        s.insert(a[i]);
    }
}