diff --git a/src/planar.jl b/src/planar.jl index 65c450db..1a6e96ec 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,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 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