#1384 - Add four point case for convex hull

docs/src/lib/utils.md

@@ -21,6 +21,7 @@ get_radius!
 inner
 isinvertible
 ispermutation
+iscounterclockwise
 issquare
 is_right_turn
 is_tighter_same_dir_2D

src/concrete_convex_hull.jl

@@ -28,8 +28,9 @@ 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`, `2` and `3` points in any dimension.
-For `4` or more points, the algorithm used depends on the dimension.
+A pre-processing step treats the cases with up to two points for one dimension
+and up to four points for two dimensions.
+For more points in one resp. two dimensions, we use more general algorithms.
 
 For the one-dimensional case we return the minimum and maximum points, in that order.
 
@@ -114,6 +115,9 @@ function convex_hull!(points::Vector{VN};
         elseif m == 3
             # three points case in 2d
             return _three_points_2d!(points)
+        elseif m == 4
+            # four points case in 2d
+            return _four_points_2d!(points)
         else
             # general case in 2d
             return _convex_hull_2d!(points, algorithm=algorithm)
@@ -186,6 +190,124 @@ function _three_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N
     return points
 end
 
+function _collinear_case!(points::Vector{VN}, A::VN, B::VN, C::VN, D::VN) where {N<:Real, VN<:AbstractVector{N}}
+    # A, B and C collinear, D is the extra point
+    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])
+            # the three points are approximately equal
+            points[1], points[2] = A, D
+            pop!(points)
+            pop!(points)
+            return _two_points_2d!(points)
+        else
+            # assign the points with max and min value in their second component to the
+            # firsts points and the extra point to the third place, then pop the point that was in the middle
+            points[1], points[2], points[3] = points[argmin([A[2], B[2], C[2]])], points[argmax([A[2], B[2], C[2]])], D
+            pop!(points)
+        end
+    else
+        # assign the points with max and min value in their first component to the
+        # firsts points and the extra point to the third place, then pop the point that was in the middle
+        points[1], points[2], points[3] = points[argmin([A[1], B[1], C[1]])], points[argmax([A[1], B[1], C[1]])], D
+        pop!(points)
+    end
+    return _three_points_2d!(points)
+end
+
+function _four_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N<:Real}
+    A, B, C, D = points[1], points[2], points[3], points[4]
+    tri_ABC = right_turn(A, B, C)
+    tri_ABD = right_turn(A, B, D)
+    tri_BCD = right_turn(B, C, D)
+    tri_CAD = right_turn(C, A, D)
+    key = 0
+    if tri_ABC > zero(N)
+        key = key + 1000
+    end
+    if tri_ABD > zero(N)
+        key = key + 100
+    end
+    if tri_BCD > zero(N)
+        key = key + 10
+    end
+    if tri_CAD > zero(N)
+        key = key + 1
+    end
+
+    if isapproxzero(tri_ABC)
+        return _collinear_case!(points, A, B, C, D)
+    end
+    if isapproxzero(tri_ABD)
+        return _collinear_case!(points, A, B, D, C)
+    end
+    if isapproxzero(tri_BCD)
+        return _collinear_case!(points, B, C, D, A)
+    end
+    if isapproxzero(tri_CAD)
+        return _collinear_case!(points, C, A, D, B)
+    end
+
+    # ABC  ABD  BCD  CAD  hull
+    # ------------------------
+    #  +    +    +    +   ABC
+    #  +    +    +    -   ABCD
+    #  +    +    -    +   ABDC
+    #  +    +    -    -   ABD
+    #  +    -    +    +   ADBC
+    #  +    -    +    -   BCD
+    #  +    -    -    +   CAD
+    #  +    -    -    -   [should not happen]
+    #  -    +    +    +   [should not happen]
+    #  -    +    +    -   ACD
+    #  -    +    -    +   DCB
+    #  -    +    -    -   DACB
+    #  -    -    +    +   ADB
+    #  -    -    +    -   ACDB
+    #  -    -    -    +   ADCB
+    #  -    -    -    -   ACB
+    if key == 1111
+        points[1], points[2], points[3] = A, B, C # +    +    +    +   ABC
+        pop!(points)
+    elseif key == 1110
+        points[1], points[2], points[3], points[4] = A, B, C, D # +    +    +    -   ABCD
+    elseif key == 1101
+        points[1], points[2], points[3], points[4] = A, B, D, C #  +    +    -    +   ABDC
+    elseif key == 1100
+        points[1], points[2], points[3] = A, B, D #  +    +    -    -   ABD
+        pop!(points)
+    elseif key == 1011
+        points[1], points[2], points[3], points[4] = A, D, B, C #  +    -    +    +   ADBC
+    elseif key == 1010
+        points[1], points[2], points[3] = B, C, D #  +    -    +    -   BCD
+        pop!(points)
+    elseif key == 1001
+        points[1], points[2], points[3] = C, A, D #  +    -    -    +    CAD
+        pop!(points)
+    elseif key == 0110
+        points[1], points[2], points[3] = A, C, D #  -    +    +    -   ACD
+        pop!(points)
+    elseif key == 0101
+        points[1], points[2], points[3] = D, C, B #  -    +    -    +   DCB
+        pop!(points)
+    elseif key == 0100
+        points[1], points[2], points[3], points[4] = D, A, C, B #  -    +    -    -   DACB
+    elseif key == 0011
+        points[1], points[2], points[3] = A, D, B #  -    -    +    +   ADB
+        pop!(points)
+    elseif key == 0010
+        points[1], points[2], points[3], points[4] = A, C, D, B #  -    -    +    -   ACDB
+    elseif key == 0001
+        points[1], points[2], points[3], points[4] = A, D, C, B #  -    -    -    +   ADCB
+    elseif key == 0000
+        points[1], points[2], points[3] = A, C, B #  -    -    -    -   ACB
+        pop!(points)
+    else
+        @assert false "unexpected case in convex_hull"
+    end
+    return points
+end
+
 function _convex_hull_1d!(points::Vector{VN})::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
     points[1:2] = [minimum(points), maximum(points)]
     return resize!(points, 2)

src/helper_functions.jl

@@ -645,3 +645,23 @@ A copy of the vector where each negative entry is replaced by zero.
 function rectify(x::AbstractVector{N}) where {N<:Real}
     return map(xi -> max(xi, zero(N)), x)
 end
+
+"""
+    iscounterclockwise(result, correct_expr) where {N<:Real}
+
+Returns a boolean, true if the elements in `result` are ordered in
+a couter-clockwise fashion and in the same order as `correct_expr`.
+
+### Input
+
+- `result` -- vector
+- `correct_expr` -- paragon vector with elements in ccw order
+
+### Output
+
+A boolean indicating if the elements of `result` are in the same order as
+`correct_expr` or any of its cyclic permutations.
+"""
+function iscounterclockwise(result, correct_expr)
+    return any([result == circshift(correct_expr, i) for i in 0:length(result)-1])
+end

test/runtests.jl

@@ -11,7 +11,7 @@ using TaylorModels: set_variables, TaylorModelN
 include("to_N.jl")
 
 # non-exported helper functions
-using LazySets: ispermutation, isinvertible, inner
+using LazySets: ispermutation, isinvertible, inner, iscounterclockwise
 using LazySets.Approximations: UnitVector
 
 global test_suite_basic = true

test/unit_convex_hull.jl

@@ -39,10 +39,6 @@ for N in [Float64, Rational{Int}]
     @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]]
@@ -55,6 +51,47 @@ for N in [Float64, Rational{Int}]
     @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
 
+    # four vertex case in 2 dimentions
+    A = N[1, 0]
+    B = N[1, 1]
+    C = N[-1, 1]
+    D = N[-1, 0]
+    expr = [A, B, C, D]
+    points = [A, B, C, D]
+    @test  iscounterclockwise(convex_hull!(points), expr) # ABCD
+    points = [A, D, C, B]
+    @test  iscounterclockwise(convex_hull!(points), expr) # ADCB
+    points = [A, B, D, C]
+    @test  iscounterclockwise(convex_hull!(points), expr) # ABDC
+    points = [A, D, B, C]
+    @test  iscounterclockwise(convex_hull!(points), expr) # ADBC
+    points = [D, A, C, B]
+    @test  iscounterclockwise(convex_hull!(points), expr) # DACB
+    points = [A, C, D, B]
+    @test  iscounterclockwise(convex_hull!(points), expr) # ACDB
+    A = N[0, 1]
+    B = N[-1, -1]
+    C = N[1, -1]
+    D = N[0, 0]
+    expr = [A, B, C]
+    points = [A, B, C, D]
+    @test iscounterclockwise(convex_hull!(points), expr) # ABC
+    points = [A, C, B, D]
+    @test iscounterclockwise(convex_hull!(points), expr) # CBA
+    points = [D, B, C, A]
+    @test iscounterclockwise(convex_hull!(points), expr) # BCD
+    points = [D, C, B, A]
+    @test iscounterclockwise(convex_hull!(points), expr) # DCB
+    points = [B, D, C, A]
+    @test iscounterclockwise(convex_hull!(points), expr) # ACD
+    points = [C, D, B, A]
+    @test iscounterclockwise(convex_hull!(points), expr) # CAD
+    points = [B, C, D, A]
+    @test iscounterclockwise(convex_hull!(points), expr) # ABD
+    points = [C, B, D, A]
+    @test iscounterclockwise(convex_hull!(points), expr) # ADB
+    @test ispermutation(convex_hull!([N[1, 1], N[2, 2], N[3, 3], N[4, 4]]), [N[1, 1], N[4, 4]]) # points aligned
+
     # five-vertices case in 2D
     points = to_N(N, [[0.9, 0.2], [0.4, 0.6], [0.2, 0.1], [0.1, 0.3], [0.3, 0.28]])
     points_copy = copy(points)