diff --git a/docs/src/lib/sets/Line.md b/docs/src/lib/sets/Line.md
index d9cdbaaabf..4a56b3b0bd 100644
--- a/docs/src/lib/sets/Line.md
+++ b/docs/src/lib/sets/Line.md
@@ -1,35 +1,45 @@
 ```@meta
-CurrentModule = LazySets
+CurrentModule = LazySets.LineModule
 ```
 
 # [Line](@id def_Line)
 
 ```@docs
 Line
-dim(::Line)
-ρ(::AbstractVector, ::Line)
-σ(::AbstractVector, ::Line)
-∈(::AbstractVector, ::Line)
+```
+
+## Operations
+
+```@docs
 an_element(::Line)
+constraints_list(::Line)
+dim(::Line)
 direction(::Line)
-rand(::Type{Line})
 isbounded(::Line)
-isuniversal(::Line; ::Bool=false)
 isempty(::Line)
-constraints_list(::Line)
-translate!(::Line, ::AbstractVector)
+isuniversal(::Line; ::Bool=false)
 normalize(::Line{N}, ::Real=N(2)) where{N}
 normalize!(::Line{N}, ::Real=N(2)) where{N}
+rand(::Type{Line})
 distance(::AbstractVector, ::Line; ::Real=2.0)
+∈(::AbstractVector, ::Line)
 linear_map(::AbstractMatrix, ::Line)
+ρ(::AbstractVector, ::Line)
+σ(::AbstractVector, ::Line)
+translate!(::Line, ::AbstractVector)
 ```
+
+```@meta
+CurrentModule = LazySets
+```
+
 Inherited from [`LazySet`](@ref):
-* [`low`](@ref low(::LazySet))
+* [`diameter`](@ref diameter(::LazySet, ::Real))
 * [`high`](@ref high(::LazySet))
+* [`high`](@ref high(::LazySet, ::Int))
+* [`low`](@ref low(::LazySet))
+* [`low`](@ref low(::LazySet, ::Int))
 * [`norm`](@ref norm(::LazySet, ::Real))
 * [`radius`](@ref radius(::LazySet, ::Real))
-* [`diameter`](@ref diameter(::LazySet, ::Real))
-* [`low`](@ref low(::LazySet, ::Int))
-* [`high`](@ref high(::LazySet, ::Int))
 * [`reflect`](@ref reflect(::LazySet))
 * [`translate`](@ref translate(::LazySet, ::AbstractVector))
diff --git a/src/LazySets.jl b/src/LazySets.jl
index 11f7479b95..cecc661914 100644
--- a/src/LazySets.jl
+++ b/src/LazySets.jl
@@ -22,7 +22,7 @@ import LinearAlgebra: ×, normalize, normalize!
 import RecipesBase: apply_recipe
 
 export Arrays
-export ×, normalize, subtypes
+export ×, normalize, normalize!, subtypes
 
 using LinearAlgebra, RecipesBase, Requires, SparseArrays
 import GLPK, JuMP, Random, ReachabilityBase
@@ -123,7 +123,10 @@ include("Sets/HPolyhedron.jl")
 include("Sets/Hyperplane.jl")
 include("Sets/Hyperrectangle.jl")
 include("Sets/Line2D.jl")
-include("Sets/Line.jl")
+
+include("Sets/Line/LineModule.jl")
+@reexport using ..LineModule: Line, direction
+
 include("Sets/Polygon.jl")
 include("Sets/RotatedHyperrectangle.jl")
 
diff --git a/src/Sets/Line.jl b/src/Sets/Line.jl
deleted file mode 100644
index b674520313..0000000000
--- a/src/Sets/Line.jl
+++ /dev/null
@@ -1,493 +0,0 @@
-export Line
-
-"""
-    Line{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
-
-Type that represents a line of the form
-
-```math
-    \\{y ∈ ℝ^n: y = p + λd, λ ∈ ℝ\\}
-```
-where ``p`` is a point on the line and ``d`` is its direction vector (not
-necessarily normalized).
-
-### Fields
-
-- `p` -- point on the line
-- `d` -- direction
-
-### Examples
-
-There are three constructors. The optional keyword argument `normalize`
-(default: `false`) can be used to normalize the direction of the resulting line
-to have norm 1 (w.r.t. the Euclidean norm).
-
-1) The default constructor takes the fields `p` and `d`:
-
-The line passing through the point ``[-1, 2, 3]`` and parallel to the vector
-``[3, 0, -1]``:
-
-```jldoctest
-julia> Line([-1.0, 2, 3], [3.0, 0, -1])
-Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [3.0, 0.0, -1.0])
-
-julia> Line([-1.0, 2, 3], [3.0, 0, -1]; normalize=true)
-Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [0.9486832980505138, 0.0, -0.31622776601683794])
-```
-
-2) The second constructor takes two points, `from` and `to`, as keyword
-arguments, and returns the line through them. See the algorithm section for
-details.
-
-```jldoctest
-julia> Line(from=[-1.0, 2, 3], to=[-4.0, 2, 4])
-Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [3.0, 0.0, -1.0])
-```
-
-3) The third constructor resembles `Line2D` and only works for two-dimensional
-lines. It takes two inputs, `a` and `b`, and constructs the line such that
-``a ⋅ x = b``.
-
-```jldoctest
-julia> Line([2.0, 0], 1.)
-Line{Float64, Vector{Float64}}([0.5, 0.0], [0.0, -1.0])
-```
-
-### Algorithm
-
-Given two points ``p ∈ ℝ^n`` and ``q ∈ ℝ^n``, the line that
-passes through these two points is
-`L: `\\{y ∈ ℝ^n: y = p + λ(q - p), λ ∈ ℝ\\}``.
-"""
-struct Line{N,VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
-    p::VN
-    d::VN
-
-    # default constructor with length constraint
-    function Line(p::VN, d::VN; check_direction::Bool=true,
-                  normalize=false) where {N,VN<:AbstractVector{N}}
-        if check_direction
-            @assert !iszero(d) "a line needs a non-zero direction vector"
-        end
-        d_n = normalize ? LinearAlgebra.normalize(d) : d
-        return new{N,VN}(p, d_n)
-    end
-end
-
-function Line(; from::AbstractVector, to::AbstractVector, normalize=false)
-    d = from - to
-    @assert !iszero(d) "points `$from` and `$to` should be distinct"
-    return Line(from, d; normalize=normalize)
-end
-
-function Line(a::AbstractVector{N}, b::N; normalize=false) where {N}
-    @assert length(a) == 2 "expected a normal vector of length two, but it " *
-                           "is $(length(a))-dimensional"
-
-    got_horizontal = iszero(a[1])
-    got_vertical = iszero(a[2])
-
-    if got_horizontal && got_vertical
-        throw(ArgumentError("the vector $a must be non-zero"))
-    end
-
-    if got_horizontal
-        α = b / a[2]
-        p = [zero(N), α]
-        q = [one(N), α]
-    elseif got_vertical
-        β = b / a[1]
-        p = [β, zero(N)]
-        q = [β, one(N)]
-    else
-        α = b / a[2]
-        μ = a[1] / a[2]
-        p = [zero(N), α]
-        q = [one(N), α - μ]
-    end
-    return Line(; from=p, to=q, normalize=normalize)
-end
-
-isoperationtype(::Type{<:Line}) = false
-
-"""
-    direction(L::Line)
-
-Return the direction of the line.
-
-### Input
-
-- `L` -- line
-
-### Output
-
-The direction of the line.
-
-### Notes
-
-The direction is not necessarily normalized.
-See [`normalize(::Line, ::Real)`](@ref) / [`normalize!(::Line, ::Real)`](@ref).
-"""
-direction(L::Line) = L.d
-
-"""
-    normalize!(L::Line{N}, p::Real=N(2)) where {N}
-
-Normalize the direction of a line storing the result in `L`.
-
-### Input
-
-- `L` -- line
-- `p` -- (optional, default: `2.0`) vector `p`-norm used in the normalization
-
-### Output
-
-A line whose direction has unit norm w.r.t. the given `p`-norm.
-"""
-function normalize!(L::Line{N}, p::Real=N(2)) where {N}
-    normalize!(L.d, p)
-    return L
-end
-
-"""
-    normalize(L::Line{N}, p::Real=N(2)) where {N}
-
-Normalize the direction of a line.
-
-### Input
-
-- `L` -- line
-- `p` -- (optional, default: `2.0`) vector `p`-norm used in the normalization
-
-### Output
-
-A line whose direction has unit norm w.r.t. the given `p`-norm.
-
-### Notes
-
-See also [`normalize!(::Line, ::Real)`](@ref) for the in-place version.
-"""
-function normalize(L::Line{N}, p::Real=N(2)) where {N}
-    return normalize!(copy(L), p)
-end
-
-"""
-    constraints_list(L::Line)
-
-Return the list of constraints of a line.
-
-### Input
-
-- `L` -- line
-
-### Output
-
-A list containing `2n-2` half-spaces whose intersection is `L`, where `n` is the
-ambient dimension of `L`.
-"""
-function constraints_list(L::Line)
-    p = L.p
-    n = length(p)
-    d = reshape(L.d, 1, n)
-    K = nullspace(d)
-    m = size(K, 2)
-    @assert m == n - 1 "expected $(n - 1) normal half-spaces, got $m"
-
-    N, VN = _parameters(L)
-    out = Vector{HalfSpace{N,VN}}(undef, 2m)
-    idx = 1
-    @inbounds for j in 1:m
-        Kj = K[:, j]
-        b = dot(Kj, p)
-        out[idx] = HalfSpace(Kj, b)
-        out[idx + 1] = HalfSpace(-Kj, -b)
-        idx += 2
-    end
-    return out
-end
-
-function _parameters(::Line{N,VN}) where {N,VN}
-    return (N, VN)
-end
-
-"""
-    dim(L::Line)
-
-Return the ambient dimension of a line.
-
-### Input
-
-- `L` -- line
-
-### Output
-
-The ambient dimension of the line.
-"""
-dim(L::Line) = length(L.p)
-
-"""
-    ρ(d::AbstractVector, L::Line)
-
-Evaluate the support function of a line in a given direction.
-
-### Input
-
-- `d` -- direction
-- `L` -- line
-
-### Output
-
-Evaluation of the support function in the given direction.
-"""
-function ρ(d::AbstractVector, L::Line)
-    if isapproxzero(dot(d, L.d))
-        return dot(d, L.p)
-    else
-        N = eltype(L)
-        return N(Inf)
-    end
-end
-
-"""
-    σ(d::AbstractVector, L::Line)
-
-Return a support vector of a line in a given direction.
-
-### Input
-
-- `d` -- direction
-- `L` -- line
-
-### Output
-
-A support vector in the given direction.
-"""
-function σ(d::AbstractVector, L::Line)
-    if isapproxzero(dot(d, L.d))
-        return L.p
-    else
-        throw(ArgumentError("the support vector is undefined because the " *
-                            "line is unbounded in the given direction"))
-    end
-end
-
-"""
-    isbounded(L::Line)
-
-Determine whether a line is bounded.
-
-### Input
-
-- `L` -- line
-
-### Output
-
-`false`.
-"""
-isbounded(::Line) = false
-
-"""
-    isuniversal(L::Line; [witness::Bool]=false)
-
-Check whether a line is universal.
-
-### Input
-
-- `P`       -- line
-- `witness` -- (optional, default: `false`) compute a witness if activated
-
-### Output
-
-* If `witness` is `false`: `true` if the ambient dimension is one, `false`
-otherwise.
-
-* If `witness` is `true`: `(true, [])` if the ambient dimension is one,
-`(false, v)` where ``v ∉ P`` otherwise.
-"""
-isuniversal(L::Line; witness::Bool=false) = isuniversal(L, Val(witness))
-
-# TODO implement case with witness
-isuniversal(L::Line, ::Val{false}) = dim(L) == 1
-
-"""
-    an_element(L::Line)
-
-Return some element of a line.
-
-### Input
-
-- `L` -- line
-
-### Output
-
-An element on the line.
-"""
-an_element(L::Line) = L.p
-
-"""
-    ∈(x::AbstractVector, L::Line)
-
-Check whether a given point is contained in a line.
-
-### Input
-
-- `x` -- point/vector
-- `L` -- line
-
-### Output
-
-`true` iff `x ∈ L`.
-
-### Algorithm
-
-The point ``x`` belongs to the line ``L : p + λd`` if and only if ``x - p`` is
-proportional to the direction ``d``.
-"""
-function ∈(x::AbstractVector, L::Line)
-    @assert length(x) == dim(L) "expected the point and the line to have the " *
-                                "same dimension, but they are $(length(x)) and $(dim(L)) respectively"
-
-    # the following check is necessary because the zero vector is a special case
-    _isapprox(x, L.p) && return true
-
-    return first(ismultiple(x - L.p, L.d))
-end
-
-"""
-    rand(::Type{Line}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
-         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)
-
-Create a random line.
-
-### Input
-
-- `Line` -- type for dispatch
-- `N`    -- (optional, default: `Float64`) numeric type
-- `dim`  -- (optional, default: 2) dimension
-- `rng`  -- (optional, default: `GLOBAL_RNG`) random number generator
-- `seed` -- (optional, default: `nothing`) seed for reseeding
-
-### Output
-
-A random line.
-
-### Algorithm
-
-All numbers are normally distributed with mean 0 and standard deviation 1.
-"""
-function rand(::Type{Line};
-              N::Type{<:Real}=Float64,
-              dim::Int=2,
-              rng::AbstractRNG=GLOBAL_RNG,
-              seed::Union{Int,Nothing}=nothing)
-    rng = reseed!(rng, seed)
-    d = randn(rng, N, dim)
-    while iszero(d)
-        d = randn(rng, N, dim)
-    end
-    p = randn(rng, N, dim)
-    return Line(p, d)
-end
-
-"""
-    isempty(L::Line)
-
-Check whether a line is empty.
-
-### Input
-
-- `L` -- line
-
-### Output
-
-`false`.
-"""
-isempty(::Line) = false
-
-"""
-    translate!(L::Line, v::AbstractVector)
-
-Translate (i.e., shift) a line by a given vector, in-place.
-
-### Input
-
-- `L` -- line
-- `v` -- translation vector
-
-### Output
-
-The line `L` translated by `v`.
-"""
-function translate!(L::Line, v::AbstractVector)
-    @assert length(v) == dim(L) "cannot translate a $(dim(L))-dimensional " *
-                                "set by a $(length(v))-dimensional vector"
-
-    L.p .+= v
-    return L
-end
-
-"""
-    distance(x::AbstractVector, L::Line; [p]::Real=2.0)
-
-Compute the distance between point `x` and the line with respect to the given
-`p`-norm.
-
-### Input
-
-- `x` -- point/vector
-- `L` -- line
-- `p` -- (optional, default: `2.0`) the `p`-norm used; `p = 2.0` corresponds to
-         the usual Euclidean norm
-
-### Output
-
-A scalar representing the distance between `x` and the line `L`.
-"""
-@commutative function distance(x::AbstractVector, L::Line; p::Real=2.0)
-    d = L.d  # direction of the line
-    t = dot(x - L.p, d) / dot(d, d)
-    return distance(x, L.p + t * d; p=p)
-end
-
-"""
-    linear_map(M::AbstractMatrix, L::Line)
-
-Concrete linear map of a line.
-
-### Input
-
-- `M` -- matrix
-- `L` -- line
-
-### Output
-
-The line obtained by applying the linear map, if that still results in a line,
-or a `Singleton` otherwise.
-
-### Algorithm
-
-We apply the linear map to the point and direction of `L`.
-If the resulting direction is zero, the result is a singleton.
-"""
-function linear_map(M::AbstractMatrix, L::Line)
-    @assert dim(L) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
-                                 "applied to a set of dimension $(dim(L))"
-
-    Mp = M * L.p
-    Md = M * L.d
-    if iszero(Md)
-        return Singleton(Mp)
-    end
-    return Line(Mp, Md)
-end
-
-function project(L::Line{N}, block::AbstractVector{Int}; kwargs...) where {N}
-    d = L.d[block]
-    if iszero(d)
-        return Singleton(L.p[block])  # projected out all nontrivial dimensions
-    elseif length(d) == 1
-        return Universe{N}(1)  # special case: 1D line is a universe
-    else
-        return Line(L.p[block], d)
-    end
-end
diff --git a/src/Sets/Line/Line.jl b/src/Sets/Line/Line.jl
new file mode 100644
index 0000000000..f3c65b3c0c
--- /dev/null
+++ b/src/Sets/Line/Line.jl
@@ -0,0 +1,107 @@
+"""
+    Line{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
+
+Type that represents a line of the form
+
+```math
+    \\{y ∈ ℝ^n: y = p + λd, λ ∈ ℝ\\}
+```
+where ``p`` is a point on the line and ``d`` is its direction vector (not
+necessarily normalized).
+
+### Fields
+
+- `p` -- point on the line
+- `d` -- direction
+
+### Examples
+
+There are three constructors. The optional keyword argument `normalize`
+(default: `false`) can be used to normalize the direction of the resulting line
+to have norm 1 (w.r.t. the Euclidean norm).
+
+1) The default constructor takes the fields `p` and `d`:
+
+The line passing through the point ``[-1, 2, 3]`` and parallel to the vector
+``[3, 0, -1]``:
+
+```jldoctest
+julia> Line([-1.0, 2, 3], [3.0, 0, -1])
+Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [3.0, 0.0, -1.0])
+
+julia> Line([-1.0, 2, 3], [3.0, 0, -1]; normalize=true)
+Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [0.9486832980505138, 0.0, -0.31622776601683794])
+```
+
+2) The second constructor takes two points, `from` and `to`, as keyword
+arguments, and returns the line through them. See the algorithm section for
+details.
+
+```jldoctest
+julia> Line(from=[-1.0, 2, 3], to=[-4.0, 2, 4])
+Line{Float64, Vector{Float64}}([-1.0, 2.0, 3.0], [3.0, 0.0, -1.0])
+```
+
+3) The third constructor resembles `Line2D` and only works for two-dimensional
+lines. It takes two inputs, `a` and `b`, and constructs the line such that
+``a ⋅ x = b``.
+
+```jldoctest
+julia> Line([2.0, 0], 1.)
+Line{Float64, Vector{Float64}}([0.5, 0.0], [0.0, -1.0])
+```
+
+### Algorithm
+
+Given two points ``p ∈ ℝ^n`` and ``q ∈ ℝ^n``, the line that
+passes through these two points is
+`L: `\\{y ∈ ℝ^n: y = p + λ(q - p), λ ∈ ℝ\\}``.
+"""
+struct Line{N,VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
+    p::VN
+    d::VN
+
+    # default constructor with length constraint
+    function Line(p::VN, d::VN; check_direction::Bool=true,
+                  normalize=false) where {N,VN<:AbstractVector{N}}
+        if check_direction
+            @assert !iszero(d) "a line needs a non-zero direction vector"
+        end
+        d_n = normalize ? LinearAlgebra.normalize(d) : d
+        return new{N,VN}(p, d_n)
+    end
+end
+
+function Line(; from::AbstractVector, to::AbstractVector, normalize=false)
+    d = from - to
+    @assert !iszero(d) "points `$from` and `$to` should be distinct"
+    return Line(from, d; normalize=normalize)
+end
+
+function Line(a::AbstractVector{N}, b::N; normalize=false) where {N}
+    @assert length(a) == 2 "expected a normal vector of length two, but it " *
+                           "is $(length(a))-dimensional"
+
+    got_horizontal = iszero(a[1])
+    got_vertical = iszero(a[2])
+
+    if got_horizontal && got_vertical
+        throw(ArgumentError("the vector $a must be non-zero"))
+    end
+
+    if got_horizontal
+        α = b / a[2]
+        p = [zero(N), α]
+        q = [one(N), α]
+    elseif got_vertical
+        β = b / a[1]
+        p = [β, zero(N)]
+        q = [β, one(N)]
+    else
+        α = b / a[2]
+        μ = a[1] / a[2]
+        p = [zero(N), α]
+        q = [one(N), α - μ]
+    end
+    return Line(; from=p, to=q, normalize=normalize)
+end
diff --git a/src/Sets/Line/LineModule.jl b/src/Sets/Line/LineModule.jl
new file mode 100644
index 0000000000..bf161e59b6
--- /dev/null
+++ b/src/Sets/Line/LineModule.jl
@@ -0,0 +1,47 @@
+module LineModule
+
+using Reexport, Requires
+
+using ..LazySets: AbstractPolyhedron, HalfSpace, @commutative
+using Random: AbstractRNG, GLOBAL_RNG
+using ReachabilityBase.Arrays: ismultiple
+using ReachabilityBase.Distribution: reseed!
+using ReachabilityBase.Comparison: _isapprox, isapproxzero
+using ReachabilityBase.Require
+import LinearAlgebra
+using LinearAlgebra: dot, nullspace
+
+@reexport import ..API: an_element, constraints_list, dim, isbounded, isempty,
+                        isoperationtype, isuniversal, project, rand, distance,
+                        ∈, linear_map, ρ, σ, translate!
+@reexport import ..LazySets: normalize
+@reexport import ..LinearAlgebra: normalize!
+@reexport using ..API
+
+export Line,
+       direction
+
+include("Line.jl")
+
+include("an_element.jl")
+include("constraints_list.jl")
+include("dim.jl")
+include("isbounded.jl")
+include("isempty.jl")
+include("isoperationtype.jl")
+include("isuniversal.jl")
+include("normalize.jl")
+include("project.jl")
+include("rand.jl")
+include("distance.jl")
+include("in.jl")
+include("linear_map.jl")
+include("support_function.jl")
+include("support_vector.jl")
+include("translate.jl")
+
+include("direction.jl")
+
+include("init.jl")
+
+end  # module
diff --git a/src/Sets/Line/an_element.jl b/src/Sets/Line/an_element.jl
new file mode 100644
index 0000000000..4cc55c93db
--- /dev/null
+++ b/src/Sets/Line/an_element.jl
@@ -0,0 +1,14 @@
+"""
+    an_element(L::Line)
+
+Return some element of a line.
+
+### Input
+
+- `L` -- line
+
+### Output
+
+An element on the line.
+"""
+an_element(L::Line) = L.p
diff --git a/src/Sets/Line/constraints_list.jl b/src/Sets/Line/constraints_list.jl
new file mode 100644
index 0000000000..0691fd6d5a
--- /dev/null
+++ b/src/Sets/Line/constraints_list.jl
@@ -0,0 +1,38 @@
+"""
+    constraints_list(L::Line)
+
+Return the list of constraints of a line.
+
+### Input
+
+- `L` -- line
+
+### Output
+
+A list containing `2n-2` half-spaces whose intersection is `L`, where `n` is the
+ambient dimension of `L`.
+"""
+function constraints_list(L::Line)
+    p = L.p
+    n = length(p)
+    d = reshape(L.d, 1, n)
+    K = nullspace(d)
+    m = size(K, 2)
+    @assert m == n - 1 "expected $(n - 1) normal half-spaces, got $m"
+
+    N, VN = _parameters(L)
+    out = Vector{HalfSpace{N,VN}}(undef, 2m)
+    idx = 1
+    @inbounds for j in 1:m
+        Kj = K[:, j]
+        b = dot(Kj, p)
+        out[idx] = HalfSpace(Kj, b)
+        out[idx + 1] = HalfSpace(-Kj, -b)
+        idx += 2
+    end
+    return out
+end
+
+function _parameters(::Line{N,VN}) where {N,VN}
+    return (N, VN)
+end
diff --git a/src/Sets/Line/dim.jl b/src/Sets/Line/dim.jl
new file mode 100644
index 0000000000..4c7841ea74
--- /dev/null
+++ b/src/Sets/Line/dim.jl
@@ -0,0 +1,14 @@
+"""
+    dim(L::Line)
+
+Return the ambient dimension of a line.
+
+### Input
+
+- `L` -- line
+
+### Output
+
+The ambient dimension of the line.
+"""
+dim(L::Line) = length(L.p)
diff --git a/src/Sets/Line/direction.jl b/src/Sets/Line/direction.jl
new file mode 100644
index 0000000000..0dc522e8ab
--- /dev/null
+++ b/src/Sets/Line/direction.jl
@@ -0,0 +1,19 @@
+"""
+    direction(L::Line)
+
+Return the direction of the line.
+
+### Input
+
+- `L` -- line
+
+### Output
+
+The direction of the line.
+
+### Notes
+
+The direction is not necessarily normalized.
+See [`normalize(::Line, ::Real)`](@ref) / [`normalize!(::Line, ::Real)`](@ref).
+"""
+direction(L::Line) = L.d
diff --git a/src/Sets/Line/distance.jl b/src/Sets/Line/distance.jl
new file mode 100644
index 0000000000..480a5e8384
--- /dev/null
+++ b/src/Sets/Line/distance.jl
@@ -0,0 +1,22 @@
+"""
+    distance(x::AbstractVector, L::Line; [p]::Real=2.0)
+
+Compute the distance between point `x` and the line with respect to the given
+`p`-norm.
+
+### Input
+
+- `x` -- point/vector
+- `L` -- line
+- `p` -- (optional, default: `2.0`) the `p`-norm used; `p = 2.0` corresponds to
+         the usual Euclidean norm
+
+### Output
+
+A scalar representing the distance between `x` and the line `L`.
+"""
+@commutative function distance(x::AbstractVector, L::Line; p::Real=2.0)
+    d = L.d  # direction of the line
+    t = dot(x - L.p, d) / dot(d, d)
+    return distance(x, L.p + t * d; p=p)
+end
diff --git a/src/Sets/Line/in.jl b/src/Sets/Line/in.jl
new file mode 100644
index 0000000000..9e1e37d6c8
--- /dev/null
+++ b/src/Sets/Line/in.jl
@@ -0,0 +1,28 @@
+"""
+    ∈(x::AbstractVector, L::Line)
+
+Check whether a given point is contained in a line.
+
+### Input
+
+- `x` -- point/vector
+- `L` -- line
+
+### Output
+
+`true` iff `x ∈ L`.
+
+### Algorithm
+
+The point ``x`` belongs to the line ``L : p + λd`` if and only if ``x - p`` is
+proportional to the direction ``d``.
+"""
+function ∈(x::AbstractVector, L::Line)
+    @assert length(x) == dim(L) "expected the point and the line to have the " *
+                                "same dimension, but they are $(length(x)) and $(dim(L)) respectively"
+
+    # the following check is necessary because the zero vector is a special case
+    _isapprox(x, L.p) && return true
+
+    return first(ismultiple(x - L.p, L.d))
+end
diff --git a/src/Sets/Line/init.jl b/src/Sets/Line/init.jl
new file mode 100644
index 0000000000..ddd845644b
--- /dev/null
+++ b/src/Sets/Line/init.jl
@@ -0,0 +1,3 @@
+function __init__()
+    @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("init_LazySets.jl")
+end
diff --git a/src/Sets/Line/init_LazySets.jl b/src/Sets/Line/init_LazySets.jl
new file mode 100644
index 0000000000..6aad752a32
--- /dev/null
+++ b/src/Sets/Line/init_LazySets.jl
@@ -0,0 +1,2 @@
+using .LazySets.SingletonModule: Singleton
+using .LazySets.UniverseModule: Universe
diff --git a/src/Sets/Line/isbounded.jl b/src/Sets/Line/isbounded.jl
new file mode 100644
index 0000000000..599f52bd7b
--- /dev/null
+++ b/src/Sets/Line/isbounded.jl
@@ -0,0 +1,14 @@
+"""
+    isbounded(L::Line)
+
+Determine whether a line is bounded.
+
+### Input
+
+- `L` -- line
+
+### Output
+
+`false`.
+"""
+isbounded(::Line) = false
diff --git a/src/Sets/Line/isempty.jl b/src/Sets/Line/isempty.jl
new file mode 100644
index 0000000000..1d591fc8ed
--- /dev/null
+++ b/src/Sets/Line/isempty.jl
@@ -0,0 +1,14 @@
+"""
+    isempty(L::Line)
+
+Check whether a line is empty.
+
+### Input
+
+- `L` -- line
+
+### Output
+
+`false`.
+"""
+isempty(::Line) = false
diff --git a/src/Sets/Line/isoperationtype.jl b/src/Sets/Line/isoperationtype.jl
new file mode 100644
index 0000000000..396b8c5e1b
--- /dev/null
+++ b/src/Sets/Line/isoperationtype.jl
@@ -0,0 +1 @@
+isoperationtype(::Type{<:Line}) = false
diff --git a/src/Sets/Line/isuniversal.jl b/src/Sets/Line/isuniversal.jl
new file mode 100644
index 0000000000..39fc646752
--- /dev/null
+++ b/src/Sets/Line/isuniversal.jl
@@ -0,0 +1,22 @@
+"""
+    isuniversal(L::Line; [witness::Bool]=false)
+
+Check whether a line is universal.
+
+### Input
+
+- `P`       -- line
+- `witness` -- (optional, default: `false`) compute a witness if activated
+
+### Output
+
+* If `witness` is `false`: `true` if the ambient dimension is one, `false`
+otherwise.
+
+* If `witness` is `true`: `(true, [])` if the ambient dimension is one,
+`(false, v)` where ``v ∉ P`` otherwise.
+"""
+isuniversal(L::Line; witness::Bool=false) = isuniversal(L, Val(witness))
+
+# TODO implement case with witness
+isuniversal(L::Line, ::Val{false}) = dim(L) == 1
diff --git a/src/Sets/Line/linear_map.jl b/src/Sets/Line/linear_map.jl
new file mode 100644
index 0000000000..a91c0c3695
--- /dev/null
+++ b/src/Sets/Line/linear_map.jl
@@ -0,0 +1,33 @@
+"""
+    linear_map(M::AbstractMatrix, L::Line)
+
+Concrete linear map of a line.
+
+### Input
+
+- `M` -- matrix
+- `L` -- line
+
+### Output
+
+The line obtained by applying the linear map, if that still results in a line,
+or a `Singleton` otherwise.
+
+### Algorithm
+
+We apply the linear map to the point and direction of `L`.
+If the resulting direction is zero, the result is a singleton.
+"""
+function linear_map(M::AbstractMatrix, L::Line)
+    @assert dim(L) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
+                                 "applied to a set of dimension $(dim(L))"
+
+    Mp = M * L.p
+    Md = M * L.d
+    if iszero(Md)
+        require(@__MODULE__, :LazySets; fun_name="linear_map")
+
+        return Singleton(Mp)
+    end
+    return Line(Mp, Md)
+end
diff --git a/src/Sets/Line/normalize.jl b/src/Sets/Line/normalize.jl
new file mode 100644
index 0000000000..5ec8ac8e18
--- /dev/null
+++ b/src/Sets/Line/normalize.jl
@@ -0,0 +1,40 @@
+"""
+    normalize(L::Line{N}, p::Real=N(2)) where {N}
+
+Normalize the direction of a line.
+
+### Input
+
+- `L` -- line
+- `p` -- (optional, default: `2.0`) vector `p`-norm used in the normalization
+
+### Output
+
+A line whose direction has unit norm w.r.t. the given `p`-norm.
+
+### Notes
+
+See also [`normalize!(::Line, ::Real)`](@ref) for the in-place version.
+"""
+function normalize(L::Line{N}, p::Real=N(2)) where {N}
+    return normalize!(copy(L), p)
+end
+
+"""
+    normalize!(L::Line{N}, p::Real=N(2)) where {N}
+
+Normalize the direction of a line storing the result in `L`.
+
+### Input
+
+- `L` -- line
+- `p` -- (optional, default: `2.0`) vector `p`-norm used in the normalization
+
+### Output
+
+A line whose direction has unit norm w.r.t. the given `p`-norm.
+"""
+function normalize!(L::Line{N}, p::Real=N(2)) where {N}
+    normalize!(L.d, p)
+    return L
+end
diff --git a/src/Sets/Line/project.jl b/src/Sets/Line/project.jl
new file mode 100644
index 0000000000..1b18d48127
--- /dev/null
+++ b/src/Sets/Line/project.jl
@@ -0,0 +1,14 @@
+function project(L::Line{N}, block::AbstractVector{Int}; kwargs...) where {N}
+    d = L.d[block]
+    if iszero(d)
+        require(@__MODULE__, :LazySets; fun_name="project")
+
+        return Singleton(L.p[block])  # projected out all nontrivial dimensions
+    elseif length(d) == 1
+        require(@__MODULE__, :LazySets; fun_name="project")
+
+        return Universe{N}(1)  # special case: 1D line is a universe
+    else
+        return Line(L.p[block], d)
+    end
+end
diff --git a/src/Sets/Line/rand.jl b/src/Sets/Line/rand.jl
new file mode 100644
index 0000000000..a660b27ae8
--- /dev/null
+++ b/src/Sets/Line/rand.jl
@@ -0,0 +1,35 @@
+"""
+    rand(::Type{Line}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
+         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)
+
+Create a random line.
+
+### Input
+
+- `Line` -- type for dispatch
+- `N`    -- (optional, default: `Float64`) numeric type
+- `dim`  -- (optional, default: 2) dimension
+- `rng`  -- (optional, default: `GLOBAL_RNG`) random number generator
+- `seed` -- (optional, default: `nothing`) seed for reseeding
+
+### Output
+
+A random line.
+
+### Algorithm
+
+All numbers are normally distributed with mean 0 and standard deviation 1.
+"""
+function rand(::Type{Line};
+              N::Type{<:Real}=Float64,
+              dim::Int=2,
+              rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int,Nothing}=nothing)
+    rng = reseed!(rng, seed)
+    d = randn(rng, N, dim)
+    while iszero(d)
+        d = randn(rng, N, dim)  # COV_EXCL_LINE
+    end  # COV_EXCL_LINE
+    p = randn(rng, N, dim)
+    return Line(p, d)
+end
diff --git a/src/Sets/Line/support_function.jl b/src/Sets/Line/support_function.jl
new file mode 100644
index 0000000000..c67c57d00d
--- /dev/null
+++ b/src/Sets/Line/support_function.jl
@@ -0,0 +1,22 @@
+"""
+    ρ(d::AbstractVector, L::Line)
+
+Evaluate the support function of a line in a given direction.
+
+### Input
+
+- `d` -- direction
+- `L` -- line
+
+### Output
+
+Evaluation of the support function in the given direction.
+"""
+function ρ(d::AbstractVector, L::Line)
+    if isapproxzero(dot(d, L.d))
+        return dot(d, L.p)
+    else
+        N = eltype(L)
+        return N(Inf)
+    end
+end
diff --git a/src/Sets/Line/support_vector.jl b/src/Sets/Line/support_vector.jl
new file mode 100644
index 0000000000..ff6f4ebc20
--- /dev/null
+++ b/src/Sets/Line/support_vector.jl
@@ -0,0 +1,22 @@
+"""
+    σ(d::AbstractVector, L::Line)
+
+Return a support vector of a line in a given direction.
+
+### Input
+
+- `d` -- direction
+- `L` -- line
+
+### Output
+
+A support vector in the given direction.
+"""
+function σ(d::AbstractVector, L::Line)
+    if isapproxzero(dot(d, L.d))
+        return L.p
+    else
+        throw(ArgumentError("the support vector is undefined because the " *
+                            "line is unbounded in the given direction"))
+    end
+end
diff --git a/src/Sets/Line/translate.jl b/src/Sets/Line/translate.jl
new file mode 100644
index 0000000000..5a36b2a988
--- /dev/null
+++ b/src/Sets/Line/translate.jl
@@ -0,0 +1,21 @@
+"""
+    translate!(L::Line, v::AbstractVector)
+
+Translate (i.e., shift) a line by a given vector, in-place.
+
+### Input
+
+- `L` -- line
+- `v` -- translation vector
+
+### Output
+
+The line `L` translated by `v`.
+"""
+function translate!(L::Line, v::AbstractVector)
+    @assert length(v) == dim(L) "cannot translate a $(dim(L))-dimensional " *
+                                "set by a $(length(v))-dimensional vector"
+
+    L.p .+= v
+    return L
+end
diff --git a/test/Sets/Line.jl b/test/Sets/Line.jl
index 7055d7fa00..42cce2c8e1 100644
--- a/test/Sets/Line.jl
+++ b/test/Sets/Line.jl
@@ -15,13 +15,26 @@ for N in [Float64, Rational{Int}, Float32]
     @test N[1, 0] ∈ ll && N[1, 1] ∈ ll
     ll = Line(N[1, -1], N(0)) # x == y
     @test N[0, 0] ∈ ll && N[1, 1] ∈ ll
+    @test_throws ArgumentError Line(zeros(N, 2), N(0))
 
     # the lines are the same modulo the sign of the normal vector
     @test l1.p ≈ l2.p && l1.d ≈ -l2.d
 
+    # direction
+    @test direction(l2) == N[1, 0]
+
+    # normalize
+    l3 = Line(N[0, 1], N[2, 0])
+    @test normalize(l3) == Line(N[0, 1], N[1, 0])
+    normalize!(l3)
+    @test l3 == Line(N[0, 1], N[1, 0])
+
     # dimension
     @test dim(l1) == 2
 
+    # isoperationtype
+    @test !isoperationtype(Line)
+
     # support function
     @test ρ(N[0, 1], l1) == N(1)
     @test ρ(N[1, 0], l1) == N(Inf)