分治算法求解凸包问题

分治算法（Divide and Conquer）是一种解决问题的算法设计策略，它将一个大问题分解成若干个规模较小且相互独立的子问题，然后将这些子问题的解合并起来，从而得到原问题的解。

分治算法通常包括以下三个步骤：

• 分解：将原问题分解为一组子问题，子问题与原问题类似，但是规模更小
• 解决：递归求解子问题，如果子问题足够小，停止递归，直接求解
• 合并：将子问题的解组合成原问题的解

相关概念

• 凸多边形：连接多边形中任意两点的线段全部在多边形内部
• 平面上某个点集的凸包：包含所有点的最小凸多边形
• 凸包的表示：最小凸多边形的顶点
  1. 集合（无序）
  2. 序列（有序）如上图按顺时针:<q,r,s,t,p>
• 点和线段的方向关系:
  给定线段$AB$和点$C$，判断$C$与$\overrightarrow{AB}$的位置关系
  计算叉积：$(x_B-x_A,y_B-y_A)\times (x_C-x_B,y_C-y_B)$
  • 大于0： $C$ 在 $\overrightarrow{AB}$ 的左侧
  • 等于0： $C$ 在 $\overrightarrow{AB}$ 所在的直线上
  • 小于0： $C$ 在 $\overrightarrow{AB}$ 的右侧

凸包问题

给定⼆维平面上的n个点集$S=(x_i,y_i)|i=1,2,...,n$ ，找到其凸包

• 假设没有两点具有相同的横坐标
• 假设没有两点具有相同的纵坐标
• 假设没有三点共线

//点结构
struct Point {
    int x, y;
};
//边结构
struct Line {
    Point start, end;
};
//判断点和线段的方向关系
int crossProduct(Point A, Point B, Point C) {
    return (B.x - A.x) * (C.y - B.y) - (B.y - A.y) * (C.x - B.x);
}

1.穷举法求凸包

点穷举

由于凸包是点集S最小的凸多边形的顶点或边，我们可以从几何的角度进行观察，对于每⼀个点，如果它在其他任意三点构成的三角形中，那么它就被凸包所围，换句话说，它不可能是构成最小凸包的点。而那些构成凸包的点则不会处在任意三角形的内部。

根据这个性质，我们可以对点进行穷举：

• 对于每⼀个点，检查其是否在任意其它三点构成的三角形中
• $O(n)$个点， $O(n^3)$个三角形，检查需要$O(1)$时间
• 总的时间复杂度：$O(n^4)$

bool isInsideTriangle(Point A, Point B, Point C, Point P) {
    int crossABP = crossProduct(A, B, P);
    int crossBCP = crossProduct(B, C, P);
    int crossCAP = crossProduct(C, A, P);

    return (crossABP >= 0 && crossBCP >= 0 && crossCAP >= 0) ||
           (crossABP <= 0 && crossBCP <= 0 && crossCAP <= 0);
}

---

边穷举

根据凸包定义，我们也可以通过两个点组成的边来判断其是否属于凸包：

• 属于：所有其它的点都在该线段的⼀侧
• $O(n^2)$条线段，测试需要$O(n)$时间
• 总的时间复杂度：$O(n^3)$
  

bool isConvexHullEdge(vector<Point>& points, Line edge) {
	//遍历点
    for (Point p : points) {
        if (p != edge.start && p != edge.end) {
            int crossProductValue = crossProduct(edge.start, edge.end, p);
            if (crossProductValue > 0) {
                return false;
            }
        }
    }
    return true;
}

---

2.分治法求凸包

对于点集的分解目前有两种分解策略：

• 1/n-1分：分为An-1个点和B1个点
• 二等分：分为A/B两个大小相等的点集

我们先考虑第一种情况，即1/n-1分治

插入凸包

• 分解:把n个点分成A、B两部分，其中A有n -1个点，B只有1个点
  1. 选取一个特别的点q放入B: 最右边的点 (横坐标最大的点)
  2. 该点一定属于新的凸包
• 递归求解A的凸包CH(A)
  • 基本情况: 三个点，直接计算
• 合并CH(A) 和 q
• 算法时间复杂度: $T(n)=T(n-1)+T_{merge}$
  • $T_{merge}=O(n^2)\Rightarrow T(n)=O(n^3)$
  • $T_{merge}=O(n)\Rightarrow T(n)=O(n^2)$
  • $T_{merge}=O(\log n)\Rightarrow T(n)=O(n\log n)$

可以看到，降低插入凸包的复杂度关键是降低Merge时的复杂度。

• 凸包的支撑线: 与凸包仅相交于一点的直线
• 过凸包外一点有且仅有两条该凸包的支撑线 ($q{p}_3$和$q{p}_5$)
• 两个交点 ($p_3$和$p_5$)将凸包边界分成两个链: 近链和远链
• 新的凸包由远链和点q决定
  • 从左支撑线开始，在凸包边界上顺时针前进直到右支撑线

如何求左右支撑线:

暴力求解：

• $O(n)$条候选支撑线$q{p}_i$;
• 检查所有其他点是否在$q{p}_i$同一侧: $O(n)$时间
• 共需$O(n^2)$时间

优化搜索测略：

• 寻找左支撑线
  • 从距离$q$最近的$p_i$开始，检查$p_{(i-1)\%n}$和$p_{(i+1)\%n}$是否都在$q{p}_i$右侧
• 右侧寻找右支撑线
  • 从距离$q$最近的$p_i$开始，检查$p_{(i-1)\%n}$和$p_{(i+1)\%n}$是否都在$q{p}_i$左侧
• $T(n)=T(n-1)+O(n)\Rightarrow T(n)=O(n^2)$

function ConvexHull(S):
    if |S| <= 3:
        return ComputeConvexHull(S)
    else:
        Divide S into A and B
        q = rightmost point in B
        CH_A = ConvexHull(A)
        Merge CH_A with q to form the new convex hull CH

function Merge(CH_A, q):
    Find left and right support lines for CH_A and q
    Traverse the boundary of CH_A from left support line to right support line and update CH

function ComputeConvexHull(S):
    // Implement a convex hull algorithm for small inputs (e.g., Gift Wrapping)

// Initial call to ConvexHull with the entire set of points S
ConvexHull(S)

---

并归凸包

• 分解: 把n个点分成A、B两部分，其中A、B各有n/2个点
• 预处理: 将所有点按横坐标排序
• 递归求解A和B的凸包，记为P和O
  • 基本情况:是不是三个点?
• 合并P和Q
• 算法时间复杂度: $T(n)=2T(n/2)+T_{merge}$
  IF $T_{merge}=O(n)\Rightarrow T(n)=O(nlogn)$

合并过程设计

• 凸包P和Q的桥: 既是P的支撑线，也是Q的支撑线 (比如$a_4{b}_1$,和$a_2{b}_1$)
  • 上桥: 如果P和O都在桥的下方
  • 下桥: 如果P和O都在桥的上方
• 找到P和O的上桥和下桥
  • 每个凸包边界被其与桥的交点分成左链和右链
• 合并后的凸包=上桥 + P的左链 + 下桥 + O的右链

寻找上下桥

• 令$a_i$是P中距离L最近的点，$b_j$是Q中距离L最近的点
• 如果$a_i{b}_j$不是$a_i$的左支撑线，找到$a_i$的左支撑线$a_i{b}_j$
  • 顺时针检查$b_j$的下一个点
• 如果$a_i{b}_j$不是$b_j$的右支撑线，找到$b_j$的右支撑线$a_i{b}_j$
  • 逆时针检查$a_i$的下一个点
• 重复上述步骤直到$a_i{b}_j$既是$a_i$的左支撑线，又是$b_j$的右支撑线，即为上桥
• 下桥类似求解
• 寻找上桥和下桥时间复杂度: O(n)
• 合并时间复杂度: O(n):
  $T(n)=2T(n/2)+O(n)\Rightarrow T(n)=O(nlogn)$

---

快速凸包

• 按照x轴排序后选择最大和最小两个极端点a和b (一定属于凸包)
• $\overrightarrow{ab}$将整个点集分成两部分: $\overrightarrow{ab}$左边的点集A、$\overrightarrow{ab}$右边的点集B

1. A中哪些点一定属于凸包?

• 距离 $\overrightarrow{ab}$ 最远的点c: O(1)时间可确定c到$\overrightarrow{ab}$的距离
• $△abc$ 中的点一定不属于凸包
• $△abc$ 外的点?
  • $\overrightarrow{bc}$右边的点: 递归求解哪些点属于凸包
  • $\overrightarrow{ca}$右边的点: 递归求解哪些点属于凸包

2. B中哪些点一定属于凸包?

quick_hull($S$) :
$(a,b)$ ← extreme_points($S$)
$A$ ← right_of $(S,\overrightarrow{ba})$
$B$ ← right_of $(S,\overrightarrow{ab})$
$Q_A$ ← quick_half_hull$(A,\overrightarrow{ba})$
$Q_B$ ← quick_half_hull$(A,\overrightarrow{ab})$
return { a } ∪ Q A ∪ { b } ∪ Q B \{a\} ∪ Q_{A} ∪ \{b\} ∪ Q_{B} {a}∪QA​∪{b}∪QB​

quick_half_hull$(S,\overrightarrow{ba})$:
if $(s=\oslash )$ return $\oslash$
c ← furthest$(S,\overrightarrow{ba})$
A ← right_of$(S,\overrightarrow{bc})$
B ← right_of$(S,\overrightarrow{ca})$
$Q_A$ ← quick_half_hull$(S,\overrightarrow{bc})$
$Q_B$ ← quick_half_hull$(S,\overrightarrow{ca})$
return Q A ∪ { c } ∪ Q B Q_{A} ∪ \{c\} ∪ Q_{B} QA​∪{c}∪QB​

时间复杂度分析
令 $|S|=n，|A|=\alpha ,|B|=\beta ,\alpha +\beta \le n-1$

• quick_half_hull
  • $T(n)=T(\alpha )+T(\beta )+O(n)$
    1.$\alpha =\beta =\frac{n}{2}\Rightarrow T(n)=O(nlogn)$
    2.$\alpha =0,\beta =n-1\Rightarrow T(n)=O(n^2)$
• quick_hull
  • 预排序计算极端点：$O(nlogn)$
  • 计算A和B：$O(n)$
  • 递归求解A和B：$<2T(n)$
  • 最坏情况下也是$O(n^2)$

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

using namespace std;

// 定义点的结构体
struct Point {
    int x, y;
};

// 计算点P到直线AB的距离
double distanceToLine(const Point& A, const Point& B, const Point& P) {
    return abs((B.x - A.x) * (A.y - P.y) - (A.x - P.x) * (B.y - A.y)) /
        sqrt(pow(B.x - A.x, 2) + pow(B.y - A.y, 2));
}

// 计算叉积，判断点C在向量AB的右侧、左侧或共线
int crossProduct(const Point& A, const Point& B, const Point& C) {
    return (B.x - A.x) * (C.y - B.y) - (B.y - A.y) * (C.x - B.x);
}

// 判断点C是否在向量AB的右侧
bool isRightOf(const Point& A, const Point& B, const Point& C) {
    return crossProduct(A, B, C) < 0;
}

// 快速凸包算法的部分实现
vector<Point> quick_half_hull(vector<Point>& S, const Point& a, const Point& b) {
    if (S.empty()) return {};

    double maxDist = -1;
    Point c;
    for (const auto& point : S) {
        double dist = distanceToLine(a, b, point);
        if (dist > maxDist) {
            maxDist = dist;
            c = point;
        }
    }

    vector<Point> A, B;

    // 筛选出右侧的点集A和B
    for (const auto& point : S) {
        if (isRightOf(b, c, point)) {
            A.push_back(point);
        }
        else if (isRightOf(c, a, point)) {
            B.push_back(point);
        }
    }

    vector<Point> QA = quick_half_hull(A, a, c);
    vector<Point> QB = quick_half_hull(B, c, b);

    vector<Point> result;
    result.insert(result.end(), QA.begin(), QA.end());
    result.push_back(c);
    result.insert(result.end(), QB.begin(), QB.end());

    return result;
}

vector<Point> quick_hull(vector<Point>& S) {
    if (S.size() < 3) return S;

    // 寻找极端点a和b
    sort(S.begin(), S.end(), [](const Point& p1, const Point& p2) {
        return (p1.x < p2.x) || (p1.x == p2.x && p1.y < p2.y);
        });

    Point a = S[0];
    Point b = S[S.size() - 1];

    vector<Point> A, B;

    // 筛选出右侧的点集A和B
    for (const auto& point : S) {
        if (isRightOf(b, a, point)) {
            A.push_back(point);
        }
        else if (isRightOf(a, b, point)) {
            B.push_back(point);
        }
    }

    vector<Point> QA = quick_half_hull(A, a, b);
    vector<Point> QB = quick_half_hull(B, b, a);

    vector<Point> result;
    result.push_back(a);
    result.insert(result.end(), QA.begin(), QA.end());
    result.push_back(b);
    result.insert(result.end(), QB.begin(), QB.end());

    return result;
}

int main() {
    // 示例点集
    vector<Point> points = { {0, 0}, {4, 4}, {4, 0}, {0, 4}, {3, 3} ,{3,1},{2,1},{1,2},{3,2},{4,3},{2,4} ,{0,2},{0,5} };

    // 调用凸包算法
    vector<Point> result = quick_hull(points);

    // 输出凸包顶点
    cout << "Convex Hull Points:\n";
    for (const auto& point : result) {
        cout << "(" << point.x << ", " << point.y << ")\n";
    }

    return 0;
}