二维凸包算法 Julia实现

问题描述：给定平面上 $n$ 个点的集合 $Q$，求其子集 $P$ 构成 $Q$ 的凸包，即 $\forall p\in Q,\exists 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 .- p0c0, c1, c2 = cross(v1, v2), cross(v0, v1), cross(v0, v2)return c0 * c1 < 0 && c0 * c2 > 0 && abs(c2 - c1) < abs(c0)endm = argmin(last.(Q))p0, Q = Q[m], vcat(Q[begin : m - 1], Q[m + 1 : end])for i = length(Q) : -1 : 1f = falsefor j = 2 : length(Q)if i == j continue endfor k = 1 : j - 1if i == k continue endif inΔ(p0, Q[j], Q[k], Q[i])f = true; breakendendif f break endendif fdeleteat!(Q, i)endendsort!(Q, lt = (p, q) -> cross(p .- p0, q .- p0) > 0)return vcat(p0, Q)
end

Jarvis March 算法

按逆时针对 $P$ 排序，将问题转化成已知 $p_0,p_1,...,p_k$ 如何求 $p_{k+1}$。

$p_0$ 固定为纵坐标最小的点。此时 Q − { p k } Q - \{p_k\} Q−{pk​} 中所有点都在向量 $p_k-p_{k-1}$ 的左侧，即向量集 { p − p k ∣ p ∈ Q − { p k } } \{ p - p_k | p \in Q - \{ p_k \} \} {p−pk​∣p∈Q−{pk​}} 可在叉积运算的正负性上构成偏序集，最右的向量即为 $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 = iendendreturn rendP = [march(Q, (p, q) -> last(p) < last(q))]while trueq = march(Q, (p, q) -> cross(p .- P[end], q .- P[end]) > 0)if P[begin] == qbreakendpush!(P, q)endreturn 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 = \{ p_0, p_1, .., p_k \} Pk​={p0​,p1​,..,pk​}，注意到 $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} \} Pk+1​=Pk​+{pk+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}$。

用栈存储中间过程的凸包计算，$Q$ 中每个点仅会进出栈一次，故扫描为线性时间，主要耗时在排序。时间复杂度 $O(n\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]) > 0pop!(P)endpush!(P, q)endreturn P
endfunction 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(n)=T(\frac{n}{2})+O(n)$，得时间复杂度 $O(n\log n)$，但常数很大。代码：

function DC(Q)function argMaxMin(A, cmp)u, v = 1, 1for i = 2 : length(A)if cmp(A[i], A[u]) > 0u = iendif cmp(A[i], A[v]) < 0v = iendendreturn u, vendfunction 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 || fcif fa && (!fb || cmp(A[ia], B[ib]) < 0) && (!fc || cmp(A[ia], C[ic]) < 0)push!(R, A[ia]); ia -= 1; fa = ia > 0elseif fb && (!fc || cmp(B[ib], C[ic]) < 0)push!(R, B[ib]); ib -= 1; fb = ib > 0elsepush!(R, C[ic]); ic -= 1; fc = ic > 0endendreturn reverse(R)endfunction CH(Q)if length(Q) <= 3m = 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) < 0Q = reverse(Q)endreturn vcat(p0, Q)endm = length(Q) ÷ 2L, R = CH(Q[begin : m]), CH(Q[m + 1 : end])if L[1][2] > R[1][2]L, R = R, Lendp0, 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))))endQ = sort(Q, by = first)return CH(Q)
end

实验测试

随机在 $\{(x,y)|0\le 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)