#3126 - Point in Polygon check

src/Sets/Polygon/PolygonModule.jl

@@ -6,7 +6,7 @@ using ..LazySets: LazySet, _plot_recipe_2d_vlist
 using ReachabilityBase.Require: require
 
 @reexport import ..API: convex_hull, dim, isconvextype, isbounded,
- isboundedtype, isempty, isoperationtype, ρ, σ
+ isboundedtype, isempty, isoperationtype, ∈, ρ, σ
 import ..LazySets: plot_recipe
 @reexport using ..API
 
@@ -21,6 +21,7 @@ include("isbounded.jl")
 include("isboundedtype.jl")
 include("isempty.jl")
 include("isoperationtype.jl")
+include("in.jl")
 include("support_function.jl")
 include("support_vector.jl")
 

src/Sets/Polygon/in.jl

@@ -0,0 +1,53 @@
+# Algorithm:
+# Choose an arbitrary ray through `x`` and count whether the number of
+# intersections with edges is odd. We choose the ray that is vertical upward.
+# Vertical line segments are ignored (after checking whether `x` is a member).
+function ∈(x::AbstractVector, P::Polygon)
+ require(@__MODULE__, :LazySets; fun_name="in")
+ @assert length(x) == dim(P) "a vector of length $(length(x)) is " *
+ "incompatible with a $(dim(P))-dimensional set"
+
+ N = promote_type(eltype(x), eltype(P))
+
+ vlist = P.vertices
+ if length(vlist) < 2
+ if isempty(vlist)
+ return false
+ elseif length(vlist) == 1
+ return @inbounds x == vlist[1]
+ end
+ end
+
+ # vertical ray as a vertical line (compare y coordinates later)
+ vline = Line2D([one(N), zero(N)], @inbounds x[1])
+
+ p = @inbounds vlist[end]
+ odd = false
+ @inbounds for q in vlist
+ if (p[1] <= x[1] && q[1] >= x[1]) || (q[1] <= x[1] && p[1] >= x[1])
+ # line segment pq intersects vertical line through x
+ if p[1] == q[1]
+ # vertical line segment
+ if (p[2] <= x[2] && q[2] >= x[2]) || (q[2] <= x[2] && p[2] >= x[2])
+ # x is on the line segment
+ return true
+ end
+ else
+ # non-vertical line segment -> intersect (Line2D intersection is used)
+ line2 = Line2D(p, q)
+ y = element(intersection(line2, vline))
+ # compare y coordinate
+ if y[2] >= x[2]
+ if y == x
+ # x is on line segment
+ return true
+ end
+ odd = !odd
+ end
+ end
+ end
+
+ p = q
+ end
+ return odd
+end

src/Sets/Polygon/init_LazySets.jl

@@ -1,2 +1,4 @@
+using .LazySets: element
+using .LazySets.Line2DModule: Line2D
 using .LazySets.VPolygonModule: VPolygon
 using .LazySets.VPolytopeModule: _ρ_vertices, _σ_vertices

test/Sets/PolygonNC.jl

@@ -27,6 +27,17 @@ for N in [Float64, Rational{Int}, Float32]
  # plot_recipe (only check for correct output type)
  x, y = LazySets.plot_recipe(P)
  @test x isa Vector{N} && y isa Vector{N}
+
+ # membership
+ P1 = Polygon([N[0, 0], N[2, 1], N[18//10, 12//10], N[1, 1], N[1//2, 3//2], N[-1//2, 2]]);
+ P2 = Polygon([N[0, 1], N[2, 1]]);
+ P3 = Polygon([N[3//2, 0], N[3//2, 11//10]]);
+ x = N[3//2, 1];
+ y = N[3//2, 3//2];
+ for P in (P1, P2, P3)
+ @test x ∈ P
+ @test y ∉ P
+ end
 end
 
 # default Float64 constructor