From 3b2b4621abaa73f100920d6235af53cf591f1572 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Wed, 14 Feb 2024 16:00:59 +0100
Subject: [PATCH] Fix edge case for planar hull

---
 src/planar.jl     | 18 ++++++++++++++++--
 test/Project.toml |  1 +
 test/planar.jl    | 33 +++++++++++++++++++++++++++++++++
 3 files changed, 50 insertions(+), 2 deletions(-)

diff --git a/src/planar.jl b/src/planar.jl
index 65c450db..0acda15c 100644
--- a/src/planar.jl
+++ b/src/planar.jl
@@ -5,6 +5,7 @@ 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
@@ -12,6 +13,7 @@ function _semi_hull(ps::Vector{PT}, sign_sense, counterclockwise, sweep_norm, yr
     # 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,18 @@ 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)
+                    sense = flat_start ? sign_sense : -sign_sense
+                    sense = sign_sense
+                    ccjsweep = sense * counterclockwise(psj_vec, sweep_norm)
+#                    if ccjsweep > 0
+#                        hull[end] = ps[j]
+#                        prev = last(hull)
+#                        break
+#                    end
+                    if sense * counterclockwise(cur_vec, sweep_norm) < ccjsweep
                         skip = true
                         break
                     end
@@ -46,6 +57,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
diff --git a/test/Project.toml b/test/Project.toml
index 3b9d50da..70379bd8 100644
--- a/test/Project.toml
+++ b/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"
diff --git a/test/planar.jl b/test/planar.jl
index 4c7c4cf3..172fec44 100644
--- a/test/planar.jl
+++ b/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