Fix edge case for planar hull (#334)

* Fix edge case for planar hull

* Fix

Project.toml

@@ -1,7 +1,7 @@
 name = "Polyhedra"
 uuid = "67491407-f73d-577b-9b50-8179a7c68029"
 repo = "https://github.com/JuliaPolyhedra/Polyhedra.jl.git"
-version = "0.7.7"
+version = "0.7.8"
 
 [deps]
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"

src/planar.jl

@@ -5,13 +5,15 @@ function _semi_hull(ps::Vector{PT}, sign_sense, counterclockwise, sweep_norm, yr
     end
     prev = sign_sense == 1 ? first(ps) : last(ps)
     cur = prev
+    flat_start = true
     # Invariant:
     # We either have:
     # * `hull` is empty and `cur == prev` or
     # * `hull` is nonempty and `prev = last(hull)`.
     # In any case, the semihull of `ps[1:(j-1)]` is given by `[hull; cur]`.
     for j in (sign_sense == 1 ? (2:length(ps)) : ((length(ps)-1):-1:1))
         skip = false
+        flat = false
         while prev != cur
             cur_vec = cur - prev
             psj_vec = ps[j] - prev
@@ -26,9 +28,30 @@ function _semi_hull(ps::Vector{PT}, sign_sense, counterclockwise, sweep_norm, yr
                 #    The one that is closer to `prev` is redundant.
                 dcur = dot(cur_vec, sweep_norm)
                 dpsj = dot(psj_vec, sweep_norm)
-                if isapproxzero(dcur) && isapproxzero(dpsj)
+                flat = isapproxzero(dcur) && isapproxzero(dpsj)
+                if flat
                     # Case 1
-                    if sign_sense * counterclockwise(cur_vec, sweep_norm) < sign_sense * counterclockwise(psj_vec, sweep_norm)
+                    # There can be two flat plateaus. The first one at the start and then one at the end
+                    # `flat_start` indicates in which case we are. We need to a different direction in both cases
+                    sense = flat_start ? sign_sense : -sign_sense
+                    ccjsweep = sense * counterclockwise(psj_vec, sweep_norm)
+                    cccursweep = sense * counterclockwise(cur_vec, sweep_norm)
+                    if cccursweep < 0 < ccjsweep
+                        # In this case, `prev` -> `cur` was in the right direction but now that
+                        # we discover `ps[j]`, we realize that `prev` is in the interval
+                        # `[ps[j], cur]`. Even if `ps[j]` -> `cur` seems to be in the right direction, we cannot
+                        # keep it because it is in the wrong order in `ps`, it will be picked up by the other
+                        # call to `_semi_hull` in which case `ps[j]` should be skipped (this is typically the case for the first flat plateau)
+                        # For the second (and last flat plateau), it will rather be the case that `cur` should be removed.
+                        # The only thing that is sure here is that `prev` should be removed so we remove `prev` and we `continue`
+                        # as the code should then automatically decide which of `ps[j]` or `cur` should be removed
+                        pop!(hull)
+                        prev = cur
+                        if !isempty(hull)
+                            prev = last(hull)
+                        end
+                        continue
+                    elseif cccursweep < ccjsweep
                         skip = true
                         break
                     end
@@ -46,6 +69,9 @@ function _semi_hull(ps::Vector{PT}, sign_sense, counterclockwise, sweep_norm, yr
                 prev = last(hull)
             end
         end
+        if prev != cur && !flat
+            flat_start = false
+        end
         if !skip
             push!(hull, cur)
             prev = cur

test/Project.toml

@@ -6,6 +6,7 @@ HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

test/planar.jl

@@ -108,6 +108,38 @@ function test_issue_271()
     @test collect(points(p)) == expected
 end
 
+function test_issue_333()
+    vs = [
+        0.5 -0.4
+        -0.5  1
+        0    1
+        0.8  1
+    ]
+
+    h = Polyhedra.planar_hull(vrep(vs))
+    @test all(points(h) .== [vs[i, :] for i in [1, 2, 4]])
+
+    vs = [-0.96     -0.4;
+        -0.5      -0.4;
+        0.0      -0.4;
+        0.5      -0.4;
+        1.0      -0.36065655;
+        -0.995    -0.05;
+        1.0      -0.05;
+        -1.0       0.3;
+        0.97      0.3;
+        -1.0       0.65;
+        0.930928  0.65;
+        -1.0       1.0;
+        -0.5       1.0;
+        0.0       1.0;
+        0.5       1.0;
+        0.865979 1.0;]
+
+    h = Polyhedra.planar_hull(vrep(vs))
+    @test all(points(h) .== [vs[i, :] for i in [1, 6, 8, 12, 16, 11, 9, 7, 5, 4]])
+end
+
 @testset "Planar" begin
     for T in [Float64, Int]
         for (VT, d) in [(Vector{T}, 2), (SVector{2, T}, StaticArrays.Size(2))]
@@ -116,4 +148,5 @@ end
     end
     test_planar_square_with_more()
     test_issue_271()
+    test_issue_333()
 end