分治算法求解凸包问题

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

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

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

相关概念

在这里插入图片描述

  • 凸多边形:连接多边形中任意两点的线段全部在多边形内部
  • 平面上某个点集的凸包:包含所有点的最小凸多边形
  • 凸包的表示:最小凸多边形的顶点
    1. 集合(无序)
    2. 序列(有序)如上图按顺时针:<q,r,s,t,p>
  • 点和线段的方向关系:
    给定线段 A B AB AB和点 C C C,判断 C C C A B → \overrightarrow{AB} AB 的位置关系
    计算叉积: ( x B − x A , y B − y A ) × ( x C − x B , y C − y B ) (x_{B}-x_{A},y_{B}-y_{A})\times(x_{C}-x_{B},y_{C}-y_{B}) (xBxA,yByA)×(xCxB,yCyB)
    • 大于0: C C C A B → \overrightarrow{AB} AB 的左侧
    • 等于0: C C C A B → \overrightarrow{AB} AB 所在的直线上
    • 小于0: C C C A B → \overrightarrow{AB} AB 的右侧

凸包问题

给定⼆维平面上的n个点集 S = ( x i , y i ) ∣ i = 1 , 2 , . . . , n S={(x_{i},y_{i})|i=1,2,...,n} S=(xi,yi)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) O(n)个点, O ( n 3 ) O(n^{3}) O(n3)个三角形,检查需要 O ( 1 ) O(1) O(1)时间
  • 总的时间复杂度: O ( n 4 ) O(n^{4}) O(n4)
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^{2}) O(n2)条线段,测试需要 O ( n ) O(n) O(n)时间
  • 总的时间复杂度: O ( n 3 ) O(n^{3}) O(n3)
    在这里插入图片描述
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 m e r g e T(n)= T( n - 1)+ T_{merge} T(n)=T(n1)+Tmerge
    • T m e r g e = O ( n 2 ) ⇒ T ( n ) = O ( n 3 ) T_{merge} = O(n^{2}) \Rightarrow T(n) = O(n^{3}) Tmerge=O(n2)T(n)=O(n3)
    • T m e r g e = O ( n ) ⇒ T ( n ) = O ( n 2 ) T_{merge} = O(n) \Rightarrow T(n) = O(n^{2}) Tmerge=O(n)T(n)=O(n2)
    • T m e r g e = O ( log ⁡ n ) ⇒ T ( n ) = O ( n log ⁡ n ) T_{merge} = O(\log{n}) \Rightarrow T(n) = O(n\log{n}) Tmerge=O(logn)T(n)=O(nlogn)

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

在这里插入图片描述

  • 凸包的支撑线: 与凸包仅相交于一点的直线
  • 过凸包外一点有且仅有两条该凸包的支撑线 ( q p 3 qp_{3} qp3 q p 5 qp_{5} qp5)
  • 两个交点 ( p 3 p_{3} p3 p 5 p_{5} p5)将凸包边界分成两个链: 近链和远链
  • 新的凸包由远链和点q决定
    • 从左支撑线开始,在凸包边界上顺时针前进直到右支撑线

如何求左右支撑线:

暴力求解:

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

优化搜索测略:

  • 寻找左支撑线
    • 从距离 q q q最近的 p i p_{i} pi开始,检查 p ( i − 1 ) % n p_{(i-1)\%n} p(i1)%n p ( i + 1 ) % n p_{(i+1)\%n} p(i+1)%n是否都在 q p i qp_{i} qpi右侧
  • 右侧寻找右支撑线
    • 从距离 q q q最近的 p i p_{i} pi开始,检查 p ( i − 1 ) % n p_{(i-1)\%n} p(i1)%n p ( i + 1 ) % n p_{(i+1)\%n} p(i+1)%n是否都在 q p i qp_{i} qpi左侧
  • T ( n ) = T ( n − 1 ) + O ( n ) ⇒ T ( n ) = O ( n 2 ) T(n) = T(n - 1) + O(n)\Rightarrow T(n) = O(n^{2}) T(n)=T(n1)+O(n)T(n)=O(n2)
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 ) = 2 T ( n / 2 ) + T m e r g e T(n) = 2T(n/2)+ T_{merge} T(n)=2T(n/2)+Tmerge
    IF T m e r g e = O ( n ) ⇒ T ( n ) = O ( n l o g n ) T_{merge} = O(n) \Rightarrow T(n) = O(nlogn) Tmerge=O(n)T(n)=O(nlogn)

合并过程设计
在这里插入图片描述

  • 凸包P和Q的桥: 既是P的支撑线,也是Q的支撑线 (比如 a 4 b 1 a_{4}b_{1} a4b1,和 a 2 b 1 a_{2}b_{1} a2b1)
    • 上桥: 如果P和O都在桥的下方
    • 下桥: 如果P和O都在桥的上方
  • 找到P和O的上桥和下桥
    • 每个凸包边界被其与桥的交点分成左链和右链
  • 合并后的凸包=上桥 + P的左链 + 下桥 + O的右链

寻找上下桥

  • a i a_{i} ai是P中距离L最近的点, b j b_{j} bj是Q中距离L最近的点
  • 如果 a i b j a_{i}b_{j} aibj不是 a i a_{i} ai的左支撑线,找到 a i a_{i} ai的左支撑线 a i b j a_{i}b_{j} aibj
    • 顺时针检查 b j b_{j} bj的下一个点
  • 如果 a i b j a_{i}b_{j} aibj不是 b j b_{j} bj的右支撑线,找到 b j b_{j} bj的右支撑线 a i b j a_{i}b_{j} aibj
    • 逆时针检查 a i a_{i} ai的下一个点
  • 重复上述步骤直到 a i b j a_{i}b_{j} aibj既是 a i a_{i} ai的左支撑线,又是 b j b_{j} bj的右支撑线,即为上桥
  • 下桥类似求解
  • 寻找上桥和下桥时间复杂度: O(n)
  • 合并时间复杂度: O(n):
    T ( n ) = 2 T ( n / 2 ) + O ( n ) ⇒ T ( n ) = O ( n l o g n ) T(n) = 2T(n/2) + O(n) \Rightarrow T(n) = O(nlogn) T(n)=2T(n/2)+O(n)T(n)=O(nlogn)

快速凸包

在这里插入图片描述

  • 按照x轴排序后选择最大和最小两个极端点a和b (一定属于凸包)
  • a b → \overrightarrow{ab} ab 将整个点集分成两部分: a b → \overrightarrow{ab} ab 左边的点集A、 a b → \overrightarrow{ab} ab 右边的点集B
  1. A中哪些点一定属于凸包?
  • 距离 a b → \overrightarrow{ab} ab 最远的点c: O(1)时间可确定c到 a b → \overrightarrow{ab} ab 的距离
  • △ a b c \bigtriangleup abc abc 中的点一定不属于凸包
  • △ a b c \bigtriangleup abc abc 外的点?
    • b c → \overrightarrow{bc} bc 右边的点: 递归求解哪些点属于凸包
    • c a → \overrightarrow{ca} ca 右边的点: 递归求解哪些点属于凸包
  1. B中哪些点一定属于凸包?

在这里插入图片描述

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

在这里插入图片描述
quick_half_hull ( S , b a → ) (S, \overrightarrow{ba}) (S,ba ):
if ( s = ⊘ ) (s=\oslash ) (s=) return ⊘ \oslash
c ← furthest ( S , b a → ) (S, \overrightarrow{ba}) (S,ba )
A ← right_of ( S , b c → ) (S, \overrightarrow{bc}) (S,bc )
B ← right_of ( S , c a → ) (S, \overrightarrow{ca}) (S,ca )
Q A Q_{A} QA ← quick_half_hull ( S , b c → ) (S, \overrightarrow{bc}) (S,bc )
Q B Q_{B} QB ← quick_half_hull ( S , c a → ) (S, \overrightarrow{ca}) (S,ca )
return Q A ∪ { c } ∪ Q B Q_{A} ∪ \{c\} ∪ Q_{B} QA{c}QB

时间复杂度分析
∣ S ∣ = n , ∣ A ∣ = α , ∣ B ∣ = β , α + β ≤ n − 1 |S| =n,|A|=\alpha,|B|=\beta,\alpha+\beta\le n-1 S=nA=α,B=β,α+βn1

  • quick_half_hull
    • T ( n ) = T ( α ) + T ( β ) + O ( n ) T(n)=T({\alpha})+T(\beta)+O(n) T(n)=T(α)+T(β)+O(n)
      1. α = β = n 2 ⇒ T ( n ) = O ( n l o g n ) \alpha=\beta=\frac{n}{2} \Rightarrow T(n)=O(nlogn) α=β=2nT(n)=O(nlogn)
      2. α = 0 , β = n − 1 ⇒ T ( n ) = O ( n 2 ) \alpha=0,\beta=n-1 \Rightarrow T(n)=O(n^{2}) α=0,β=n1T(n)=O(n2)
  • quick_hull
    • 预排序计算极端点: O ( n l o g n ) O(nlogn) O(nlogn)
    • 计算A和B: O ( n ) O(n) O(n)
    • 递归求解A和B: < 2 T ( n ) <2T(n) <2T(n)
    • 最坏情况下也是 O ( n 2 ) O(n^{2}) O(n2)
#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;
}
  • 0
    点赞
  • 31
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 5
    评论
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

恭仔さん

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值