问题描述

有n个三维地理信息数据：。
要求：距离最近的局部高点。
试采用递归分治的思想设计算法，并以不同规模的数据输入，与非递归分治算法比较，作算法时间复杂度分析。

子问题1. 求所有局部高点

（1）求所有局部高点
局部高点的坐标输入：把平面坐标 规定在规则连续的网格区域上（如图）， 坐标存在差异。但注意网格点数据输入时仍是无序的，因此，需要先对 和 坐标的排序。
局部高点的条件：在上述网格规定下，一般考虑某一点与周围的八个点之间的关系（四个相邻角点、或 向紧相邻四点）。可以只通过其中四点来定义高点条件（四个角点，或紧相邻x-y向的四个点，或者都考虑，说明即可）。

解决方案

定义局部高点

定义局部高点为一个点的高度大于他上下左右四个点的高度

离散化处理

注意到数据是以网格区域的形式给出, 那么根据网格的特征, 必定会有很多点的 坐标和 坐标处于同一水平线上, 由于给出的数据是浮点数, 为了方便处理, 我们可以对数据先进行离散化处理, 把其从浮点数的值域映射到整数值域, 然后找局部高点.
具体离散化的方法是:

• 对 坐标进行排序, 去重后在原数组上二分搜索后赋值
• 对 坐标进行排序, 去重后在原数组上二分搜索后赋值

这样我们得到了一个网格图, 表示 行 列的点的高度(即 坐标第 小, 坐标第 小)

查找局部高点

此时, 只需要遍历一遍点集(网格图), 对于每一个点, 找到它周围的四个点, 并判断它们的高度关系即可

复杂度分析

因为在建立网格图的过程中需要进行排序和离散化处理， 令 为点集的数目, 那么时间复杂度是 , 此外, 对于每一个点集都需要遍历一次并判断周围四个点, 时间复杂度是 , 于是总体的时间复杂度是 的, 空间复杂度是

Code

/*
  autor: Kurisu
  2020年4月28日15:48:12
*/
#include<bits/stdc++.h>
using namespace std;

const long long inf = 1e18;
const int N = 5e5 + 5;
const double eps = 1e-10;
const int mod = 1e9 + 7;
typedef long long ll;

struct Point { // 定义点的结构体
  double x, y, z;
  Point(double _x = 0, double _y = 0, double _z = 0): x(_x), y(_y), z(_z) {}
}point[N];
vector<pair<int, int> > ans; // 存储局部高点的向量
double p1[N], p2[N];
int up[N]; // 离散化的辅助数组
int Point_num; // 点数
int n, m;
int dist[4][2] = {0, 1, 0, -1, 1, 0, -1, 0}; // 4个方向

double solve(vector<vector<double>> &Graph) {
  for(int i = 1; i <= n; i++) {
    for(int j = 1; j <= up[i]; j++) { // 遍历每一个点
      int flag = 0;
      for(int k = 0; k < 4; k++) { // 四个方向
        int xx = i + dist[k][0];
        int yy = j + dist[k][1];
        if(xx < 0 || xx > n || yy < 0 || yy > up[xx] || Graph[xx][yy] == -1e9) continue; // 超出这个范围或者这个点不存在
        double zz = Graph[xx][yy];
        if(zz >= Graph[i][j]) { // 出现一个点比它高 跳出
          flag = 1;
          break;
        }
      }
      if(!flag) { // 周围四个点没有一个比它高的, 是局部高点
        ans.push_back({i, j});
      }
    }
  }
  return 0;
}
vector<vector<double>> Discretization() { // 离散化
  sort(p1 + 1, p1 + 1 + Point_num); // x 坐标排序
  sort(p2 + 1, p2 + 1 + Point_num); // y 坐标排序
  int len1 = unique(p1 + 1, p1 + 1 + Point_num) - p1 - 1; // x 坐标去重
  int len2 = unique(p2 + 1, p2 + 1 + Point_num) - p2 - 1; // y 坐标去重
  for(int i = 1; i <= Point_num; i++) {
    point[i].x = lower_bound(p1 + 1, p1 + 1 + len1, point[i].x) - p1; // 对 x 坐标离散化
    point[i].y = lower_bound(p2 + 1, p2 + 1 + len2, point[i].y) - p2; // 对 y 坐标离散化
  }
 for(int i = 1; i <= Point_num; i++) n = max(n, int(point[i].x)); // 离散化后 x 坐标最大值
 for(int i = 1; i <= Point_num; i++) up[int(point[i].x)] = max(up[int(point[i].x)], int(point[i].y)); // 离散化后 y 坐标最大值
  vector<vector<double> > Graph(n + 1); // 二维向量存储图
  for(int i = 1; i <= n; i++) {
    Graph[i].resize(up[i] + 1); // x 坐标为 i 的向量容量扩充到这一行 y坐标最大值
  }
  for(int i = 1; i <= n; i++) {
    for(int j = 1; j <= up[i]; j++) {
      Graph[i][j] = -1e9; // 初始化 不存在的点高度为无穷小
    }
  }
  for(int i = 1; i <= Point_num; i++) {
    Graph[int(point[i].x)][int(point[i].y)] = point[i].z; // 赋值
  }
  return Graph;
}
void input() { // 输入
  cin >> Point_num;
  for(int i = 1; i <= Point_num; i++) {
    cin >> point[i].x >> point[i].y >> point[i].z;
    p1[i] = point[i].x, p2[i] = point[i].y;
  }
}
int main () {
  freopen("5.in", "r", stdin);
  input();
  vector<vector<double> > Graph = Discretization(); // 对图进行离散化
  solve(Graph);
  //cout << "局部高点的坐标为:\n";
  for(int i = 0; i < ans.size(); i++) {
    int x = ans[i].first, y = ans[i].second; // x 坐标和 y 坐标
    cout << fixed << setprecision(2) << p1[x] << ' ' << p2[y] << ' ' << Graph[x][y] << "\n";
  }
  return 0;
}

子问题2. 求距离最近点对

（与局部高点问题无关)
求如图不规则无序点集中的最近点对。
给两种选项：
A. 最近距离只需按照平面距离确定（只考虑xy坐标）。（相当于实践标准的最近点对算法）
B. 考虑三维空间点的距离（计算距离时需考虑z坐标）。（相当于对标准算法增加了一维）
要求完成A、B之一。最起码按A来实施。

问题A的解决方案

朴素做法

这个问题是经典的求平面最近点对问题, 对于这个问题, 我们有朴素的做法:

double solve(){
  double ans = inf; // 无穷大
  for(int i = 1; i <= Point_num; i++){
    for(int j = i + 1; j <= Point_num; j++){
      ans = min(ans, distance(point[i], point[j]));
    }
  }
  return ans;
}

分析上述算法, 不妨视distance函数和min函数为常数级别的运算, 则设总的点数为
总的计算次数是
那么时间复杂度是级别的

分治做法

注意到我们是求 区间的平面最近点对, 那么我们从分治的思想去考虑,
能不能转化为求 和 的最近点对呢?
但是这样做有一个问题, 最近点对不一定单独出现在某个区间
也就是说, 可能存在一个最近点对, 它的一个点在, 一个点在
所以, 我们在分治的时候还需要对这一步进行检验, 怎么检验呢？
我们先对点集按 坐标升序排序
我们求得 和 的最近点对距离 和 后
取， 在中找到跟 的 距离小于 的点集
再把 中的点按 坐标排序, 对于当前点集, 我们直接统计两两点对的最小距离即可

Code

/*
author: Kurisu
2020年4月18日10:27:40
*/
#include <bits/stdc++.h>
using namespace std;

const int N = 2e5 + 5;
const long long inf = 1e18; // 无限大
struct Point{ // 定义点的结构体
  double x, y;
  Point(double _x = 0, double _y = 0): x(_x), y(_y){}
}point[N], tmp[N];
int Point_num; // 点的数量
void input(){
  cin >> Point_num; // 输入点的数量
  for(int i = 1; i <= Point_num; i++){
    cin >> point[i].x >> point[i].y; // 按顺序输入每个点的 x, y坐标
  }
}
double distance(Point p1, Point p2){
  return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)); // 计算两点的距离
}
double solve(int l, int r){
  double d = inf;
  if(l == r) return d; // 只有一个点
  else if(r - l == 1){ // 两个点结束分治
    return distance(point[l], point[r]);
  }
  int mid = l + r >> 1;
  d = min(solve(l, mid), solve(mid + 1, r)); // 分治
  int num = 0;
  for(int i = l; i <= r; i++){
    if(fabs(point[i].x - point[mid].x) < d){
      tmp[++num] = point[i]; // 找到离中线距离小于 d 的点
    }
  }
  sort(tmp + 1, tmp + 1 + num, [](const Point &x, const Point &y){return x.y < y.y;});
  // STL 内置排序算法 和 lamda 表达式重载比较运算符, 按y坐标排序
  for(int i = 1; i <= num; i++){
    for(int j = i + 1; j <= num && fabs(tmp[j].y - tmp[i].y) < d; j++){
      d = min(d, distance(tmp[j], tmp[i])); // 剩余的点里找最小
    }
  }
  return d;
}
int main(){
  input(); // 输入点
  sort(point + 1, point + 1 + Point_num, [](const Point &x, const Point &y){return x.x < y.x;});
  // STL 内置排序算法 和 lamda 表达式重载比较运算符
  cout << fixed << setprecision(4) << solve(1, Point_num) << "\n"; // 输出最近距离, 保留四位小数
  return 0;
}

在上面代码中有一步按 坐标排序后找最小值的过程
可以看到它有两层循环, 但其实实际的运算次数很少, 可以看做是一个常数
这里按照老师发的算法导论简单证明下:


如上图, 我们在某一递归过程中出现最近点对的两个端点分别处于不同区域, 即一个端点属于 , 另一个属于
那么我们沿着中线 画出一个长 , 宽 的小矩形区域, 对于这个矩形区域, 最多会有八个点, 这个是可以证明的
因为我们知道当前的最小距离是 , 那么对于左边的边长为 的正方形, 任意两点的距离都不会小于
显然, 最多的情况不超过四个点, 怎么证明呢?
这里我们可以用反证法:
在上述的 边长为 的小正方形中, 我们划分成四个 的小正方形, 若这个区域存在至少五个点, 由鸽巢原理
得至少有两个点位于同一个小正方形, 而两者最远的距离为对角线距离为 矛盾, 因此这个区域的点不大于四个
同理, 对于右边的小正方形, 最多也只有四个点
那么总的这个矩形区域的点最多不会超过八个点
因此, 我们在上述的二重循环中最多只会做每个点之后的七个点, 那么这个操作可以看做是在常数时间内完成的
于是, 该算法的分割步骤和合并步骤总共耗时是的, 排序的耗时是, 而对 排序耗时为常数级别
因此算法耗费时间 满足递归方程
可得, 故总体时间复杂度是

问题B的解决方案

我们还是从朴素的做法开始入手

朴素做法

最简单的做法是直接对每两个点求距离, 在距离中取最小值

double solve() {
  double ans = 1e18;
  for(int i = 1; i <= Point_num; i++) {
    for(int j = i + 1; j <= Point_num; j++) {
      ans = min(ans, distance(point[i], point[j]));
    }
  }
  return ans;
}

把 和 函数看作是常数级别的运算, 则设总的点数为
总的计算次数是
那么时间复杂度是级别的

分治做法——从二维扩展到三维[^1]

我们都知道二维平面最近点对的分治做法, 那么我们能否将这个算法推广到三维的情况呢?
还是那个思路, 在求 1 到 的最近点对问题, 可以转化为求 到 的最近点对
同样的思路, 我们先对点按 坐标排序, 截取一个平面 ， 对于平面的左右部分分治
找到左右两端点集中最近距离
此时同样会出现最近点对在分界面两边的情况, 比较难处理的是合并部分
我们考虑位于平面 左右两端的点， 找出左端距离平面小于 的点加入到点集
找出右端距离平面小于的点加入到点集
那么点集 到 的距离都需要我们去验证是否为最近点对
最坏情况下, 和 的点集规模都是 的, 直接找复杂度是
但是, 考虑遍历 的每一个点, 都可以对它建立一个 的长方体, 所有可选方案都将包含在长方体中
因此可选点集具有稀疏的特点, 可以证明, 最多不超过 个点存在于这个长方体中
证明思路依然是利用鸽巢原理：
把长方体分割成 个大小相等的 的小长方体,
若长方体中有多于 个点, 那么至少有一个小长方体里有两个点,
它们最远的距离是
这与 是最近点对距离矛盾, 因此大长方体里最多有 个点
于是我们只需要对 和 按 坐标为第一关键字, 坐标为第二关键字,
利用尺取法的思想, 每一次取一段满足条件的区域, 如果遍历超过 个点就跳出

复杂度分析

算法耗费时间 满足递归方程
可得, 故总体时间复杂度是

[^1]: 参考论文 分治法实现最接近点对问题的三维推广算法——张晓红, 三维空间的最接近点对问题——刘志辉,

Code

/*
  autor: Kurisu
  2020年4月23日10:27:06
*/
#include<bits/stdc++.h>
using namespace std;

const long long inf = 1e18;
const int N = 1e5 + 5;

struct Point { // 定义点集
  double x, y, z;
}point[N], tmp1[N], tmp2[N];
int Point_num;
void Input() { // 输入数据
  cin >> Point_num; // 点的数目
  for(int i = 1; i <= Point_num; i++) {
    cin >> point[i].x >> point[i].y >> point[i].z; // 输入点
  }
}
double distance(Point x, Point y) { // 计算距离
  double ans = (x.x - y.x) * (x.x - y.x) + (x.y - y.y) * (x.y - y.y) + (x.z - y.z) * (x.z - y.z);
  return sqrt(ans);
}
double solve(int l, int r){ // 分治
  double d = inf;
  if(l == r) return d; // 只有一个点
  else if(r - l == 1) { // 只有两个点, 直接计算两点距离
    return distance(point[l], point[r]);
  }
  int mid = l + r >> 1; // 找到中间的 mid 平面
  d = min(solve(l, mid), solve(mid + 1, r)); // 分治
  int num1 = 0, num2 = 0;
  for(int i = l; i <= mid; i++) {
    if(fabs(point[i].x - point[mid].x) < d) { // mid 平面的左边 离平面距离小于 d 的点
      tmp1[++num1] = point[i];
    }
  }
  for(int i = mid + 1; i <= r; i++) {
    if(fabs(point[i].x - point[mid].x) < d) { // mid 平面的右边 离平面距离小于 d 的点
      tmp2[++num2] = point[i];
    }
  }
  sort(tmp1 + 1, tmp1 + 1 + num1, [](const Point &x, const Point &y){return x.y == y.y ? x.z < y.z : x.y < y.y;}); // 对tmp1数组以 y 为第一关键字 z 为第二关键字升序排序
  sort(tmp2 + 1, tmp2 + 1 + num2, [](const Point &x, const Point &y){return x.y == y.y ? x.z < y.z : x.y < y.y;}); // 对tmp2数组以 y 为第一关键字 z 为第二关键字升序排序
  int pl = 1, pr = 1; // 滑动窗口
  for(int i = 1; i <= num1; i++) {
    if(tmp1[i].y < tmp2[pl].y) {  // 点在窗口的额左边, 此时只需要与左端点比较
      int cur = pl;
      while(cur <= num2 && fabs(tmp2[cur].y - tmp1[cur].y) < d) {
        d = min(d, distance(tmp1[i], tmp2[cur]));
        cur++;
      }
    }
    while(pr <= num2 && pr < tmp1[i].y) pr++; // 滑动右端点
    while(pl <= num1 && tmp2[pl].y < tmp1[i].y && fabs(tmp2[pl].y - tmp1[i].y) > d) pl++; // 滑动左端点
    while(pr <= num2 && fabs(tmp2[pr].y - tmp1[i].y) < d) pr++; // 滑动右端点
    if(pr > pl) {
      int cur = 0;
      for(int j = pl; j < pr; j++) {
        d = min(d, distance(tmp1[i], tmp2[j]));
        ++cur;
        if(cur == 24) break; // 只找24个点
      }
    }
  }
  return d;
}
int main() {
  //clock_t start,endz;
  //start = clock();
  //freopen("4.in", "r", stdin);
  Input(); // 输入数据
  sort(point + 1, point + 1 + Point_num, [](const Point &x, const Point &y){return x.x < y.x;}); // 按 x 坐标升序排序
  cout << fixed << setprecision(4) << solve(1, Point_num) << "\n";
  //endz = clock();
  //cout << (double)(endz - start) / CLOCKS_PER_SEC * 1000 << "ms\n";
  return 0;
}