diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl
index 5fd08339e4..f32b1ac9b6 100644
--- a/src/HPolyhedron.jl
+++ b/src/HPolyhedron.jl
@@ -499,6 +499,9 @@ The `HPolyhedron` (resp. `HPolytope`) obtained by the concrete convex hull of
 
 ### Notes
 
+For performance reasons, it is suggested to use the `CDDLib.Library()` backend
+for the `convex_hull`.
+
 For further information on the supported backends see
 [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/).
 """
@@ -506,8 +509,9 @@ function convex_hull(P1::HPoly{N},
                      P2::HPoly{N};
                      backend=default_polyhedra_backend(P1, N)) where {N}
     @assert isdefined(@__MODULE__, :Polyhedra) "the function `convex_hull` needs " *
-                                        "the package 'Polyhedra' to be loaded"
+                                               "the package 'Polyhedra' to be loaded"
     Pch = convexhull(polyhedron(P1; backend=backend), polyhedron(P2; backend=backend))
+    removehredundancy!(Pch)
     return convert(typeof(P1), Pch)
 end
 
diff --git a/src/VPolytope.jl b/src/VPolytope.jl
index 52f2c2d47d..4b3b08a8a0 100644
--- a/src/VPolytope.jl
+++ b/src/VPolytope.jl
@@ -3,7 +3,8 @@ using MathProgBase, GLPKMathProgInterface
 import Base.rand
 
 export VPolytope,
-       vertices_list
+       vertices_list,
+       convex_hull
 
 """
     VPolytope{N<:Real} <: AbstractPolytope{N}
@@ -187,16 +188,48 @@ function constraints_list(P::VPolytope{N})::Vector{LinearConstraint{N}} where {N
     return constraints_list(tohrep(P))
 end
 
+"""
+    convex_hull(P1::VPolytope{N}, P2::VPolytope{N};
+                [backend]=default_polyhedra_backend(P1, N)) where {N}
+
+Compute the convex hull of the set union of two polytopes in V-representation.
+
+### Input
+
+- `P1`         -- polytope
+- `P2`         -- another polytope
+- `backend`    -- (optional, default: `default_polyhedra_backend(P1, N)`) the polyhedral
+                  computations backend, see [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1)
+                  for further information
+
+### Output
 
-# --- functions that use Polyhedra.jl ---
+The `VPolytope` obtained by the concrete convex hull of `P1` and `P2`.
+
+### Notes
+
+For performance reasons, it is suggested to use the `CDDLib.Library()` backend
+for the `convex_hull`.
+"""
+function convex_hull(P1::VPolytope{N}, P2::VPolytope{N};
+                     backend=default_polyhedra_backend(P1, N)) where {N}
+    @assert isdefined(@__MODULE__, :Polyhedra) "the function `convex_hull` needs " *
+                                               "the package 'Polyhedra' to be loaded"
+    Pch = convexhull(polyhedron(P1; backend=backend),
+                     polyhedron(P2; backend=backend))
+    removevredundancy!(Pch)
+    return VPolytope(Pch)
+end
 
+# ==========================================
+# Lower level methods that use Polyhedra.jl
+# ==========================================
 
 function load_polyhedra_vpolytope() # function to be loaded by Requires
 return quote
 # see the interface file AbstractPolytope.jl for the imports
 
-export convex_hull,
-       cartesian_product,
+export cartesian_product,
        vertices_list,
        tohrep,
        tovrep
@@ -247,31 +280,6 @@ function polyhedron(P::VPolytope{N};
     return polyhedron(Polyhedra.vrep(V), backend)
 end
 
-"""
-    convex_hull(P1::VPolytope{N}, P2::VPolytope{N};
-                [backend]=default_polyhedra_backend(P1, N)) where {N}
-
-Compute the convex hull of the set union of two polytopes in V-representation.
-
-### Input
-
-- `P1`         -- polytope
-- `P2`         -- another polytope
-- `backend`    -- (optional, default: `default_polyhedra_backend(P1, N)`) the polyhedral
-                  computations backend, see [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1)
-                  for further information
-
-### Output
-
-The `VPolytope` obtained by the concrete convex hull of `P1` and `P2`.
-"""
-function convex_hull(P1::VPolytope{N}, P2::VPolytope{N};
-                     backend=default_polyhedra_backend(P1, N)) where {N}
-    Pch = convexhull(polyhedron(P1; backend=backend),
-                     polyhedron(P2; backend=backend))
-    return VPolytope(Pch)
-end
-
 """
     cartesian_product(P1::VPolytope{N}, P2::VPolytope{N};
                       [backend]=default_polyhedra_backend(P1, N)) where {N}
diff --git a/test/unit_Polytope.jl b/test/unit_Polytope.jl
index 3beba210a7..9b83d94be6 100644
--- a/test/unit_Polytope.jl
+++ b/test/unit_Polytope.jl
@@ -229,11 +229,7 @@ if test_suite_polyhedra
         p2 = VPolytope([v3, v4])
         ch = convex_hull(p1, p2)
         vl = vertices_list(ch)
-        @test v1 ∈ vl && v2 ∈ vl && v3 ∈ vl && v4 ∈ vl
-        # Note: The redundant vertex v5 is not removed (see #561).
-        # This test can be removed (and the length above should be corrected)
-        # when that issue is resolved.
-        @test length(vl) == 5 && v5 ∈ vl
+        @test ispermutation(vl, [v1, v2, v3, v4])
 
         # Cartesian product
         p1 = VPolytope([N[0, 0], N[1, 1]])