Generalize check if point is in hyperplane (#3304)

* generalize check if point is in hyperplane

* Update test/Sets/Tetrahedron.jl

Co-authored-by: Christian Schilling <git@christianschilling.net>

---------

Co-authored-by: Marcelo Forets <marcelo@planting.space>
Co-authored-by: Christian Schilling <git@christianschilling.net>

Author: 3 people

src/Sets/Tetrahedron.jl

@@ -12,7 +12,7 @@ Type that represents a tetrahedron in vertex representation.
 ### Examples
 
 A tetrahedron can be constructed by passing the list of vertices.
-The following builds the tetrahedron with edge length 2 from the [wikipedia page Tetrahedron](https://en.wikipedia.org/wiki/Tetrahedron): 
+The following builds the tetrahedron with edge length 2 from the [wikipedia page Tetrahedron](https://en.wikipedia.org/wiki/Tetrahedron):
 
 ```jldoctest
 julia> vertices = [[1, 0, -1/sqrt(2)], [-1, 0, -1/sqrt(2)], [0, 1, 1/sqrt(2)], [0, -1, 1/sqrt(2)]];
@@ -110,20 +110,22 @@ function ∈(x::AbstractVector, T::Tetrahedron)
            same_side(v[2], v[3], v[4], v[1], x)
 end
 
-# Return `true` iff point `x` lies in the same half-space as `v4` with respect to the hyperplane determined by `v1`, `v2` and `v3`.
+# Return `true` iff point `x` lies in the same half-space as `v4` with respect to the hyperplane `H` determined by `v1`, `v2` and `v3`.
 function same_side(v1, v2, v3, v4, x)
     normal = cross(v2 - v1, v3 - v1)
     dotv4 = dot(normal, v4 - v1)
     dotx = dot(normal, x - v1)
-    return signbit(dotv4) == signbit(dotx)
+    # If the point `x` lies in `H` then `dotx` is zero.
+    return signbit(dotv4) == signbit(dotx) || isapproxzero(dotx)
 end
 
-function rand(::Type{Tetrahedron}; N::Type{<:Real}=Float64, rng::AbstractRNG=GLOBAL_RNG, seed::Union{Int,Nothing}=nothing)
+function rand(::Type{Tetrahedron}; N::Type{<:Real}=Float64, rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int,Nothing}=nothing)
     P = rand(VPolytope; N=N, dim=3, rng=rng, seed=seed, num_vertices=4)
     vertices = P.vertices
     return Tetrahedron(vertices)
 end
 
 function constraints_list(T::Tetrahedron)
-    constraints_list(convert(VPolytope, T))
+    return constraints_list(convert(VPolytope, T))
 end

test/Sets/Tetrahedron.jl

@@ -1,6 +1,7 @@
 for N in [Float64, Float32, Rational{Int}]
     # constructor
-    vertices = [N[1, 0, -1/sqrt(2)], N[-1, 0, -1/sqrt(2)], N[0, 1, 1/sqrt(2)], N[0, -1, 1/sqrt(2)]]
+    vertices = [N[1, 0, -1 / sqrt(2)], N[-1, 0, -1 / sqrt(2)], N[0, 1, 1 / sqrt(2)],
+                N[0, -1, 1 / sqrt(2)]]
     T = Tetrahedron(vertices)
     @test eltype(T) == N
     @test dim(T) == 3
@@ -10,11 +11,15 @@ for N in [Float64, Float32, Rational{Int}]
     @test Tmat == T
 
     # suport vector
-    @test σ(ones(3), T) == N[0, 1, 1/sqrt(2)]
+    @test σ(ones(3), T) == N[0, 1, 1 / sqrt(2)]
 
     # membership
     @test zeros(N, 3) ∈ T
 
     # is_polyhedral
     @test is_polyhedral(T)
+
+    # LazySets#3303
+    T = Tetrahedron([N[0, 0, 0], N[0, 1, 0], N[1, 0, 0], N[1, 0, 1]])
+    @test zeros(N, 3) ∈ T
 end