Fix planar hull on corner case (#320)

* Fix planar hull on corner case

* Fix first and last of semi hulls

Author: blegat

src/planar.jl

@@ -1,6 +1,6 @@
-function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray::Nothing = nothing) where PT
+function _semi_hull(ps::Vector{PT}, sign_sense, counterclockwise, sweep_norm, yray::Nothing = nothing) where PT
     hull = PT[]
-    if length(ps) == 0
+    if isempty(ps)
         return hull
     end
     prev = sign_sense == 1 ? first(ps) : last(ps)
@@ -17,12 +17,22 @@ function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray::Nothing
             psj_vec = ps[j] - prev
             cc = counterclockwise(cur_vec, psj_vec)
             if isapproxzero(cc)
-                # `cur` and `ps[j]` are on the same ray starting from `prev`.
-                # The one that is closer to `prev` is redundant.
-                # If `j` is the last index and redundant (it may happen if this
-                # ray is perpendicular to the direction of sorting) then we should
-                # also avoid adding `ps[j]` to `hull` so we set `skip` to `true`.
-                if norm(cur_vec, 1) > norm(psj_vec, 1)
+                # `prev`, `cur` and `ps[j]` are on the same line
+                # Two cases here:
+                # 1) `prev`, `cur`, `ps[j]` are on a line perpendicular to `sweep_norm`
+                #    The one that is not clockwise is redundant.
+                # 2) `cur`, `ps[j]` and on the same ray starting from `prev` of direction
+                #    `sweep_norm * sign_sense`
+                #    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)
+                    # Case 1
+                    if sign_sense * counterclockwise(cur_vec, sweep_norm) < sign_sense * counterclockwise(psj_vec, sweep_norm)
+                        skip = true
+                        break
+                    end
+                elseif sign_sense * dcur > sign_sense * dpsj
                     skip = true
                     break
                 end
@@ -46,14 +56,16 @@ function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray::Nothing
     return hull
 end
 
-function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray) where PT
-    hull = getsemihull(ps, sign_sense, counterclockwise)
+function _semi_hull(ps::Vector{PT}, sign_sense, counterclockwise, sweep_norm, yray) where PT
+    hull = _semi_hull(ps, sign_sense, counterclockwise, sweep_norm)
     while length(hull) >= 2 && counterclockwise(hull[end] - hull[end - 1], yray) >= 0
         pop!(hull)
     end
     return hull
 end
 
+_colinear(counterclockwise, a, b, c) = isapproxzero(counterclockwise(b - a, c - a))
+
 function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate)
     line = nothing
     lineleft = false
@@ -124,32 +136,31 @@ function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate)
     else
         sweep_norm = rotate(coord(line))
     end
-    sort!(points, by = x -> dot(x, sweep_norm))
-
-    if !isempty(points)
-        # `getsemihull` fails if `points` starts or end with several points on the same sweep line that are not ordered clockwise
-        start_line = dot(points[1], sweep_norm)
-        # `isapprox` won't work well with `atol=0` (which is the default) if `start_line` is zero, so we set a nonzero `atol`.
-        # TODO We should also multiply it with a scaling.
-        end_start = something(findfirst(x -> !isapprox(dot(x, sweep_norm), start_line, atol=Base.rtoldefault(typeof(start_line))), points), length(points) + 1) - 1
-        if end_start > 1
-            sort!(view(points, 1:end_start), by = x -> counterclockwise(x, sweep_norm), rev=true)
-        end
-        end_line = dot(points[end], sweep_norm)
-        start_end = something(findlast(x -> !isapprox(dot(x, sweep_norm), end_line, atol=Base.rtoldefault(typeof(start_line))), points), 0) + 1
-        if start_end < length(points)
-            sort!(view(points, start_end:length(points)), by = x -> counterclockwise(x, sweep_norm))
-        end
-    end
+    sort!(points, by = Base.Fix2(dot, sweep_norm))
 
     _points = eltype(points)[]
     _lines = eltype(lines)[]
     _rays = eltype(rays)[]
     if line === nothing
-        append!(_points, getsemihull(points, 1, counterclockwise, yray))
+        half_points = _semi_hull(points, 1, counterclockwise, sweep_norm, yray)
         if yray === nothing
-            append!(_points, getsemihull(points, -1, counterclockwise, yray)[2:end-1])
+            other_half = _semi_hull(points, -1, counterclockwise, sweep_norm, yray)[2:end-1]
+            if !isempty(other_half) && length(half_points) > 1 &&
+                _colinear(counterclockwise, half_points[1], half_points[2], other_half[end])
+                start = 2
+            else
+                start = 1
+            end
+            if !isempty(other_half) && length(half_points) > 1 &&
+                _colinear(counterclockwise, half_points[end], half_points[end - 1], other_half[1])
+                stop = length(half_points) - 1
+            else
+                stop = length(half_points)
+            end
+            append!(_points, half_points[start:stop])
+            append!(_points, other_half)
         else
+            append!(_points, half_points)
             push!(_rays, Ray(xray))
             if !(Ray(xray) ≈ Ray(yray))
                 push!(_rays, Ray(yray))

src/recipe.jl

@@ -18,8 +18,9 @@ function planar_contour(p::Polyhedron)
     end
     sort!(ps, by = x -> x[1])
     counterclockwise(p1, p2) = dot(cross([p1; 0], [p2; 0]), [0, 0, 1])
-    top = getsemihull(ps,  1, counterclockwise)
-    bot = getsemihull(ps, -1, counterclockwise)
+    sweep_norm = basis(eltype(ps), fulldim(p), 1)
+    top = _semi_hull(ps,  1, counterclockwise, sweep_norm)
+    bot = _semi_hull(ps, -1, counterclockwise, sweep_norm)
     if !isempty(top) && !isempty(bot)
         @assert top[end] == bot[1]
         pop!(top)

test/issue311.jl

@@ -0,0 +1,26 @@
+function issue311test(lib::Polyhedra.Library)
+    W = [
+        [4, 8],
+        [3, 9],
+        [9, 3],
+        [8, 4],
+        [9, 8],
+        [8, 9],
+        [4, 8],
+        [8, 4],
+        [8, 8],
+        [3, 8],
+        [8, 3],
+        [8, 8],
+    ]
+    V = [
+        [3, 9],
+        [3, 8],
+        [9, 8],
+        [8, 9],
+        [8, 3],
+        [9, 3],
+    ]
+    p = polyhedron(vrep(W), lib)
+    generator_fulltest(p, V)
+end

test/misc.jl

@@ -6,6 +6,7 @@ include("docexample.jl")
 include("issue48.jl")
 include("issue224.jl")
 include("issue301.jl")
+include("issue311.jl")
 include("empty.jl")
 include("sparse.jl")
 include("sparserect.jl")
@@ -25,6 +26,7 @@ const misctests = Dict(
     "issue48" => issue48test,
     "issue224" => issue224test,
     "issue301" => issue301test,
+    "issue311" => issue311test,
     "empty" => emptytest,
     "sparse" => sparsetest,
     "sparserect" => sparserecttest,

test/planar.jl

@@ -8,7 +8,7 @@ function test_planar_square(VT, d)
     no = -ei[1] - ei[2]
     ne = ei[1] - ei[2]
     se = ei[1] + ei[2]
-    expected = [ne, no, so, se]
+    expected = [no, so, se, ne]
     for ps in [
         [no, so, ne, se],
         [no, se, so, ne],
@@ -61,7 +61,7 @@ end
 
 # Square + redundant points + [+-1.5, -+0.25]
 function test_planar_square_with_more()
-    expected = [[1.0, -1.0], [-1.0, -1.0], [-1.0, 1.0], [1.0, 1.0]]
+    expected = [[-1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [1.0, -1.0]]
     v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0])
     v1 = Polyhedra.planar_hull(v0)
     @test collect(points(v1)) == expected
@@ -72,10 +72,10 @@ function test_planar_square_with_more()
     v1 = Polyhedra.planar_hull(v0)
     @test collect(points(v1)) == expected
     expected = [[1.0, -1.0], [-1.0, -1.0], [-1.0, 1.0], [0.0, 1.0 + 1e-4], [1.0, 1.0]]
-    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0 + 1e-4])
+    v0 = convexhull([1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0 + 1e-4], [-1.0, -1.0])
     v1 = Polyhedra.planar_hull(v0)
     @test collect(points(v1)) == expected
-    expected = [[1.0, -1.0], [-1.0, -1.0], [-1.5, 0.25], [-1.0, 1.0], [1.0, 1.0], [1.5, -0.25]]
+    expected = [[-1.0, -1.0], [-1.5, 0.25], [-1.0, 1.0], [1.0, 1.0], [1.5, -0.25], [1.0, -1.0]]
     v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0], [-1.5, 0.25], [1.5, -0.25])
     v1 = Polyhedra.planar_hull(v0)
     @test collect(points(v1)) == expected
@@ -91,11 +91,11 @@ function test_planar_square_with_more()
 end
 
 function test_issue_271()
-    expected = [[1.0, 0.0], [-0.5, 0.0], [0.5, 0.5]]
+    expected = [[-0.5, 0.0], [0.5, 0.5], [1.0, 0.0]]
     v = convexhull([0.0, 0.0], [1, 0], [0.5, 0.0], [-0.5, 0.0], [0.5, 0.5])
     p = Polyhedra.planar_hull(v)
     @test collect(points(p)) == expected
-    v = convexhull([1, 0], [0.0, 0.0], [0.5, 0.0], [-0.5, 0.0], [0.5, 0.5])
+    v = convexhull([0.0, 0.0], [0.5, 0.0], [-0.5, 0.0], [0.5, 0.5], [1, 0])
     p = Polyhedra.planar_hull(v)
     @test collect(points(p)) == expected
 

test/utils.jl

@@ -1,4 +1,4 @@
-using Test
+using LinearAlgebra, Test
 using Polyhedra
 
 _isapprox(x::Real, y::Real) = _isapprox(promote(x, y)...)