Merge pull request #251 from JuliaReach/schillic/231

#231 - LineSegment

docs/src/lib/binary_functions.md

@@ -41,4 +41,5 @@ is_intersection_empty(::AbstractSingleton{Float64}, ::AbstractSingleton{Float64}
 is_intersection_empty(::Zonotope{Float64}, ::Hyperplane{Float64})
 is_intersection_empty(::Hyperplane{Float64}, ::Zonotope{Float64})
 is_intersection_empty(::Ball2{Float64}, ::Ball2{Float64})
+is_intersection_empty(::LineSegment{Float64}, ::LineSegment{Float64})
 ```

docs/src/lib/representations.md

@@ -135,9 +135,20 @@ low(::Hyperrectangle{Float64})
 
 ```@docs
 Line
+dim(::Line{Float64})
+σ(::AbstractVector{Float64}, ::Line{Float64})
 intersection(::Line{Float64}, ::Line{Float64})
 ```
 
+## Line segment
+
+```@docs
+LineSegment
+dim(::LineSegment{Float64})
+σ(::AbstractVector{Float64}, ::LineSegment{Float64})
+∈(::AbstractVector{Float64}, ::LineSegment{Float64})
+```
+
 ## Polygons
 
 ### Constraint representation

src/LazySet.jl

@@ -30,7 +30,7 @@ Every concrete `LazySet` must define the following functions:
 
 ```jldoctest
 julia> subtypes(LazySet)
-17-element Array{Union{DataType, UnionAll},1}:
+18-element Array{Union{DataType, UnionAll},1}:
  LazySets.AbstractPointSymmetric
  LazySets.AbstractPolytope
  LazySets.CartesianProduct
@@ -45,6 +45,7 @@ julia> subtypes(LazySet)
  LazySets.Intersection
  LazySets.Line
  LazySets.LinearMap
+ LazySets.LineSegment
  LazySets.MinkowskiSum
  LazySets.MinkowskiSumArray
  LazySets.PolynomialZonotope

src/LazySets.jl

@@ -42,6 +42,7 @@ include("HPolytope.jl")
 include("Hyperplane.jl")
 include("Hyperrectangle.jl")
 include("Line.jl")
+include("LineSegment.jl")
 include("Singleton.jl")
 include("VPolygon.jl")
 include("VPolytope.jl")

src/LineSegment.jl

@@ -0,0 +1,146 @@
+import Base.∈
+
+export LineSegment
+
+"""
+    LineSegment{N<:Real} <: LazySet{N}
+
+Type that represents a line segment in 2D between two points ``p`` and ``q``.
+
+### Fields
+
+- `p` -- first point
+- `q` -- second point
+
+### Examples
+
+A line segment along the ``x = y`` diagonal:
+
+```jldoctest linesegment_constructor
+julia> s = LineSegment([0., 0], [1., 1.])
+LazySets.LineSegment{Float64}([0.0, 0.0], [1.0, 1.0])
+julia> dim(s)
+2
+```
+
+Use ``plot(s)`` to plot the extreme points of `s` and the line segment joining
+them.
+Membership test is computed with ∈ (`in`):
+
+```jldoctest linesegment_constructor
+julia> [0., 0] ∈ s && [.25, .25] ∈ s && [1., 1] ∈ s && !([.5, .25] ∈ s)
+true
+```
+
+We can check the intersection with another line segment, and optionally compute
+a witness (which is just the common point in this case):
+
+```jldoctest linesegment_constructor
+julia> sn = LineSegment([1., 0], [0., 1.])
+LazySets.LineSegment{Float64}([1.0, 0.0], [0.0, 1.0])
+julia> isempty(s ∩ sn)
+false
+julia> is_intersection_empty(s, sn, true)
+(false, [0.5, 0.5])
+```
+"""
+struct LineSegment{N<:Real} <: LazySet{N}
+    p::AbstractVector{N}
+    q::AbstractVector{N}
+
+    # default constructor with length constraint
+    LineSegment{N}(p::AbstractVector{N}, q::AbstractVector{N}) where {N<:Real} =
+        (length(p) == length(q) == 2 ? new{N}(p, q) : throw(DimensionMismatch))
+end
+
+# type-less convenience constructor
+LineSegment(p::AbstractVector{N}, q::AbstractVector{N}) where {N<:Real} =
+    LineSegment{N}(p, q)
+
+
+# --- LazySet interface functions ---
+
+
+"""
+    dim(L::LineSegment)::Int
+
+Return the ambient dimension of a line segment.
+
+### Input
+
+- `L` -- line segment
+
+### Output
+
+The ambient dimension of the line segment, which is 2.
+"""
+function dim(L::LineSegment)::Int
+    return 2
+end
+
+"""
+    σ(d::AbstractVector{<:Real}, L::LineSegment)::AbstractVector{<:Real}
+
+Return the support vector of a line segment in a given direction.
+
+### Input
+
+- `d` -- direction
+- `L` -- line segment
+
+### Output
+
+The support vector in the given direction.
+
+### Algorithm
+
+If the angle between the vector ``q - p`` and ``d`` is bigger than 90° and less
+than 270° (measured in counter-clockwise order), the result is ``p``, otherwise
+it is ``q``.
+If the angle is exactly 90° or 270°, or if the direction has norm zero, this
+implementation returns ``q``.
+"""
+function σ(d::AbstractVector{<:Real}, L::LineSegment)::AbstractVector{<:Real}
+    return sign(dot(L.q - L.p, d)) >= 0 ? L.q : L.p
+end
+
+"""
+    ∈(x::AbstractVector{N}, L::LineSegment{N})::Bool where {N<:Real}
+
+Check whether a given point is contained in a line segment.
+
+### Input
+
+- `x` -- point/vector
+- `L` -- line segment
+
+### Output
+
+`true` iff ``x ∈ L``.
+
+### Algorithm
+
+Let ``L = (p, q)`` be the line segment with extremes ``p`` and ``q``, and let
+``x`` be the given point.
+
+1. A necessary condition for ``x ∈ (p, q)`` is that the three points are aligned,
+   thus their cross product should be zero.
+2. It remains to check that ``x`` belongs to the box approximation of ``L``.
+   This amounts to comparing each coordinate with those of the extremes ``p``
+   and ``q``.
+
+### Notes
+
+The algorithm is inspired from [here](https://stackoverflow.com/a/328110).
+"""
+function ∈(x::AbstractVector{N}, L::LineSegment{N})::Bool where {N<:Real}
+    @assert length(x) == dim(L)
+    # check if the point is on the line through the line segment
+    if (x[2] - L.p[2]) * (L.q[1] - L.p[1]) -
+            (x[1] - L.p[1]) * (L.q[2] - L.p[2]) != 0
+        return false
+    end
+    # check if the point is inside the box approximation of the line segment
+    return min(L.p[1], L.q[1]) <= x[1] <= max(L.p[1], L.q[1]) &&
+           min(L.p[2], L.q[2]) <= x[2] <= max(L.p[2], L.q[2])
+end

src/is_intersection_empty.jl

@@ -418,3 +418,124 @@ function is_intersection_empty(H::Hyperplane{N},
                               )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
     return is_intersection_empty(Z, H, witness)
 end
+
+"""
+    is_intersection_empty(ls1::LineSegment{N},
+                          ls2::LineSegment{N},
+                          witness::Bool=false
+                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+
+Check whether two line segments do not intersect, and otherwise optionally
+compute a witness.
+
+### Input
+
+- `ls1` -- first line segment
+- `ls2` -- second line segment
+- `witness` -- (optional, default: `false`) compute a witness if activated
+
+### Output
+
+* If `witness` option is deactivated: `true` iff ``ls1 ∩ ls2 = ∅``
+* If `witness` option is activated:
+  * `(true, [])` iff ``ls1 ∩ ls2 = ∅``
+  * `(false, v)` iff ``ls1 ∩ ls2 ≠ ∅`` and ``v ∈ ls1 ∩ ls2``
+
+### Algorithm
+
+The algorithm is inspired from [here](https://stackoverflow.com/a/565282), which
+again is the special 2D case of a 3D algorithm by Ronald Goldman's article on the
+*Intersection of two lines in three-space* in Graphics Gems, Andrew S. (ed.), 1990.
+
+We first check if the two line segments are parallel, and if so, if they are
+collinear.
+In the latter case, we check containment of any of the end points in the other
+line segment.
+Otherwise the lines are not parallel, so we can solve an equation of the
+intersection point, if it exists.
+"""
+function is_intersection_empty(ls1::LineSegment{N},
+                               ls2::LineSegment{N},
+                               witness::Bool=false
+                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    function cross(ls1::Vector{N}, ls2::Vector{N})::N where {N<:Real}
+        return ls1[1] * ls2[2] - ls1[2] * ls2[1]
+    end
+
+    r = ls1.q - ls1.p
+    if iszero(r)
+        # first line segment is a point
+        empty_intersection = !∈(ls1.q, ls2)
+        if witness
+            return (empty_intersection, empty_intersection ? N[] : ls1.q)
+        else
+            return empty_intersection
+        end
+    end
+
+    s = ls2.q - ls2.p
+    if iszero(s)
+        # second line segment is a point
+        empty_intersection = !∈(ls2.q, ls1)
+        if witness
+            return (empty_intersection, empty_intersection ? N[] : ls2.q)
+        else
+            return empty_intersection
+        end
+    end
+
+    p1p2 = ls2.p - ls1.p
+    u_numerator = cross(p1p2, r)
+    u_denominator = cross(r, s)
+
+    if u_denominator == 0
+        # line segments are parallel
+        if u_numerator == 0
+            # line segments are collinear
+            if ∈(ls1.p, ls2)
+                empty_intersection = false
+                if witness
+                    v = ls1.p
+                end
+            elseif ∈(ls1.q, ls2)
+                empty_intersection = false
+                if witness
+                    v = ls1.q
+                end
+            elseif ∈(ls2.p, ls1)
+                empty_intersection = false
+                if witness
+                    v = ls2.p
+                end
+            elseif ∈(ls2.q, ls1)
+                empty_intersection = false
+                if witness
+                    v = ls2.q
+                end
+            else
+                empty_intersection = true
+            end
+        else
+            # line segments are parallel and not collinear
+            empty_intersection = true
+        end
+    else
+        # line segments are not parallel
+        u = u_numerator / u_denominator
+        if u < 0 || u > 1
+            empty_intersection = true
+        else
+            t = cross(p1p2, s) / u_denominator
+            empty_intersection = t < 0 || t > 1
+            if witness
+                v = ls1.p + t * r
+            end
+        end
+    end
+
+    if witness
+        return (empty_intersection, empty_intersection ? N[] : v)
+    else
+        return empty_intersection
+    end
+end

src/is_subset.jl

@@ -365,8 +365,8 @@ end
               witness::Bool=false
              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
 
-Check whether a ball in the 2-norm is contained in a set with a single value,
-and if not, optionally compute a witness.
+Check whether a ball in the 2-norm or p-norm is contained in a set with a single
+value, and if not, optionally compute a witness.
 
 ### Input
 
@@ -405,6 +405,52 @@ function is_subset(B::Union{Ball2{N}, Ballp{N}},
 end
 
 
+# --- LineSegment ---
+
+# TODO commented until AbstractConvexSet is implemented
+# """
+#     is_subset(L::LineSegment{N},
+#               S::AbstractConvexSet{N},
+#               witness::Bool=false
+#              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+#
+# Check whether a line segment is contained in a convex set, and if not,
+# optionally compute a witness.
+#
+# ### Input
+#
+# - `L` -- inner line segment
+# - `S` -- outer convex set
+# - `witness` -- (optional, default: `false`) compute a witness if activated
+#
+# ### Output
+#
+# * If `witness` option is deactivated: `true` iff ``L ⊆ S``
+# * If `witness` option is activated:
+#   * `(true, [])` iff ``L ⊆ S``
+#   * `(false, v)` iff ``L \\not\\subseteq S`` and ``v ∈ L \\setminus S``
+#
+# ### Algorithm
+#
+# Since ``S`` is convex, ``L ⊆ S`` iff ``p ∈ S`` and ``q ∈ S``, where ``p, q`` are
+# the end points of ``L``.
+# """
+# function is_subset(L::LineSegment{N},
+#                    S::AbstractConvexSet{N},
+#                    witness::Bool=false
+#                   )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+#     p_in_S = ∈(L.p, S)
+#     result = p_in_S && ∈(L.q, S)
+#     if !witness
+#         return result
+#     elseif result
+#         return (result, N[])
+#     else
+#         return (result, p_in_S ? L.q : L.p)
+#     end
+# end
+
+
 # --- alias ---