Merge pull request #376 from JuliaReach/schillic/5

#5 - Implement hrep from vrep for VPolygon

docs/src/lib/representations.md

@@ -98,6 +98,8 @@ dim(::HalfSpace{Float64})
 σ(::AbstractVector{Float64}, ::HalfSpace{Float64})
 an_element(::HalfSpace{Float64})
 ∈(::AbstractVector{Float64}, ::HalfSpace{Float64})
+LazySets.halfspace_left(::AbstractVector{Float64}, ::AbstractVector{Float64})
+LazySets.halfspace_right(::AbstractVector{Float64}, ::AbstractVector{Float64})
 ```
 
 ## Hyperplane
@@ -164,6 +166,8 @@ LineSegment
 dim(::LineSegment{Float64})
 σ(::AbstractVector{Float64}, ::LineSegment{Float64})
 ∈(::AbstractVector{Float64}, ::LineSegment{Float64})
+LazySets.halfspace_left(::LineSegment)
+LazySets.halfspace_right(::LineSegment)
 ```
 
 ## Polygons

src/HalfSpace.jl

@@ -116,3 +116,94 @@ We just check if ``x`` satisfies ``a⋅x ≤ b``.
 function ∈(x::AbstractVector{N}, hs::HalfSpace{N})::Bool where {N<:Real}
     return dot(x, hs.a) <= hs.b
 end
+
+"""
+    halfspace_left(p::AbstractVector{N},
+                   q::AbstractVector{N})::HalfSpace{N} where {N<:Real}
+
+Return a half-space describing the 'left' of a two-dimensional line segment through
+two points.
+
+### Input
+
+- `p` -- first point
+- `q` -- second point
+
+### Output
+
+The half-space whose boundary goes through the two points `p` and `q` and which
+describes the left-hand side of the directed line segment `pq`.
+
+### Algorithm
+
+The implementation is simple: the half-space ``a⋅x ≤ b`` is calculated as
+`a = [dy, -dx]`, where ``d = (dx, dy)`` denotes the line segment
+`pq`, that is, ``\\vec{d} = \\vec{p} - \\vec{q}``, and `b = dot(p, a)`.
+
+### Examples
+
+The left half-space of the "east" and "west" directions in two-dimensions are the
+upper and lower half-spaces:
+
+```jldoctest halfspace_left
+julia> import LazySets.halfspace_left
+
+julia> halfspace_left([0.0, 0.0], [1.0, 0.0])
+LazySets.HalfSpace{Float64}([0.0, -1.0], 0.0)
+
+julia> halfspace_left([0.0, 0.0], [-1.0, 0.0])
+LazySets.HalfSpace{Float64}([0.0, 1.0], 0.0)
+```
+
+We create a box from the sequence of line segments that describe its edges:
+
+```jldoctest halfspace_left
+julia> H1 = halfspace_left([-1.0, -1.0], [1.0, -1.0]);
+
+julia> H2 = halfspace_left([1.0, -1.0], [1.0, 1.0]);
+
+julia> H3 = halfspace_left([1.0, 1.0], [-1.0, 1.0]);
+
+julia> H4 = halfspace_left([-1.0, 1.0], [-1.0, -1.0]);
+
+julia> H = HPolygon([H1, H2, H3, H4]);
+
+julia> B = BallInf([0.0, 0.0], 1.0);
+
+julia> B ⊆ H && H ⊆ B
+true
+```
+"""
+function halfspace_left(p::AbstractVector{N},
+                        q::AbstractVector{N})::HalfSpace{N} where {N<:Real}
+    @assert length(p) == length(q) == 2 "the points must be two-dimensional"
+    @assert p != q "the points must not be equal"
+    a = [q[2] - p[2], p[1] - q[1]]
+    return HalfSpace(a, dot(p, a))
+end
+
+"""
+    halfspace_right(p::AbstractVector{N},
+                    q::AbstractVector{N})::HalfSpace{N} where {N<:Real}
+
+Return a half-space describing the 'right' of a two-dimensional line segment through
+two points.
+
+### Input
+
+- `p` -- first point
+- `q` -- second point
+
+### Output
+
+The half-space whose boundary goes through the two points `p` and `q` and which
+describes the right-hand side of the directed line segment `pq`.
+
+### Algorithm
+
+See the documentation of `halfspace_left`. 
+"""
+function halfspace_right(p::AbstractVector{N},
+                         q::AbstractVector{N})::HalfSpace{N} where {N<:Real}
+    return halfspace_right(q, p)
+end

src/LazySets.jl

@@ -26,8 +26,8 @@ include("AbstractPointSymmetric.jl")
 include("AbstractPointSymmetricPolytope.jl")
 include("AbstractHyperrectangle.jl")
 include("AbstractPolygon.jl")
-include("AbstractHPolygon.jl")
 include("AbstractSingleton.jl")
+include("AbstractHPolygon.jl")
 
 # =============================
 # Types representing basic sets

src/LineSegment.jl

@@ -144,3 +144,37 @@ function ∈(x::AbstractVector{N}, L::LineSegment{N})::Bool where {N<:Real}
     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
+
+"""
+    halfspace_left(L::LineSegment)
+
+Return a half-space describing the 'left' of a two-dimensional line segment through
+two points.
+
+### Input
+
+ - `L` -- line segment
+
+### Output
+
+The half-space whose boundary goes through the two points `p` and `q` and which
+describes the left-hand side of the directed line segment `pq`.
+"""
+halfspace_left(L::LineSegment) = halfspace_left(L.p, L.q)
+
+"""
+    halfspace_right(L::LineSegment)
+
+Return a half-space describing the 'right' of a two-dimensional line segment through
+two points.
+
+### Input
+
+ - `L` -- line segment
+
+### Output
+
+The half-space whose boundary goes through the two points `p` and `q` and which
+describes the right-hand side of the directed line segment `pq`.
+"""
+halfspace_right(L::LineSegment) = halfspace_right(L.p, L.q)

src/VPolygon.jl

@@ -59,20 +59,62 @@ function tovrep(P::VPolygon{N})::VPolygon{N} where {N<:Real}
 end
 
 """
-    tohrep(P::VPolygon{N})::AbstractHPolygon{N} where {N<:Real}
+    tohrep(P::VPolygon{N}, ::Type{HPOLYGON}=HPolygon
+          )::AbstractHPolygon{N} where {N<:Real, HPOLYGON<:AbstractHPolygon}
 
 Build a constraint representation of the given polygon.
 
 ### Input
 
-- `P` -- polygon in vertex representation
+- `P`        -- polygon in vertex representation
+- `HPOLYGON` -- (optional, default: `HPolygon`) type of target polygon
 
 ### Output
 
 The same polygon but in constraint representation, an `AbstractHPolygon`.
+
+### Algorithm
+
+The algorithms consists of adding an edge for each consecutive pair of vertices.
+Since the vertices are already ordered in counter-clockwise fashion (CWW), the
+constraints will be sorted automatically (CCW) if we start with the first edge
+between the first and second vertex.
 """
-function tohrep(P::VPolygon{N})::AbstractHPolygon{N} where {N<:Real}
-    error("this function is not implemented yet, see issue #5")
+function tohrep(P::VPolygon{N}, ::Type{HPOLYGON}=HPolygon
+               )::AbstractHPolygon{N} where {N<:Real, HPOLYGON<:AbstractHPolygon}
+    vl = vertices_list(P)
+    n = length(vl)
+    if n == 0
+        # no vertex -> no constraint
+        constraints_list = Vector{LinearConstraint{N}}(0)
+    elseif n == 1
+        # only one vertex -> use function for singletons
+        return convert(HPOLYGON, Singleton(vl[1]))
+    elseif n == 2
+        # only two vertices -> use function for line segments
+        return convert(HPOLYGON, LineSegment(vl[1], vl[2]))
+    else
+        # find right-most vertex
+        i = div(n, 2)
+        x = vl[i][1]
+        while i > 1 && vl[i-1][1] > x
+            # search forward in list
+            i = i - 1
+            x = vl[i][1]
+        end
+        while i < n && vl[i+1][1] > x
+            # search backward in list
+            i = i + 1
+            x = vl[i][1]
+        end
+
+        # create constraints ordered in CCW starting at the right-most index
+        upper_hull = [halfspace_left(vl[j], vl[j+1]) for j in i:length(vl)-1]
+        mid_hull = [halfspace_left(vl[end], vl[1])]
+        lower_hull = [halfspace_left(vl[j], vl[j+1]) for j in 1:i-1]
+        constraints_list = vcat(upper_hull, mid_hull, lower_hull)
+    end
+    return HPOLYGON{N}(constraints_list)
 end
 
 

src/convert.jl

@@ -89,6 +89,64 @@ function convert(::Type{Zonotope}, H::AbstractHyperrectangle{N}) where {N}
     return Zonotope{N}(center(H), diagm(radius_hyperrectangle(H)))
 end
 
+"""
+    convert(::Type{HPOLYGON}, S::AbstractSingleton{N}
+           ) where {N, HPOLYGON<:AbstractHPolygon}
+
+Convert from singleton to polygon in H-representation.
+
+### Input
+
+- `type` -- target type
+- `S`    -- singleton
+
+### Output
+
+A polygon in constraint representation with the minimal number of constraints
+(three).
+"""
+function convert(::Type{HPOLYGON}, S::AbstractSingleton{N}
+                ) where {N, HPOLYGON<:AbstractHPolygon}
+    constraints_list = Vector{LinearConstraint{N}}(3)
+    o = one(N)
+    z = zero(N)
+    v = element(S)
+    constraints_list[1] = LinearConstraint([o, o], v[1] + v[2])
+    constraints_list[2] = LinearConstraint([-o, z], -v[1])
+    constraints_list[3] = LinearConstraint([z, -o], -v[2])
+    return HPOLYGON{N}(constraints_list)
+end
+
+"""
+    convert(::Type{HPOLYGON}, L::LineSegment{N}
+          ) where {N, HPOLYGON<:AbstractHPolygon}
+
+Convert from line segment to polygon in H-representation.
+
+### Input
+
+- `type` -- target type
+- `L`    -- line segment
+
+### Output
+
+A flat polygon in constraint representation with the minimal number of
+constraints (four).
+"""
+function convert(::Type{HPOLYGON}, L::LineSegment{N}
+                ) where {N, HPOLYGON<:AbstractHPolygon}
+    H = HPOLYGON{N}()
+    c = halfspace_left(L.p, L.q)
+    addconstraint!(H, c)
+    addconstraint!(H, LinearConstraint(-c.a, -c.b))
+    line_dir = L.q - L.p
+    c = LinearConstraint(line_dir, dot(L.q, line_dir))
+    addconstraint!(H, c)
+    line_dir = -line_dir
+    addconstraint!(H, LinearConstraint(line_dir, dot(L.p, line_dir)))
+    return H
+end
+
 import IntervalArithmetic.AbstractInterval
 
 """

test/unit_Polygon.jl

@@ -33,9 +33,9 @@ for N in [Float64, Float32, Rational{Int}]
     HPolytope(po)
 
     # support vector of empty polygon
-    @test_throws ErrorException σ([0.], HPolygon())
-    @test_throws ErrorException σ([0.], HPolygonOpt(HPolygon()))
-    @test_throws ErrorException σ([0.], HPolytope())
+    @test_throws ErrorException σ(N[0.], HPolygon{N}())
+    @test_throws ErrorException σ(N[0.], HPolygonOpt(HPolygon{N}()))
+    @test_throws ErrorException σ(N[0.], HPolytope{N}())
 
     # HPolygon/HPolygonOpt tests
     for p in [p, po]
@@ -145,6 +145,45 @@ for N in [Float64, Float32, Rational{Int}]
                 @test vi <= v[j]
         end
     end
+
+    # hrep conversion
+    v1 = to_N(N, [0.9, 0.2])
+    v2 = to_N(N, [0.4, 0.6])
+    v3 = to_N(N, [0.2, 0.1])
+    v4 = to_N(N, [0.1, 0.3])
+    points5 = [v1, v2, v3, v4]
+    for i in [0, 1, 2, 4]
+        points = i == 0 ? Vector{Vector{N}}() : points5[1:i]
+        vp = VPolygon(points, apply_convex_hull=i > 0)
+        h1 = tohrep(vp)
+        if i == 0
+            @test isempty(h1.constraints)
+        elseif i == 1
+            @test v1 ∈ h1
+        elseif i == 2
+            c = h1.constraints[1]
+            @test c.a ≈ to_N(N, [0.4, 0.5])&& c.b ≈ to_N(N, (0.46))
+            c = h1.constraints[3]
+            @test c.a ≈ to_N(N, [-0.4, -0.5])&& c.b ≈ to_N(N, (-0.46))
+        elseif i == 4
+            @test length(h1.constraints) == 4
+            c = h1.constraints[1]
+            @test c.a ≈ to_N(N, [0.4, 0.5])&& c.b ≈ to_N(N, (0.46))
+            c = h1.constraints[2]
+            @test c.a ≈ to_N(N, [-0.3, 0.3])&& c.b ≈ to_N(N, (0.06))
+            c = h1.constraints[3]
+            @test c.a ≈ to_N(N, [-0.2, -0.1])&& c.b ≈ to_N(N, (-0.05))
+            c = h1.constraints[4]
+            @test c.a ≈ to_N(N, [0.1, -0.7])&& c.b ≈ to_N(N, (-0.05))
+        end
+
+        # check that constraints are sorted correctly
+        h2 = HPolygon{N}()
+        for c in h1.constraints
+            addconstraint!(h2, c)
+        end
+        @test h1.constraints == h2.constraints
+    end
 end
 
 # default Float64 constructors