diff --git a/src/PlanarConvexHulls.jl b/src/PlanarConvexHulls.jl index ffd73d3..2deef27 100644 --- a/src/PlanarConvexHulls.jl +++ b/src/PlanarConvexHulls.jl @@ -1,33 +1,3 @@ module PlanarConvexHulls -export - ConvexHull, - CW, - CCW, - vertices, - num_vertices, - area, - centroid, - is_ordered_and_convex, - closest_point, - jarvis_march!, - hrep, - hrep! - -using LinearAlgebra -using StaticArrays -using StaticArrays: arithmetic_closure - -include("order.jl") -include("util.jl") -include("exceptions.jl") -include("core_types.jl") -include("area.jl") -include("convexity_test.jl") -include("centroid.jl") -include("membership.jl") -include("closest_point.jl") -include("jarvis_march.jl") -include("hrep.jl") - end # module diff --git a/src/area.jl b/src/area.jl deleted file mode 100644 index 9a3e350..0000000 --- a/src/area.jl +++ /dev/null @@ -1,14 +0,0 @@ -function area(hull::ConvexHull) - # https://en.wikipedia.org/wiki/Shoelace_formula - T = eltype(hull) - vertices = hull.vertices - n = length(vertices) - n <= 2 && return zero(arithmetic_closure(T)) - @inbounds begin - ret = cross2(vertices[n], vertices[1]) - @simd for i in Base.OneTo(n - 1) - ret += cross2(vertices[i], vertices[i + 1]) - end - return abs(ret) / 2 - end -end diff --git a/src/centroid.jl b/src/centroid.jl deleted file mode 100644 index 5dca538..0000000 --- a/src/centroid.jl +++ /dev/null @@ -1,28 +0,0 @@ -function centroid(hull::ConvexHull) - # https://en.wikipedia.org/wiki/Centroid#Of_a_polygon - T = eltype(hull) - R = arithmetic_closure(T) - vertices = hull.vertices - n = length(vertices) - @inbounds begin - if n === 0 - error() - elseif n === 1 - return R.(vertices[1]) - elseif n === 2 - return (vertices[1] + vertices[2]) / 2 - else - c = cross2(vertices[n], vertices[1]) - centroid = (vertices[n] + vertices[1]) * c - double_area = c - @simd for i in Base.OneTo(n - 1) - c = cross2(vertices[i], vertices[i + 1]) - centroid += (vertices[i] + vertices[i + 1]) * c - double_area += c - end - double_area = abs(double_area) - centroid /= 3 * double_area - return centroid - end - end -end diff --git a/src/closest_point.jl b/src/closest_point.jl deleted file mode 100644 index 75e24ad..0000000 --- a/src/closest_point.jl +++ /dev/null @@ -1,37 +0,0 @@ -function closest_point(p::PointLike, hull::ConvexHull) - op = orientation_comparator(hull) - vertices = hull.vertices - n = length(vertices) - @inbounds begin - if n === 0 - throw(ArgumentError()) - elseif n === 1 - return vertices[1] - else - closest = p - i = 1 - v1 = vertices[n] - v2 = vertices[1] - while true - δ = v2 - v1 - p′ = p - v1 - strictly_outside_edge = op(cross2(p′, δ), 0) - if strictly_outside_edge - # Find closest point to current edge. - λ = clamp(p′ ⋅ δ / (δ ⋅ δ), false, true) - closest_to_edge = v1 + λ * δ - - # Accept if previous candidate was p itself or if it's closer. - if closest == p || normsquared(closest_to_edge - p) < normsquared(closest - p) - closest = closest_to_edge - end - end - i == n && break - v1 = vertices[i] - v2 = vertices[i + 1] - i += 1 - end - return closest - end - end -end diff --git a/src/convexity_test.jl b/src/convexity_test.jl deleted file mode 100644 index 6e2ea13..0000000 --- a/src/convexity_test.jl +++ /dev/null @@ -1,15 +0,0 @@ -function is_ordered_and_convex(vertices::AbstractVector{<:PointLike}, ::Type{O}) where {O<:VertexOrder} - op = orientation_comparator(O) - n = length(vertices) - n <= 2 && return true - @inbounds begin - δprev = vertices[n] - vertices[n - 1] - δnext = vertices[1] - vertices[n] - for i in Base.OneTo(n - 1) - op(cross2(δprev, δnext), 0) || return false - δprev = δnext - δnext = vertices[i + 1] - vertices[i] - end - return op(cross2(δprev, δnext), 0) - end -end diff --git a/src/core_types.jl b/src/core_types.jl deleted file mode 100644 index 862ae0d..0000000 --- a/src/core_types.jl +++ /dev/null @@ -1,27 +0,0 @@ -const PointLike{T} = StaticVector{2, T} - -struct ConvexHull{O<:VertexOrder, T, P<:PointLike{T}, V<:AbstractVector{P}} - vertices::V - - # TODO: use constructor to compute convex hull? - function ConvexHull{O}(vertices::V; check=true) where {O<:VertexOrder, T, P<:PointLike{T}, V<:AbstractVector{P}} - if check - is_ordered_and_convex(vertices, O) || throw(OrderedStronglyConvexError()) - end - new{O, T, P, V}(vertices) - end -end - -ConvexHull{O, T}() where {O<:VertexOrder, T} = ConvexHull{O}(SVector{2, T}[], check=false) - -Base.eltype(::Type{<:ConvexHull{<:Any, T}}) where {T} = T -vertex_order(::Type{<:ConvexHull{O}}) where {O} = O -vertex_order(hull::ConvexHull) = vertex_order(typeof(hull)) -edge_normal_sign_operator(hull::ConvexHull) = edge_normal_sign_operator(vertex_order(hull)) - -orientation_comparator(hull::ConvexHull) = orientation_comparator(vertex_order(hull)) -vertices(hull::ConvexHull) = hull.vertices -num_vertices(hull::ConvexHull) = length(vertices(hull)) -Base.isempty(hull::ConvexHull) = num_vertices(hull) > 0 -Base.empty!(hull::ConvexHull) = (empty!(hull.vertices); hull) -Base.sizehint!(hull::ConvexHull, n) = (sizehint!(hull.vertices, n); hull) diff --git a/src/exceptions.jl b/src/exceptions.jl deleted file mode 100644 index c204b24..0000000 --- a/src/exceptions.jl +++ /dev/null @@ -1,6 +0,0 @@ -struct OrderedStronglyConvexError <: Exception -end - -function Base.showerror(io::IO, e::OrderedStronglyConvexError) - print(io, "Points are not ordered or do not represent a strongly convex set") -end diff --git a/src/hrep.jl b/src/hrep.jl deleted file mode 100644 index ad5bfce..0000000 --- a/src/hrep.jl +++ /dev/null @@ -1,58 +0,0 @@ -""" -Return the equivalent halfspace representation of the convex hull, i.e. -matrix ``A`` and vector ``b`` such that the set of points inside the hull -is - -```math -\\left\\{ x \\mid A x \\le b \\right\\} -``` -""" -hrep(hull::ConvexHull) = _hrep(hull, Length(vertices(hull))) - -function _hrep(hull::ConvexHull, ::Length{N}) where N - T = eltype(hull) - R = arithmetic_closure(T) - vertices = hull.vertices - if N === StaticArrays.Dynamic() - n = length(vertices) - A = similar(vertices, R, n, 2) - b = similar(vertices, R, n) - hrep!(A, b, hull) - return A, b - else - Amut = similar(vertices, R, Size(N, 2)) - bmut = similar(vertices, R, Size(N)) - hrep!(Amut, bmut, hull) - A = convert(similar_type(vertices, R, Size(Amut)), Amut) - b = convert(similar_type(vertices, R, Size(bmut)), bmut) - return A, b - end -end - -@inline function hrep!(A::AbstractMatrix, b::AbstractVector, hull::ConvexHull) - signop = edge_normal_sign_operator(hull) - vertices = hull.vertices - n = length(vertices) - @boundscheck begin - size(A) == (n, 2) || throw(DimensionMismatch()) - length(b) == n || throw(DimensionMismatch()) - end - @inbounds @simd ivdep for i = Base.OneTo(n) - Ai, bi = hrep_kernel(vertices, signop, i) - A[i, 1] = Ai[1] - A[i, 2] = Ai[2] - b[i] = bi - end -end - -Base.@propagate_inbounds function hrep_kernel(vertices, signop, i) - n = length(vertices) - v1 = vertices[i] - v2 = vertices[ifelse(i == n, 1, i + 1)] - δ = v2 - v1 - δx, δy = unpack(δ) - outward_normal = signop(SVector(δy, -δx)) - Ai = transpose(outward_normal) - bi = Ai * v1 - Ai, bi -end diff --git a/src/jarvis_march.jl b/src/jarvis_march.jl deleted file mode 100644 index d1c2c74..0000000 --- a/src/jarvis_march.jl +++ /dev/null @@ -1,55 +0,0 @@ -function jarvis_march!(hull::ConvexHull, points::AbstractVector{<:PointLike}) - # Adapted from https://www.algorithm-archive.org/contents/jarvis_march/jarvis_march.html. - op = orientation_comparator(hull) - n = length(points) - vertices = hull.vertices - @inbounds begin - if n <= 2 - resize!(vertices, n) - vertices .= points - else - # Preallocate - resize!(vertices, n) - - # Find an initial hull vertex using lexicographic ordering. - start = last(points) - for i in Base.OneTo(n - 1) - p = points[i] - if Tuple(p) < Tuple(start) - start = p - end - end - - i = 1 - current = start - while true - # Add point - vertices[i] = current - - # Next point is the one with extremal internal angle. - next = last(points) - δnext = next - current - for i in Base.OneTo(n - 1) - p = points[i] - δ = p - current - c = cross2(δnext, δ) - - # Note the last clause here, which ensures strong convexity in the presence of - # collinear points by accepting `p` if it's farther away from `current` than - # `next`. - if next == current || op(0, c) || (c == 0 && δ ⋅ δ > δnext ⋅ δnext) - next = p - δnext = δ - end - end - current = next - current == first(vertices) && break - i += 1 - end - - # Shrink to computed number of vertices. - resize!(vertices, i) - end - end - return hull -end diff --git a/src/membership.jl b/src/membership.jl deleted file mode 100644 index 41e5f53..0000000 --- a/src/membership.jl +++ /dev/null @@ -1,23 +0,0 @@ -function Base.in(point::PointLike, hull::ConvexHull) - op = orientation_comparator(hull) - vertices = hull.vertices - n = length(vertices) - @inbounds begin - if n === 0 - return false - elseif n === 1 - return point == hull.vertices[1] - elseif n === 2 - p′ = point - vertices[1] - δ = vertices[2] - vertices[1] - cross2(p′, δ) == 0 && 0 <= p′ ⋅ δ <= δ ⋅ δ - else - δ = vertices[1] - vertices[n] - for i in Base.OneTo(n - 1) - op(cross2(point - vertices[i], δ), 0) && return false - δ = vertices[i + 1] - vertices[i] - end - return !op(cross2(point - vertices[n], δ), 0) - end - end -end diff --git a/src/order.jl b/src/order.jl deleted file mode 100644 index c00c9d4..0000000 --- a/src/order.jl +++ /dev/null @@ -1,17 +0,0 @@ -abstract type VertexOrder end -orientation_comparator(o::VertexOrder) = orientation_comparator(typeof(o)) -edge_normal_sign_operator(o::VertexOrder) = orientation_comparator(typeof(o)) - -""" -Counterclockwise order. -""" -struct CCW <: VertexOrder end -orientation_comparator(::Type{CCW}) = > -edge_normal_sign_operator(::Type{CCW}) = + - -""" -Clockwise order. -""" -struct CW <: VertexOrder end -orientation_comparator(::Type{CW}) = < -edge_normal_sign_operator(::Type{CW}) = - diff --git a/src/util.jl b/src/util.jl deleted file mode 100644 index 137d7b8..0000000 --- a/src/util.jl +++ /dev/null @@ -1,9 +0,0 @@ -unpack(v::StaticVector{2}) = @inbounds return v[1], v[2] - -@inline function cross2(v1::StaticVector{2}, v2::StaticVector{2}) - x1, y1 = unpack(v1) - x2, y2 = unpack(v2) - x1 * y2 - x2 * y1 -end - -normsquared(x) = sum(y -> y^2, x) diff --git a/test/runtests.jl b/test/runtests.jl index fdca216..547eb03 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,235 +1,94 @@ module PlanarConvexHullsTest -using PlanarConvexHulls using StaticArrays -using Test -using Random using LinearAlgebra -using Statistics - -const Point{T} = SVector{2, T} - -@testset "is_ordered_and_convex" begin - @test is_ordered_and_convex([Point(1, 2)], CCW) - @test is_ordered_and_convex([Point(1, 2), Point(3, 4)], CCW) - @test is_ordered_and_convex([Point(1, 2), Point(3, 4), Point(2, 4)], CCW) - @test !is_ordered_and_convex([Point(1, 2), Point(3, 4), Point(2, 3)], CCW) # on a line - - v1 = Point(0, 0) - v2 = Point(1, 0) - v3 = Point(1, 1) - v4 = Point(0, 1) - vertices = [v1, v2, v3, v4] - for i in 0 : 4 - shifted = circshift(vertices, i) - @test is_ordered_and_convex(shifted, CCW) - @test !is_ordered_and_convex(reverse(shifted), CCW) - @test is_ordered_and_convex(reverse(shifted), CW) - end - - for i in eachindex(vertices) - for j in eachindex(vertices) - if i != j - vertices′ = copy(vertices) - vertices′[i] = vertices[j] - vertices′[j] = vertices[i] - @test !is_ordered_and_convex(vertices′, CCW) - end - end - end -end -@testset "area" begin - @test area(ConvexHull{CCW}(SVector((SVector(1, 1),)))) == 0 +abstract type VertexOrder end +orientation_comparator(o::VertexOrder) = orientation_comparator(typeof(o)) - @test area(ConvexHull{CCW}(SVector(SVector(1, 1), SVector(2, 3)))) == 0 +struct CCW <: VertexOrder end +orientation_comparator(::Type{CCW}) = > - triangle = ConvexHull{CCW}(SVector(Point(1, 1), Point(2, 1), Point(3, 3))) - @test area(triangle) == 1.0 +unpack(v::StaticVector{2}) = @inbounds return v[1], v[2] - square = ConvexHull{CCW}(SVector(Point(1, 1), Point(4, 1), Point(4, 3), Point(1, 3))) - @test area(square) == 3 * 2 +@inline function cross2(v1::StaticVector{2}, v2::StaticVector{2}) + x1, y1 = v1[1], v1[2] + x2, y2 = v2[1], v2[2] + x1 * y2 - x2 * y1 end -@testset "in" begin - @testset "point" begin - rng = MersenneTwister(1) - p = SVector(1, 1) - C = ConvexHull{CCW}(SVector((p,))) - @test p ∈ C - @test Float64.(p) ∈ C - for i in 1 : 10 - @test p + SVector(randn(rng), randn(rng)) ∉ C - end - end - - @testset "line segment" begin - p1 = SVector(1, 1) - p2 = SVector(3, 5) - linesegment = ConvexHull{CCW}([p1, p2]) - @test p1 ∈ linesegment - @test p2 ∈ linesegment - @test div.(p1 + p2, 2) ∈ linesegment - @test p1 + 2 * (p2 - p1) ∉ linesegment - end - - @testset "triangle" begin - rng = MersenneTwister(1) - triangle = ConvexHull{CCW}(SVector(Point(1, 1), Point(2, 1), Point(3, 3))) - for p in vertices(triangle) - @test p ∈ triangle - end - for i = 1 : 100_000 - weights = normalize(rand(rng, SVector{3}), 1) - p = reduce(+, vertices(triangle) .* weights) - @test p ∈ triangle - end - end +struct ConvexHull{O<:VertexOrder, T, P<:StaticVector{2, T}, V<:AbstractVector{P}} + vertices::V - @testset "rectangle" begin - rng = MersenneTwister(1) - width = 4 - height = 3 - origin = Point(2, 4) - rectangle = ConvexHull{CCW}(map(x -> x + origin, SVector(Point(0, 0), Point(width, 0), SVector(width, height), SVector(0, height)))) - for p in vertices(rectangle) - @test p ∈ rectangle - end - for i = 1 : 100_000 - p = origin + SVector(width * rand(rng), height * rand(rng)) - @test p ∈ rectangle - end - for i = 1 : 10 - p = origin + SVector(width * rand(rng), height * rand(rng)) - @test setindex(p, origin[1] + width + rand(rng), 1) ∉ rectangle - @test setindex(p, origin[1] - rand(rng), 1) ∉ rectangle - @test setindex(p, origin[2] + height + rand(rng), 2) ∉ rectangle - @test setindex(p, origin[2] - rand(rng), 2) ∉ rectangle - end + function ConvexHull{O}(vertices::V; check=true) where {O<:VertexOrder, T, P<:StaticVector{2, T}, V<:AbstractVector{P}} + new{O, T, P, V}(vertices) end end -@testset "centroid" begin - @testset "point" begin - p = Point(1, 2) - @test centroid(ConvexHull{CCW}([p])) === Float64.(p) - end - - @testset "line segment" begin - p1 = SVector(1, 1) - p2 = SVector(3, 5) - linesegment = ConvexHull{CCW}([p1, p2]) - @test centroid(linesegment) == Point(2.0, 3.0) - end - - @testset "triangle" begin - triangle = ConvexHull{CCW}(SVector(Point(1, 1), Point(2, 1), Point(3, 3))) - @test centroid(triangle) ≈ mean(vertices(triangle)) atol=1e-15 - end -end - -function convex_hull_alg_test(hull_alg!) - @testset "random $order" for order in [CCW, CW] - hull = ConvexHull{order, Float64}() - rng = MersenneTwister(2) - for n = 1 : 10 - sizehint!(hull, n) # just to get code coverage - for _ = 1 : 10_000 - points = [rand(rng, Point{Float64}) for i = 1 : n] - hull_alg!(hull, points) - @test is_ordered_and_convex(vertices(hull), order) +ConvexHull{O, T}() where {O<:VertexOrder, T} = ConvexHull{O}(SVector{2, T}[], check=false) +vertex_order(::Type{<:ConvexHull{O}}) where {O} = O +vertex_order(hull::ConvexHull) = vertex_order(typeof(hull)) +orientation_comparator(hull::ConvexHull) = orientation_comparator(vertex_order(hull)) + +function jarvis_march!(hull::ConvexHull, points::AbstractVector{<:StaticVector{2}}) + # Adapted from https://www.algorithm-archive.org/contents/jarvis_march/jarvis_march.html. + op = orientation_comparator(hull) + # @show op # Uncommenting this makes tests pass with code coverage on! + n = length(points) + vertices = hull.vertices + if n <= 2 + error() + else + # Preallocate + resize!(vertices, n) + + # Find an initial hull vertex using lexicographic ordering. + start = last(points) + for i in Base.OneTo(n - 1) + p = points[i] + if Tuple(p) < Tuple(start) + start = p end end - end - - @testset "collinear input $order" for order in [CCW, CW] - hull = ConvexHull{order, Float64}() - points = [Point(0., 0.), Point(0., 1.), Point(0., 2.), Point(1., 0.), Point(1., 1.), Point(1., 2.)] - for i = 1 : 10 - shuffle!(points) - hull_alg!(hull, points) - @test is_ordered_and_convex(vertices(hull), order) - @test isempty(symdiff(vertices(hull), [Point(0., 0.), Point(1., 0.), Point(1., 2.), Point(0., 2.)])) - end - end -end - -@testset "jarvis_march!" begin - convex_hull_alg_test(jarvis_march!) -end - -@testset "closest_point $order" for order in [CCW, CW] - hull = ConvexHull{order, Float64}() - rng = MersenneTwister(3) - for n = 1 : 10 - for _ = 1 : 100 - points = [rand(rng, Point{Float64}) for i = 1 : n] - jarvis_march!(hull, points) - - for point in vertices(hull) - closest = closest_point(point, hull) - @test closest == point - end - for _ = 1 : 100 - point = rand(rng, Point{Float64}) - closest = closest_point(point, hull) - if point ∈ hull - @test closest == point - else - closest_to_closest = closest_point(closest, hull) - @test closest_to_closest ≈ closest atol=1e-14 + i = 1 + current = start + while true + # Add point + vertices[i] = current + + # Next point is the one with extremal internal angle. + next = last(points) + δnext = next - current + for i in Base.OneTo(n - 1) + p = points[i] + δ = p - current + c = cross2(δnext, δ) + + # Note the last clause here, which ensures strong convexity in the presence of + # collinear points by accepting `p` if it's farther away from `current` than + # `next`. + if next == current || op(0, c) || (c == 0 && δ ⋅ δ > δnext ⋅ δnext) + next = p + δnext = δ end end + current = next + current == first(vertices) && break + i += 1 + if i > n + error("Should never happen") + end end - end - - # visualize = true - # if visualize - # flush(stdout) - # hull = ConvexHull{CCW, Float64}() - # rng = MersenneTwister(3) - # n = 5 - # points = [rand(rng, Point{Float64}) for i = 1 : n] - # jarvis_march!(hull, points) - # projections = map(1 : 1_000_000) do _ - # closest_point(rand(rng, Point{Float64}), hull) - # end - # plt = scatterplot(getindex.(projections, 1), getindex.(projections, 2)) - # scatterplot!(plt, getindex.(vertices(hull), 1), getindex.(vertices(hull), 2), color = :red) - # display(plt) - # end -end - -function hreptest(hull::H, rng) where {H<:ConvexHull} - A, b = hrep(hull) - for _ = 1 : 100 - testpoint = rand(rng, Point{Float64}) - @test (testpoint ∈ hull) == all(A * testpoint .<= b) - end - if Length(vertices(hull)) !== Length{StaticArrays.Dynamic()}() - @test(@allocated(hrep(hull)) == 0) - end -end -@testset "hrep $order" for order in [CCW, CW] - dynamichull = ConvexHull{order, Float64}() - rng = MersenneTwister(4) - for n = 2 : 10 - for _ = 1 : 100 - points = [rand(rng, Point{Float64}) for i = 1 : n] - jarvis_march!(dynamichull, points) - hreptest(dynamichull, rng) - - h = num_vertices(dynamichull) - statichull = ConvexHull{order}(SVector{h}(vertices(dynamichull))) - hreptest(statichull, rng) - end + # Shrink to computed number of vertices. + resize!(vertices, i) end + return hull end -@testset "benchmarks" begin - include(joinpath("..", "perf", "runbenchmarks.jl")) -end +hull = ConvexHull{CCW, Float64}() +points = StaticArrays.SArray{Tuple{2},Float64,1,2}[[0.0271377, 0.785383], [0.689278, 0.242258], [0.711323, 0.403459]] +jarvis_march!(hull, points) end # module