From de45252a1ff8b5ca9d413d8705b11ba21d842953 Mon Sep 17 00:00:00 2001 From: Sebastian Guadalupe Date: Wed, 29 May 2019 17:09:27 -0300 Subject: [PATCH] Add three point case for concrete convex hull --- src/concrete_convex_hull.jl | 112 ++++++++++++++++++++++++++---------- test/runtests.jl | 1 + test/unit_convex_hull.jl | 25 ++++++-- 3 files changed, 103 insertions(+), 35 deletions(-) diff --git a/src/concrete_convex_hull.jl b/src/concrete_convex_hull.jl index f022fe08f5..3ac54df7d0 100644 --- a/src/concrete_convex_hull.jl +++ b/src/concrete_convex_hull.jl @@ -28,8 +28,8 @@ The convex hull as a list of vectors with the coordinates of the points. ### Algorithm -A pre-processing step treats the cases with `0`, `1` and `2` points in any dimension. -For more than `3` points, the algorithm used depends on the dimension. +A pre-processing step treats the cases with `0`, `1`, `2` and `3` points in any dimension. +For `4` or more points, the algorithm used depends on the dimension. For the one-dimensional case we return the minimum and maximum points, in that order. @@ -99,41 +99,91 @@ function convex_hull!(points::Vector{VN}; n = length(first(points)) - # two points - if m == 2 - if n == 1 - p1, p2 = points[1], points[2] - if p1 == p2 # check for redundancy - pop!(points) - elseif p1[1] > p2[1] - points[1], points[2] = p2, p1 - end - elseif n == 2 - # special case, see #876 - p1, p2 = points[1], points[2] - if p1 == p2 # check for redundancy - pop!(points) - elseif p1 <= p2 - nothing - else - points[1], points[2] = p2, p1 - end + if n == 1 # dimensional check + if m == 2 + # two points case in 1d + return _two_points_1d!(points) + else + # general case in 1d + return _convex_hull_1d!(points) end - return points + elseif n == 2 + if m == 2 + # two points case in 2d + return _two_points_2d!(points) + elseif m == 3 + # three points case in 2d + return _three_points_2d!(points) + else + # general case in 2d + return _convex_hull_2d!(points, algorithm=algorithm) + end + else + # general case in nd + return _convex_hull_nd!(points, backend=backend) + end +end + +function _two_points_1d!(points) + p1, p2 = points[1], points[2] + if p1 == p2 + # check for redundancy + pop!(points) + elseif p1[1] > p2[1] + points[1], points[2] = p2, p1 end + return points +end - # =========================== - # Dispatch for general cases - # =========================== +function _two_points_2d!(points) - # _convex_hull_x can assume that there are at least three points - if n == 1 - return _convex_hull_1d!(points) - elseif n == 2 - return _convex_hull_2d!(points, algorithm=algorithm) + # special case, see #876 + p1, p2 = points[1], points[2] + if p1 == p2 + # check for redundancy + pop!(points) + elseif p1 <= p2 + nothing else - return _convex_hull_nd!(points, backend=backend) + points[1], points[2] = p2, p1 + end + return points +end + +function _three_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N<:Real} + # Algorithm: the function takes three points and uses the formula + # from here: https://stackoverflow.com/questions/2122305/convex-hull-of-4-points/2122620#2122620 + # to decide if the points are ordered in a counter-clockwise fashion or not, the result is saved + # in the 'turn' boolean, then returns the points in ordered ccw acting according to 'turn'. For the + # cases where the points are collinear we pass the points with the minimum and maximum first + # component to the function for two points in 2d(_two_points_2d), if those are equal, we do the same + # but with the second component. + A, B, C = points[1], points[2], points[3] + turn = (A[2] - B[2]) * C[1] + (B[1] - A[1]) * C[2] + (A[1] * B[2] - B[1] * A[2]) + + if isapproxzero(turn) + # ABC are collinear + if isapprox(A[1], B[1]) && isapprox(B[1], C[1]) && isapprox(C[1], A[1]) + # points are approximately equal in their first component + if isapprox(A[2], B[2]) && isapprox(B[2], C[2]) && isapprox(C[2], A[2]) + # all points are approximately equal + pop!(points) + pop!(points) + return points + else + points[1], points[2] = points[argmin([A[2], B[2], C[2]])], points[argmax([A[2], B[2], C[2]])] + end + else + points[1], points[2] = points[argmin([A[1], B[1], C[1]])], points[argmax([A[1], B[1], C[1]])] + end + pop!(points) + _two_points_2d!(points) + elseif turn < zero(N) + # ABC is CW + points[1], points[2], points[3] = C, B, A end + # else ABC is CCW => nothing to do + return points end function _convex_hull_1d!(points::Vector{VN})::Vector{VN} where {N<:Real, VN<:AbstractVector{N}} diff --git a/test/runtests.jl b/test/runtests.jl index 5ee20e222f..9324b4b993 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -104,6 +104,7 @@ if test_suite_basic @time @testset "LazySets.CartesianProduct" begin include("unit_CartesianProduct.jl") end @time @testset "LazySets.ResetMap" begin include("unit_ResetMap.jl") end @time @testset "LazySets.SymmetricIntervalHull" begin include("unit_SymmetricIntervalHull.jl") end + @time @testset "LazySets.concrete_convex_hull" begin include("unit_convex_hull.jl") end # ====================== # Testing set interfaces diff --git a/test/unit_convex_hull.jl b/test/unit_convex_hull.jl index e9be3b3530..3916e30d62 100644 --- a/test/unit_convex_hull.jl +++ b/test/unit_convex_hull.jl @@ -33,16 +33,33 @@ for N in [Float64, Rational{Int}] convex_hull!(points_2D) # check in-place version @test points_2D == [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]] + # three vertex case in 2 dimensions + function iscounterclockwise(result, correct_expr) + # this function checks if the result equals any of the correct answer cyclic permutation + result == correct_expr || result == circshift(correct_expr, 1) || result == circshift(correct_expr, 2) + end + ccw_points = [N[1, 1], N[-1, 1], N[-1, 0]] + ccw_p = convex_hull!(ccw_points) + ccw_expr = [N[1, 1], N[-1, 1], N[-1, 0]] + @test iscounterclockwise(ccw_p, ccw_expr) # points sorted in a counter-clockwise fashion + cw_points = [N[-1, 1], N[1, 1], N[-1, 0]] + cw_p = convex_hull!(cw_points) + cw_expr = [N[1, 1], N[-1, 1], N[-1, 0]] + @test iscounterclockwise(cw_p, cw_expr) # points sorted in clockwise fashion + @test ispermutation(convex_hull!([N[1, 1], N[2, 2], N[3, 3]]), [N[1, 1], N[3, 3]]) # points aligned + @test ispermutation(convex_hull!([N[0, 1], N[0, 2], N[0, 3]]), [N[0, 1], N[0, 3]]) # three points on a vertical line + @test convex_hull!([N[0, 1], N[0, 1], N[0, 1]]) == [N[0, 1]] # three equal points + # higher dimension if test_suite_polyhedra && N != Float32 # no backend supporting Float32 points_3D = [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(0), N(0), N(0)], [N(1), N(2), N(3)]] - @test convex_hull(points_3D) == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)], + @test ispermutation(convex_hull(points_3D), [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]) convex_hull!(points_3D) # check in-place version - @test points_3D == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)], + @test ispermutation(points_3D, [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]) end @@ -53,7 +70,7 @@ for N in [Float64, Rational{Int}] V1 = VPolygon([[N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]) V2 = VPolygon([[N(-1), N(-1)], [N(1/2), N(1/2)]]) ch = convex_hull(V1, V2) - @test vertices_list(ch) == ispermutation([[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]) + @test ispermutation(vertices_list(ch), [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]) if test_suite_polyhedra && N != Float32 # no backend supporting Float32 V1 = VPolytope([[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)]]) @@ -61,7 +78,7 @@ for N in [Float64, Rational{Int}] [N(1), N(2), N(3)]]) ch = convex_hull(V1, V2) @test ch isa VPolytope - @test vertices_list(ch) == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)], + @test ispermutation(vertices_list(ch), [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]) end