二维凸包算法 Julia实现

Foliciatarier

已于 2024-08-07 07:48:11 修改

阅读量899

点赞数 28

文章标签：算法 julia

于 2024-08-04 15:07:27 首次发布

本文链接：https://blog.csdn.net/Foliciatarier/article/details/140905490

版权

问题描述：给定平面上 $n$ 个点的集合 $Q$ ，求其子集 $P$ 构成 $Q$ 的凸包，即 $\forall p \in Q, \exist p_0, p_1, p_2 \in P$ 使得点 $p$ 在以点 $p_0, p_1, p_2$ 构成的三角形中或边上。

朴素枚举

初始设 $P := Q$ ，枚举 $p_0, p_1, p_2, p$ ，判断 $p$ 是否在以点 $p_0, p_1, p_2$ 构成的三角形中，是则删去 $p$ ，最后剩下的即为凸包。注意到纵坐标最小的点必在结果中，可将 $p_0$ 固定为之。

算法时间复杂度 $O(n^3)$ 。代码：

cross(p, q) = p[1] * q[2] - q[1] * p[2]

function NE(Q)
    function inΔ(p0, p1, p2, p)
        v0, v1, v2 = p .- p0, p1 .- p0, p2 .- p0
        c0, c1, c2 = cross(v1, v2), cross(v0, v1), cross(v0, v2)
        return c0 * c1 < 0 && c0 * c2 > 0 && abs(c2 - c1) < abs(c0)
    end
    m = argmin(last.(Q))
    p0, Q = Q[m], vcat(Q[begin : m - 1], Q[m + 1 : end])
    for i = length(Q) : -1 : 1
        f = false
        for j = 2 : length(Q)
            if i == j continue end
            for k = 1 : j - 1
                if i == k continue end
                if inΔ(p0, Q[j], Q[k], Q[i])
                    f = true; break
                end
            end
            if f break end
        end
        if f
            deleteat!(Q, i)
        end
    end
    sort!(Q, lt = (p, q) -> cross(p .- p0, q .- p0) > 0)
    return vcat(p0, Q)
end

Jarvis March 算法

按逆时针对 $P$ 排序， $p_0$ 固定为纵坐标最小的点，将问题转化成已知 $p_0, p_1, ..., p_k$ 如何求 $p_{k+1}$ 。

此时 $Q - \{p_k\}$ 中所有点都在向量 $p_k - p_{k-1}$ 的左侧，即向量集 $\{ p - p_k | p \in Q - \{ p_k \} \}$ 可在叉积运算的正负性上构成偏序集，最右的向量即为 $p_{k+1} - p_k$ 。直白点说，为了使得 $Q - \{p_{k+1}\}$ 中所有点都在向量 $p_{k+1} - p_k$ 的左侧，就要从向量集 $\{ p - p_k | p \in Q - \{ p_k \} \}$ 中选最右的向量作为 $p_{k+1} - p_k$ 。当 $p_{k+1} = p_0$ 时算法终止。

时间复杂度 $O (nh)$ ， $h$ 为凸包的点数。代码：

function Jarvis(Q)
    function march(A, cmp)
        r = A[begin]
        for i in A[2 : end]
            if cmp(i, r)
                r = i
            end
        end
        return r
    end
    P = [march(Q, (p, q) -> last(p) < last(q))]
    while true
        q = march(Q, (p, q) -> cross(p .- P[end], q .- P[end]) > 0)
        if P[begin] == q
            break
        end
        push!(P, q)
    end
    return P
end

Graham Scan 算法

$p_0$ 固定为纵坐标最小的点。按与 $p_0$ 构成向量的极角序对 $Q$ 排序，将问题转化为已知 $Q$ 中前 $k$ 个点的凸包 $P_k$ ，如何求 $Q$ 中前 $k + 1$ 个点的凸包 $P_{k+1}$ 。

不妨设 $P_k = \{ p_0, p_1, .., p_k \}$ ，注意到 $p_{k+1}$ 的在 $p_k - p_0$ 的左边，而 $P_k$ 中其它点在 $p_k - p_0$ 的右边，故一定有 $p_{k+1} \in P_{k+1}$ ，设 $P_{k+1} = P_k + \{ p_{k+1} \}$ 。判断 $p_{k+1}$ 在 $p_k - p_{k-1}$ 的哪边，若在右边则从凸包 $P_{k+1}$ 中抛出 $p_k$ ，继续迭代判断 $p_{k+1}$ 在 $p_{k-1} - p_{k-2}$ 的哪边，直到在左边则结束迭代，最后剩下的即为凸包 $P_{k+1}$ 。 $P_{\left| Q \right|}$ 即为所求的 $P$ 。

用栈存储中间过程的凸包计算， $Q$ 中每个点仅会进出栈一次，故扫描为线性时间，主要耗时在排序。时间复杂度 $\log n)$ 。代码：

function GrahamScan(Q)
    P = Q[1 : 3]
    for q in Q[4 : end]
        while cross(q .- P[end], P[end] .- P[end - 1]) > 0
            pop!(P)
        end
        push!(P, q)
    end
    return P
end

function Graham(Q)
    m = argmin(last.(Q))
    p0, Q = Q[m], vcat(Q[begin : m - 1], Q[m + 1 : end])
    Q = sort(Q, lt = (p, q) -> cross(p .- p0, q .- p0) > 0)
    return GrahamScan(vcat(p0, Q))
end

分治算法

将当前点集均分为左右两部分，分别用分治算法得到各自的凸包，再将两个凸包合并成一个凸包。

合并过程是，取一边凸包上某一点作为基点，求另一边凸包上极角最大和最小的点，则一边凸包的点按基点的极角序有序，另一边凸包的点分为两条从极角最小到极角最大的极角序有序点集，将三个极角序有序点集按三路归并成一个极角序有序的点集，并进行 Graham Scan 算法流程。

时间函数 $T(\frac{n}{2}) + O(n)$ ，得时间复杂度 $\log n)$ ，但常数很大。代码：

function DC(Q)
    function argMaxMin(A, cmp)
        u, v = 1, 1
        for i = 2 : length(A)
            if cmp(A[i], A[u]) > 0
                u = i
            end
            if cmp(A[i], A[v]) < 0
                v = i
            end
        end
        return u, v
    end
    function merge3(A, B, C, cmp)
        ia, ib, ic = length(A), length(B), length(C)
        fa, fb, fc, R::typeof(A) = ia > 0, ib > 0, ic > 0, []
        while fa || fb || fc
            if fa && (!fb || cmp(A[ia], B[ib]) < 0) && (!fc || cmp(A[ia], C[ic]) < 0)
                push!(R, A[ia]); ia -= 1; fa = ia > 0
            elseif fb && (!fc || cmp(B[ib], C[ic]) < 0)
                push!(R, B[ib]); ib -= 1; fb = ib > 0
            else
                push!(R, C[ic]); ic -= 1; fc = ic > 0
            end
        end
        return reverse(R)
    end
    function CH(Q)
        if length(Q) <= 3
            m = argmin(last.(Q))
            p0, Q = Q[m], vcat(Q[begin : m - 1], Q[m + 1 : end])
            if length(Q) == 2 && cross(Q[1] .- p0, Q[2] .- p0) < 0
                Q = reverse(Q)
            end
            return vcat(p0, Q)
        end
        m = length(Q) ÷ 2
        L, R = CH(Q[begin : m]), CH(Q[m + 1 : end])
        if L[1][2] > R[1][2]
            L, R = R, L
        end
        p0, L = L[1], L[2 : end]
        u, v = argMaxMin(R, (p, q) -> cross(p .- p0, q .- p0))
        slice(A, l, r) = l <= r ? A[l : r] : vcat(A[l : end], A[begin : r])
        R0, R1 = slice(R, u, v), slice(R, v, u)[end - 1 : -1 : begin + 1]
        return GrahamScan(vcat(p0, merge3(L, R0, R1, (p, q) -> cross(p .- p0, q .- p0))))
    end
    Q = sort(Q, by = first)
    return CH(Q)
end

实验测试

随机在 ${ (x, y) | 0 ≤ x, y < 100 \}$ 区域生成 3000 个点，测试各算法凸包结果和效率。代码：

using Random
function genData(n)
    Random.seed!(2024)
    return [q for q in zip(100 * rand(n), 100 * rand(n))]
end
Q = genData(3000)
@time P1 = NE(Q)
@time P2 = Jarvis(Q)
@time P3 = Graham(Q)
@time P4 = DC(Q)
@show length(P1)
@show length(P2), P2 == P1
@show length(P3), P3 == P1
@show length(P4), P4 == P1

结果：

  0.684359 seconds (585.47 k allocations: 40.423 MiB, 8.73% gc time, 53.56% compilation time)
  0.047316 seconds (27.23 k allocations: 2.719 MiB, 97.94% compilation time)
  0.175105 seconds (142.39 k allocations: 9.483 MiB, 3.61% gc time, 96.17% compilation time)
  0.412407 seconds (423.40 k allocations: 25.049 MiB, 98.20% compilation time)
length(P1) = 19
(length(P2), P2 == P1) = (19, true)
(length(P3), P3 == P1) = (19, true)
(length(P4), P4 == P1) = (19, true)