diff --git a/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp b/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
index 6b4da26dbf96a9a54026d9732e528ca6d7651da3..061ee9a05fb0291ca7e8c4bfeefa35b12c725616 100644
--- a/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
@@ -6,11 +6,8 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CuboidVolumeBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
+
 #include <cmath>
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
@@ -84,11 +81,9 @@ class CuboidVolumeBounds : public VolumeBounds {
   bool inside(const Vector3D& pos, double tol = 0.) const override;
 
   /// Method to decompose the Bounds into boundarySurfaces
-  /// @note this method is a pure factory the volume is resposible
-  /// for the memory management
   ///
   /// @param transformPtr is the transfrom of the volume
-  std::vector<std::shared_ptr<const Surface>> decomposeToSurfaces(
+  SurfacePtrVector decomposeToSurfaces(
       const Transform3D* transformPtr) const override;
 
   /// Construct bounding box for this shape
diff --git a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
index b4738c311e25075f53c888ec7a01198a6d573b84..90c9aff5611f13812463de1e612b539b4296b1c6 100644
--- a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
@@ -8,10 +8,11 @@
 
 #pragma once
 
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
 #include "Acts/Surfaces/DiscSurface.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
@@ -19,116 +20,87 @@ namespace Acts {
 
 class IVisualization;
 
-/**
- * Class which implements a cutout cylinder. This shape is bascially a cylinder,
- * with another, smaller cylinder subtracted from the center.
- * --------------------- rmax
- * |                   |
- * |    |---------|    | rmed
- * |    |         |    |
- * ------         ------ rmin
- *       -- dz2 --
- * -------- dz1 -------
- *
- */
+/// Class which implements a cutout cylinder. This shape is bascially a
+/// cylinder, with another, smaller cylinder subtracted from the center.
+/// --------------------- rmax
+/// |                   |
+/// |    |---------|    | rmed
+/// |    |         |    |
+/// ------         ------ rmin
+///       -- dz2 --
+/// -------- dz1 -------
+///
+///
 class CutoutCylinderVolumeBounds : public VolumeBounds {
  public:
-  /**
-   * Constructor from defining parameters
-   * @param rmin Minimum radius at the "choke points"
-   * @param rmed The medium radius (outer radius of the cutout)
-   * @param rmax The outer radius of the overall shape
-   * @param dz1 The longer halflength of the shape
-   * @param dz2 The shorter halflength of the shape
-   */
+  /// Constructor from defining parameters
+  ///
+  /// @param rmin Minimum radius at the "choke points"
+  /// @param rmed The medium radius (outer radius of the cutout)
+  /// @param rmax The outer radius of the overall shape
+  /// @param dz1 The longer halflength of the shape
+  /// @param dz2 The shorter halflength of the shape
   CutoutCylinderVolumeBounds(double rmin, double rmed, double rmax, double dz1,
                              double dz2)
       : m_rmin(rmin), m_rmed(rmed), m_rmax(rmax), m_dz1(dz1), m_dz2(dz2) {}
 
-  /**
-   * Virtual default constructor
-   */
+  /// Virtual default constructor
   ~CutoutCylinderVolumeBounds() override = default;
 
-  /**
-   * Clone method.
-   * @return Pointer to a copy of the shape
-   */
+  /// Clone method.
+  /// @return Pointer to a copy of the shape
   VolumeBounds* clone() const override;
 
-  /**
-   * Inside method to test whether a point is inside the shape
-   * @param gpos The point to test
-   * @param tol The tolerance to test with
-   * @return Whether the point is inside or not.
-   */
+  /// Inside method to test whether a point is inside the shape
+  ///
+  /// @param gpos The point to test
+  /// @param tol The tolerance to test with
+  /// @return Whether the point is inside or not.
   bool inside(const Vector3D& gpos, double tol = 0) const override;
 
-  /**
-   * Method to decompose the Bounds into Surfaces
-   *
-   * @param transform is the transform to position the surfaces in 3D space
-   * @note this is a factory method
-   *
-   * @return vector of surfaces from the decopmosition
-   */
+  /// Method to decompose the Bounds into Surfaces
+  ///
+  /// @param transform is the transform to position the surfaces in 3D space
+  ///
+  /// @return vector of surfaces from the decopmosition
+  ///
   std::vector<std::shared_ptr<const Surface>> decomposeToSurfaces(
-      const Transform3D* transform) const override;
-
-  /**
-   * Construct bounding box for this shape
-   * @param trf Optional transform
-   * @param envelope Optional envelope to add / subtract from min/max
-   * @param entity Entity to associate this bounding box with
-   * @return Constructed bounding box
-   */
+      const Transform3D* transform = nullptr) const override;
+
+  /// Construct bounding box for this shape
+  ///
+  /// @param trf Optional transform
+  /// @param envelope Optional envelope to add / subtract from min/max
+  /// @param entity Entity to associate this bounding box with
+  /// @return Constructed bounding box
   Volume::BoundingBox boundingBox(const Transform3D* trf = nullptr,
                                   const Vector3D& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
-  /**
-   * Write information about this instance to an outstream
-   * @param sl The outstream
-   * @return The outstream
-   */
+  /// Write information about this instance to an outstream
+  ///
+  /// @param sl The outstream
+  /// @return The outstream
   std::ostream& toStream(std::ostream& sl) const override;
 
-  /**
-   * Draw this shape using a visualization helper
-   * @param helper The visualizatin helper
-   * @param transform Optional transformation matrix
-   */
-  void draw(IVisualization& helper,
-            const Transform3D& transform = Transform3D::Identity()) const;
-
-  /**
-   * Return the minimum radius
-   * @return The minimum radius
-   */
+  /// Return the minimum radius
+  /// @return The minimum radius
   double rMin() const { return m_rmin; }
 
-  /**
-   * Return the medium radius
-   * @return The medium radius
-   */
+  /// Return the medium radius
+  /// @return The medium radius
   double rMed() const { return m_rmed; }
 
-  /**
-   * Return the maximum radius
-   * @return The maximum radius
-   */
+  /// Return the maximum radius
+  /// @return The maximum radius
   double rMax() const { return m_rmax; }
 
-  /**
-   * Return the longer halflength in z.
-   * @return The halflength
-   */
+  /// Return the longer halflength in z.
+  /// @return The halflength
   double dZ1() const { return m_dz1; }
 
-  /**
-   * Return the shorter halflength in z.
-   * @return The halflength
-   */
+  /// Return the shorter halflength in z.
+  /// @return The halflength
   double dZ2() const { return m_dz2; }
 
  private:
diff --git a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
index 291f394cdf5a5f884d8d223f1c658610174723b5..4239e632aa35a140994a8a05b5cfbb5a97a4c29a 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
@@ -6,12 +6,9 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CylinderVolumeBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
+#include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
@@ -143,9 +140,8 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// Method to decompose the Bounds into boundarySurfaces
   /// @param transformPtr is the transform where the boundary surfaces are
   /// situated
-  /// @note this surface is a factory method, the volume handles the memory
   std::vector<std::shared_ptr<const Surface>> decomposeToSurfaces(
-      const Transform3D* transformPtr) const override;
+      const Transform3D* transformPtr = nullptr) const override;
 
   /// Construct bounding box for this shape
   /// @param trf Optional transform
@@ -187,12 +183,6 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// Output Method for std::ostream
   std::ostream& toStream(std::ostream& sl) const override;
 
-  /// Draw this cylinder to a given helper
-  /// @param helper The helper instance
-  /// @param transform An additional transform, default is identity
-  void draw(IVisualization& helper,
-            const Transform3D& transform = Transform3D::Identity()) const;
-
  private:
   /// templated dumpT method
   template <class T>
diff --git a/Core/include/Acts/Geometry/Extent.hpp b/Core/include/Acts/Geometry/Extent.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e8826cc47870946d54ce6f894bbb1c4b8b8c88e5
--- /dev/null
+++ b/Core/include/Acts/Geometry/Extent.hpp
@@ -0,0 +1,93 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+
+#include <iosfwd>
+#include <utility>
+#include <vector>
+
+namespace Acts {
+
+using Range = std::pair<double, double>;
+
+// @brief Extent in space
+///
+/// This is a nested struct to the GeometryObject representation
+/// which can be retrieved and used for surface parsing and will
+/// give you the maximal extent in 3D space/
+struct Extent {
+  /// Possible maximal value
+  static constexpr double maxval = std::numeric_limits<double>::max();
+
+  /// Start value
+  static constexpr Range maxrange = {maxval, -maxval};
+
+  // The different ranges
+  std::vector<Range> ranges = std::vector<Range>(9, maxrange);
+
+  // Constructor
+  Extent() = default;
+
+  /// Extend with another extent
+  /// @param other is the source Extent
+  void extend(const Extent& other) {
+    for (auto ir = 0; ir < other.ranges.size(); ++ir) {
+      ranges[ir].first = std::min(ranges[ir].first, other.ranges[ir].first);
+      ranges[ir].second = std::max(ranges[ir].second, other.ranges[ir].second);
+    }
+  }
+
+  /// Convert to output stream for screen output
+  /// @param sl [in,out] The output stream
+  std::ostream& toStream(std::ostream& sl) const;
+
+  /// Access the minimum parameter
+  /// @param bval the binning identification
+  double min(BinningValue bval) const { return ranges[bval].first; }
+
+  /// Access the max parameter
+  /// @param bval the binning identification
+  double max(BinningValue bval) const { return ranges[bval].second; }
+
+  /// Access the medium parameter
+  /// @param bval the binning identification
+  double medium(BinningValue bval) const {
+    return 0.5 * (ranges[bval].first + ranges[bval].second);
+  }
+
+  /// Access the range - always positive
+  /// @param bval the binning identification
+  double range(BinningValue bval) const {
+    return std::abs(ranges[bval].second - ranges[bval].first);
+  }
+
+  /// Check the vertex
+  /// @param vtx the Vertex to be checked
+  void check(const Vector3D& vtx) {
+    // min/max value check
+    auto minMax = [&](BinningValue bval, double value) -> void {
+      ranges[bval].first = std::min(value, ranges[bval].first);
+      ranges[bval].second = std::max(value, ranges[bval].second);
+    };
+    // Walk through the binning parameters
+    for (int bval = 0; bval < binValues; ++bval) {
+      BinningValue bValue = static_cast<BinningValue>(bval);
+      minMax(bValue, VectorHelpers::cast(vtx, bValue));
+    }
+  }
+};
+
+/// Overload of << operator for std::ostream for debug output
+std::ostream& operator<<(std::ostream& sl, const Extent& ext);
+
+}  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Geometry/GeometryObject.hpp b/Core/include/Acts/Geometry/GeometryObject.hpp
index cb33dc613bb006e2fb5c1345f9bb39e605ac78b3..4faa288adffcb0debf3bab653d18b7292ed58040 100644
--- a/Core/include/Acts/Geometry/GeometryObject.hpp
+++ b/Core/include/Acts/Geometry/GeometryObject.hpp
@@ -1,19 +1,16 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// GeometryObject.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryID.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/Helpers.hpp"
@@ -30,15 +27,18 @@ namespace Acts {
 ///
 class GeometryObject {
  public:
+  /// Defaulted construrctor
   GeometryObject() = default;
+
+  /// Defaulted copy constructor
   GeometryObject(const GeometryObject&) = default;
 
-  /// constructor from a ready-made value
+  /// Constructor from a value
   ///
   /// @param geoID the geometry identifier of the object
   GeometryObject(const GeometryID& geoID) : m_geoID(geoID) {}
 
-  /// assignment operator
+  /// Assignment operator
   ///
   /// @param geoID the source geoID
   GeometryObject& operator=(const GeometryObject& geoID) {
@@ -48,7 +48,6 @@ class GeometryObject {
     return *this;
   }
 
-  /// Return the value
   /// @return the geometry id by reference
   const GeometryID& geoID() const;
 
@@ -89,28 +88,6 @@ inline void GeometryObject::assignGeoID(const GeometryID& geoID) {
 
 inline double GeometryObject::binningPositionValue(const GeometryContext& gctx,
                                                    BinningValue bValue) const {
-  using VectorHelpers::perp;
-  // now switch
-  switch (bValue) {
-    // case x
-    case Acts::binX: {
-      return binningPosition(gctx, bValue).x();
-    } break;
-    // case y
-    case Acts::binY: {
-      return binningPosition(gctx, bValue).y();
-    } break;
-    // case z
-    case Acts::binZ: {
-      return binningPosition(gctx, bValue).z();
-    } break;
-    // case - to be overwritten by disks
-    case Acts::binR: {
-      return perp(binningPosition(gctx, bValue));
-    } break;
-    // do nothing for the default
-    default:
-      return 0.;
-  }
+  return VectorHelpers::cast(binningPosition(gctx, bValue), bValue);
 }
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/Polyhedron.hpp b/Core/include/Acts/Geometry/Polyhedron.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..28b823d8ce09abec10bbc5b3cb4d50e3fb8054f3
--- /dev/null
+++ b/Core/include/Acts/Geometry/Polyhedron.hpp
@@ -0,0 +1,104 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2018-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Geometry/Extent.hpp"
+#include "Acts/Surfaces/BoundaryCheck.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+
+#include <limits>
+#include <vector>
+
+namespace Acts {
+
+/// @class Polyhedron
+///
+/// Struct which contains a cartesian approximation for any surface type.
+/// It contains a list of cartesian vertices in the global frame, and
+/// additionally
+/// a list of lists of indices which indicate which vertices form a face.
+/// Each entry in @c faces is a face, which is in turn a list of vertices
+/// that need to be connected to form a face.
+/// This allows the @c objString method to produce a ready-to-go obj output.
+struct Polyhedron {
+  using Face = std::vector<size_t>;
+
+  /// Default constructor
+  Polyhedron() = default;
+
+  /// Default constructor from a vector of vertices and a vector of faces
+  /// @param verticesIn The 3D global vertices that make up the object
+  /// @param facesIn List of lists of indices for faces.
+  /// @param triangularMeshIn List of lists of indices for a triangular mesh
+  /// @param radialCheck A dedicated check for radial extent done
+  /// @note This creates copies of the input vectors
+  Polyhedron(const std::vector<Vector3D>& verticesIn,
+             const std::vector<Face>& facesIn,
+             const std::vector<Face>& triangularMeshIn,
+             bool radialCheck = false)
+      : vertices(verticesIn),
+        faces(facesIn),
+        triangularMesh(triangularMeshIn) {}
+
+  /// List of 3D vertices as vectors
+  std::vector<Vector3D> vertices;
+
+  /// List of faces connecting the vertices.
+  /// each face is a list of vertices v
+  /// corresponding to the vertex vector above
+  std::vector<Face> faces;
+
+  /// List of faces connecting the vertices.
+  /// each face is a list of vertices v
+  /// - in this case restricted to a triangular representation
+  std::vector<Face> triangularMesh;
+
+  /// Merge another Polyhedron into this one
+  ///
+  /// @param other is the source representation
+  void merge(const Polyhedron& other);
+
+  /// Draw method for polyhedrons
+  ///
+  /// @tparam helper_t The draw helper
+  ///
+  /// @param helper The draw helper object (visitor pattern)
+  /// @param triangulate Force the faces to be a triangular mesh
+  /// @param decompose Boolean that forces a decomposition into
+  /// individual faces with unique vertices.
+  template <typename helper_t>
+  void draw(helper_t& helper, bool triangulate = false,
+            bool decompose = false) const {
+    // vertices and faces are
+    if (not decompose and not triangulate) {
+      helper.faces(vertices, faces);
+    } else if (triangulate) {
+      helper.faces(vertices, triangularMesh);
+    } else {
+      for (const auto& face : faces) {
+        std::vector<Vector3D> face_vtx;
+        for (size_t i : face) {
+          face_vtx.push_back(vertices[i]);
+        }
+        helper.face(face_vtx);
+      }
+    }
+  }
+
+  /// Maximum extent of the polyhedron in space
+  ///
+  /// @param transform An (optional) transform
+  /// to apply to the vertices for estimation the extent
+  /// with respect to a given coordinate frame
+  ///
+  /// @return ranges that describe the space taken by this surface
+  Extent extent(const Transform3D& transform = Transform3D::Identity()) const;
+};
+}  // namespace Acts
diff --git a/Core/include/Acts/Geometry/VolumeBounds.hpp b/Core/include/Acts/Geometry/VolumeBounds.hpp
index b141b5b70d10de8a95d7100d3a187e20cfce6084..e39dcedffc70301babc1bf5d1272ac01a13271dd 100644
--- a/Core/include/Acts/Geometry/VolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/VolumeBounds.hpp
@@ -6,10 +6,6 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// VolumeBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <iomanip>
 #include <iostream>
@@ -19,6 +15,20 @@
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
+#ifndef VOLUMEBOUNDS_VALUESTORE_FILL
+#define VOLUMEBOUNDS_VALUESTORE_FILL(val) m_valueStore[bv_##val] = val
+#endif
+
+#ifndef VOLUMEBOUNDS_VALUESTORE_ACCESS
+#define VOLUMEBOUNDS_VALUESTORE_ACCESS(val) \
+  double val() const { return m_valueStore[bv_##val]; }
+#endif
+
+#ifndef VOLUMEBOUNDS_DERIVED_ACCESS
+#define VOLUMEBOUNDS_DERIVED_ACCESS(derived) \
+  double derived() const { return m_##derived; }
+#endif
+
 namespace Acts {
 
 class Surface;
@@ -27,6 +37,8 @@ class Volume;
 class VolumeBounds;
 using VolumeBoundsPtr = std::shared_ptr<const VolumeBounds>;
 
+using SurfacePtr = std::shared_ptr<const Surface>;
+using SurfacePtrVector = std::vector<SurfacePtr>;
 /// @class VolumeBounds
 ///
 /// Pure Absract Base Class for Volume bounds.
@@ -45,8 +57,10 @@ class VolumeBounds {
  public:
   /// Default Constructor*/
   VolumeBounds() = default;
+
   /// Destructor
   virtual ~VolumeBounds() = default;
+
   ///  clone() method to make deep copy in Volume copy constructor and for
   /// assigment operator  of the Surface class.
   virtual VolumeBounds* clone() const = 0;
@@ -64,10 +78,9 @@ class VolumeBounds {
   ///
   /// @param transform is the 3D transform to be applied to the boundary
   /// surfaces to position them in 3D space
-  /// @note this is factory method
   ///
   /// @return a vector of surfaces bounding this volume
-  virtual std::vector<std::shared_ptr<const Surface>> decomposeToSurfaces(
+  virtual SurfacePtrVector decomposeToSurfaces(
       const Transform3D* transform) const = 0;
 
   /// Construct bounding box for this shape
diff --git a/Core/include/Acts/Surfaces/AnnulusBounds.hpp b/Core/include/Acts/Surfaces/AnnulusBounds.hpp
index 6136fee114cc0af3c8943d06e0bfbfab827ba97b..6b302da3773c732f88160b8e29838bbb3f1520ef 100644
--- a/Core/include/Acts/Surfaces/AnnulusBounds.hpp
+++ b/Core/include/Acts/Surfaces/AnnulusBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// AnnulusBounds.hpp
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include "Acts/Surfaces/DiscBounds.hpp"
@@ -133,11 +129,20 @@ class AnnulusBounds : public DiscBounds {
   std::vector<Vector2D> corners() const;
 
   /// This method returns the xy coordinates of the four corners of the
-  /// bounds in module coorindates
+  /// bounds in module coorindates (in x/y)
   /// Starting from the upper right (max R, pos locX) and proceding clock-wise
   /// i.e. (max R; pos locX), (min R; pos locX), (min R; neg loc X), (max R: neg
   /// locX)
-  std::vector<Vector2D> vertices() const;
+  ///
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note that the extremas are given, which may slightly alter the
+  /// number of segments returned
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg) const;
 
  private:
   double m_rMin;
diff --git a/Core/include/Acts/Surfaces/BoundaryCheck.hpp b/Core/include/Acts/Surfaces/BoundaryCheck.hpp
index 5a9fc9bc2c03d19be35ddcfe075023c6db8d4af2..754fbc956d44e7a5abe1ae355d30ea4259570d36 100644
--- a/Core/include/Acts/Surfaces/BoundaryCheck.hpp
+++ b/Core/include/Acts/Surfaces/BoundaryCheck.hpp
@@ -1,21 +1,18 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// BoundaryCheck.hpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cfloat>
 #include <cmath>
 #include <iterator>
 #include <vector>
 
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/ParameterDefinitions.hpp"
 
@@ -136,15 +133,6 @@ class BoundaryCheck {
   ///          for the tolerance based check.
   BoundaryCheck transformed(const ActsMatrixD<2, 2>& jacobian) const;
 
-  /// Check if the point is inside the polygon w/o any tolerances.
-  template <typename Vector2DContainer>
-  bool isInsidePolygon(const Vector2D& point,
-                       const Vector2DContainer& vertices) const;
-
-  /// Check if the point is inside the aligned box
-  bool isInsideRectangle(const Vector2D& point, const Vector2D& lowerLeft,
-                         const Vector2D& upperRight) const;
-
   /// Check if the distance vector is within the absolute or relative limits.
   bool isTolerated(const Vector2D& delta) const;
 
@@ -173,7 +161,7 @@ class BoundaryCheck {
   friend class RectangleBounds;
   // To be able to use `transformed`
   friend class CylinderBounds;
-  friend class DiscTrapezoidalBounds;
+  friend class DiscTrapezoidBounds;
   // EllipseBounds needs a custom implementation
   friend class EllipseBounds;
 };
@@ -230,7 +218,7 @@ inline bool Acts::BoundaryCheck::isInside(
   if (m_type == Type::eNone) {
     // The null boundary check always succeeds
     return true;
-  } else if (isInsidePolygon(point, vertices)) {
+  } else if (detail::VerticesHelper::isInsidePolygon(point, vertices)) {
     // If the point falls inside the polygon, the check always succeeds
     return true;
   } else if (m_tolerance == Vector2D(0., 0.)) {
@@ -253,53 +241,10 @@ inline bool Acts::BoundaryCheck::isInside(
   }
 }
 
-template <typename Vector2DContainer>
-inline bool Acts::BoundaryCheck::isInsidePolygon(
-    const Vector2D& point, const Vector2DContainer& vertices) const {
-  // when we move along the edges of a convex polygon, a point on the inside of
-  // the polygon will always appear on the same side of each edge.
-  // a point on the outside will switch sides at least once.
-
-  // returns which side of the connecting line between `ll0` and `ll1` the point
-  // `p` is on. computes the sign of the z-component of the cross-product
-  // between the line normal vector and the vector from `ll0` to `p`.
-  auto lineSide = [&](auto&& ll0, auto&& ll1) {
-    auto normal = ll1 - ll0;
-    auto delta = point - ll0;
-    return std::signbit((normal[0] * delta[1]) - (normal[1] * delta[0]));
-  };
-
-  auto iv = std::begin(vertices);
-  Vector2D l0 = *iv;
-  Vector2D l1 = *(++iv);
-  // use vertex0 to vertex1 to define reference sign and compare w/ all edges
-  auto reference = lineSide(l0, l1);
-  for (++iv; iv != std::end(vertices); ++iv) {
-    l0 = l1;
-    l1 = *iv;
-    if (lineSide(l0, l1) != reference) {
-      return false;
-    }
-  }
-  // manual check for last edge from last vertex back to the first vertex
-  if (lineSide(l1, *std::begin(vertices)) != reference) {
-    return false;
-  }
-  // point was always on the same side. point must be inside.
-  return true;
-}
-
-inline bool Acts::BoundaryCheck::isInsideRectangle(
-    const Vector2D& point, const Vector2D& lowerLeft,
-    const Vector2D& upperRight) const {
-  return (lowerLeft[0] <= point[0]) && (point[0] < upperRight[0]) &&
-         (lowerLeft[1] <= point[1]) && (point[1] < upperRight[1]);
-}
-
 inline bool Acts::BoundaryCheck::isInside(const Vector2D& point,
                                           const Vector2D& lowerLeft,
                                           const Vector2D& upperRight) const {
-  if (isInsideRectangle(point, lowerLeft, upperRight)) {
+  if (detail::VerticesHelper::isInsideRectangle(point, lowerLeft, upperRight)) {
     return true;
   } else {
     Vector2D closestPoint;
@@ -328,7 +273,7 @@ inline double Acts::BoundaryCheck::distance(
   // TODO 2017-04-06 msmk: this should be calculable directly
   double d = squaredNorm(point - computeClosestPointOnPolygon(point, vertices));
   d = std::sqrt(d);
-  return isInsidePolygon(point, vertices) ? -d : d;
+  return detail::VerticesHelper::isInsidePolygon(point, vertices) ? -d : d;
 }
 
 inline double Acts::BoundaryCheck::distance(const Acts::Vector2D& point,
@@ -339,7 +284,10 @@ inline double Acts::BoundaryCheck::distance(const Acts::Vector2D& point,
     double d = (point - computeEuclideanClosestPointOnRectangle(
                             point, lowerLeft, upperRight))
                    .norm();
-    return isInsideRectangle(point, lowerLeft, upperRight) ? -d : d;
+    return detail::VerticesHelper::isInsideRectangle(point, lowerLeft,
+                                                     upperRight)
+               ? -d
+               : d;
 
   } else /* Type::eChi2 */ {
     Vector2D vertices[] = {{lowerLeft[0], lowerLeft[1]},
diff --git a/Core/include/Acts/Surfaces/ConeBounds.hpp b/Core/include/Acts/Surfaces/ConeBounds.hpp
index 877e0a89d7d6642460423eb4299932cedd43d2c1..9f8d4a964edc2af54ed10da04a22229a2b32a044 100644
--- a/Core/include/Acts/Surfaces/ConeBounds.hpp
+++ b/Core/include/Acts/Surfaces/ConeBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// ConeBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cfloat>
 
@@ -41,6 +37,7 @@ class ConeBounds : public SurfaceBounds {
     bv_length = 5
   };
 
+  // Deleted default constructor
   ConeBounds() = delete;
 
   /// Constructor - open cone with alpha, by default a full cone
@@ -65,12 +62,16 @@ class ConeBounds : public SurfaceBounds {
   ConeBounds(double alpha, double zmin, double zmax, double halfphi = M_PI,
              double avphi = 0.);
 
-  ~ConeBounds() override;
+  /// Defaulted destructor
+  ~ConeBounds() override = default;
 
+  /// Virtual constructor
   ConeBounds* clone() const final;
 
+  /// The type enumeration
   BoundsType type() const final;
 
+  /// The value store for persistency
   std::vector<TDD_real_t> valueStore() const final;
 
   /// inside method for local position
diff --git a/Core/include/Acts/Surfaces/ConeSurface.hpp b/Core/include/Acts/Surfaces/ConeSurface.hpp
index b8c93cd549241dff882088fd80fc42358dc9df48..fd2b36f3a9c287634dac1c05d65df32792a912cc 100644
--- a/Core/include/Acts/Surfaces/ConeSurface.hpp
+++ b/Core/include/Acts/Surfaces/ConeSurface.hpp
@@ -1,18 +1,15 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// ConeSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/ConeBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Definitions.hpp"
@@ -196,7 +193,7 @@ class ConeSurface : public Surface {
                                 const Vector3D& direction,
                                 const BoundaryCheck& bcheck) const final;
 
-  /// the pathCorrection for derived classes with thickness
+  /// The pathCorrection for derived classes with thickness
   ///
   /// @param gctx The current geometry context object, e.g. alignment
   /// @param position is the global potion at the correction point
@@ -205,6 +202,19 @@ class ConeSurface : public Surface {
   double pathCorrection(const GeometryContext& gctx, const Vector3D& position,
                         const Vector3D& direction) const final;
 
+  /// Return a Polyhedron for the surfaces
+  ///
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param lseg Number of segments along curved lines, it represents
+  /// the full 2*M_PI coverange, if lseg is set to 1 only the extrema
+  /// are given
+  /// @note that a surface transform can invalidate the extrema
+  /// in the transformed space
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                      size_t lseg) const override;
+
   /// Return properly formatted class name for screen output
   std::string name() const override;
 
diff --git a/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp b/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp
index 3ebfffcb5c61e7711f44d45b9d62a06d691dfb40..d40b967b19a3cbe9242f423e4a8e06a359500dc9 100644
--- a/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp
+++ b/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -17,136 +17,109 @@
 
 namespace Acts {
 
-/**
- * This class serves as a base class for the actual bounds class.
- * The only deriving type is the templated `ConvexPolygonBounds`.
- */
+/// @brief base class for convex polygon bounds
+///
+/// This class serves as a base class for the actual bounds class.
+/// The only deriving type is the templated `ConvexPolygonBounds`.
+///
 class ConvexPolygonBoundsBase : public PlanarBounds {
  public:
-  /**
-   * Output Method for std::ostream
-   * @param sl is the ostream to be written into
-   */
+  /// Output Method for std::ostream
+  /// @param sl is the ostream to be written into
   std::ostream& toStream(std::ostream& sl) const final;
 
-  /**
-   * Return vector containing defining parameters
-   * @return the parameters
-   */
+  /// Return vector containing defining parameters
+  /// @return the parameters
   std::vector<TDD_real_t> valueStore() const final;
 
  protected:
-  /**
-   * Return a rectangle bounds instance that encloses a set of vertices.
-   * @param vertices A collection of vertices to enclose.
-   * @return Enclosing rectangle.
-   */
+  /// Return a rectangle bounds instance that encloses a set of vertices.
+  /// @param vertices A collection of vertices to enclose.
+  /// @return Enclosing rectangle.
   template <typename coll_t>
   static RectangleBounds makeBoundingBox(const coll_t& vertices);
 
-  /**
-   * Calculates whether a set of vertices forms a convex polygon. This is
-   * generic over the number of vertices, so it's factored out of the concrete
-   * classes and into this base class.
-   * @param vertices A collection of vertices.
-   * @return Whether the vertices form a convex polygon.
-   */
+  /// Calculates whether a set of vertices forms a convex polygon. This is
+  /// generic over the number of vertices, so it's factored out of the concrete
+  /// classes and into this base class.
+  /// @param vertices A collection of vertices.
+  /// @return Whether the vertices form a convex polygon.
   template <typename coll_t>
   static bool convex_impl(const coll_t& vertices);
 };
 
-/**
- * This is the actual implementation of the bounds.
- * It is templated on the number of vertices, but there is a specialization for
- * *dynamic* number of vertices, where the underlying storage is then a vector.
- *
- * @tparam N Number of vertices
- */
+/// This is the actual implementation of the bounds.
+/// It is templated on the number of vertices, but there is a specialization for
+/// *dynamic* number of vertices, where the underlying storage is then a vector.
+///
+/// @tparam N Number of vertices
 template <int N>
 class ConvexPolygonBounds : public ConvexPolygonBoundsBase {
-  /**
-   * Expose number of vertices given as template parameter.
-   */
+  /// Expose number of vertices given as template parameter.
+  ///
   static constexpr size_t num_vertices = N;
-  /**
-   * Type that's used to store the vertices, in this case a fixed size array.
-   */
+  /// Type that's used to store the vertices, in this case a fixed size array.
+  ///
   using vertex_array = std::array<Vector2D, num_vertices>;
 
  public:
   static_assert(N >= 3, "ConvexPolygonBounds needs at least 3 sides.");
 
-  /**
-   * Default constructor, deleted
-   */
+  /// Default constructor, deleted
   ConvexPolygonBounds() = delete;
 
-  /**
-   * Constructor from a vector of vertices, to facilitate construction.
-   * This will throw if the vector size does not match `num_vertices`.
-   * This will throw if the vertices do not form a convex polygon.
-   * @param vertices The list of vertices.
-   */
+  /// Constructor from a vector of vertices, to facilitate construction.
+  /// This will throw if the vector size does not match `num_vertices`.
+  /// This will throw if the vertices do not form a convex polygon.
+  /// @param vertices The list of vertices.
   ConvexPolygonBounds(const std::vector<Vector2D>& vertices);
 
-  /**
-   * Constructor from a fixed size array of vertices.
-   * This will throw if the vertices do not form a convex polygon.
-   * @param vertices The vertices
-   */
+  /// Constructor from a fixed size array of vertices.
+  /// This will throw if the vertices do not form a convex polygon.
+  /// @param vertices The vertices
   ConvexPolygonBounds(const vertex_array& vertices);
 
-  /**
-   * Defaulted destructor
-   */
+  /// Defaulted destructor
   ~ConvexPolygonBounds() override = default;
 
-  /**
-   * Return a copy of this bounds object.
-   * @return The cloned instance
-   */
+  /// Return a copy of this bounds object.
+  /// @return The cloned instance
   ConvexPolygonBounds<N>* clone() const final;
 
-  /**
-   * Return the bounds type of this bounds object.
-   * @return The bounds type
-   */
+  /// Return the bounds type of this bounds object.
+  /// @return The bounds type
   BoundsType type() const final;
 
-  /**
-   * Return whether a local 2D point lies inside of the bounds defined by this
-   * object.
-   * @param lposition The local position to check
-   * @param bcheck The `BoundaryCheck` object handling tolerances.
-   * @return Whether the points is inside
-   */
+  /// Return whether a local 2D point lies inside of the bounds defined by this
+  /// object.
+  /// @param lposition The local position to check
+  /// @param bcheck The `BoundaryCheck` object handling tolerances.
+  /// @return Whether the points is inside
   bool inside(const Vector2D& lposition,
               const BoundaryCheck& bcheck) const final;
 
-  /**
-   * Return the smallest distance to any point on the boundary of this bounds
-   * object.
-   * @param lposition The local position to get the distance to
-   * @return The smallest distance to the boundary.
-   */
+  /// Return the smallest distance to any point on the boundary of this bounds
+  /// object.
+  /// @param lposition The local position to get the distance to
+  /// @return The smallest distance to the boundary.
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /**
-   * Returns a vector containing the vertices making up the bounds.
-   * @return Vector of vertices
-   */
-  std::vector<Vector2D> vertices() const final;
-
-  /**
-   * Return a rectangle bounds object that encloses this polygon.
-   * @return The rectangular bounds
-   */
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note the number of segements is ignored in this representation
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg = 1) const final;
+
+  /// Return a rectangle bounds object that encloses this polygon.
+  /// @return The rectangular bounds
   const RectangleBounds& boundingBox() const final;
 
-  /**
-   * Return whether this bounds class is in fact convex
-   * @return Whether the bounds are convex.
-   */
+  /// Return whether this bounds class is in fact convex
+  /// @return Whether the bounds are convex.
   bool convex() const;
 
  private:
@@ -154,81 +127,68 @@ class ConvexPolygonBounds : public ConvexPolygonBoundsBase {
   RectangleBounds m_boundingBox;
 };
 
-/**
- * Tag to trigger specialization of a dynamic polygon
- */
+/// Tag to trigger specialization of a dynamic polygon
 constexpr int PolygonDynamic = -1;
 
-/**
- * This is the specialization handling a polygon with a dynamic number of
- * points. It can accept any number of points.
- */
+/// This is the specialization handling a polygon with a dynamic number of
+/// points. It can accept any number of points.
+///
 template <>
 class ConvexPolygonBounds<PolygonDynamic> : public ConvexPolygonBoundsBase {
  public:
-  /**
-   * Default constructor, deleted
-   */
+  /// Default constructor, deleted
   ConvexPolygonBounds() = delete;
 
-  /**
-   * Defaulted destructor
-   */
+  /// Defaulted destructor
   ~ConvexPolygonBounds() override = default;
 
-  /**
-   * Constructor from a vector of vertices, to facilitate construction.
-   * This will throw if the vertices do not form a convex polygon.
-   * @param vertices The list of vertices.
-   */
+  /// Constructor from a vector of vertices, to facilitate construction.
+  /// This will throw if the vertices do not form a convex polygon.
+  /// @param vertices The list of vertices.
   ConvexPolygonBounds(const std::vector<Vector2D>& vertices);
 
-  /**
-   * Return a copy of this bounds object.
-   * @return The cloned instance
-   */
+  /// Return a copy of this bounds object.
+  /// @return The cloned instance
   ConvexPolygonBounds<PolygonDynamic>* clone() const final;
 
-  /**
-   * Return the bounds type of this bounds object.
-   * @return The bounds type
-   */
+  /// Return the bounds type of this bounds object.
+  /// @return The bounds type
   BoundsType type() const final;
 
-  /**
-   * Return whether a local 2D point lies inside of the bounds defined by this
-   * object.
-   * @param lposition The local position to check
-   * @param bcheck The `BoundaryCheck` object handling tolerances.
-   * @return Whether the points is inside
-   */
+  /// Return whether a local 2D point lies inside of the bounds defined by this
+  /// object.
+  /// @param lposition The local position to check
+  /// @param bcheck The `BoundaryCheck` object handling tolerances.
+  /// @return Whether the points is inside
   bool inside(const Vector2D& lposition,
               const BoundaryCheck& bcheck) const final;
 
-  /**
-   * Return the smallest distance to any point on the boundary of this bounds
-   * object.
-   * @param lpos The lposition position to get the distance to
-   * @return The smallest distance to the boundary.
-   */
+  /// Return the smallest distance to any point on the boundary of this bounds
+  /// object.
+  /// @param lpos The lposition position to get the distance to
+  /// @return The smallest distance to the boundary.
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /**
-   * Returns a vector containing the vertices making up the bounds.
-   * @return Vector of vertices
-   */
-  std::vector<Vector2D> vertices() const final;
-
-  /**
-   * Return a rectangle bounds object that encloses this polygon.
-   * @return The rectangular bounds
-   */
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note the number of segements is ignored in this representation
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg = 1) const final;
+
+  ///
+  /// Return a rectangle bounds object that encloses this polygon.
+  /// @return The rectangular bounds
+  ///
   const RectangleBounds& boundingBox() const final;
 
-  /**
-   * Return whether this bounds class is in fact convex
-   * @return Whether the bounds are convex.
-   */
+  ///
+  /// Return whether this bounds class is in fact convex
+  /// @return Whether the bounds are convex.
+  ///
   bool convex() const;
 
  private:
diff --git a/Core/include/Acts/Surfaces/ConvexPolygonBounds.ipp b/Core/include/Acts/Surfaces/ConvexPolygonBounds.ipp
index 0ca3cda2c3638760d43bb0f0d2363eb38f14f559..b7c8eb7c1774fd21fdd4891d47d4d2c2bf1872b6 100644
--- a/Core/include/Acts/Surfaces/ConvexPolygonBounds.ipp
+++ b/Core/include/Acts/Surfaces/ConvexPolygonBounds.ipp
@@ -126,7 +126,8 @@ double Acts::ConvexPolygonBounds<N>::distanceToBoundary(
 }
 
 template <int N>
-std::vector<Acts::Vector2D> Acts::ConvexPolygonBounds<N>::vertices() const {
+std::vector<Acts::Vector2D> Acts::ConvexPolygonBounds<N>::vertices(
+    unsigned int /*lseg*/) const {
   return {m_vertices.begin(), m_vertices.end()};
 }
 
@@ -164,8 +165,8 @@ double Acts::ConvexPolygonBounds<Acts::PolygonDynamic>::distanceToBoundary(
   return BoundaryCheck(true).distance(lposition, m_vertices);
 }
 
-std::vector<Acts::Vector2D>
-Acts::ConvexPolygonBounds<Acts::PolygonDynamic>::vertices() const {
+std::vector<Acts::Vector2D> Acts::ConvexPolygonBounds<
+    Acts::PolygonDynamic>::vertices(unsigned int /*lseg*/) const {
   return {m_vertices.begin(), m_vertices.end()};
 }
 
diff --git a/Core/include/Acts/Surfaces/CylinderBounds.hpp b/Core/include/Acts/Surfaces/CylinderBounds.hpp
index 2f3e77d4169e448be560ee541b8337ddb584d6ab..039c04c77a419d796f615e96e2ae9f2387620897 100644
--- a/Core/include/Acts/Surfaces/CylinderBounds.hpp
+++ b/Core/include/Acts/Surfaces/CylinderBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CylinderBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
 
@@ -72,12 +68,16 @@ class CylinderBounds : public SurfaceBounds {
   CylinderBounds(double radius, double averagePhi, double halfPhi,
                  double halfZ);
 
-  ~CylinderBounds() override;
+  /// Defaulted destructor
+  ~CylinderBounds() override = default;
 
+  /// Virtual constructor
   CylinderBounds* clone() const final;
 
+  /// Type enumeration
   BoundsType type() const final;
 
+  /// The value store
   std::vector<TDD_real_t> valueStore() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
diff --git a/Core/include/Acts/Surfaces/CylinderSurface.hpp b/Core/include/Acts/Surfaces/CylinderSurface.hpp
index 3d9381a9341e967f5ddaec5ebca2b0d3b97e0136..c91769264915ca510e12a67b5752796b9619db25 100644
--- a/Core/include/Acts/Surfaces/CylinderSurface.hpp
+++ b/Core/include/Acts/Surfaces/CylinderSurface.hpp
@@ -1,23 +1,19 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CylinderSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include <cmath>
 
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryStatics.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/CylinderBounds.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/detail/RealQuadraticEquation.hpp"
@@ -228,13 +224,16 @@ class CylinderSurface : public Surface {
   /// Return method for properly formatted output string
   std::string name() const override;
 
-  /// Return a PolyhedronRepresentation for this object
+  /// Return a Polyhedron for a cylinder
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param l0div Number of divisions along l0 (phi)
-  /// @param l1div Number of divisions along l1 (z)
-  virtual PolyhedronRepresentation polyhedronRepresentation(
-      const GeometryContext& gctx, size_t l0div = 10, size_t l1div = 1) const;
+  /// @param lseg Number of segments along curved lines, it represents
+  /// the full 2*M_PI coverange, if lseg is set to 1 only the extrema
+  /// are given
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                      size_t lseg) const override;
 
  protected:
   std::shared_ptr<const CylinderBounds> m_bounds;  //!< bounds (shared)
diff --git a/Core/include/Acts/Surfaces/DiamondBounds.hpp b/Core/include/Acts/Surfaces/DiamondBounds.hpp
index bdafab730702f7b9150a5c08155a704c5b206d0d..b0078dbfc0f542143e91fbdf806942036d2be489 100644
--- a/Core/include/Acts/Surfaces/DiamondBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiamondBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiamondBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
 
@@ -48,12 +44,16 @@ class DiamondBounds : public PlanarBounds {
   /// @image html DiamondBounds.svg
   DiamondBounds(double x1, double x2, double x3, double y1, double y2);
 
-  ~DiamondBounds() override;
+  /// Defaulted desctructor
+  ~DiamondBounds() override = default;
 
+  /// Virtual constructor
   DiamondBounds* clone() const final;
 
+  /// Enumeration type
   BoundsType type() const final;
 
+  /// The value store for persistency
   std::vector<TDD_real_t> valueStore() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
@@ -72,8 +72,15 @@ class DiamondBounds : public PlanarBounds {
   /// @return is a signed distance parameter
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /// Return the vertices - or, the points of the extremas
-  std::vector<Vector2D> vertices() const final;
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note the number of segements is ignored for this representation
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg = 1) const final;
 
   // Bounding box representation
   const RectangleBounds& boundingBox() const final;
diff --git a/Core/include/Acts/Surfaces/DiscBounds.hpp b/Core/include/Acts/Surfaces/DiscBounds.hpp
index e3cc0f86e3fa292a1949a9e5aab36af7ff1664ad..599574fe8e4979e2749379c8899806f146c9f333 100644
--- a/Core/include/Acts/Surfaces/DiscBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiscBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiscBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include "Acts/Surfaces/SurfaceBounds.hpp"
 
@@ -22,6 +18,23 @@ namespace Acts {
 
 class DiscBounds : public SurfaceBounds {
  public:
+  /// Return method for inner Radius
+  virtual double rMin() const = 0;
+
+  /// Return method for outer Radius
+  virtual double rMax() const = 0;
+
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line, the number referrs to full 2*PI
+  ///
+  /// @note that the extremas are given, which may slightly alter the
+  /// number of segments returned
+  ///
+  /// @return vector for vertices in 2D
+  virtual std::vector<Vector2D> vertices(unsigned int lseg) const = 0;
+
   /// Returns a reference radius for binning
   virtual double binningValueR() const = 0;
 
diff --git a/Core/include/Acts/Surfaces/DiscSurface.hpp b/Core/include/Acts/Surfaces/DiscSurface.hpp
index 1996329040b75f85af3ee9cb98271c116a8662a9..f105139d85a69ccd2c17a0b5751bfa0a41a8effb 100644
--- a/Core/include/Acts/Surfaces/DiscSurface.hpp
+++ b/Core/include/Acts/Surfaces/DiscSurface.hpp
@@ -1,22 +1,18 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiscSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryStatics.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/DiscBounds.hpp"
 #include "Acts/Surfaces/InfiniteBounds.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
 #include "Acts/Utilities/Definitions.hpp"
@@ -63,7 +59,7 @@ class DiscSurface : public Surface {
 
   /// Constructor for Discs from Transform3D, \f$ r_{min}, r_{max}, hx_{min},
   /// hx_{max} \f$
-  /// This is n this case you have DiscTrapezoidalBounds
+  /// This is n this case you have DiscTrapezoidBounds
   ///
   /// @param htrans is transform that places the disc in the global 3D space
   /// (can be nullptr)
@@ -310,12 +306,16 @@ class DiscSurface : public Surface {
   /// Return properly formatted class name for screen output
   std::string name() const override;
 
-  /// Return a PolyhedronRepresentation for this object
+  /// Return a Polyhedron for the surfaces
+  ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param l0div Number of divisions along l0 (phi)
-  /// @param l1div Number of divisions along l1 (r)
-  virtual PolyhedronRepresentation polyhedronRepresentation(
-      const GeometryContext& gctx, size_t l0div = 10, size_t l1div = 1) const;
+  /// @param lseg Number of segments along curved lines, it represents
+  /// the full 2*M_PI coverange, if lseg is set to 1 only the extrema
+  /// are given
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                      size_t lseg) const override;
 
  protected:
   std::shared_ptr<const DiscBounds> m_bounds;  ///< bounds (shared)
diff --git a/Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp b/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp
similarity index 72%
rename from Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp
rename to Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp
index 21adcd9aef57e50bf3b7e282f86204103f63d1fc..1d9057977c16c648d5fa3f0259d3ae1c2d296629 100644
--- a/Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiscTrapezoidalBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
 
@@ -20,14 +16,14 @@
 namespace Acts {
 
 ///
-/// @class DiscTrapezoidalBounds
+/// @class DiscTrapezoidBounds
 ///
 /// Class to describe the bounds for a planar DiscSurface.
 /// By providing an argument for hphisec, the bounds can
 /// be restricted to a phi-range around the center position.
 ///
 
-class DiscTrapezoidalBounds : public DiscBounds {
+class DiscTrapezoidBounds : public DiscBounds {
  public:
   /// @enum BoundValues
   /// enumeration for readability
@@ -41,26 +37,29 @@ class DiscTrapezoidalBounds : public DiscBounds {
     bv_length = 6
   };
 
-  DiscTrapezoidalBounds() = delete;
+  DiscTrapezoidBounds() = delete;
 
   /// Constructor for a symmetric Trapezoid giving min X length, max X length,
   /// Rmin and R max
   /// @param minhalfx half length in X at min radius
   /// @param maxhalfx half length in X at maximum radius
-  /// @param maxR outer radius
   /// @param minR inner radius
+  /// @param maxR outer radius
   /// @param avephi average phi value
   /// @param stereo optional stero angle applied
-  DiscTrapezoidalBounds(double minhalfx, double maxhalfx, double maxR,
-                        double minR, double avephi = M_PI_2,
-                        double stereo = 0.);
+  DiscTrapezoidBounds(double minhalfx, double maxhalfx, double minR,
+                      double maxR, double avephi = M_PI_2, double stereo = 0.);
 
-  ~DiscTrapezoidalBounds() override;
+  /// Defaulted Destructor
+  ~DiscTrapezoidBounds() override = default;
 
-  DiscTrapezoidalBounds* clone() const final;
+  /// Overloaded clone method
+  DiscTrapezoidBounds* clone() const final;
 
+  /// Type identifier
   SurfaceBounds::BoundsType type() const final;
 
+  /// Value store for persistency
   std::vector<TDD_real_t> valueStore() const final;
 
   ///  This method cheks if the radius given in the LocalPosition is inside
@@ -121,6 +120,17 @@ class DiscTrapezoidalBounds : public DiscBounds {
   /// Return a reference phi for binning
   double binningValuePhi() const final;
 
+  /// This method returns the xy coordinates of the four corners of the
+  /// bounds in module coorindates (in xy)
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note that the number of segments are ignored for this surface
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg) const;
+
  private:
   double m_rMin, m_rMax, m_minHalfX, m_maxHalfX, m_avgPhi;
   double m_stereo;  // TODO 2017-04-09 msmk: what is this good for?
@@ -138,62 +148,62 @@ class DiscTrapezoidalBounds : public DiscBounds {
   ActsMatrixD<2, 2> jacobianToLocalCartesian(const Vector2D& lposition) const;
 };
 
-inline double DiscTrapezoidalBounds::rMin() const {
+inline double DiscTrapezoidBounds::rMin() const {
   return m_rMin;
 }
 
-inline double DiscTrapezoidalBounds::rMax() const {
+inline double DiscTrapezoidBounds::rMax() const {
   return m_rMax;
 }
 
-inline double DiscTrapezoidalBounds::minHalflengthX() const {
+inline double DiscTrapezoidBounds::minHalflengthX() const {
   return m_minHalfX;
 }
 
-inline double DiscTrapezoidalBounds::maxHalflengthX() const {
+inline double DiscTrapezoidBounds::maxHalflengthX() const {
   return m_maxHalfX;
 }
 
-inline double DiscTrapezoidalBounds::averagePhi() const {
+inline double DiscTrapezoidBounds::averagePhi() const {
   return m_avgPhi;
 }
 
-inline double DiscTrapezoidalBounds::stereo() const {
+inline double DiscTrapezoidBounds::stereo() const {
   return m_stereo;
 }
 
-inline double DiscTrapezoidalBounds::halfPhiSector() const {
+inline double DiscTrapezoidBounds::halfPhiSector() const {
   auto minHalfPhi = std::asin(m_minHalfX / m_rMin);
   auto maxHalfPhi = std::asin(m_maxHalfX / m_rMax);
   return std::max(minHalfPhi, maxHalfPhi);
 }
 
-inline double DiscTrapezoidalBounds::rCenter() const {
+inline double DiscTrapezoidBounds::rCenter() const {
   auto hmin = std::sqrt(m_rMin * m_rMin - m_minHalfX * m_minHalfX);
   auto hmax = std::sqrt(m_rMax * m_rMax - m_maxHalfX * m_maxHalfX);
   return (hmin + hmax) / 2.0;
 }
 
-inline double DiscTrapezoidalBounds::halflengthY() const {
+inline double DiscTrapezoidBounds::halflengthY() const {
   auto hmin = std::sqrt(m_rMin * m_rMin - m_minHalfX * m_minHalfX);
   auto hmax = std::sqrt(m_rMax * m_rMax - m_maxHalfX * m_maxHalfX);
   return (hmax - hmin) / 2.0;
 }
 
-inline bool DiscTrapezoidalBounds::coversFullAzimuth() const {
+inline bool DiscTrapezoidBounds::coversFullAzimuth() const {
   return false;
 }
 
-inline bool DiscTrapezoidalBounds::insideRadialBounds(double R,
-                                                      double tolerance) const {
+inline bool DiscTrapezoidBounds::insideRadialBounds(double R,
+                                                    double tolerance) const {
   return (R + tolerance > m_rMin and R - tolerance < m_rMax);
 }
 
-inline double DiscTrapezoidalBounds::binningValueR() const {
+inline double DiscTrapezoidBounds::binningValueR() const {
   return 0.5 * (m_rMin + m_rMax);
 }
 
-inline double DiscTrapezoidalBounds::binningValuePhi() const {
+inline double DiscTrapezoidBounds::binningValuePhi() const {
   return m_avgPhi;
 }
 
diff --git a/Core/include/Acts/Surfaces/EllipseBounds.hpp b/Core/include/Acts/Surfaces/EllipseBounds.hpp
index 5f6e51c86133a5fcdb7220d5efd932b18629ebe6..be5999e0288ead8efbda4c7cf8bfcbb0629d20ec 100644
--- a/Core/include/Acts/Surfaces/EllipseBounds.hpp
+++ b/Core/include/Acts/Surfaces/EllipseBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// EllipseBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
 #include <cstdlib>
@@ -22,8 +18,8 @@ namespace Acts {
 
 /// @class EllipseBounds
 ///
-/// Class to describe the bounds for a planar EllipseSurface,
-/// i.e. the surface between two ellipses.
+/// Class to describe the bounds for a planar ellispoid
+/// surface.
 /// By providing an argument for hphisec, the bounds can
 /// be restricted to a phi-range around the center position.
 ///
@@ -42,6 +38,7 @@ class EllipseBounds : public PlanarBounds {
     bv_length = 6
   };
 
+  /// Deleted default constructor
   EllipseBounds() = delete;
 
   /// Constructor for full of an ellipsoid disc
@@ -56,12 +53,16 @@ class EllipseBounds : public PlanarBounds {
                 double maxRadius1, double averagePhi = 0.,
                 double halfPhi = M_PI);
 
-  ~EllipseBounds() override;
+  /// Defaulted destructor
+  ~EllipseBounds() override = default;
 
+  /// Clone method for surface cloning
   EllipseBounds* clone() const final;
 
+  /// Type enumeration
   BoundsType type() const final;
 
+  /// Complete value store for persistency
   std::vector<TDD_real_t> valueStore() const final;
 
   /// This method checks if the point given in the local coordinates is between
@@ -80,8 +81,16 @@ class EllipseBounds : public PlanarBounds {
   /// @return is a signed distance parameter
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /// Return the vertices - or, the points of the extremas
-  std::vector<Vector2D> vertices() const final;
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line, here it refers to the full 2PI Ellipse
+  ///
+  /// @note the number of segements to may be altered by also providing
+  /// the extremas in all direction
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg) const final;
 
   // Bounding box representation
   const RectangleBounds& boundingBox() const final;
diff --git a/Core/include/Acts/Surfaces/InfiniteBounds.hpp b/Core/include/Acts/Surfaces/InfiniteBounds.hpp
index 2a17de9462f1ef3dfaddb592605591baa119e6be..cb3dde3a4cd7254fbb8197aea1bfa051fa3f3235 100644
--- a/Core/include/Acts/Surfaces/InfiniteBounds.hpp
+++ b/Core/include/Acts/Surfaces/InfiniteBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// InfiniteBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include "Acts/Surfaces/SurfaceBounds.hpp"
 #include "Acts/Utilities/Definitions.hpp"
diff --git a/Core/include/Acts/Surfaces/LineBounds.hpp b/Core/include/Acts/Surfaces/LineBounds.hpp
index 8b82880decc45865fb5609adf2f273111e440000..3636317887819ff45178575dcb5cd98db9f3db71 100644
--- a/Core/include/Acts/Surfaces/LineBounds.hpp
+++ b/Core/include/Acts/Surfaces/LineBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// LineBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include "Acts/Surfaces/SurfaceBounds.hpp"
 #include "Acts/Utilities/Definitions.hpp"
@@ -33,12 +29,16 @@ class LineBounds : public SurfaceBounds {
   /// @param halez is the half length in z, defualt = 0.
   LineBounds(double radius = 0., double halez = 0.);
 
-  ~LineBounds() override;
+  /// Defaulted destructor
+  ~LineBounds() override = default;
 
+  /// Virtual constructor
   LineBounds* clone() const final;
 
+  /// Type enumeration
   BoundsType type() const final;
 
+  /// Return the intrinsic values
   std::vector<TDD_real_t> valueStore() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
diff --git a/Core/include/Acts/Surfaces/LineSurface.hpp b/Core/include/Acts/Surfaces/LineSurface.hpp
index fbe2a9c94f013bae014e4acdd20cd1e0ea2fcace..514d98ac052f4a130b5b1b203472a52484b5c6a9 100644
--- a/Core/include/Acts/Surfaces/LineSurface.hpp
+++ b/Core/include/Acts/Surfaces/LineSurface.hpp
@@ -1,18 +1,15 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-/////////////////////////////////////////////////////////////////
-// LineSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryStatics.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/InfiniteBounds.hpp"
 #include "Acts/Surfaces/LineBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
diff --git a/Core/include/Acts/Surfaces/PerigeeSurface.hpp b/Core/include/Acts/Surfaces/PerigeeSurface.hpp
index 70100e57de1e2309a773da8c2f0064452511dce4..90bca3c7fd681d682a8d1f054e694cabfb9728cd 100644
--- a/Core/include/Acts/Surfaces/PerigeeSurface.hpp
+++ b/Core/include/Acts/Surfaces/PerigeeSurface.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-/////////////////////////////////////////////////////////////////
-// PerigeeSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include "Acts/Geometry/GeometryStatics.hpp"
 #include "Acts/Surfaces/InfiniteBounds.hpp"
@@ -86,6 +82,15 @@ class PerigeeSurface : public LineSurface {
   std::ostream& toStream(const GeometryContext& gctx,
                          std::ostream& sl) const final;
 
+  /// Return a Polyhedron for the surfaces
+  ///
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param lseg is ignored for a perigee @note ignored
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                      size_t /*ignored*/) const final;
+
  private:
   /// Clone method implementation
   ///
diff --git a/Core/include/Acts/Surfaces/PlanarBounds.hpp b/Core/include/Acts/Surfaces/PlanarBounds.hpp
index b510044db90dab511ffb999ebba3bc9d9ea2ae6d..8fafc816ce46e672d29453dd5ea45bd3106af9f1 100644
--- a/Core/include/Acts/Surfaces/PlanarBounds.hpp
+++ b/Core/include/Acts/Surfaces/PlanarBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// PlanarBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <vector>
 
@@ -26,10 +22,20 @@ class RectangleBounds;
 ///
 class PlanarBounds : public SurfaceBounds {
  public:
-  /// Return the vertices - or, the points of the extremas
-  virtual std::vector<Vector2D> vertices() const = 0;
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note that the extremas are given, which may slightly alter the
+  /// number of segments returned
+  ///
+  /// @return vector for vertices in 2D
+  virtual std::vector<Vector2D> vertices(unsigned int lseg = 1) const = 0;
 
-  // Bounding box parameters
+  /// Bounding box parameters
+  ///
+  /// @return rectangle bounds for a bounding box
   virtual const RectangleBounds& boundingBox() const = 0;
 };
 
diff --git a/Core/include/Acts/Surfaces/PlaneSurface.hpp b/Core/include/Acts/Surfaces/PlaneSurface.hpp
index e8069f8b8c7e52471e45d6d686772a027dce6801..d13584beeaa61689606af29cf025934cdd4d64a6 100644
--- a/Core/include/Acts/Surfaces/PlaneSurface.hpp
+++ b/Core/include/Acts/Surfaces/PlaneSurface.hpp
@@ -1,20 +1,17 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2019 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// PlaneSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include <limits>
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryStatics.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/InfiniteBounds.hpp"
 #include "Acts/Surfaces/PlanarBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
@@ -191,6 +188,17 @@ class PlaneSurface : public Surface {
       const Vector3D& direction,
       const BoundaryCheck& bcheck = false) const final;
 
+  /// Return a Polyhedron for the surfaces
+  ///
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param lseg Number of segments along curved lines, it represents
+  /// the full 2*M_PI coverange, if lseg is set to 1 only the extrema
+  /// are given
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                      size_t lseg) const override;
+
   /// Return properly formatted class name for screen output
   std::string name() const override;
 
diff --git a/Core/include/Acts/Surfaces/PolyhedronRepresentation.hpp b/Core/include/Acts/Surfaces/PolyhedronRepresentation.hpp
deleted file mode 100644
index 607ad7f573c8c94bec01158ea8c77139c3a9b647..0000000000000000000000000000000000000000
--- a/Core/include/Acts/Surfaces/PolyhedronRepresentation.hpp
+++ /dev/null
@@ -1,60 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2018 CERN for the benefit of the Acts project
-//
-// This Source Code Form is subject to the terms of the Mozilla Public
-// License, v. 2.0. If a copy of the MPL was not distributed with this
-// file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#pragma once
-
-#include <vector>
-
-#include "Acts/Utilities/Definitions.hpp"
-
-namespace Acts {
-
-/// @class PolyhedronRepresentation
-///
-/// Struct which contains a cartesian approximation for any surface type.
-/// It contains a list of cartesian vertices in the global frame, and
-/// additionally
-/// a list of lists of indices which indicate which vertices form a face.
-/// Each entry in @c faces is a face, which is in turn a list of vertices
-/// that need to be connected to form a face.
-/// This allows the @c objString method to produce a ready-to-go obj output.
-struct PolyhedronRepresentation {
-  /// Default constructor from a vector of vertices and a vector of faces
-  /// @param verticesIn The 3D global vertices that make up the object
-  /// @param facesIn List of lists of indices for faces.
-  /// @note This creates copies of the input vectors
-  PolyhedronRepresentation(const std::vector<Vector3D>& verticesIn,
-                           const std::vector<std::vector<size_t>>& facesIn)
-      : vertices(verticesIn), faces(facesIn) {}
-
-  /// list of 3D vertices as vectors
-  std::vector<Vector3D> vertices;
-
-  /// list of faces connecting the vertices.
-  /// each face is a list of vertices v
-  /// corresponding to the vertex vector above
-  std::vector<std::vector<size_t>> faces;
-
-  /// Method to return an obj string.
-  /// @param vtxOffset Optional obj vertex enumeration offset
-  /// @note Vertices in obj are enumerated globally. The offset is required
-  ///       to allow for multiple output objects in one obj file.
-  std::string objString(size_t vtxOffset = 0) const;
-
-  template <typename helper_t>
-  void draw(helper_t& helper) const {
-    for (const auto& face : faces) {
-      std::vector<Vector3D> face_vtx;
-      for (size_t i : face) {
-        face_vtx.push_back(vertices[i]);
-      }
-      helper.face(face_vtx);
-    }
-  }
-};
-}  // namespace Acts
diff --git a/Core/include/Acts/Surfaces/RadialBounds.hpp b/Core/include/Acts/Surfaces/RadialBounds.hpp
index 157e4bb41bbf63de151360efc7db275634e8f93f..c5192452ab859ee447a67e3206c2811b4a9b4f1c 100644
--- a/Core/include/Acts/Surfaces/RadialBounds.hpp
+++ b/Core/include/Acts/Surfaces/RadialBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// RadialBounds.h, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
 
@@ -59,12 +55,16 @@ class RadialBounds : public DiscBounds {
   /// coverage)
   RadialBounds(double minrad, double maxrad, double avephi, double hphisec);
 
-  ~RadialBounds() override;
+  /// Defaulted Destructor
+  ~RadialBounds() override = default;
 
+  /// Virtual constructor
   RadialBounds* clone() const final;
 
+  /// Return the type enumerator
   SurfaceBounds::BoundsType type() const final;
 
+  /// Value store to be written out
   std::vector<TDD_real_t> valueStore() const final;
 
   /// For disc surfaces the local position in (r,phi) is checked
@@ -122,6 +122,18 @@ class RadialBounds : public DiscBounds {
   ///
   /// @param lposition The local position in polar coordinates
   Vector2D shifted(const Vector2D& lposition) const;
+
+  /// This method returns the xy coordinates of vertices along
+  /// the radial bounds
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note that the extremas are given, which may slightly alter the
+  /// number of segments returned
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg) const;
 };
 
 inline double RadialBounds::rMin() const {
diff --git a/Core/include/Acts/Surfaces/RectangleBounds.hpp b/Core/include/Acts/Surfaces/RectangleBounds.hpp
index 3efe6174d88193203fac953939c8bf9dbc69532e..95b13be8757fe943e0de58389f92248d3b9aa037 100644
--- a/Core/include/Acts/Surfaces/RectangleBounds.hpp
+++ b/Core/include/Acts/Surfaces/RectangleBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// RectangleBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include "Acts/Surfaces/PlanarBounds.hpp"
 #include "Acts/Utilities/Definitions.hpp"
@@ -69,8 +65,15 @@ class RectangleBounds : public PlanarBounds {
   /// @return is a signed distance parameter
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /// Return the vertices - or, the points of the extremas
-  std::vector<Vector2D> vertices() const final;
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note the number of segements is ignored in this representation
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg = 1) const final;
 
   // Bounding box representation
   const RectangleBounds& boundingBox() const final;
diff --git a/Core/include/Acts/Surfaces/StrawSurface.hpp b/Core/include/Acts/Surfaces/StrawSurface.hpp
index cf55599bec38f8496bc194dc894bc2d72d701913..90da119bd7e328e7632160bd30eb7819db06f582 100644
--- a/Core/include/Acts/Surfaces/StrawSurface.hpp
+++ b/Core/include/Acts/Surfaces/StrawSurface.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-/////////////////////////////////////////////////////////////////
-// StrawSurface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include "Acts/Geometry/GeometryContext.hpp"
@@ -20,7 +16,7 @@
 namespace Acts {
 
 class DetectorElementBase;
-struct PolyhedronRepresentation;
+struct Polyhedron;
 
 ///  @class StrawSurface
 ///
@@ -97,12 +93,16 @@ class StrawSurface : public LineSurface {
   /// Return properly formatted class name for screen output */
   std::string name() const final;
 
-  /// Return a PolyhedronRepresentation for this object
+  /// Return a Polyhedron for the surfaces
+  ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param l0div Number of divisions along l0 (phi)
-  /// @param l1div Number of divisions along l1 (z)
-  virtual PolyhedronRepresentation polyhedronRepresentation(
-      const GeometryContext& gctx, size_t l0div = 10, size_t l1div = 1) const;
+  /// @param lseg Number of segments along curved lines, it represents
+  /// the full 2*M_PI coverange, if lseg is set to 1 only the extrema
+  /// are given @note if lseg is set to 1 then only the straw is created
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                      size_t lseg) const final;
 
  private:
   /// Clone method implementation
diff --git a/Core/include/Acts/Surfaces/Surface.hpp b/Core/include/Acts/Surfaces/Surface.hpp
index 9d67852edf38c71d8f6dee9c481f35901c51d91c..6287c0cf04557951d47421299df441ec6aa9bac1 100644
--- a/Core/include/Acts/Surfaces/Surface.hpp
+++ b/Core/include/Acts/Surfaces/Surface.hpp
@@ -1,21 +1,18 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// Surface.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 
 #include "Acts/Geometry/DetectorElementBase.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryObject.hpp"
 #include "Acts/Geometry/GeometryStatics.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/BoundaryCheck.hpp"
 #include "Acts/Surfaces/SurfaceBounds.hpp"
 #include "Acts/Utilities/BinnedArray.hpp"
@@ -453,6 +450,21 @@ class Surface : public virtual GeometryObject,
   /// Return properly formatted class name
   virtual std::string name() const = 0;
 
+  /// Return a Polyhedron for this object
+  ///
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param lseg Number of segments along curved lines, if the lseg
+  /// is set to one, only the corners and the extrema are given,
+  /// otherwise it represents the number of segments for a full 2*M_PI
+  /// circle and is scaled to the relevant sector
+  ///
+  /// @note An internal surface transform can invalidate the extrema
+  /// in the transformed space
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  virtual Polyhedron polyhedronRepresentation(const GeometryContext& gctx,
+                                              size_t lseg) const = 0;
+
  protected:
   /// Transform3D definition that positions
   /// (translation, rotation) the surface in global space
diff --git a/Core/include/Acts/Surfaces/SurfaceArray.hpp b/Core/include/Acts/Surfaces/SurfaceArray.hpp
index 5e1ac370a9cb0319774c410a5667e6510506e12b..0b8c72506800add0e70aa800080f895e0f9f1cbb 100644
--- a/Core/include/Acts/Surfaces/SurfaceArray.hpp
+++ b/Core/include/Acts/Surfaces/SurfaceArray.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2017-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
diff --git a/Core/include/Acts/Surfaces/SurfaceBounds.hpp b/Core/include/Acts/Surfaces/SurfaceBounds.hpp
index f3e5805b410d16a765e47b781ac1b350ba4c30c2..7a58c5aab8ad032bd650d9ef1b2a68c770a91ad4 100644
--- a/Core/include/Acts/Surfaces/SurfaceBounds.hpp
+++ b/Core/include/Acts/Surfaces/SurfaceBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// SurfaceBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <ostream>
 
diff --git a/Core/include/Acts/Surfaces/TrapezoidBounds.hpp b/Core/include/Acts/Surfaces/TrapezoidBounds.hpp
index d948f616f6005bb62fc5018c9c900c4fa9b99ee4..42e87c1ece3e0462387502f9402ff522b53a4b4a 100644
--- a/Core/include/Acts/Surfaces/TrapezoidBounds.hpp
+++ b/Core/include/Acts/Surfaces/TrapezoidBounds.hpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// TrapezoidBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
 #pragma once
 #include <cmath>
 
@@ -24,7 +20,7 @@ namespace Acts {
 ///
 /// Bounds for a trapezoidal, planar Surface.
 ///
-/// @image html TrapezoidalBounds.gif
+/// @image html TrapezoidBounds.gif
 ///
 /// @todo can be speed optimized by calculating kappa/delta and caching it
 
@@ -109,8 +105,15 @@ class TrapezoidBounds : public PlanarBounds {
   /// @return is a signed distance parameter
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /// Return the vertices - or, the points of the extremas
-  std::vector<Vector2D> vertices() const final;
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note the number of segements is ignored in this representation
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg = 1) const final;
 
   // Bounding box representation
   const RectangleBounds& boundingBox() const final;
diff --git a/Core/include/Acts/Surfaces/TriangleBounds.hpp b/Core/include/Acts/Surfaces/TriangleBounds.hpp
index 000f4c0a610fca92527b4a76a51bad053028dd97..204ea0d8697fd770a604fd785142aab54dc762bc 100644
--- a/Core/include/Acts/Surfaces/TriangleBounds.hpp
+++ b/Core/include/Acts/Surfaces/TriangleBounds.hpp
@@ -70,8 +70,15 @@ class TriangleBounds : public PlanarBounds {
   /// @return is a signed distance parameter
   double distanceToBoundary(const Vector2D& lposition) const final;
 
-  /// This method returns the coordinates of vertices
-  std::vector<Vector2D> vertices() const final;
+  /// Return the vertices
+  ///
+  /// @param lseg the number of segments used to approximate
+  /// and eventually curved line
+  ///
+  /// @note the number of segements is ignored in this representation
+  ///
+  /// @return vector for vertices in 2D
+  std::vector<Vector2D> vertices(unsigned int lseg = 1) const final;
 
   // Bounding box representation
   const RectangleBounds& boundingBox() const final;
diff --git a/Core/include/Acts/Surfaces/detail/FacesHelper.hpp b/Core/include/Acts/Surfaces/detail/FacesHelper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..47f1932fc6a80a8f6e5004c7ad630b766654652c
--- /dev/null
+++ b/Core/include/Acts/Surfaces/detail/FacesHelper.hpp
@@ -0,0 +1,87 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Geometry/Polyhedron.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+
+#include <numeric>
+#include <utility>
+#include <vector>
+
+namespace Acts {
+
+namespace detail {
+
+/// @brief Helper for writing out faces for polyhedron representation
+struct FacesHelper {
+  using FaceVector = std::vector<Polyhedron::Face>;
+
+  /// @brief This method words for all convex type surface setups
+  /// It includes:
+  ///
+  /// Rectangle / Triangle / Polygon
+  /// Full disc, ellipse, half ring
+  /// @param vertices The vector of vertices
+  /// @param centerLast Boolean indicator if the center is given for
+  /// a better triangulation method as last element of the vector
+  static std::pair<FaceVector, FaceVector> convexFaceMesh(
+      const std::vector<Vector3D>& vertices, bool centerLast = false) {
+    FaceVector faces;
+    FaceVector triangularMesh;
+    // Write the face
+    unsigned int offset = centerLast ? 1 : 0;
+    std::vector<size_t> face(vertices.size() - offset);
+    std::iota(face.begin(), face.end(), 0);
+    faces.push_back(face);
+    /// Triangular mesh construction
+    unsigned int anker = centerLast ? vertices.size() - 1 : 0;
+    for (unsigned int it = 2 - offset; it < vertices.size() - offset; ++it) {
+      triangularMesh.push_back({anker, it - 1, it});
+    }
+    // Close for centered reference point
+    if (centerLast) {
+      triangularMesh.push_back({anker, vertices.size() - 2, 0});
+    }
+    return {faces, triangularMesh};
+  }
+
+  /// @brief This method works for all concentric type surface setups
+  /// It includes :
+  ///
+  /// - Cylinder (concentric bows on each side, separated by z)
+  /// - Cut-off cone
+  ///
+  /// The single requirement is that the #vertices are equal and the input
+  /// vector is splittable in half into the two bows.
+  ///
+  /// @param vertices The vector of vertices
+  /// @param fullTwoPi The indicator if the concentric face is closed
+  static std::pair<FaceVector, FaceVector> cylindricalFaceMesh(
+      const std::vector<Vector3D>& vertices, bool fullTwoPi = true) {
+    FaceVector faces;
+    FaceVector triangularMesh;
+    size_t nqfaces = 0.5 * vertices.size();
+    size_t reduce = (not fullTwoPi) ? 1 : 0;
+    for (size_t iface = 0; iface < nqfaces - reduce; ++iface) {
+      size_t p2 = (iface + 1 == nqfaces) ? 0 : iface + 1;
+      std::vector<size_t> face = {iface, p2, p2 + nqfaces, nqfaces + iface};
+      faces.push_back(face);
+      std::vector<size_t> triA = {iface, p2, p2 + nqfaces};
+      triangularMesh.push_back(triA);
+      std::vector<size_t> triB = {p2 + nqfaces, nqfaces + iface, iface};
+      triangularMesh.push_back(triB);
+    }
+    return {faces, triangularMesh};
+  }
+};
+
+}  // namespace detail
+
+}  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp b/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c4e55f7e799dc751ae383039e9453b420e2caba2
--- /dev/null
+++ b/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp
@@ -0,0 +1,247 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <utility>
+#include <vector>
+#include "Acts/Utilities/Definitions.hpp"
+
+namespace Acts {
+
+namespace detail {
+
+/// @brief Helper method for set of vertices for polyhedrons
+/// drawing and inside/outside checks
+namespace VerticesHelper {
+
+/// A method that inserts the cartesian extrema points and segments
+/// a curved segment into sub segments
+///
+/// @param phiMin the minimum Phi of the bounds object
+/// @param phiMax the maximum Phi of the bounds object
+/// @param phiRef is a vector of reference phi values to be included as well
+/// @param phiTolerance is the tolerance for reference phi insertion
+///
+/// @return a vector
+static std::vector<double> phiSegments(double phiMin = -M_PI,
+                                       double phiMax = M_PI,
+                                       std::vector<double> phiRefs = {},
+                                       double phiTolerance = 1e-6) {
+  // This is to ensure that the extrema are built regardless of number
+  // of segments
+  std::vector<double> phiSegments;
+  std::vector<double> quarters = {-M_PI, -0.5 * M_PI, 0., 0.5 * M_PI, M_PI};
+  // It does not cover the full azimuth
+  if (phiMin != -M_PI or phiMax != M_PI) {
+    phiSegments.push_back(phiMin);
+    for (unsigned int iq = 1; iq < 4; ++iq) {
+      if (phiMin < quarters[iq] and phiMax > quarters[iq]) {
+        phiSegments.push_back(quarters[iq]);
+      }
+    }
+    phiSegments.push_back(phiMax);
+  } else {
+    phiSegments = quarters;
+  }
+  // Insert the reference phis if
+  if (not phiRefs.empty()) {
+    for (const auto& phiRef : phiRefs) {
+      // Trying to find the right patch
+      auto match = std::find_if(
+          phiSegments.begin(), phiSegments.end(), [&](double phiSeg) {
+            return std::abs(phiSeg - phiRef) < phiTolerance;
+          });
+      if (match == phiSegments.end()) {
+        phiSegments.push_back(phiRef);
+      }
+    }
+    std::sort(phiSegments.begin(), phiSegments.end());
+  }
+  return phiSegments;
+}
+
+/// Helper method to create a regular 2 or 3 D segment
+///  between two phi values
+///
+/// @tparam vertex_t Type of vertex to be applied
+/// @tparam transform_t Optional transform
+///
+/// @param vertices [in,out] The 3D vertices to be filled
+/// @param rxy The radius description if first +/= second: ellipse
+/// @param phi1 The first phi value
+/// @param phi2 The second phi value
+/// @param lseg The number of segments for full 2*PI
+/// @param addon The additional segments to be built
+/// @param offset The out of plane offset position of the bow
+/// @param transform The transform applied (optional)
+template <typename vertex_t, typename transform_t>
+void createSegment(std::vector<vertex_t>& vertices,
+                   std::pair<double, double> rxy, double phi1, double phi2,
+                   unsigned int lseg, int addon = 0,
+                   const vertex_t& offset = vertex_t::Zero(),
+                   const transform_t& transform = transform_t::Identity()) {
+  // Calculate the number of segments - 1 is the minimum
+  unsigned int segs = std::abs(phi2 - phi1) / (2 * M_PI) * lseg;
+  segs = segs > 0 ? segs : 1;
+  double phistep = (phi2 - phi1) / segs;
+  // Create the segments
+  for (unsigned int iphi = 0; iphi < segs + addon; ++iphi) {
+    double phi = phi1 + iphi * phistep;
+    vertex_t vertex = vertex_t::Zero();
+    vertex(0) = rxy.first * std::cos(phi);
+    vertex(1) = rxy.second * std::sin(phi);
+    vertex = vertex + offset;
+    vertices.push_back(transform * vertex);
+  }
+}
+
+/// Vertices on an ellipse-like bound object
+///
+/// @param innerRx The radius of the inner ellipse (in x), 0 if sector
+/// @param innerRy The radius of the inner ellipse (in y), 0 if sector
+/// @param outerRx The radius of the outer ellipse (in x)
+/// @param outerRy The radius of the outer ellipse (in y)
+/// @param avgPhi The phi direction of the center if sector
+/// @param halfPhi The half phi sector if sector
+/// @param lseg The number of segments for for a full 2*pi segment
+///
+/// @return a vector of 2d-vectors
+static std::vector<Vector2D> ellispoidVertices(double innerRx, double innerRy,
+                                               double outerRx, double outerRy,
+                                               double avgPhi = 0.,
+                                               double halfPhi = M_PI,
+                                               unsigned int lseg = 1) {
+  // List of vertices counter-clockwise starting at smallest phi w.r.t center,
+  // for both inner/outer ring/segment
+  std::vector<Acts::Vector2D> rvertices;  // return vertices
+  std::vector<Acts::Vector2D> ivertices;  // inner vertices
+  std::vector<Acts::Vector2D> overtices;  // outer verices
+
+  bool innerExists = (innerRx > 0. and innerRy > 0.);
+  bool closed = std::abs(halfPhi - M_PI) < s_onSurfaceTolerance;
+
+  // Get the phi segments from the helper method
+  auto phiSegs = detail::VerticesHelper::phiSegments(
+      avgPhi - halfPhi, avgPhi + halfPhi, {avgPhi});
+
+  // The inner (if exists) and outer bow
+  for (unsigned int iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
+    int addon = (iseg == phiSegs.size() - 2 and not closed) ? 1 : 0;
+    if (innerExists) {
+      detail::VerticesHelper::createSegment<Vector2D, Eigen::Affine2d>(
+          ivertices, {innerRx, innerRy}, phiSegs[iseg], phiSegs[iseg + 1], lseg,
+          addon);
+    }
+    detail::VerticesHelper::createSegment<Vector2D, Eigen::Affine2d>(
+        overtices, {outerRx, outerRy}, phiSegs[iseg], phiSegs[iseg + 1], lseg,
+        addon);
+  }
+
+  // We want to keep the same counter-clockwise orientation for displaying
+  if (not innerExists) {
+    if (not closed) {
+      // Add the center case we have a sector
+      rvertices.push_back(Vector2D(0., 0.));
+    }
+    rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
+  } else if (not closed) {
+    rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
+    rvertices.insert(rvertices.end(), ivertices.rbegin(), ivertices.rend());
+  } else {
+    rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
+    rvertices.insert(rvertices.end(), ivertices.begin(), ivertices.end());
+  }
+  return rvertices;
+}
+
+/// Vertices on an disc/wheel-like bound object
+///
+/// @param innerR The radius of the inner circle (sector)
+/// @param outerR The radius of the outer circle (sector)
+/// @param avgPhi The phi direction of the center if sector
+/// @param halfPhi The half phi sector if sector
+/// @param lseg The number of segments for for a full 2*pi segment
+///
+/// @return a vector of 2d-vectors
+static std::vector<Vector2D> circularVertices(double innerR, double outerR,
+                                              double avgPhi = 0.,
+                                              double halfPhi = M_PI,
+                                              unsigned int lseg = 1) {
+  return ellispoidVertices(innerR, innerR, outerR, outerR, avgPhi, halfPhi,
+                           lseg);
+}
+
+/// Check if the point is inside the polygon w/o any tolerances
+///
+/// @tparam vertex_container_t is an iterable container
+///
+/// @param point is the Vector2DType to check
+/// @param vertices Forward iterable container of convex polygon vertices.
+///                 Calling `std::begin`/ `std::end` on the container must
+///                 return an iterator where `*it` must be convertible to
+///                 an `Vector2DType`.
+/// @return bool for inside/outside
+template <typename vertex_t, typename vertex_container_t>
+bool isInsidePolygon(const vertex_t& point,
+                     const vertex_container_t& vertices) {
+  // when we move along the edges of a convex polygon, a point on the inside of
+  // the polygon will always appear on the same side of each edge.
+  // a point on the outside will switch sides at least once.
+
+  // returns which side of the connecting line between `ll0` and `ll1` the point
+  // `p` is on. computes the sign of the z-component of the cross-product
+  // between the line normal vector and the vector from `ll0` to `p`.
+  auto lineSide = [&](auto&& ll0, auto&& ll1) {
+    auto normal = ll1 - ll0;
+    auto delta = point - ll0;
+    return std::signbit((normal[0] * delta[1]) - (normal[1] * delta[0]));
+  };
+
+  auto iv = std::begin(vertices);
+  auto l0 = *iv;
+  auto l1 = *(++iv);
+  // use vertex0 to vertex1 to define reference sign and compare w/ all edges
+  auto reference = lineSide(l0, l1);
+  for (++iv; iv != std::end(vertices); ++iv) {
+    l0 = l1;
+    l1 = *iv;
+    if (lineSide(l0, l1) != reference) {
+      return false;
+    }
+  }
+  // manual check for last edge from last vertex back to the first vertex
+  if (lineSide(l1, *std::begin(vertices)) != reference) {
+    return false;
+  }
+  // point was always on the same side. point must be inside.
+  return true;
+}
+
+/// Check if the point is inside the rectangle
+///
+/// @tparam vertex_t is vector with [0],[1] access
+///
+/// @param point is the Vector2DType to check
+/// @param vertices Forward iterable container of convex polygon vertices.
+///                 Calling `std::begin`/ `std::end` on the container must
+///                 return an iterator where `*it` must be convertible to
+///                 an `Vector2DType`.
+/// @return bool for inside/outside
+template <typename vertex_t>
+bool isInsideRectangle(const vertex_t& point, const vertex_t& lowerLeft,
+                       const vertex_t& upperRight) {
+  return (lowerLeft[0] <= point[0]) && (point[0] < upperRight[0]) &&
+         (lowerLeft[1] <= point[1]) && (point[1] < upperRight[1]);
+}
+
+};  // namespace VerticesHelper
+
+}  // namespace detail
+
+}  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Utilities/BinningType.hpp b/Core/include/Acts/Utilities/BinningType.hpp
index 3014f05e7ab39791ada53509a6e2dae88794a09b..42923c0e4138090d3d99d314a8c5db440ce5804a 100644
--- a/Core/include/Acts/Utilities/BinningType.hpp
+++ b/Core/include/Acts/Utilities/BinningType.hpp
@@ -36,16 +36,17 @@ enum BinningType { equidistant, arbitrary };
 enum BinningOption { open, closed };
 
 /// @enum BinningValue how to take the global / local position
-enum BinningValue {
-  binX,
-  binY,
-  binZ,
-  binR,
-  binPhi,
-  binRPhi,
-  binH,
-  binEta,
-  binMag
+enum BinningValue : int {
+  binX = 0,
+  binY = 1,
+  binZ = 2,
+  binR = 3,
+  binPhi = 4,
+  binRPhi = 5,
+  binH = 6,
+  binEta = 7,
+  binMag = 8,
+  binValues = 9
 };
 
 /// @brief screen output option
diff --git a/Core/include/Acts/Utilities/Helpers.hpp b/Core/include/Acts/Utilities/Helpers.hpp
index 71e310c7a05f18492012abb9490b1a1f8beba2c5..359b2ff1bb09114d8f16a130ed62652332150632 100644
--- a/Core/include/Acts/Utilities/Helpers.hpp
+++ b/Core/include/Acts/Utilities/Helpers.hpp
@@ -21,6 +21,7 @@
 #include <string>
 #include <vector>
 
+#include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/TypeTraits.hpp"
 #include "Acts/Utilities/detail/DefaultParameterDefinitions.hpp"
@@ -162,6 +163,34 @@ double eta(const Eigen::MatrixBase<Derived>& v) noexcept {
   return std::atanh(v[2] / v.norm());
 }
 
+/// Helper method to cast out the binning value from a 3D Vector
+///
+/// For this method a 3D vector is required to guarantee all potential
+/// binning values
+///
+static double cast(const Vector3D& position, BinningValue bval) {
+  if (bval < 3)
+    return position[bval];
+  switch (bval) {
+    case binR:
+      return perp(position);
+      break;
+    case binPhi:
+      return phi(position);
+      break;
+    case binH:
+      return theta(position);
+      break;
+    case binEta:
+      return eta(position);
+      break;
+    case binMag:
+      return position.norm();
+      break;
+  }
+  return 0.;
+}
+
 /// @brief Calculates column-wise cross products of a matrix and a vector and
 /// stores the result column-wise in a matrix.
 ///
diff --git a/Core/include/Acts/Utilities/IVisualization.hpp b/Core/include/Acts/Utilities/IVisualization.hpp
index 68b0bc3a8162b6370d5738ed291944723f6dc628..9f840c5a87cbc929423522c6ad2460cf8bf24464 100644
--- a/Core/include/Acts/Utilities/IVisualization.hpp
+++ b/Core/include/Acts/Utilities/IVisualization.hpp
@@ -13,75 +13,80 @@
 
 namespace Acts {
 
-/**
- * Partially abstract base class which provides an interface to visualization
- * helper classes. It provides a number of methods that all the helpers need to
- * conform to. It also provides a `color_type` typedef, but not all of the
- * helpers actually support that.
- */
+/// Partially abstract base class which provides an interface to visualization
+/// helper classes. It provides a number of methods that all the helpers need to
+/// conform to. It also provides a `color_type` typedef, but not all of the
+/// helpers actually support that.
+///
 class IVisualization {
  public:
-  /**
-   * The color typedef. It's an array of three numbers [0, 255] indicating RGB
-   * color values.
-   */
+  /// The color typedef. It's an array of three numbers [0, 255] indicating RGB
+  /// color values.
+  ///
   using color_type = std::array<int, 3>;
 
-  /**
-   * Draw a vertex at a given location and a color.
-   * @param vtx The vertex position
-   * @param color The color
-   */
+  /// The face type
+  using face_type = std::vector<size_t>;
+
+  /// Draw a vertex at a given location and a color.
+  /// @param vtx The vertex position
+  /// @param color The color
+  ///
   virtual void vertex(const Vector3D& vtx,
                       color_type color = {120, 120, 120}) = 0;
 
-  /**
-   * Draw a face that connects a list of vertices.
-   * @note Depending on the helper implementation, out of plane vertices might
-   * be handled differently.
-   * @param vtxs The vertices that make up the face
-   * @param color The color of the face
-   */
+  /// Draw a face that connects a list of vertices.
+  /// @note Depending on the helper implementation, out of plane vertices might
+  /// be handled differently.
+  /// @param vtxs The vertices that make up the face
+  /// @param color The color of the face
+  ///
   virtual void face(const std::vector<Vector3D>& vtxs,
                     color_type color = {120, 120, 120}) = 0;
 
-  /**
-   * Draw a line from a vertex to another
-   * @param a The start vertex
-   * @param b The end vertex
-   * @param color The color of the line
-   */
+  /// Draw a faces that connects a list of vertices - expert only
+  ///
+  /// @note Depending on the helper implementation, out of plane vertices might
+  /// be handled differently.
+  /// @param vtxs The vertices that make up the faceS
+  /// @param faces The face presectiotions (i.e. connecting vertices)
+  /// @param color The color of the face
+  ///
+  virtual void faces(const std::vector<Vector3D>& vtxs,
+                     const std::vector<face_type>& faces,
+                     color_type color = {120, 120, 120}) = 0;
+
+  /// Draw a line from a vertex to another
+  /// @param a The start vertex
+  /// @param b The end vertex
+  /// @param color The color of the line
+  ///
   virtual void line(const Vector3D& a, const Vector3D& b,
                     color_type color = {120, 120, 120}) = 0;
 
-  /**
-   * Write the content of the helper to an outstream.
-   * @param os The output stream
-   */
+  /// Write the content of the helper to an outstream.
+  /// @param os The output stream
+  ///
   virtual void write(std::ostream& os) const = 0;
 
-  /**
-   * Remove all contents of this helper
-   */
+  /// Remove all contents of this helper
+  ///
   virtual void clear() = 0;
 
-  /**
-   * Below are helper functions, which share the same interface as the ones
-   * above, but explicitly accept float values (instead of double), converts
-   * them and calls the above methods.
-   */
+  /// Below are helper functions, which share the same interface as the ones
+  /// above, but explicitly accept float values (instead of double), converts
+  /// them and calls the above methods.
+  ///
 
-  /**
-   * @copydoc Acts::IVisualization::vertex(const Vector3D&, color_type)
-   */
+  /// @copydoc Acts::IVisualization::vertex(const Vector3D&, color_type)
+  ///
   void vertex(const Vector3F& vtx, color_type color = {120, 120, 120}) {
     Vector3D vtxd = vtx.template cast<double>();
     vertex(vtxd, color);
   }
 
-  /**
-   * @copydoc Acts::IVisualization::face(std::vector<Vector3F>&, color_type)
-   */
+  /// @copydoc Acts::IVisualization::face(std::vector<Vector3F>&, color_type)
+  ///
   void face(const std::vector<Vector3F>& vtxs,
             color_type color = {120, 120, 120}) {
     std::vector<Vector3D> vtxsd;
@@ -90,10 +95,9 @@ class IVisualization {
     face(vtxsd, color);
   }
 
-  /**
-   * @copydoc Acts::IVisualization::line(const Vector3F&, const Vector3F&,
-   * color_type)
-   */
+  ///  @copydoc Acts::IVisualization::line(const Vector3F&, const Vector3F&,
+  /// color_type)
+  ///
   void line(const Vector3F& a, const Vector3F& b,
             color_type color = {120, 120, 120}) {
     Vector3D ad = a.template cast<double>();
diff --git a/Core/include/Acts/Utilities/ObjHelper.hpp b/Core/include/Acts/Utilities/ObjHelper.hpp
index c45abf648472720953eb13b9a1352ba8c7f8be7a..f26b940ae571c0482f9fa393eecb8f0e2d48dc6a 100644
--- a/Core/include/Acts/Utilities/ObjHelper.hpp
+++ b/Core/include/Acts/Utilities/ObjHelper.hpp
@@ -11,56 +11,47 @@
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/IVisualization.hpp"
 
+#include <utility>
+#include <vector>
+
 namespace Acts {
 
-/**
- * This helper produces output in the OBJ format. Note that colors are not
- * supported in this implementation.
- */
+/// This helper produces output in the OBJ format. Note that colors are not
+/// supported in this implementation.
+///
 template <typename T = double>
 class ObjHelper : public IVisualization {
  public:
   static_assert(std::is_same_v<T, double> || std::is_same_v<T, float>,
                 "Use either double or float");
 
-  /**
-   * Stored value type, should be double or float
-   */
+  /// Stored value type, should be double or float
   using value_type = T;
 
-  /**
-   * Type of a vertex based on the value type
-   */
+  /// Type of a vertex based on the value type
   using vertex_type = ActsVector<value_type, 3>;
 
-  /**
-   * Typedef of what a face is.
-   */
-  using face_type = std::vector<size_t>;
+  /// Type of a line
+  using line_type = std::pair<size_t, size_t>;
 
-  /**
-   * @copydoc Acts::IVisualization::vertex()
-   */
+  /// @copydoc Acts::IVisualization::vertex()
   void vertex(const Vector3D& vtx,
               IVisualization::color_type color = {0, 0, 0}) override {
     (void)color;  // suppress unused warning
     m_vertices.push_back(vtx.template cast<value_type>());
   }
 
-  /**
-   * @copydoc Acts::IVisualization::line()
-   * @note This function ist not implemented for the OBJ format.
-   */
-  void line(const Vector3D& /*a*/, const Vector3D& /*b*/,
+  /// @copydoc Acts::IVisualization::line()
+  void line(const Vector3D& a, const Vector3D& b,
             IVisualization::color_type color = {0, 0, 0}) override {
     (void)color;  // suppress unused warning
     // not implemented
-    throw std::runtime_error("Line not implemented for OBJ");
+    vertex(a);
+    vertex(b);
+    m_lines.push_back({m_vertices.size() - 2, m_vertices.size() - 1});
   }
 
-  /**
-   * @copydoc Acts::IVisualization::face()
-   */
+  /// @copydoc Acts::IVisualization::face()
   void face(const std::vector<Vector3D>& vtxs,
             IVisualization::color_type color = {0, 0, 0}) override {
     (void)color;  // suppress unused warning
@@ -73,14 +64,39 @@ class ObjHelper : public IVisualization {
     m_faces.push_back(std::move(idxs));
   }
 
-  /**
-   * @copydoc Acts::IVisualization::write()
-   */
+  /// @copydoc Acts::IVisualization::faces()
+  void faces(const std::vector<Vector3D>& vtxs,
+             const std::vector<face_type>& faces,
+             color_type color = {120, 120, 120}) {
+    // No faces given - call the face() method
+    if (faces.empty()) {
+      face(vtxs, color);
+    } else {
+      auto vtxoffs = m_vertices.size();
+      m_vertices.insert(m_vertices.end(), vtxs.begin(), vtxs.end());
+      for (const auto& face : faces) {
+        if (face.size() == 2) {
+          m_lines.push_back({face[0] + vtxoffs, face[2] + vtxoffs});
+        } else {
+          face_type rawFace = face;
+          std::transform(rawFace.begin(), rawFace.end(), rawFace.begin(),
+                         [&](size_t& iv) { return (iv + vtxoffs); });
+          m_faces.push_back(rawFace);
+        }
+      }
+    }
+  }
+
+  /// @copydoc Acts::IVisualization::write()
   void write(std::ostream& os) const override {
     for (const vertex_type& vtx : m_vertices) {
       os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
     }
 
+    for (const line_type& ln : m_lines) {
+      os << "l " << ln.first + 1 << " " << ln.second + 1 << "\n";
+    }
+
     for (const face_type& fc : m_faces) {
       os << "f";
       for (size_t i = 0; i < fc.size(); i++) {
@@ -90,16 +106,16 @@ class ObjHelper : public IVisualization {
     }
   }
 
-  /**
-   * @copydoc Acts::IVisualization::clear()
-   */
+  ///  @copydoc Acts::IVisualization::clear()
   void clear() override {
     m_vertices.clear();
     m_faces.clear();
+    m_lines.clear();
   }
 
  private:
   std::vector<vertex_type> m_vertices;
   std::vector<face_type> m_faces;
+  std::vector<line_type> m_lines;
 };
 }  // namespace Acts
diff --git a/Core/include/Acts/Utilities/PlyHelper.hpp b/Core/include/Acts/Utilities/PlyHelper.hpp
index 2b75bb687ccae82144c2eb2b6eb31da51974c717..39e62cad5acf29b48927fd8c07923383344f0d4c 100644
--- a/Core/include/Acts/Utilities/PlyHelper.hpp
+++ b/Core/include/Acts/Utilities/PlyHelper.hpp
@@ -14,37 +14,26 @@
 namespace Acts {
 
 template <typename T = double>
+
+/// @brief Helper to write out Ply visualization format
 class PlyHelper : public IVisualization {
  public:
   static_assert(std::is_same_v<T, double> || std::is_same_v<T, float>,
                 "Use either double or float");
 
-  /**
-   * Stored value type, should be double or float
-   */
+  /// Stored value type, should be double or float
   using value_type = T;
 
-  /**
-   * Type of a vertex based on the value type
-   */
+  /// Type of a vertex based on the value type
   using vertex_type = ActsVector<value_type, 3>;
 
-  /**
-   * Typedef of what a face is.
-   */
-  using face_type = std::vector<size_t>;
-
-  /**
-   * @copydoc Acts::IVisualization::vertex()
-   */
+  /// @copydoc Acts::IVisualization::vertex()
   void vertex(const Vector3D& vtx,
               IVisualization::color_type color = {120, 120, 120}) override {
     m_vertices.emplace_back(vtx.template cast<value_type>(), color);
   }
 
-  /**
-   * @copydoc Acts::IVisualization::line()
-   */
+  /// @copydoc Acts::IVisualization::line()
   void face(const std::vector<Vector3D>& vtxs,
             IVisualization::color_type color = {120, 120, 120}) override {
     face_type idxs;
@@ -56,9 +45,13 @@ class PlyHelper : public IVisualization {
     m_faces.push_back(std::move(idxs));
   }
 
-  /**
-   * @copydoc Acts::IVisualization::face()
-   */
+  /// @copydoc Acts::IVisualization::faces()
+  void faces(const std::vector<Vector3D>& vtxs, const std::vector<face_type>&,
+             color_type color = {120, 120, 120}) {
+    face(vtxs, color);
+  }
+
+  /// @copydoc Acts::IVisualization::face()
   void line(const Vector3D& a, const Vector3D& b,
             IVisualization::color_type color = {120, 120, 120}) override {
     vertex(a, color);
@@ -68,9 +61,7 @@ class PlyHelper : public IVisualization {
     m_edges.emplace_back(std::make_pair(std::make_pair(idx_a, idx_b), color));
   }
 
-  /**
-   * @copydoc Acts::IVisualization::write()
-   */
+  /// @copydoc Acts::IVisualization::write()
   void write(std::ostream& os) const override {
     os << "ply\n";
     os << "format ascii 1.0\n";
@@ -116,9 +107,7 @@ class PlyHelper : public IVisualization {
     }
   }
 
-  /**
-   * @copydoc Acts::IVisualization::clear()
-   */
+  /// @copydoc Acts::IVisualization::clear()
   void clear() override {
     m_vertices.clear();
     m_faces.clear();
diff --git a/Core/src/Geometry/CMakeLists.txt b/Core/src/Geometry/CMakeLists.txt
index 16654ca692e676036315618c7087dd5ca170df6d..f4669b6d20d9680efcd5c20c5d3fe1dafd3bb685 100644
--- a/Core/src/Geometry/CMakeLists.txt
+++ b/Core/src/Geometry/CMakeLists.txt
@@ -11,6 +11,7 @@ target_sources_local(
     CylinderVolumeBounds.cpp
     CylinderVolumeBuilder.cpp
     CylinderVolumeHelper.cpp
+    Extent.cpp
     DiscLayer.cpp
     DoubleTrapezoidVolumeBounds.cpp
     GenericApproachDescriptor.cpp
@@ -23,6 +24,7 @@ target_sources_local(
     NavigationLayer.cpp
     PassiveLayerBuilder.cpp
     PlaneLayer.cpp
+    Polyhedron.cpp
     ProtoLayer.cpp
     SurfaceArrayCreator.cpp
     TrackingGeometry.cpp
diff --git a/Core/src/Geometry/CuboidVolumeBounds.cpp b/Core/src/Geometry/CuboidVolumeBounds.cpp
index bf657093390d7175638ece78045c693b37b83069..125b3f4d31d4c0c7ce237804339ba4a630979dda 100644
--- a/Core/src/Geometry/CuboidVolumeBounds.cpp
+++ b/Core/src/Geometry/CuboidVolumeBounds.cpp
@@ -6,10 +6,6 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CuboidVolumeBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Geometry/CuboidVolumeBounds.hpp"
 #include <cmath>
 #include <iostream>
@@ -57,15 +53,14 @@ Acts::CuboidVolumeBounds& Acts::CuboidVolumeBounds::operator=(
   return *this;
 }
 
-std::vector<std::shared_ptr<const Acts::Surface>>
-Acts::CuboidVolumeBounds::decomposeToSurfaces(
+Acts::SurfacePtrVector Acts::CuboidVolumeBounds::decomposeToSurfaces(
     const Transform3D* transformPtr) const {
   // the transform - apply when given
   Transform3D transform =
       (transformPtr == nullptr) ? Transform3D::Identity() : (*transformPtr);
   const Transform3D* tTransform = nullptr;
 
-  std::vector<std::shared_ptr<const Acts::Surface>> rSurfaces;
+  SurfacePtrVector rSurfaces;
   rSurfaces.reserve(6);
   // face surfaces xy -------------------------------------
   //   (1) - at negative local z
diff --git a/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp b/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
index bd9032232f5316d7dc1c73ca29082e9a69d6fb3a..e777b74ea06e960ac3f4a0fb0646b8c055982a4c 100644
--- a/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
@@ -145,18 +145,3 @@ std::ostream& Acts::CutoutCylinderVolumeBounds::toStream(
   sl << "dz1 = " << m_dz1 << " dz2 = " << m_dz2;
   return sl;
 }
-
-void Acts::CutoutCylinderVolumeBounds::draw(
-    IVisualization& helper, const Transform3D& transform) const {
-  std::vector<std::shared_ptr<const Acts::Surface>> surfaces =
-      decomposeToSurfaces(&transform);
-  for (const auto& srf : surfaces) {
-    auto cyl = dynamic_cast<const CylinderSurface*>(srf.get());
-    auto disc = dynamic_cast<const DiscSurface*>(srf.get());
-    if (cyl != nullptr) {
-      cyl->polyhedronRepresentation(50).draw(helper);
-    } else {
-      disc->polyhedronRepresentation(50).draw(helper);
-    }
-  }
-}
diff --git a/Core/src/Geometry/CylinderVolumeBounds.cpp b/Core/src/Geometry/CylinderVolumeBounds.cpp
index 58373138d83389cec99031073ec063712751ac4a..c90f71a4a4862b5296a7cc84f2ac9bffc9879c93 100644
--- a/Core/src/Geometry/CylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CylinderVolumeBounds.cpp
@@ -202,18 +202,3 @@ Acts::Volume::BoundingBox Acts::CylinderVolumeBounds::boundingBox(
   Volume::BoundingBox box{entity, vmin - envelope, vmax + envelope};
   return trf == nullptr ? box : box.transformed(*trf);
 }
-
-void Acts::CylinderVolumeBounds::draw(IVisualization& helper,
-                                      const Transform3D& transform) const {
-  std::vector<std::shared_ptr<const Acts::Surface>> surfaces =
-      decomposeToSurfaces(&transform);
-  for (const auto& srf : surfaces) {
-    auto cyl = dynamic_cast<const CylinderSurface*>(srf.get());
-    auto disc = dynamic_cast<const DiscSurface*>(srf.get());
-    if (cyl != nullptr) {
-      cyl->polyhedronRepresentation(50).draw(helper);
-    } else {
-      disc->polyhedronRepresentation(50).draw(helper);
-    }
-  }
-}
diff --git a/Core/src/Geometry/Extent.cpp b/Core/src/Geometry/Extent.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..11c9c9fa7ad887eb6086dffcae6eaea79655b1fa
--- /dev/null
+++ b/Core/src/Geometry/Extent.cpp
@@ -0,0 +1,26 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/Geometry/Extent.hpp"
+
+#include <ostream>
+
+std::ostream& Acts::Extent::toStream(std::ostream& sl) const {
+  sl << "Extent in space : " << std::endl;
+  for (size_t ib = 0; ib < static_cast<size_t>(binValues); ++ib) {
+    sl << "  - value :" << std::setw(10) << binningValueNames[ib]
+       << " | range = [" << ranges[ib].first << ", " << ranges[ib].second << "]"
+       << std::endl;
+  }
+  return sl;
+}
+
+// Overload of << operator for std::ostream for debug output
+std::ostream& Acts::operator<<(std::ostream& sl, const Extent& ext) {
+  return ext.toStream(sl);
+}
\ No newline at end of file
diff --git a/Core/src/Geometry/Polyhedron.cpp b/Core/src/Geometry/Polyhedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6be1ce44be5b8094458292c74a33ba9ced211786
--- /dev/null
+++ b/Core/src/Geometry/Polyhedron.cpp
@@ -0,0 +1,85 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/Geometry/Polyhedron.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
+#include "Acts/Utilities/BinningType.hpp"
+
+void Acts::Polyhedron::merge(const Acts::Polyhedron& other) {
+  size_t cvert = vertices.size();
+  vertices.insert(vertices.end(), other.vertices.begin(), other.vertices.end());
+  /// Add the new faces with offsets
+  auto join = [&](std::vector<Face>& existing,
+                  const std::vector<Face>& additional) -> void {
+    for (const auto& aface : additional) {
+      Face nface = aface;
+      std::transform(nface.begin(), nface.end(), nface.begin(),
+                     [&](size_t x) { return (x + cvert); });
+      existing.push_back(nface);
+    }
+  };
+  // For faces and triangular mesh
+  join(faces, other.faces);
+  join(triangularMesh, other.triangularMesh);
+}
+
+Acts::Extent Acts::Polyhedron::extent(const Transform3D& transform) const {
+  Extent extent;
+  for (const auto& vtx : vertices) {
+    extent.check(transform * vtx);
+  }
+  // Special checks for binR:
+  if (std::abs(extent.range(binZ)) < s_onSurfaceTolerance) {
+    // Check inclusion of origin (i.e. convex around origin)
+    Vector3D origin = transform * Vector3D(0., 0., extent.medium(binZ));
+    for (const auto& face : faces) {
+      std::vector<Vector3D> tface;
+      tface.reserve(face.size());
+      for (auto f : face) {
+        tface.push_back(transform * vertices[f]);
+      }
+      if (detail::VerticesHelper::isInsidePolygon(origin, tface)) {
+        extent.ranges[binR].first = 0.;
+        break;
+      }
+    }
+    // Check for radial extent
+    auto radialDistance = [&](const Vector3D& pos1,
+                              const Vector3D& pos2) -> double {
+      Vector2D p1(pos1.x(), pos1.y());
+      Vector2D p2(pos2.x(), pos2.y());
+
+      Vector2D O(0, 0);
+      Vector2D p1p2 = (p2 - p1);
+      double L = p1p2.norm();
+      Vector2D p1O = (O - p1);
+
+      // don't do division if L is very small
+      if (L < 1e-7) {
+        return std::numeric_limits<double>::max();
+      }
+      double f = p1p2.dot(p1O) / L;
+
+      // clamp to [0, |p1p2|]
+      f = std::min(L, std::max(0., f));
+
+      Vector2D closest = f * p1p2.normalized() + p1;
+      double dist = (closest - O).norm();
+
+      return dist;
+    };
+
+    for (size_t iv = 1; iv < vertices.size() + 1; ++iv) {
+      size_t fpoint = iv < vertices.size() ? iv : 0;
+      double testR = radialDistance(transform * vertices[fpoint],
+                                    transform * vertices[iv - 1]);
+      extent.ranges[binR].first = std::min(extent.ranges[binR].first, testR);
+    }
+  }
+  return extent;
+}
diff --git a/Core/src/Geometry/ProtoLayer.cpp b/Core/src/Geometry/ProtoLayer.cpp
index 643e5b47ba8953e444aa9895df03e15c132da12f..9694dfab4f48a0302a31358b4a81e6abe2909c76 100644
--- a/Core/src/Geometry/ProtoLayer.cpp
+++ b/Core/src/Geometry/ProtoLayer.cpp
@@ -9,11 +9,11 @@
 #include <cmath>
 
 #include <algorithm>
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Geometry/ProtoLayer.hpp"
 #include "Acts/Surfaces/AnnulusBounds.hpp"
 #include "Acts/Surfaces/CylinderBounds.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 
 using Acts::VectorHelpers::perp;
@@ -136,12 +136,12 @@ void ProtoLayer::measure(const GeometryContext& gctx,
       }
     } else if (cylSurface != nullptr) {
       // this is an explicit cast and if right now.
-      // It should work with all PolyhedronRepresentations
+      // It should work with all Polyhedrons
       // @TODO: Remove the cast and if as soon as ::polyhedronRepresentation()
       //        makes it into the Surface base class
       //        The envelopes might need special treatments though
 
-      PolyhedronRepresentation ph = cylSurface->polyhedronRepresentation(gctx);
+      Polyhedron ph = cylSurface->polyhedronRepresentation(gctx, 1);
       // evaluate at all vertices
       for (const auto& vtx : ph.vertices) {
         maxX = std::max(maxX, vtx.x());
diff --git a/Core/src/Surfaces/AnnulusBounds.cpp b/Core/src/Surfaces/AnnulusBounds.cpp
index 881d2ce9cc0806b65286b7f24e065befb2208b25..1ef13638fc0042d2928628016e47837e1c98ec7a 100644
--- a/Core/src/Surfaces/AnnulusBounds.cpp
+++ b/Core/src/Surfaces/AnnulusBounds.cpp
@@ -1,16 +1,13 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// AnnulusBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/AnnulusBounds.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/detail/periodic.hpp"
 
@@ -111,9 +108,37 @@ std::vector<Acts::Vector2D> Acts::AnnulusBounds::corners() const {
           rot * m_inLeftStripPC, rot * m_inRightStripPC};
 }
 
-std::vector<Acts::Vector2D> Acts::AnnulusBounds::vertices() const {
-  return {m_outRightStripXY, m_outLeftStripXY, m_inLeftStripXY,
-          m_inRightStripXY};
+std::vector<Acts::Vector2D> Acts::AnnulusBounds::vertices(
+    unsigned int lseg) const {
+  // List of vertices counter-clockwise starting with left inner
+  std::vector<Acts::Vector2D> rvertices;
+
+  double phiMinInner = VectorHelpers::phi(m_inLeftStripXY);
+  double phiMaxInner = VectorHelpers::phi(m_inRightStripXY);
+  double phiMinOuter = VectorHelpers::phi(m_outRightStripXY);
+  double phiMaxOuter = VectorHelpers::phi(m_outLeftStripXY);
+
+  std::vector<double> phisInner =
+      detail::VerticesHelper::phiSegments(phiMinInner, phiMaxInner);
+  std::vector<double> phisOuter =
+      detail::VerticesHelper::phiSegments(phiMinOuter, phiMaxOuter);
+
+  // Inner bow from phi_min -> phi_max
+  for (unsigned int iseg = 0; iseg < phisInner.size() - 1; ++iseg) {
+    int addon = (iseg == phisInner.size() - 2) ? 1 : 0;
+    detail::VerticesHelper::createSegment<Vector2D, Eigen::Affine2d>(
+        rvertices, {rMin(), rMin()}, phisInner[iseg], phisInner[iseg + 1], lseg,
+        addon);
+  }
+  // Upper bow from phi_min -> phi_max
+  for (unsigned int iseg = 0; iseg < phisOuter.size() - 1; ++iseg) {
+    int addon = (iseg == phisOuter.size() - 2) ? 1 : 0;
+    detail::VerticesHelper::createSegment<Vector2D, Eigen::Affine2d>(
+        rvertices, {rMax(), rMax()}, phisOuter[iseg], phisOuter[iseg + 1], lseg,
+        addon);
+  }
+
+  return rvertices;
 }
 
 bool Acts::AnnulusBounds::inside(const Vector2D& lposition, double tolR,
diff --git a/Core/src/Surfaces/CMakeLists.txt b/Core/src/Surfaces/CMakeLists.txt
index cf6fe1b02b87749a42d52201130f902d41e0a788..8a1a0822d93533b4053f2a2aba22d1253e9fe721 100644
--- a/Core/src/Surfaces/CMakeLists.txt
+++ b/Core/src/Surfaces/CMakeLists.txt
@@ -8,13 +8,12 @@ target_sources_local(
     CylinderSurface.cpp
     DiamondBounds.cpp
     DiscSurface.cpp
-    DiscTrapezoidalBounds.cpp
+    DiscTrapezoidBounds.cpp
     EllipseBounds.cpp
     LineBounds.cpp
     LineSurface.cpp
     PerigeeSurface.cpp
     PlaneSurface.cpp
-    PolyhedronRepresentation.cpp
     RadialBounds.cpp
     RectangleBounds.cpp
     StrawSurface.cpp
diff --git a/Core/src/Surfaces/ConeBounds.cpp b/Core/src/Surfaces/ConeBounds.cpp
index 5e2550a5a517cdd013a546820ebeefa828cc4273..2b523d7e260daca920a3c6878fcf02fe77b4314f 100644
--- a/Core/src/Surfaces/ConeBounds.cpp
+++ b/Core/src/Surfaces/ConeBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// ConeBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/ConeBounds.hpp"
 
 #include <cmath>
@@ -33,8 +29,6 @@ Acts::ConeBounds::ConeBounds(double alpha, double zmin, double zmax,
       m_avgPhi(detail::radian_sym(avphi)),
       m_halfPhi(std::abs(halfphi)) {}
 
-Acts::ConeBounds::~ConeBounds() = default;
-
 Acts::ConeBounds* Acts::ConeBounds::clone() const {
   return new ConeBounds(*this);
 }
diff --git a/Core/src/Surfaces/ConeSurface.cpp b/Core/src/Surfaces/ConeSurface.cpp
index a40235e6680bea391f2509e20c9e4c302dd4460b..3da3084424b34231afc97ae307abbc959e434dce 100644
--- a/Core/src/Surfaces/ConeSurface.cpp
+++ b/Core/src/Surfaces/ConeSurface.cpp
@@ -1,23 +1,20 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// ConeSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
-#include "Acts/Surfaces/ConeSurface.hpp"
-
 #include <cassert>
 #include <cmath>
 #include <iomanip>
 #include <iostream>
 #include <utility>
 
+#include "Acts/Surfaces/ConeSurface.hpp"
+#include "Acts/Surfaces/detail/FacesHelper.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 #include "Acts/Utilities/detail/RealQuadraticEquation.hpp"
 
@@ -196,3 +193,87 @@ const Acts::ConeBounds& Acts::ConeSurface::bounds() const {
   // is safe because no constructor w/o bounds exists
   return (*m_bounds.get());
 }
+
+Acts::Polyhedron Acts::ConeSurface::polyhedronRepresentation(
+    const GeometryContext& gctx, size_t lseg) const {
+  // Prepare vertices and faces
+  std::vector<Vector3D> vertices;
+  std::vector<Polyhedron::Face> faces;
+  std::vector<Polyhedron::Face> triangularMesh;
+
+  if (bounds().minZ() == -std::numeric_limits<double>::infinity() or
+      bounds().maxZ() == std::numeric_limits<double>::infinity()) {
+    throw std::domain_error(
+        "Polyhedron repr of boundless surface not possible");
+  }
+
+  auto ctransform = transform(gctx);
+
+  // The tip - created only once and only
+  // if the we don't have a cut-off cone
+  bool tipExists = false;
+  if (bounds().minZ() * bounds().maxZ() <= s_onSurfaceTolerance) {
+    vertices.push_back(ctransform * Vector3D(0., 0., 0.));
+    tipExists = true;
+  }
+
+  // Cone parameters
+  double hPhiSec = bounds().halfPhiSector();
+  double avgPhi = bounds().averagePhi();
+  bool fullCone = (hPhiSec == M_PI);
+
+  // Get the phi segments from the helper
+  auto phiSegs = fullCone ? detail::VerticesHelper::phiSegments()
+                          : detail::VerticesHelper::phiSegments(
+                                avgPhi - hPhiSec, avgPhi + hPhiSec, {avgPhi});
+
+  // Negative cone if exists
+  std::vector<double> coneSides;
+  if (std::abs(bounds().minZ()) > s_onSurfaceTolerance) {
+    coneSides.push_back(bounds().minZ());
+  }
+  if (std::abs(bounds().maxZ()) > s_onSurfaceTolerance) {
+    coneSides.push_back(bounds().maxZ());
+  }
+  for (auto& z : coneSides) {
+    // Remember the first vertex
+    size_t firstIv = vertices.size();
+    // Radius and z offset
+    double r = std::abs(z) * bounds().tanAlpha();
+    Vector3D zoffset(0., 0., z);
+    for (unsigned int iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
+      int addon = (iseg == phiSegs.size() - 2 and not fullCone) ? 1 : 0;
+      detail::VerticesHelper::createSegment(vertices, {r, r}, phiSegs[iseg],
+                                            phiSegs[iseg + 1], lseg, addon,
+                                            zoffset, ctransform);
+    }
+    // Create the faces
+    if (tipExists) {
+      for (size_t iv = firstIv + 2; iv < vertices.size() + 1; ++iv) {
+        size_t one = 0, two = iv - 1, three = iv - 2;
+        if (z < 0.) {
+          std::swap(two, three);
+        }
+        faces.push_back({one, two, three});
+      }
+      // Complete cone if necessary
+      if (fullCone) {
+        if (z > 0.) {
+          faces.push_back({0, firstIv, vertices.size() - 1});
+        } else {
+          faces.push_back({0, vertices.size() - 1, firstIv});
+        }
+      }
+    }
+  }
+  // if no tip exists, connect the two bows
+  if (tipExists) {
+    triangularMesh = faces;
+  } else {
+    auto facesMesh =
+        detail::FacesHelper::cylindricalFaceMesh(vertices, fullCone);
+    faces = facesMesh.first;
+    triangularMesh = facesMesh.second;
+  }
+  return Polyhedron(vertices, faces, triangularMesh);
+}
\ No newline at end of file
diff --git a/Core/src/Surfaces/CylinderBounds.cpp b/Core/src/Surfaces/CylinderBounds.cpp
index fcc4a4d49d83c4afe3f5324df811c2ef28985989..063c145c9c0af4992f48f61bd171e2ed500ac478 100644
--- a/Core/src/Surfaces/CylinderBounds.cpp
+++ b/Core/src/Surfaces/CylinderBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CylinderBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/CylinderBounds.hpp"
 
 #include <cmath>
@@ -43,8 +39,6 @@ Acts::CylinderBounds::CylinderBounds(double radius, double averagePhi,
   }
 }
 
-Acts::CylinderBounds::~CylinderBounds() = default;
-
 Acts::CylinderBounds* Acts::CylinderBounds::clone() const {
   return new CylinderBounds(*this);
 }
diff --git a/Core/src/Surfaces/CylinderSurface.cpp b/Core/src/Surfaces/CylinderSurface.cpp
index 3ce4d42c5528760902d254a8f439e13089155fec..cd7ebd3e2cad3a681c239531a6104a60f38789fc 100644
--- a/Core/src/Surfaces/CylinderSurface.cpp
+++ b/Core/src/Surfaces/CylinderSurface.cpp
@@ -1,24 +1,20 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// CylinderSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/CylinderSurface.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
-
 #include <cassert>
 #include <cmath>
 #include <iomanip>
 #include <iostream>
 #include <utility>
 
+#include "Acts/Surfaces/detail/FacesHelper.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 
 using Acts::VectorHelpers::perp;
@@ -189,37 +185,37 @@ const Acts::CylinderBounds& Acts::CylinderSurface::bounds() const {
   return (*m_bounds.get());
 }
 
-Acts::PolyhedronRepresentation Acts::CylinderSurface::polyhedronRepresentation(
-    const GeometryContext& gctx, size_t l0div, size_t /*unused*/) const {
+Acts::Polyhedron Acts::CylinderSurface::polyhedronRepresentation(
+    const GeometryContext& gctx, size_t lseg) const {
+  // Prepare vertices and faces
   std::vector<Vector3D> vertices;
-  std::vector<std::vector<size_t>> faces;
-
-  if (l0div <= 1) {
-    throw std::domain_error(
-        "Polyhedron repr of cylinder with 1 div is undefined");
-  }
-
-  double phistep = 2 * M_PI / l0div;
-  double hlZ = bounds().halflengthZ();
-  double r = bounds().r();
-
-  Vector3D left(r, 0, -hlZ);
-  Vector3D right(r, 0, hlZ);
-
-  const Transform3D& sfTransform = transform(gctx);
-
-  for (size_t i = 0; i < l0div; i++) {
-    Transform3D rot(AngleAxis3D(i * phistep, Vector3D::UnitZ()));
-    vertices.push_back(sfTransform * rot * left);
-    vertices.push_back(sfTransform * rot * right);
+  std::vector<Polyhedron::Face> faces;
+  std::vector<Polyhedron::Face> triangularMesh;
+
+  auto ctrans = transform(gctx);
+  bool fullCylinder = bounds().coversFullAzimuth();
+
+  // Get the phi segments from the helper - ensures extra points
+  auto phiSegs = fullCylinder
+                     ? detail::VerticesHelper::phiSegments()
+                     : detail::VerticesHelper::phiSegments(
+                           bounds().averagePhi() - bounds().halfPhiSector(),
+                           bounds().averagePhi() + bounds().halfPhiSector(),
+                           {bounds().averagePhi()});
+
+  // Write the two bows/circles on either side
+  std::vector<int> sides = {-1, 1};
+  for (auto& side : sides) {
+    for (size_t iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
+      int addon = (iseg == phiSegs.size() - 2 and not fullCylinder) ? 1 : 0;
+      /// Helper method to create the segment
+      detail::VerticesHelper::createSegment(
+          vertices, {bounds().r(), bounds().r()}, phiSegs[iseg],
+          phiSegs[iseg + 1], lseg, addon,
+          Vector3D(0., 0., side * bounds().halflengthZ()), ctrans);
+    }
   }
-
-  for (size_t v = 0; v < vertices.size() - 2; v = v + 2) {
-    faces.push_back({v, v + 1, v + 3, v + 2});
-  }
-  if (l0div > 2) {
-    faces.push_back({vertices.size() - 2, vertices.size() - 1, 1, 0});
-  }
-
-  return PolyhedronRepresentation(vertices, faces);
+  auto facesMesh =
+      detail::FacesHelper::cylindricalFaceMesh(vertices, fullCylinder);
+  return Polyhedron(vertices, facesMesh.first, facesMesh.second);
 }
diff --git a/Core/src/Surfaces/DiamondBounds.cpp b/Core/src/Surfaces/DiamondBounds.cpp
index bcf8b729b1d04628416cf23bd79bfecc09db7982..b110452e0fd3a2cc12db957b91ad71cdc384e3ef 100644
--- a/Core/src/Surfaces/DiamondBounds.cpp
+++ b/Core/src/Surfaces/DiamondBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiamondBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/DiamondBounds.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 
@@ -29,8 +25,6 @@ Acts::DiamondBounds::DiamondBounds(double x1, double x2, double x3, double y1,
   throw_assert((x3 <= x2), "Hexagon must be a convex polygon");
 }
 
-Acts::DiamondBounds::~DiamondBounds() = default;
-
 Acts::DiamondBounds* Acts::DiamondBounds::clone() const {
   return new DiamondBounds(*this);
 }
@@ -59,10 +53,12 @@ double Acts::DiamondBounds::distanceToBoundary(
   return BoundaryCheck(true).distance(lposition, vertices());
 }
 
-std::vector<Acts::Vector2D> Acts::DiamondBounds::vertices() const {
-  // vertices starting from lower right in clock-wise order
-  return {{x1(), -y1()}, {x2(), 0},  {x3(), y2()},
-          {-x3(), y2()}, {-x2(), 0}, {-x1(), -y1()}};
+std::vector<Acts::Vector2D> Acts::DiamondBounds::vertices(
+    unsigned int /*lseg*/) const {
+  // Vertices starting at lower left (min rel. phi)
+  // counter-clockwise
+  return {{-x1(), -y1()}, {x1(), -y1()}, {x2(), 0.},
+          {x3(), y2()},   {-x3(), y2()}, {-x2(), 0.}};
 }
 
 const Acts::RectangleBounds& Acts::DiamondBounds::boundingBox() const {
diff --git a/Core/src/Surfaces/DiscSurface.cpp b/Core/src/Surfaces/DiscSurface.cpp
index 90f566441416f390c1116fd118f465676e61e37c..9c1ea9e3900d4601470335fbc875e794a6e4e85e 100644
--- a/Core/src/Surfaces/DiscSurface.cpp
+++ b/Core/src/Surfaces/DiscSurface.cpp
@@ -1,26 +1,25 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiscSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/DiscSurface.hpp"
 
+#include <algorithm>
 #include <cmath>
 #include <iomanip>
 #include <iostream>
+#include <numeric>
 #include <utility>
+#include <vector>
 
-#include "Acts/Surfaces/DiscTrapezoidalBounds.hpp"
+#include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
 #include "Acts/Surfaces/InfiniteBounds.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
 #include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/detail/FacesHelper.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 
@@ -48,7 +47,7 @@ Acts::DiscSurface::DiscSurface(std::shared_ptr<const Transform3D> htrans,
                                double minR, double avephi, double stereo)
     : GeometryObject(),
       Surface(std::move(htrans)),
-      m_bounds(std::make_shared<const DiscTrapezoidalBounds>(
+      m_bounds(std::make_shared<const DiscTrapezoidBounds>(
           minhalfx, maxhalfx, maxR, minR, avephi, stereo)) {}
 
 Acts::DiscSurface::DiscSurface(std::shared_ptr<const Transform3D> htrans,
@@ -99,8 +98,8 @@ bool Acts::DiscSurface::globalToLocal(const GeometryContext& gctx,
 
 const Acts::Vector2D Acts::DiscSurface::localPolarToLocalCartesian(
     const Vector2D& locpol) const {
-  const DiscTrapezoidalBounds* dtbo =
-      dynamic_cast<const Acts::DiscTrapezoidalBounds*>(&(bounds()));
+  const DiscTrapezoidBounds* dtbo =
+      dynamic_cast<const Acts::DiscTrapezoidBounds*>(&(bounds()));
   if (dtbo != nullptr) {
     double rMedium = dtbo->rCenter();
     double phi = dtbo->averagePhi();
@@ -153,42 +152,46 @@ const Acts::SurfaceBounds& Acts::DiscSurface::bounds() const {
   return s_noBounds;
 }
 
-Acts::PolyhedronRepresentation Acts::DiscSurface::polyhedronRepresentation(
-    const GeometryContext& gctx, size_t l0div, size_t /*unused*/) const {
+Acts::Polyhedron Acts::DiscSurface::polyhedronRepresentation(
+    const GeometryContext& gctx, size_t lseg) const {
+  // Prepare vertices and faces
   std::vector<Vector3D> vertices;
-  std::vector<std::vector<size_t>> faces;
-
-  if (l0div < 3) {
-    throw std::domain_error("Polyhedron repr of disk with <3 div is undefined");
-  }
+  std::vector<Polyhedron::Face> faces;
+  std::vector<Polyhedron::Face> triangularMesh;
 
-  auto bounds = std::dynamic_pointer_cast<const RadialBounds>(m_bounds);
-  if (!bounds) {
+  // Understand the disc
+  bool fullDisc = m_bounds->coversFullAzimuth();
+  bool toCenter = m_bounds->rMin() < s_onSurfaceTolerance;
+  // If you have bounds you can create a polyhedron representation
+  if (m_bounds) {
+    auto vertices2D = m_bounds->vertices(lseg);
+    vertices.reserve(vertices2D.size() + 1);
+    Vector3D wCenter(0., 0., 0);
+    for (const auto& v2D : vertices2D) {
+      vertices.push_back(transform(gctx) * Vector3D(v2D.x(), v2D.y(), 0.));
+      wCenter += (*vertices.rbegin());
+    }
+    // These are convex shapes, use the helper method
+    // For rings there's a sweet spot when this stops working
+    if (m_bounds->type() == SurfaceBounds::DiscTrapezoidal or toCenter or
+        not fullDisc) {
+      // Transformt hem into the vertex frame
+      wCenter *= 1. / vertices.size();
+      vertices.push_back(wCenter);
+      auto facesMesh = detail::FacesHelper::convexFaceMesh(vertices, true);
+      faces = facesMesh.first;
+      triangularMesh = facesMesh.second;
+    } else {
+      // Two concentric rings, we use the pure concentric method momentarily,
+      // but that creates too  many unneccesarry faces, when only two
+      // are needed to descibe the mesh, @todo investigate merging flag
+      auto facesMesh = detail::FacesHelper::cylindricalFaceMesh(vertices, true);
+      faces = facesMesh.first;
+      triangularMesh = facesMesh.second;
+    }
+  } else {
     throw std::domain_error(
-        "Polyhedron repr of disk with non RadialBounds currently unsupported");
+        "Polyhedron repr of boundless surface not possible.");
   }
-
-  double phistep = 2 * M_PI / l0div;
-  double rMin = bounds->rMin();
-  double rMax = bounds->rMax();
-
-  Vector3D inner(rMin, 0, 0);
-  Vector3D outer(rMax, 0, 0);
-
-  const Transform3D& sfTransform = transform(gctx);
-
-  for (size_t i = 0; i < l0div; i++) {
-    Transform3D rot(AngleAxis3D(i * phistep, Vector3D::UnitZ()));
-    vertices.push_back(sfTransform * rot * inner);
-    vertices.push_back(sfTransform * rot * outer);
-  }
-
-  for (size_t v = 0; v < vertices.size() - 2; v = v + 2) {
-    faces.push_back({v, v + 1, v + 3, v + 2});
-  }
-  if (l0div > 2) {
-    faces.push_back({vertices.size() - 2, vertices.size() - 1, 1, 0});
-  }
-
-  return PolyhedronRepresentation(vertices, faces);
+  return Polyhedron(vertices, faces, triangularMesh);
 }
diff --git a/Core/src/Surfaces/DiscTrapezoidalBounds.cpp b/Core/src/Surfaces/DiscTrapezoidBounds.cpp
similarity index 66%
rename from Core/src/Surfaces/DiscTrapezoidalBounds.cpp
rename to Core/src/Surfaces/DiscTrapezoidBounds.cpp
index 40d2b1c279dee0730377cab42f85b1c87b6dc2cb..f6197d0a6afc814cf84063265d28908a75fa2936 100644
--- a/Core/src/Surfaces/DiscTrapezoidalBounds.cpp
+++ b/Core/src/Surfaces/DiscTrapezoidBounds.cpp
@@ -1,16 +1,12 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// DiscTrapezoidalBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
-#include "Acts/Surfaces/DiscTrapezoidalBounds.hpp"
+#include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
 
 #include <cmath>
 #include <iomanip>
@@ -18,10 +14,9 @@
 
 #include "Acts/Utilities/detail/periodic.hpp"
 
-Acts::DiscTrapezoidalBounds::DiscTrapezoidalBounds(double minhalfx,
-                                                   double maxhalfx, double maxR,
-                                                   double minR, double avephi,
-                                                   double stereo)
+Acts::DiscTrapezoidBounds::DiscTrapezoidBounds(double minhalfx, double maxhalfx,
+                                               double minR, double maxR,
+                                               double avephi, double stereo)
     : m_rMin(std::min(std::abs(minR), std::abs(maxR))),
       m_rMax(std::max(std::abs(minR), std::abs(maxR))),
       m_minHalfX(std::abs(minhalfx)),
@@ -29,18 +24,16 @@ Acts::DiscTrapezoidalBounds::DiscTrapezoidalBounds(double minhalfx,
       m_avgPhi(detail::radian_sym(avephi)),
       m_stereo(stereo) {}
 
-Acts::DiscTrapezoidalBounds::~DiscTrapezoidalBounds() = default;
-
-Acts::DiscTrapezoidalBounds* Acts::DiscTrapezoidalBounds::clone() const {
-  return new DiscTrapezoidalBounds(*this);
+Acts::DiscTrapezoidBounds* Acts::DiscTrapezoidBounds::clone() const {
+  return new DiscTrapezoidBounds(*this);
 }
 
-Acts::SurfaceBounds::BoundsType Acts::DiscTrapezoidalBounds::type() const {
+Acts::SurfaceBounds::BoundsType Acts::DiscTrapezoidBounds::type() const {
   return SurfaceBounds::DiscTrapezoidal;
 }
 
-std::vector<TDD_real_t> Acts::DiscTrapezoidalBounds::valueStore() const {
-  std::vector<TDD_real_t> values(DiscTrapezoidalBounds::bv_length);
+std::vector<TDD_real_t> Acts::DiscTrapezoidBounds::valueStore() const {
+  std::vector<TDD_real_t> values(DiscTrapezoidBounds::bv_length);
   values[bv_rMin] = rMin();
   values[bv_rMax] = rMax();
   values[bv_minHalfX] = minHalflengthX();
@@ -50,13 +43,13 @@ std::vector<TDD_real_t> Acts::DiscTrapezoidalBounds::valueStore() const {
   return values;
 }
 
-Acts::Vector2D Acts::DiscTrapezoidalBounds::toLocalCartesian(
+Acts::Vector2D Acts::DiscTrapezoidBounds::toLocalCartesian(
     const Acts::Vector2D& lposition) const {
   return {lposition[eLOC_R] * std::sin(lposition[eLOC_PHI] - m_avgPhi),
           lposition[eLOC_R] * std::cos(lposition[eLOC_PHI] - m_avgPhi)};
 }
 
-Acts::ActsMatrixD<2, 2> Acts::DiscTrapezoidalBounds::jacobianToLocalCartesian(
+Acts::ActsMatrixD<2, 2> Acts::DiscTrapezoidBounds::jacobianToLocalCartesian(
     const Acts::Vector2D& lposition) const {
   ActsMatrixD<2, 2> jacobian;
   jacobian(0, eLOC_R) = std::sin(lposition[eLOC_PHI] - m_avgPhi);
@@ -66,7 +59,7 @@ Acts::ActsMatrixD<2, 2> Acts::DiscTrapezoidalBounds::jacobianToLocalCartesian(
   return jacobian;
 }
 
-bool Acts::DiscTrapezoidalBounds::inside(
+bool Acts::DiscTrapezoidBounds::inside(
     const Acts::Vector2D& lposition, const Acts::BoundaryCheck& bcheck) const {
   Vector2D vertices[] = {{minHalflengthX(), rMin()},
                          {maxHalflengthX(), rMax()},
@@ -77,7 +70,7 @@ bool Acts::DiscTrapezoidalBounds::inside(
                                                vertices);
 }
 
-double Acts::DiscTrapezoidalBounds::distanceToBoundary(
+double Acts::DiscTrapezoidBounds::distanceToBoundary(
     const Acts::Vector2D& lposition) const {
   Vector2D vertices[] = {{minHalflengthX(), rMin()},
                          {maxHalflengthX(), rMax()},
@@ -86,11 +79,20 @@ double Acts::DiscTrapezoidalBounds::distanceToBoundary(
   return BoundaryCheck(true).distance(toLocalCartesian(lposition), vertices);
 }
 
+std::vector<Acts::Vector2D> Acts::DiscTrapezoidBounds::vertices(
+    unsigned int /*lseg*/) const {
+  Vector2D cAxis(std::cos(m_avgPhi), std::sin(m_avgPhi));
+  Vector2D nAxis(cAxis.y(), -cAxis.x());
+  return {
+      m_rMin * cAxis - m_minHalfX * nAxis, m_rMin * cAxis + m_minHalfX * nAxis,
+      m_rMax * cAxis + m_maxHalfX * nAxis, m_rMax * cAxis - m_maxHalfX * nAxis};
+}
+
 // ostream operator overload
-std::ostream& Acts::DiscTrapezoidalBounds::toStream(std::ostream& sl) const {
+std::ostream& Acts::DiscTrapezoidBounds::toStream(std::ostream& sl) const {
   sl << std::setiosflags(std::ios::fixed);
   sl << std::setprecision(7);
-  sl << "Acts::DiscTrapezoidalBounds:  (innerRadius, outerRadius, hMinX, "
+  sl << "Acts::DiscTrapezoidBounds:  (innerRadius, outerRadius, hMinX, "
         "hMaxX, hlengthY, hPhiSector, averagePhi, rCenter, stereo) = ";
   sl << "(" << rMin() << ", " << rMax() << ", " << minHalflengthX() << ", "
      << maxHalflengthX() << ", " << halflengthY() << ", " << halfPhiSector()
diff --git a/Core/src/Surfaces/EllipseBounds.cpp b/Core/src/Surfaces/EllipseBounds.cpp
index e38aebde109f58502fb2e2c0e428b90fd26d5c16..7d4f06d691f68cb91d8b89ffeb4993ba506690b5 100644
--- a/Core/src/Surfaces/EllipseBounds.cpp
+++ b/Core/src/Surfaces/EllipseBounds.cpp
@@ -1,16 +1,13 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// EllipseBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/EllipseBounds.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 
 #include <cmath>
 #include <iomanip>
@@ -34,8 +31,6 @@ Acts::EllipseBounds::EllipseBounds(double minRadius0, double minRadius1,
       m_boundingBox(std::max(minRadius0, maxRadius0),
                     std::max(minRadius1, maxRadius1)) {}
 
-Acts::EllipseBounds::~EllipseBounds() = default;
-
 Acts::EllipseBounds* Acts::EllipseBounds::clone() const {
   return new EllipseBounds(*this);
 }
@@ -145,9 +140,10 @@ double Acts::EllipseBounds::distanceToBoundary(
   return sf;
 }
 
-std::vector<Acts::Vector2D> Acts::EllipseBounds::vertices() const {
-  // 2017-04-08 msmk: this is definitely too coarse
-  return {{rMaxX(), 0}, {0, rMaxY()}, {-rMaxX(), 0}, {0, -rMaxY()}};
+std::vector<Acts::Vector2D> Acts::EllipseBounds::vertices(
+    unsigned int lseg) const {
+  return detail::VerticesHelper::ellispoidVertices(
+      m_rMinX, m_rMinY, m_rMaxX, m_rMaxY, m_avgPhi, m_halfPhi, lseg);
 }
 
 const Acts::RectangleBounds& Acts::EllipseBounds::boundingBox() const {
diff --git a/Core/src/Surfaces/LineBounds.cpp b/Core/src/Surfaces/LineBounds.cpp
index e5bdabbee0bef31b347feb339ddbcb9fe044f141..fe29188e66412ee27b2a29ca4d146ae74ad1a05c 100644
--- a/Core/src/Surfaces/LineBounds.cpp
+++ b/Core/src/Surfaces/LineBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// LineBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/LineBounds.hpp"
 
 #include <iomanip>
@@ -18,8 +14,6 @@
 Acts::LineBounds::LineBounds(double radius, double halez)
     : m_radius(std::abs(radius)), m_halfZ(std::abs(halez)) {}
 
-Acts::LineBounds::~LineBounds() = default;
-
 Acts::LineBounds* Acts::LineBounds::clone() const {
   return new LineBounds(*this);
 }
diff --git a/Core/src/Surfaces/LineSurface.cpp b/Core/src/Surfaces/LineSurface.cpp
index e165be96306e7ba326ad18967685347703cc70b9..023f2c05006847758e2e3a6b22d873ccd239717b 100644
--- a/Core/src/Surfaces/LineSurface.cpp
+++ b/Core/src/Surfaces/LineSurface.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// LineSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/LineSurface.hpp"
 
 #include <cmath>
diff --git a/Core/src/Surfaces/PerigeeSurface.cpp b/Core/src/Surfaces/PerigeeSurface.cpp
index 18e6dd5803eb4377cdea3688cc819775dd3d4537..582fb5432e048ffffea740cbc6165cb18a9b86b4 100644
--- a/Core/src/Surfaces/PerigeeSurface.cpp
+++ b/Core/src/Surfaces/PerigeeSurface.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-/////////////////////////////////////////////////////////////////
-// PerigeeSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 
 #include <iomanip>
@@ -71,3 +67,24 @@ std::ostream& Acts::PerigeeSurface::toStream(const GeometryContext& gctx,
   sl << std::setprecision(-1);
   return sl;
 }
+
+Acts::Polyhedron Acts::PerigeeSurface::polyhedronRepresentation(
+    const GeometryContext& gctx, size_t /*lseg*/) const {
+  // Prepare vertices and faces
+  std::vector<Vector3D> vertices;
+  std::vector<Polyhedron::Face> faces;
+  std::vector<Polyhedron::Face> triangularMesh;
+
+  const Transform3D& ctransform = transform(gctx);
+  Vector3D left(0, 0, -100.);
+  Vector3D right(0, 0, 100.);
+
+  // The central wire/straw
+  vertices.push_back(ctransform * left);
+  vertices.push_back(ctransform * right);
+  faces.push_back({0, 1});
+  vertices.push_back(ctransform * Vector3D(0., 0., 0.));
+  triangularMesh.push_back({0, 2, 1});
+
+  return Polyhedron(vertices, faces, triangularMesh);
+}
diff --git a/Core/src/Surfaces/PlaneSurface.cpp b/Core/src/Surfaces/PlaneSurface.cpp
index 435d1fe4d18580622f55a13824059e41b97f913d..c55d3fa8ae96b37f26b518afad128a1a3ef98159 100644
--- a/Core/src/Surfaces/PlaneSurface.cpp
+++ b/Core/src/Surfaces/PlaneSurface.cpp
@@ -1,23 +1,22 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// PlaneSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/PlaneSurface.hpp"
 
 #include <cmath>
 #include <iomanip>
 #include <iostream>
+#include <numeric>
 
+#include "Acts/Surfaces/EllipseBounds.hpp"
 #include "Acts/Surfaces/InfiniteBounds.hpp"
 #include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Surfaces/detail/FacesHelper.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 
 Acts::PlaneSurface::PlaneSurface(const PlaneSurface& other)
@@ -118,3 +117,49 @@ const Acts::SurfaceBounds& Acts::PlaneSurface::bounds() const {
   }
   return s_noBounds;
 }
+
+Acts::Polyhedron Acts::PlaneSurface::polyhedronRepresentation(
+    const GeometryContext& gctx, size_t lseg) const {
+  // Prepare vertices and faces
+  std::vector<Vector3D> vertices;
+  std::vector<Polyhedron::Face> faces;
+  std::vector<Polyhedron::Face> triangularMesh;
+
+  // If you have bounds you can create a polyhedron representation
+  if (m_bounds) {
+    auto vertices2D = m_bounds->vertices(lseg);
+    vertices.reserve(vertices2D.size() + 1);
+    for (const auto& v2D : vertices2D) {
+      vertices.push_back(transform(gctx) * Vector3D(v2D.x(), v2D.y(), 0.));
+    }
+    bool isEllipse = bounds().type() == SurfaceBounds::Ellipse;
+    bool innerExists = false, coversFull = false;
+    if (isEllipse) {
+      auto vStore = bounds().valueStore();
+      innerExists = std::abs(vStore[EllipseBounds::BoundValues::bv_rMinX]) <
+                    s_onSurfaceTolerance;
+      coversFull =
+          std::abs(vStore[EllipseBounds::BoundValues::bv_halfPhiSector]) <
+          M_PI - s_onSurfaceTolerance;
+    }
+    // All of those can be described as convex
+    // @todo same as for Discs: coversFull is not the right criterium
+    // for triangulation
+    if (not isEllipse or not innerExists or not coversFull) {
+      auto facesMesh = detail::FacesHelper::convexFaceMesh(vertices);
+      faces = facesMesh.first;
+      triangularMesh = facesMesh.second;
+    } else {
+      // Two concentric rings, we use the pure concentric method momentarily,
+      // but that creates too  many unneccesarry faces, when only two
+      // are needed to descibe the mesh, @todo investigate merging flag
+      auto facesMesh = detail::FacesHelper::cylindricalFaceMesh(vertices, true);
+      faces = facesMesh.first;
+      triangularMesh = facesMesh.second;
+    }
+  } else {
+    throw std::domain_error(
+        "Polyhedron repr of boundless surface not possible.");
+  }
+  return Polyhedron(vertices, faces, triangularMesh);
+}
diff --git a/Core/src/Surfaces/PolyhedronRepresentation.cpp b/Core/src/Surfaces/PolyhedronRepresentation.cpp
deleted file mode 100644
index 6d37aae75983bf15455402eb59e6acee7a3d0664..0000000000000000000000000000000000000000
--- a/Core/src/Surfaces/PolyhedronRepresentation.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2018 CERN for the benefit of the Acts project
-//
-// This Source Code Form is subject to the terms of the Mozilla Public
-// License, v. 2.0. If a copy of the MPL was not distributed with this
-// file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include <sstream>
-
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
-
-std::string Acts::PolyhedronRepresentation::objString(size_t vtxOffset) const {
-  std::stringstream sstr;
-
-  for (const auto& vtx : vertices) {
-    sstr << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << std::endl;
-  }
-  for (const auto& face : faces) {
-    sstr << "f";
-    for (const auto& idx : face) {
-      sstr << " " << (vtxOffset + idx + 1);
-    }
-    sstr << std::endl;
-  }
-  return sstr.str();
-}
diff --git a/Core/src/Surfaces/RadialBounds.cpp b/Core/src/Surfaces/RadialBounds.cpp
index 577c97a7e0b17e99f01dbfdca44086300983c8eb..c4952f24e04d9e276826e1f534c194e72d4e6edd 100644
--- a/Core/src/Surfaces/RadialBounds.cpp
+++ b/Core/src/Surfaces/RadialBounds.cpp
@@ -1,16 +1,13 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// RadialBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 #include "Acts/Utilities/detail/periodic.hpp"
 
 #include <cmath>
@@ -27,8 +24,6 @@ Acts::RadialBounds::RadialBounds(double minrad, double maxrad, double avephi,
       m_avgPhi(detail::radian_sym(avephi)),
       m_halfPhi(std::abs(hphisec)) {}
 
-Acts::RadialBounds::~RadialBounds() = default;
-
 Acts::RadialBounds* Acts::RadialBounds::clone() const {
   return new RadialBounds(*this);
 }
@@ -67,6 +62,12 @@ double Acts::RadialBounds::distanceToBoundary(
                                       Vector2D(rMax(), halfPhiSector()));
 }
 
+std::vector<Acts::Vector2D> Acts::RadialBounds::vertices(
+    unsigned int lseg) const {
+  return detail::VerticesHelper::circularVertices(m_rMin, m_rMax, m_avgPhi,
+                                                  m_halfPhi, lseg);
+}
+
 // ostream operator overload
 std::ostream& Acts::RadialBounds::toStream(std::ostream& sl) const {
   sl << std::setiosflags(std::ios::fixed);
diff --git a/Core/src/Surfaces/RectangleBounds.cpp b/Core/src/Surfaces/RectangleBounds.cpp
index d5f154092d23cada67460dc3f1149443b4f8d345..abca4abcbd00d26b9712640f8c6eaa32d3bc3a4a 100644
--- a/Core/src/Surfaces/RectangleBounds.cpp
+++ b/Core/src/Surfaces/RectangleBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// RectangleBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
 
@@ -44,14 +40,10 @@ double Acts::RectangleBounds::distanceToBoundary(
   return BoundaryCheck(true).distance(lposition, m_min, m_max);
 }
 
-std::vector<Acts::Vector2D> Acts::RectangleBounds::vertices() const {
-  // counter-clockwise starting from bottom-right corner
-  return {
-      {m_max.x(), m_min.y()},
-      m_max,
-      {m_min.x(), m_max.y()},
-      m_min,
-  };
+std::vector<Acts::Vector2D> Acts::RectangleBounds::vertices(
+    unsigned int /*lseg*/) const {
+  // counter-clockwise starting from bottom-left corner
+  return {m_min, {m_max.x(), m_min.y()}, m_max, {m_min.x(), m_max.y()}};
 }
 
 const Acts::RectangleBounds& Acts::RectangleBounds::boundingBox() const {
diff --git a/Core/src/Surfaces/StrawSurface.cpp b/Core/src/Surfaces/StrawSurface.cpp
index f8d55714f2dcf859212a8b34aecb4809bf46db14..91a3dce097184de378f3ba18bdc7a51ae98f9211 100644
--- a/Core/src/Surfaces/StrawSurface.cpp
+++ b/Core/src/Surfaces/StrawSurface.cpp
@@ -1,21 +1,18 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// StrawSurface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/StrawSurface.hpp"
-#include "Acts/Surfaces/PolyhedronRepresentation.hpp"
-
 #include <iomanip>
 #include <iostream>
 #include <utility>
+#include "Acts/Geometry/Polyhedron.hpp"
+#include "Acts/Surfaces/detail/FacesHelper.hpp"
+#include "Acts/Surfaces/detail/VerticesHelper.hpp"
 
 #include "Acts/Surfaces/InfiniteBounds.hpp"
 
@@ -58,36 +55,43 @@ Acts::StrawSurface* Acts::StrawSurface::clone_impl(
   return new StrawSurface(gctx, *this, shift);
 }
 
-Acts::PolyhedronRepresentation Acts::StrawSurface::polyhedronRepresentation(
-    const GeometryContext& gctx, size_t l0div, size_t /*l1div*/) const {
+Acts::Polyhedron Acts::StrawSurface::polyhedronRepresentation(
+    const GeometryContext& gctx, size_t lseg) const {
+  // Prepare vertices and faces
   std::vector<Vector3D> vertices;
-  std::vector<std::vector<size_t>> faces;
-
-  if (l0div == 1) {
-    throw std::domain_error("Polyhedron repr of straw with 1 div is undefined");
-  }
-
-  double phistep = 2 * M_PI / l0div;
-  double hlZ = m_bounds->halflengthZ();
-  double r = m_bounds->r();
-
-  Vector3D left(r, 0, -hlZ);
-  Vector3D right(r, 0, hlZ);
-
-  const Transform3D& sfTransform = transform(gctx);
-
-  for (size_t i = 0; i < l0div; i++) {
-    Transform3D rot(AngleAxis3D(i * phistep, Vector3D::UnitZ()));
-    vertices.push_back(sfTransform * rot * left);
-    vertices.push_back(sfTransform * rot * right);
-  }
-
-  for (size_t v = 0; v < vertices.size() - 2; v = v + 2) {
-    faces.push_back({v, v + 1, v + 3, v + 2});
-  }
-  if (l0div > 2) {
-    faces.push_back({vertices.size() - 2, vertices.size() - 1, 1, 0});
+  std::vector<Polyhedron::Face> faces;
+  std::vector<Polyhedron::Face> triangularMesh;
+
+  const Transform3D& ctransform = transform(gctx);
+  // Draw the bounds if more than one segment are chosen
+  if (lseg > 1) {
+    auto phiSegs = detail::VerticesHelper::phiSegments();
+    // Write the two bows/circles on either side
+    std::vector<int> sides = {-1, 1};
+    for (auto& side : sides) {
+      for (size_t iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
+        int addon = (iseg == phiSegs.size() - 2) ? 1 : 0;
+        /// Helper method to create the segment
+        detail::VerticesHelper::createSegment(
+            vertices, {m_bounds->r(), m_bounds->r()}, phiSegs[iseg],
+            phiSegs[iseg + 1], lseg, addon,
+            Vector3D(0., 0., side * m_bounds->halflengthZ()), ctransform);
+      }
+    }
+    auto facesMesh = detail::FacesHelper::cylindricalFaceMesh(vertices);
+    faces = facesMesh.first;
+    triangularMesh = facesMesh.second;
   }
 
-  return PolyhedronRepresentation(vertices, faces);
+  size_t bvertices = vertices.size();
+  Vector3D left(0, 0, -m_bounds->halflengthZ());
+  Vector3D right(0, 0, m_bounds->halflengthZ());
+  // The central wire/straw
+  vertices.push_back(ctransform * left);
+  vertices.push_back(ctransform * right);
+  faces.push_back({bvertices, bvertices + 1});
+  vertices.push_back(ctransform * Vector3D(0., 0., 0.));
+  triangularMesh.push_back({bvertices, bvertices + 2, bvertices + 1});
+
+  return Polyhedron(vertices, faces, triangularMesh);
 }
diff --git a/Core/src/Surfaces/Surface.cpp b/Core/src/Surfaces/Surface.cpp
index a7055a673b8fcec479c12c924a7568bb5eb2b88f..c35580c6ae9a9ea6f70375c27a131716f7d5142f 100644
--- a/Core/src/Surfaces/Surface.cpp
+++ b/Core/src/Surfaces/Surface.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// Surface.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/Surface.hpp"
 
 #include <iomanip>
diff --git a/Core/src/Surfaces/TrapezoidBounds.cpp b/Core/src/Surfaces/TrapezoidBounds.cpp
index a96d9698b404baa773e68283e7b292caee4bf805..8f6e1cb23e95a166f2c3cbf62052adcf56f62ad2 100644
--- a/Core/src/Surfaces/TrapezoidBounds.cpp
+++ b/Core/src/Surfaces/TrapezoidBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// TrapezoidBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/TrapezoidBounds.hpp"
 
 #include <cmath>
@@ -51,7 +47,8 @@ double Acts::TrapezoidBounds::distanceToBoundary(
   return BoundaryCheck(true).distance(lposition, vertices());
 }
 
-std::vector<Acts::Vector2D> Acts::TrapezoidBounds::vertices() const {
+std::vector<Acts::Vector2D> Acts::TrapezoidBounds::vertices(
+    unsigned int /*lseg*/) const {
   // counter-clockwise from bottom-right corner
   return {{minHalflengthX(), -halflengthY()},
           {maxHalflengthX(), halflengthY()},
diff --git a/Core/src/Surfaces/TriangleBounds.cpp b/Core/src/Surfaces/TriangleBounds.cpp
index 5ffa21564ab9e25d4691e92ce363beea65cc2fec..a9e9e09013bbf7d66ca03141f5a618454326dba2 100644
--- a/Core/src/Surfaces/TriangleBounds.cpp
+++ b/Core/src/Surfaces/TriangleBounds.cpp
@@ -1,15 +1,11 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-///////////////////////////////////////////////////////////////////
-// TriangleBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
 #include "Acts/Surfaces/TriangleBounds.hpp"
 
 #include <iomanip>
@@ -57,7 +53,8 @@ double Acts::TriangleBounds::distanceToBoundary(
   return BoundaryCheck(true).distance(lposition, m_vertices);
 }
 
-std::vector<Acts::Vector2D> Acts::TriangleBounds::vertices() const {
+std::vector<Acts::Vector2D> Acts::TriangleBounds::vertices(
+    unsigned int /*lseg*/) const {
   return {std::begin(m_vertices), std::end(m_vertices)};
 }
 
diff --git a/Tests/CommonHelpers/Acts/Tests/CommonHelpers/LineSurfaceStub.hpp b/Tests/CommonHelpers/Acts/Tests/CommonHelpers/LineSurfaceStub.hpp
index 56e5f06aff1a90042ac1b42951e85b31765fad9d..41ed0a6b803ef3e0e222e16a13a5b9429b310aab 100644
--- a/Tests/CommonHelpers/Acts/Tests/CommonHelpers/LineSurfaceStub.hpp
+++ b/Tests/CommonHelpers/Acts/Tests/CommonHelpers/LineSurfaceStub.hpp
@@ -57,6 +57,17 @@ class LineSurfaceStub : public LineSurface {
 
   using Surface::normal;
 
+  /// Return a Polyhedron for the surfaces
+  ///
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param lseg is ignored for a perigee @note ignored
+  ///
+  /// @return A list of vertices and a face/facett description of it
+  Polyhedron polyhedronRepresentation(const GeometryContext& /*gctx*/,
+                                      size_t /*lseg*/) const final {
+    return Polyhedron({}, {}, {});
+  }
+
  private:
   Surface* clone_impl(const GeometryContext& /*gctx*/,
                       const Transform3D& /*unused*/) const {
diff --git a/Tests/CommonHelpers/Acts/Tests/CommonHelpers/ObjTestWriter.hpp b/Tests/CommonHelpers/Acts/Tests/CommonHelpers/ObjTestWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..506b324efcb9785cb308b40300204bed8dd6c66c
--- /dev/null
+++ b/Tests/CommonHelpers/Acts/Tests/CommonHelpers/ObjTestWriter.hpp
@@ -0,0 +1,118 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <fstream>
+#include <vector>
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
+#include "Acts/Geometry/Volume.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Utilities/BoundingBox.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/ObjHelper.hpp"
+
+namespace Acts {
+
+using IdentifiedPolyderon = std::tuple<std::string, bool, Polyhedron>;
+
+namespace Test {
+
+struct ObjTestWriter {
+  /// Helper method to write sector lines for 2D views
+  /// @param name The output file name
+  /// @param linaA The first line
+  /// @param lineB The second line
+  static void writeSectorLinesObj(const std::string& name,
+                                  const std::pair<Vector3D, Vector3D>& lineA,
+                                  const std::pair<Vector3D, Vector3D>& lineB) {
+    std::ofstream ostream;
+    ostream.open(name + ".obj");
+    ObjHelper objH;
+    objH.line(lineA.first, lineA.second);
+    objH.line(lineB.first, lineB.second);
+    objH.write(ostream);
+    ostream.close();
+  }
+
+  /// Helper method to write sector planes for two dimensional sectors
+  /// (symmetric)
+  /// @param name The output file name
+  /// @param phiSec The opening angle
+  /// @param phiAvg The average phi (= center phi position)
+  /// @param hX The half length in X of the sector plane
+  /// @param hY the half length in Y of the sector plane
+  static void writeSectorPlanesObj(const std::string& name, double phiSec,
+                                   double phiAvg, double hX, double hY) {
+    // Construct the helper planes for sectoral building
+    auto sectorBounds = std::make_shared<RectangleBounds>(hX, hY);
+
+    Vector3D helperColX(0., 0., 1.);
+    Vector3D helperColY(1., 0., 0.);
+    Vector3D helperColZ(0., 1., 0.);
+    RotationMatrix3D helperRotation;
+    helperRotation.col(0) = helperColX;
+    helperRotation.col(1) = helperColY;
+    helperRotation.col(2) = helperColZ;
+    // curvilinear surfaces are boundless
+    Transform3D helperTransform{helperRotation};
+
+    auto sectorTransformM = std::make_shared<Transform3D>(helperTransform);
+    sectorTransformM->prerotate(AngleAxis3D(phiAvg - phiSec, helperColX));
+
+    auto sectorTransformP = std::make_shared<Transform3D>(helperTransform);
+    sectorTransformP->prerotate(AngleAxis3D(phiAvg + phiSec, helperColX));
+
+    auto sectorPlaneM =
+        Surface::makeShared<PlaneSurface>(sectorTransformM, sectorBounds);
+
+    auto sectorPlaneP =
+        Surface::makeShared<PlaneSurface>(sectorTransformP, sectorBounds);
+
+    std::ofstream ostream;
+    ostream.open(name + ".obj");
+    ObjHelper objH;
+    sectorPlaneM->polyhedronRepresentation(GeometryContext(), 1).draw(objH);
+    sectorPlaneP->polyhedronRepresentation(GeometryContext(), 1).draw(objH);
+    objH.write(ostream);
+    ostream.close();
+  }
+
+  /// Helper method to be called from sub tests
+  /// It will draw the polyhedron and create a file, the boolean
+  /// steers whether the obj should be triangulated
+  /// @param iphs The Identified Polyhedrons (= with name and boolean)
+  static void writeObj(const std::vector<IdentifiedPolyderon>& iphs) {
+    for (const auto& iph : iphs) {
+      std::ofstream ostream;
+      ostream.open(std::get<std::string>(iph) + ".obj");
+      ObjHelper objH;
+      std::get<Polyhedron>(iph).draw(objH, std::get<bool>(iph));
+      objH.write(ostream);
+      ostream.close();
+    }
+  }
+
+  /// Helper method to be called from sub tests to write bounding boxes
+  /// It will draw the polyhedron and create a file, the boolean
+  /// steers whether the obj should be triangulated
+  /// @param iphs The Identified Polyhedrons (= with name and boolean)
+  static void writeObj(const std::string& name,
+                       const Volume::BoundingBox& bBox) {
+    std::ofstream ostream;
+    ostream.open(name + std::string(".obj"));
+    ObjHelper objH;
+    bBox.draw(objH);
+    objH.write(ostream);
+    ostream.close();
+  }
+};
+
+}  // namespace Test
+
+}  // namespace Acts
\ No newline at end of file
diff --git a/Tests/UnitTests/Core/Geometry/CMakeLists.txt b/Tests/UnitTests/Core/Geometry/CMakeLists.txt
index 3f7c02fd9b046f3c0cbe689d412193dabc532571..7535c5ce500bcb67d3660857861ac11ba2b7f08f 100644
--- a/Tests/UnitTests/Core/Geometry/CMakeLists.txt
+++ b/Tests/UnitTests/Core/Geometry/CMakeLists.txt
@@ -1,24 +1,25 @@
 add_unittest(AlignmentContextTests AlignmentContextTests.cpp)
-add_unittest(ConeLayerTests ConeLayerTests.cpp)
 add_unittest(CuboidVolumeBuilderTests CuboidVolumeBuilderTests.cpp)
-add_unittest(CutoutCylinderVolumeBoundsTest CutoutCylinderVolumeBoundsTest.cpp)
+add_unittest(CutoutCylinderVolumeBoundsTests CutoutCylinderVolumeBoundsTests.cpp)
 add_unittest(CylinderLayerTests CylinderLayerTests.cpp)
 add_unittest(CylinderVolumeBoundsTests CylinderVolumeBoundsTests.cpp)
 add_unittest(CylinderVolumeBuilderTests CylinderVolumeBuilderTests.cpp)
 add_unittest(DiscLayerTests DiscLayerTests.cpp)
-add_unittest(DoubleTrapezoidVolumeBoundsTest DoubleTrapezoidVolumeBoundsTest.cpp)
+add_unittest(DoubleTrapezoidVolumeBoundsTests DoubleTrapezoidVolumeBoundsTests.cpp)
+add_unittest(ExtentTests ExtentTests.cpp)
 add_unittest(GenericApproachDescriptorTests GenericApproachDescriptorTests.cpp)
-add_unittest(GenericCuboidVolumeBoundsTest GenericCuboidVolumeBoundsTest.cpp)
+add_unittest(GenericCuboidVolumeBoundsTests GenericCuboidVolumeBoundsTests.cpp)
 add_unittest(GeometryIDTests GeometryIDTests.cpp)
 add_unittest(LayerCreatorTests LayerCreatorTests.cpp)
 add_unittest(LayerTests LayerTests.cpp)
 add_unittest(NavigationLayerTests NavigationLayerTests.cpp)
 add_unittest(PlaneLayerTests PlaneLayerTests.cpp)
 add_unittest(ProtoLayerTests ProtoLayerTests.cpp)
+add_unittest(PolyhedronTests PolyhedronTests.cpp)
 add_unittest(SimpleGeometryTests SimpleGeometryTests.cpp)
 add_unittest(SurfaceArrayCreatorTests SurfaceArrayCreatorTests.cpp)
 add_unittest(TrackingGeometryClosureGeometryTests TrackingGeometryClosureTests.cpp)
 add_unittest(TrackingGeometryCreationTests TrackingGeometryCreationTests.cpp)
 add_unittest(TrackingGeometryGeoIDTests TrackingGeometryGeoIDTests.cpp)
 add_unittest(TrackingVolumeTests TrackingVolumeTests.cpp)
-add_unittest(TrapezoidVolumeBoundsTest TrapezoidVolumeBoundsTest.cpp)
+add_unittest(TrapezoidVolumeBoundsTests TrapezoidVolumeBoundsTests.cpp)
diff --git a/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTest.cpp b/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp
similarity index 78%
rename from Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTest.cpp
rename to Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp
index fb5cf36888a3c5be99f04563f9179c65cef8dcc9..0fe43355cbf588a10e34f1ea2fde5034f66c7cbc 100644
--- a/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTest.cpp
+++ b/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp
@@ -13,8 +13,11 @@
 #include <memory>
 
 #include "Acts/Geometry/CutoutCylinderVolumeBounds.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Tests/CommonHelpers/ObjTestWriter.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/PlyHelper.hpp"
 
@@ -28,15 +31,6 @@ BOOST_AUTO_TEST_CASE(construction_test) {
   ccvb.toStream(std::cout);
 }
 
-BOOST_AUTO_TEST_CASE(decomposeToSurfaces_test) {
-  CutoutCylinderVolumeBounds ccvb(5, 10, 15, 30, 25);
-  PlyHelper<double> ply;
-  ccvb.draw(ply);
-
-  std::ofstream os("ccvb.ply");
-  os << ply;
-}
-
 BOOST_AUTO_TEST_CASE(inside_test) {
   CutoutCylinderVolumeBounds ccvb(5, 10, 15, 30, 25);
 
@@ -108,10 +102,34 @@ BOOST_AUTO_TEST_CASE(inside_test) {
 }
 
 BOOST_AUTO_TEST_CASE(boundingbox_test) {
+  GeometryContext tgContext = GeometryContext();
+  std::vector<IdentifiedPolyderon> tPolyhedrons;
+
+  auto combineAndDecompose = [&](const SurfacePtrVector& surfaces,
+                                 const std::string& name) -> void {
+    std::string writeBase = std::string("CutoutCylinderVolumeBounds") + name;
+
+    Polyhedron phCombined;
+    size_t is = 0;
+    for (const auto& sf : surfaces) {
+      Polyhedron phComponent = sf->polyhedronRepresentation(tgContext, 72);
+      phCombined.merge(phComponent);
+      tPolyhedrons.push_back(
+          {writeBase + std::string("_comp_") + std::to_string(is++), false,
+           phComponent});
+    }
+    tPolyhedrons.push_back({writeBase, false, phCombined});
+  };
+
   CutoutCylinderVolumeBounds ccvb(5, 10, 15, 30, 25);
   auto box = ccvb.boundingBox();
   CHECK_CLOSE_ABS(box.min(), Vector3D(-15, -15, -30), 1e-6);
   CHECK_CLOSE_ABS(box.max(), Vector3D(15, 15, 30), 1e-6);
+
+  auto ccvbSurfaces = ccvb.decomposeToSurfaces();
+  combineAndDecompose(ccvbSurfaces, "");
+  ObjTestWriter::writeObj("CutoutCylinderVolumeBounds_BB", box);
+  ObjTestWriter::writeObj(tPolyhedrons);
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp b/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp
index b96d90655f9306de47aa1fb7dcad0ba21a5e123a..90c48854bcec73b6338ddd7211e0c581617f7a51 100644
--- a/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp
@@ -10,8 +10,10 @@
 #include <boost/test/unit_test.hpp>
 
 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Tests/CommonHelpers/ObjTestWriter.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
@@ -21,6 +23,7 @@ namespace tt = boost::test_tools;
 namespace Acts {
 
 namespace Test {
+
 BOOST_AUTO_TEST_SUITE(Volumes)
 
 /// Unit test for testing the decomposeToSurfaces() function
@@ -95,10 +98,32 @@ BOOST_DATA_TEST_CASE(CylinderVolumeBounds_decomposeToSurfaces,
 }
 
 BOOST_AUTO_TEST_CASE(bounding_box_creation) {
+  GeometryContext tgContext = GeometryContext();
+  std::vector<IdentifiedPolyderon> tPolyhedrons;
+
+  auto combineAndDecompose = [&](const SurfacePtrVector& surfaces,
+                                 const std::string& name) -> void {
+    std::string writeBase = std::string("CylinderVolumeBounds_") + name;
+
+    Polyhedron phCombined;
+    size_t is = 0;
+    for (const auto& sf : surfaces) {
+      Polyhedron phComponent = sf->polyhedronRepresentation(tgContext, 72);
+      phCombined.merge(phComponent);
+      tPolyhedrons.push_back(
+          {writeBase + std::string("_comp_") + std::to_string(is++), false,
+           phComponent});
+    }
+    tPolyhedrons.push_back({writeBase, false, phCombined});
+  };
+
   float tol = 1e-4;
 
   CylinderVolumeBounds cvb(5, 10);
+  auto cvbSurfaces = cvb.decomposeToSurfaces();
+  combineAndDecompose(cvbSurfaces, "Solid");
   auto bb = cvb.boundingBox();
+  ObjTestWriter::writeObj("CylinderVolumeBounds_Solid_BB", bb);
 
   Transform3D rot;
   rot = AngleAxis3D(M_PI / 2., Vector3D::UnitX());
@@ -117,6 +142,9 @@ BOOST_AUTO_TEST_CASE(bounding_box_creation) {
   BOOST_CHECK_EQUAL(bb.entity(), nullptr);
   BOOST_CHECK_EQUAL(bb.max(), Vector3D(8, 8, 12));
   BOOST_CHECK_EQUAL(bb.min(), Vector3D(-8, -8, -12));
+  cvbSurfaces = cvb.decomposeToSurfaces();
+  combineAndDecompose(cvbSurfaces, "Tube");
+  ObjTestWriter::writeObj("CylinderVolumeBounds_Tube_BB", bb);
 
   double angle = M_PI / 8.;
   cvb = CylinderVolumeBounds(5, 8, angle, 13);
@@ -125,6 +153,9 @@ BOOST_AUTO_TEST_CASE(bounding_box_creation) {
   CHECK_CLOSE_ABS(bb.max(), Vector3D(8, 8 * std::sin(angle), 13), tol);
   CHECK_CLOSE_ABS(
       bb.min(), Vector3D(5 * std::cos(angle), -8 * std::sin(angle), -13), tol);
+  cvbSurfaces = cvb.decomposeToSurfaces();
+  combineAndDecompose(cvbSurfaces, "TubeSector");
+  ObjTestWriter::writeObj("CylinderVolumeBounds_TubeSector_BB", bb);
 
   rot = AngleAxis3D(M_PI / 2., Vector3D::UnitZ());
   bb = cvb.boundingBox(&rot);
@@ -138,6 +169,7 @@ BOOST_AUTO_TEST_CASE(bounding_box_creation) {
   BOOST_CHECK_EQUAL(bb.entity(), nullptr);
   CHECK_CLOSE_ABS(bb.max(), Vector3D(8.40007, 15.2828, 3.88911), tol);
   CHECK_CLOSE_ABS(bb.min(), Vector3D(-7.27834, -8.12028, -14.2182), tol);
+  ObjTestWriter::writeObj(tPolyhedrons);
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTest.cpp b/Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTests.cpp
similarity index 100%
rename from Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTest.cpp
rename to Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTests.cpp
diff --git a/Tests/UnitTests/Core/Geometry/ExtentTests.cpp b/Tests/UnitTests/Core/Geometry/ExtentTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f488ec088a47e285e1a098ed7445c2c9dad6a2d
--- /dev/null
+++ b/Tests/UnitTests/Core/Geometry/ExtentTests.cpp
@@ -0,0 +1,77 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/data/test_case.hpp>
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+// Helper
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+
+// The class to test
+#include "Acts/Geometry/Extent.hpp"
+
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/Units.hpp"
+
+namespace Acts {
+
+using namespace UnitLiterals;
+
+namespace Test {
+
+BOOST_AUTO_TEST_SUITE(Geometry)
+
+/// Unit tests for Polyderon construction & operator +=
+BOOST_AUTO_TEST_CASE(ExtentTest) {
+  std::vector<Vector3D> vertices = {
+      Vector3D(15_mm, -3_mm, -10_mm), Vector3D(18_mm, 0_mm, -10_mm),
+      Vector3D(15_mm, 3_mm, -10_mm),  Vector3D(15_mm, -3_mm, 10_mm),
+      Vector3D(18_mm, 0_mm, 10_mm),   Vector3D(15_mm, 3_mm, 10_mm)};
+
+  // Create an Extent
+  Extent gExt;
+  for (const auto& v : vertices) {
+    gExt.check(v);
+  }
+
+  double phiMin = std::atan2(-3_mm, 15_mm);
+  double phiMax = std::atan2(3_mm, 15_mm);
+  double rMin = std::sqrt(15_mm * 15_mm + 3_mm * 3_mm);
+
+  CHECK_CLOSE_ABS(gExt.min(binX), 15_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.max(binX), 18_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.min(binY), -3_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.max(binY), 3_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.min(binZ), -10_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.max(binZ), 10_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.min(binR), rMin, 1e-6);
+  CHECK_CLOSE_ABS(gExt.max(binR), 18_mm, 1e-6);
+  CHECK_CLOSE_ABS(gExt.min(binPhi), phiMin, 1e-6);
+  CHECK_CLOSE_ABS(gExt.max(binPhi), phiMax, 1e-6);
+
+  // Create a second Extent
+  Extent otherExt;
+  otherExt.extend(gExt);
+
+  CHECK_CLOSE_ABS(otherExt.min(binX), 15_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.max(binX), 18_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.min(binY), -3_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.max(binY), 3_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.min(binZ), -10_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.max(binZ), 10_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.min(binR), rMin, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.max(binR), 18_mm, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.min(binPhi), phiMin, 1e-6);
+  CHECK_CLOSE_ABS(otherExt.max(binPhi), phiMax, 1e-6);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}  // namespace Test
+}  // namespace Acts
\ No newline at end of file
diff --git a/Tests/UnitTests/Core/Geometry/GenericCuboidVolumeBoundsTest.cpp b/Tests/UnitTests/Core/Geometry/GenericCuboidVolumeBoundsTests.cpp
similarity index 100%
rename from Tests/UnitTests/Core/Geometry/GenericCuboidVolumeBoundsTest.cpp
rename to Tests/UnitTests/Core/Geometry/GenericCuboidVolumeBoundsTests.cpp
diff --git a/Tests/UnitTests/Core/Geometry/PolyhedronTests.cpp b/Tests/UnitTests/Core/Geometry/PolyhedronTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..334727a1b60e46431006b6fd476e70be6c3bac58
--- /dev/null
+++ b/Tests/UnitTests/Core/Geometry/PolyhedronTests.cpp
@@ -0,0 +1,122 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/data/test_case.hpp>
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+// Helper
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+
+// The class to test
+#include <fstream>
+#include "Acts/Geometry/Polyhedron.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/ObjHelper.hpp"
+#include "Acts/Utilities/Units.hpp"
+
+namespace Acts {
+
+using namespace UnitLiterals;
+
+namespace Test {
+
+BOOST_AUTO_TEST_SUITE(Geometry)
+
+/// Unit tests for Polyderon construction & operator +=
+BOOST_AUTO_TEST_CASE(PolyhedronTest) {
+  std::vector<Vector3D> tvertices = {Vector3D(-1, -1, 0.), Vector3D(1., -1, 0.),
+                                     Vector3D(0., 1., 0.)};
+  std::vector<std::vector<size_t>> tfaces = {{0, 1, 2}};
+
+  Polyhedron triangle(tvertices, tfaces, tfaces);
+  BOOST_CHECK(tvertices == triangle.vertices);
+  BOOST_CHECK(tfaces == triangle.faces);
+  BOOST_CHECK(tfaces == triangle.triangularMesh);
+
+  std::ofstream tStream;
+  tStream.open("PolyhedronTriangle.obj");
+  ObjHelper objtH;
+  triangle.draw(objtH);
+  objtH.write(tStream);
+  tStream.close();
+
+  std::vector<Vector3D> rvertices = {Vector3D(-1, -2, 0.), Vector3D(1., -2, 0.),
+                                     Vector3D(1., -1., 0.),
+                                     Vector3D(-1., -1., 0.)};
+  std::vector<std::vector<size_t>> rfaces = {{0, 1, 2, 3}};
+  std::vector<std::vector<size_t>> rmesh = {{0, 1, 2}, {2, 3, 0}};
+  Polyhedron rectangle(rvertices, rfaces, rmesh);
+  BOOST_CHECK(rvertices == rectangle.vertices);
+  BOOST_CHECK(rfaces == rectangle.faces);
+  BOOST_CHECK(rmesh == rectangle.triangularMesh);
+
+  std::ofstream rStream;
+  rStream.open("PolyhedronRectangle.obj");
+  ObjHelper objrH;
+  rectangle.draw(objrH);
+  objrH.write(rStream);
+  rStream.close();
+
+  // Now add them
+  Polyhedron tr;
+  tr.merge(triangle);
+  BOOST_CHECK(tr.vertices == triangle.vertices);
+  BOOST_CHECK(tr.faces == triangle.faces);
+  BOOST_CHECK(tr.triangularMesh == triangle.triangularMesh);
+
+  tr.merge(rectangle);
+
+  std::ofstream trStream;
+  trStream.open("PolyhedronTriangleRectangle.obj");
+  ObjHelper objtrH;
+  tr.draw(objtrH);
+  objtrH.write(trStream);
+  trStream.close();
+}
+
+/// Unit tests for Polyderon construction & operator +=
+BOOST_AUTO_TEST_CASE(PolyhedronExtent) {
+  std::vector<Vector3D> rvertices = {Vector3D(-1, -2, 0.), Vector3D(1., -2, 0.),
+                                     Vector3D(1., -1., 0.),
+                                     Vector3D(-1., -1., 0.)};
+
+  std::vector<std::vector<size_t>> rfaces = {{0, 1, 2, 3}};
+  std::vector<std::vector<size_t>> rmesh = {{0, 1, 2}, {2, 3, 0}};
+  Polyhedron rectangle(rvertices, rfaces, rmesh);
+
+  auto rExtent = rectangle.extent();
+  CHECK_CLOSE_ABS(rExtent.min(binX), -1., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binX), 1., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.min(binY), -2., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binY), -1., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.min(binZ), 0., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binZ), 0., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.min(binR), 1., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binR), VectorHelpers::perp(rvertices[0]), 1e-6);
+  CHECK_CLOSE_ABS(rExtent.min(binPhi), VectorHelpers::phi(rvertices[3]), 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binPhi), VectorHelpers::phi(rvertices[2]), 1e-6);
+
+  // Now shift the Extent
+  Vector3D shift(-1., 0., 1.);
+  Transform3D shiftedTransform = Transform3D::Identity();
+  shiftedTransform.pretranslate(shift);
+  rExtent = rectangle.extent(shiftedTransform);
+  CHECK_CLOSE_ABS(rExtent.min(binX), -2., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binX), 0., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.min(binY), -2., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binY), -1., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.min(binZ), 1., 1e-6);
+  CHECK_CLOSE_ABS(rExtent.max(binZ), 1., 1e-6);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}  // namespace Test
+
+}  // namespace Acts
\ No newline at end of file
diff --git a/Tests/UnitTests/Core/Geometry/TrapezoidVolumeBoundsTest.cpp b/Tests/UnitTests/Core/Geometry/TrapezoidVolumeBoundsTests.cpp
similarity index 100%
rename from Tests/UnitTests/Core/Geometry/TrapezoidVolumeBoundsTest.cpp
rename to Tests/UnitTests/Core/Geometry/TrapezoidVolumeBoundsTests.cpp
diff --git a/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp
index 6d14368f96d6ae75d6e6348c2881bc83aa335904..6f5193eec60311412010e505d73723bfcfc00dc5 100644
--- a/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp
@@ -32,80 +32,6 @@ double maxPhi = 1.33970;
 
 Vector2D offset(-2., 2.);
 
-bool writeObj = false;
-
-static void writeAnnulusDiscObj(
-    std::ofstream& stream, double scalor, unsigned int nSegments,
-    const Acts::Transform3D& transform,
-    const std::vector<Acts::Vector2D>& vertices,
-    const Acts::Vector3D& os = Acts::Vector3D(0., 0., 0.)) {
-  unsigned int cvc = 0;
-
-  std::vector<Acts::Vector3D> gVertices;
-  // Write the vertices
-  for (auto& v : vertices) {
-    gVertices.push_back(transform * Acts::Vector3D(v.x(), v.y(), 0.));
-  }
-
-  // Complete if there are segments defined for a bow
-  if (nSegments > 1) {
-    /// Fill the vertices in betwen
-    auto fillBow =
-        [&](const Acts::Vector3D& first,
-            const Acts::Vector3D& second) -> std::vector<Acts::Vector3D> {
-      // The completed list of vertices
-      std::vector<Acts::Vector3D> completed3D;
-
-      double oseg = 1. / nSegments;
-      double phif = Acts::VectorHelpers::phi(first);
-      double phis = Acts::VectorHelpers::phi(second);
-      double phiD = (phis - phif);
-
-      if (std::abs(phiD) > M_PI and phif * phis < 0.) {
-        phiD += phiD < 0. ? 2 * M_PI : -2 * M_PI;
-      }
-
-      double rf = Acts::VectorHelpers::perp(first);
-      double rs = Acts::VectorHelpers::perp(second);
-      double zf = first.z();
-      double zs = second.z();
-
-      phiD *= oseg;
-      double rD = (rs - rf) * oseg;
-      double zD = (zs - zf) * oseg;
-
-      for (unsigned int is = 0; is < nSegments + 1; ++is) {
-        double r = rf + is * rD;
-        double phi = phif + is * phiD;
-        double z = zf + is * zD;
-        completed3D.push_back(Acts::Vector3D(r * cos(phi), r * sin(phi), z) -
-                              os);
-      }
-      // Reassing the global points
-      return completed3D;
-    };
-    // Fill the bows
-    auto completedBow1 = fillBow(gVertices[0], gVertices[1]);
-    auto completedBow2 = fillBow(gVertices[2], gVertices[3]);
-    // Clear the vertices
-    gVertices = completedBow1;
-    gVertices.insert(gVertices.end(), completedBow2.begin(),
-                     completedBow2.end());
-  }
-
-  // Write the vertices
-  for (const auto& gv : gVertices) {
-    stream << "v " << scalor * gv.x() << " " << scalor * gv.y() << " "
-           << scalor * gv.z() << std::endl;
-  }
-  // Write the faces
-  stream << "f";
-  for (unsigned int iv = 0; iv < gVertices.size(); ++iv) {
-    stream << " " << cvc + iv + 1;
-  }
-  stream << std::endl;
-}
-
 /// Unit tests for AnnulusBounds constrcuctors
 BOOST_AUTO_TEST_CASE(AnnulusBoundsConstruction) {
   //
@@ -115,56 +41,6 @@ BOOST_AUTO_TEST_CASE(AnnulusBoundsConstruction) {
 
   AnnulusBounds copied(original);
   BOOST_CHECK_EQUAL(copied.type(), SurfaceBounds::Annulus);
-
-  if (writeObj) {
-    std::ofstream rOut;
-    rOut.open("AnnulusDisc.obj");
-    unsigned int nSegments = 72;
-    double phiStep = 2 * M_PI / nSegments;
-    for (unsigned int is = 0; is < nSegments; ++is) {
-      double phi = is * phiStep - M_PI;
-      double cphi = std::cos(phi);
-      double sphi = std::sin(phi);
-      rOut << "v " << minRadius * cphi << " " << minRadius * sphi << " 0."
-           << std::endl;
-      rOut << "v " << maxRadius * cphi << " " << maxRadius * sphi << " 0."
-           << std::endl;
-    }
-    for (unsigned int il = 1; il < 2 * nSegments; ++il) {
-      rOut << "l " << il << " " << il + 2 << std::endl;
-    }
-    rOut.close();
-
-    std::ofstream abOut;
-    abOut.open("AnnulusBoundTests.obj");
-    writeAnnulusDiscObj(abOut, 1., 12, Transform3D::Identity(),
-                        original.vertices(),
-                        Vector3D(offset.x(), offset.y(), 0.));
-    abOut << std::endl;
-    abOut.close();
-
-    std::ofstream linesOut;
-    linesOut.open("AnnulusLines.obj");
-    auto writeVertex = [&](const Vector2D& v,
-                           const Vector2D& o = Vector2D(0., 0.)) -> void {
-      linesOut << "v " << v.x() - o.x() << " " << v.y() - o.y() << " 0"
-               << std::endl;
-    };
-
-    // write extra lines
-    auto vertices = original.vertices();
-    Vector2D sideA = vertices[0] - vertices[3];
-    Vector2D sideB = vertices[1] - vertices[2];
-
-    writeVertex(-offset);
-    writeVertex(vertices[3] - 10 * sideA, offset);
-    writeVertex(vertices[3] + 11 * sideA, offset);
-    writeVertex(vertices[2] - 10 * sideB, offset);
-    writeVertex(vertices[2] + 11 * sideB, offset);
-    linesOut << "l 2 3" << std::endl;
-    linesOut << "l 4 5" << std::endl;
-    linesOut.close();
-  }
 }
 
 /// Unit tests for AnnulusBounds properties
@@ -208,12 +84,6 @@ BOOST_AUTO_TEST_CASE(AnnulusBoundsProperties) {
   BOOST_CHECK(aBounds.insideRadialBounds(9.));
   BOOST_CHECK(!aBounds.insideRadialBounds(18.));
 
-  /// Check binning value
-  // double binningValueR();
-
-  /// Check the value in phi
-  // CHECK_CLOSE_ABS(0.5*(minPhi+maxPhi), aBounds.binningValuePhi(), 1e-6);
-
   /// Test rMin
   BOOST_CHECK_EQUAL(aBounds.rMin(), minRadius);
   //
@@ -224,15 +94,6 @@ BOOST_AUTO_TEST_CASE(AnnulusBoundsProperties) {
   //
   /// Test phiMax
   BOOST_CHECK_EQUAL(aBounds.rMax(), maxRadius);
-
-  if (writeObj) {
-    std::ofstream testOut;
-    testOut.open("AnnulusTestPoints.obj");
-    for (const auto& vc : testPoints) {
-      testOut << "v " << vc.x() << " " << vc.y() << " " << 0 << std::endl;
-    }
-    testOut.close();
-  }
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Surfaces/CMakeLists.txt b/Tests/UnitTests/Core/Surfaces/CMakeLists.txt
index 90fa4a48077486b8cf5c2c1d99e7cef8b7f47df2..fcc253d910a9f041da85ad5332e485d8de9d9af9 100644
--- a/Tests/UnitTests/Core/Surfaces/CMakeLists.txt
+++ b/Tests/UnitTests/Core/Surfaces/CMakeLists.txt
@@ -2,12 +2,12 @@ add_unittest(BoundaryCheckTests BoundaryCheckTests.cpp)
 add_unittest(AnnulusBoundsTests AnnulusBoundsTests.cpp)
 add_unittest(ConeBoundsTests ConeBoundsTests.cpp)
 add_unittest(ConeSurfaceTests ConeSurfaceTests.cpp)
-add_unittest(ConvexPolygonBoundsTest ConvexPolygonBoundsTest.cpp)
+add_unittest(ConvexPolygonBoundsTests ConvexPolygonBoundsTests.cpp)
 add_unittest(CylinderBoundsTests CylinderBoundsTests.cpp)
 add_unittest(CylinderSurfaceTests CylinderSurfaceTests.cpp)
 add_unittest(DiamondBoundsTests DiamondBoundsTests.cpp)
 add_unittest(DiscSurfaceTests DiscSurfaceTests.cpp)
-add_unittest(DiscTrapezoidalBoundsTests DiscTrapezoidalBoundsTests.cpp)
+add_unittest(DiscTrapezoidBoundsTests DiscTrapezoidBoundsTests.cpp)
 add_unittest(EllipseBoundsTests EllipseBoundsTests.cpp)
 add_unittest(InfiniteBoundsTests InfiniteBoundsTests.cpp)
 add_unittest(LineBoundsTests LineBoundsTests.cpp)
diff --git a/Tests/UnitTests/Core/Surfaces/ConvexPolygonBoundsTest.cpp b/Tests/UnitTests/Core/Surfaces/ConvexPolygonBoundsTests.cpp
similarity index 100%
rename from Tests/UnitTests/Core/Surfaces/ConvexPolygonBoundsTest.cpp
rename to Tests/UnitTests/Core/Surfaces/ConvexPolygonBoundsTests.cpp
diff --git a/Tests/UnitTests/Core/Surfaces/DiamondBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/DiamondBoundsTests.cpp
index 8b0d9a5591e408357786e84de3286d8e011daf42..20cdcdc8249f2dc1d29dbe3026a6f18e134148f2 100644
--- a/Tests/UnitTests/Core/Surfaces/DiamondBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/DiamondBoundsTests.cpp
@@ -110,8 +110,8 @@ BOOST_AUTO_TEST_CASE(DiamondBoundsProperties) {
   /// Test vertices (does this need to be implemented in this class??
   // auto v=diamondBoundsObject.vertices();
   std::vector<Vector2D> referenceVertices{
-      {minHalfX, -halfY1}, {midHalfX, 0.},  {maxHalfX, halfY2},
-      {-maxHalfX, halfY2}, {-midHalfX, 0.}, {-minHalfX, -halfY1}};
+      {-minHalfX, -halfY1}, {minHalfX, -halfY1}, {midHalfX, 0.},
+      {maxHalfX, halfY2},   {-maxHalfX, halfY2}, {-midHalfX, 0.}};
   const auto& actualVertices = diamondBoundsObject.vertices();
   BOOST_CHECK_EQUAL_COLLECTIONS(actualVertices.cbegin(), actualVertices.cend(),
                                 referenceVertices.cbegin(),
diff --git a/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp b/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp
index 459938f914d0ca4cdb2703c10538950e4bce9227..4557d176e29938b7bf5e23a880fa2e89068f9f13 100644
--- a/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp
@@ -14,7 +14,7 @@
 
 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
 #include "Acts/Surfaces/DiscSurface.hpp"
-#include "Acts/Surfaces/DiscTrapezoidalBounds.hpp"
+#include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
 #include "Acts/Surfaces/RadialBounds.hpp"
 #include "Acts/Tests/CommonHelpers/DetectorElementStub.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
diff --git a/Tests/UnitTests/Core/Surfaces/DiscTrapezoidBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/DiscTrapezoidBoundsTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fec8a40f19eb02760ef368deb0f040c504eb163b
--- /dev/null
+++ b/Tests/UnitTests/Core/Surfaces/DiscTrapezoidBoundsTests.cpp
@@ -0,0 +1,129 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/data/test_case.hpp>
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <limits>
+
+#include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+
+namespace Acts {
+
+namespace Test {
+BOOST_AUTO_TEST_SUITE(Surfaces)
+/// Unit tests for DiscTrapezoidBounds constrcuctors
+BOOST_AUTO_TEST_CASE(DiscTrapezoidBoundsConstruction) {
+  double minHalfX(1.0), maxHalfX(5.0), rMin(2.0), rMax(6.0), averagePhi(0.0),
+      stereo(0.1);
+  // test default construction
+  // DiscTrapezoidBounds defaultConstructedDiscTrapezoidBounds;  should be
+  // deleted
+  //
+  /// Test construction with dimensions and default stereo
+  BOOST_CHECK_EQUAL(
+      DiscTrapezoidBounds(minHalfX, maxHalfX, rMin, rMax, averagePhi).type(),
+      SurfaceBounds::DiscTrapezoidal);
+  //
+  /// Test construction with all dimensions
+  BOOST_CHECK_EQUAL(
+      DiscTrapezoidBounds(minHalfX, maxHalfX, rMin, rMax, averagePhi, stereo)
+          .type(),
+      SurfaceBounds::DiscTrapezoidal);
+  //
+  /// Copy constructor
+  DiscTrapezoidBounds original(minHalfX, maxHalfX, rMin, rMax, averagePhi);
+  DiscTrapezoidBounds copied(original);
+  BOOST_CHECK_EQUAL(copied.type(), SurfaceBounds::DiscTrapezoidal);
+}
+
+/// Unit tests for DiscTrapezoidBounds properties
+BOOST_AUTO_TEST_CASE(DiscTrapezoidBoundsProperties) {
+  double minHalfX(1.0), maxHalfX(5.0), rMin(2.0), rMax(6.0),
+      averagePhi(0.0) /*, stereo(0.1)*/;
+  /// Test clone
+  DiscTrapezoidBounds DiscTrapezoidBoundsObject(minHalfX, maxHalfX, rMin, rMax,
+                                                averagePhi);
+  auto pClonedDiscTrapezoidBounds = DiscTrapezoidBoundsObject.clone();
+  BOOST_CHECK_NE(pClonedDiscTrapezoidBounds, nullptr);
+  delete pClonedDiscTrapezoidBounds;
+  //
+  /// Test type() (redundant; already used in constructor confirmation)
+  BOOST_CHECK_EQUAL(DiscTrapezoidBoundsObject.type(),
+                    SurfaceBounds::DiscTrapezoidal);
+  //
+  /// Test distanceToBoundary
+  Vector2D origin(0., 0.);
+  Vector2D outside(30., 0.);
+  Vector2D inSurface(2., 0.0);
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.distanceToBoundary(origin), 2.0,
+                  1e-6);
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.distanceToBoundary(outside), 24.0,
+                  1e-6);
+  //
+  /// Test dump
+  boost::test_tools::output_test_stream dumpOuput;
+  DiscTrapezoidBoundsObject.toStream(dumpOuput);
+  BOOST_CHECK(dumpOuput.is_equal(
+      "Acts::DiscTrapezoidBounds:  (innerRadius, outerRadius, hMinX, "
+      "hMaxX, hlengthY, hPhiSector, averagePhi, rCenter, stereo) = "
+      "(2.0000000, 6.0000000, 1.0000000, 5.0000000, 0.7922870, 0.9851108, "
+      "0.0000000, 2.5243378, 0.0000000)"));
+  //
+  /// Test inside
+  BOOST_CHECK(DiscTrapezoidBoundsObject.inside(inSurface, BoundaryCheck(true)));
+  BOOST_CHECK(!DiscTrapezoidBoundsObject.inside(outside, BoundaryCheck(true)));
+  //
+  /// Test rMin
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.rMin(), rMin, 1e-6);
+  //
+  /// Test rMax
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.rMax(), rMax, 1e-6);
+  //
+  /// Test averagePhi
+  CHECK_SMALL(DiscTrapezoidBoundsObject.averagePhi(), 1e-9);
+  //
+  /// Test rCenter (redundant; not configurable)
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.rCenter(), 2.524337798, 1e-6);
+  //
+  /// Test halfPhiSector (redundant; not configurable)
+  CHECK_SMALL(DiscTrapezoidBoundsObject.stereo(), 1e-6);
+  //
+  /// Test minHalflengthX
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.minHalflengthX(), minHalfX, 1e-6);
+  //
+  /// Test maxHalflengthX
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.maxHalflengthX(), maxHalfX, 1e-6);
+  //
+  /// Test halflengthY
+  CHECK_CLOSE_REL(DiscTrapezoidBoundsObject.halflengthY(), 0.792286991, 1e-6);
+}
+/// Unit test for testing DiscTrapezoidBounds assignment
+BOOST_AUTO_TEST_CASE(DiscTrapezoidBoundsAssignment) {
+  double minHalfX(1.0), maxHalfX(5.0), rMin(2.0), rMax(6.0), averagePhi(0.0),
+      stereo(0.1);
+  DiscTrapezoidBounds DiscTrapezoidBoundsObject(minHalfX, maxHalfX, rMin, rMax,
+                                                averagePhi, stereo);
+  // operator == not implemented in this class
+  //
+  /// Test assignment
+  DiscTrapezoidBounds assignedDiscTrapezoidBoundsObject(2.1, 6.6, 3.4, 4.2,
+                                                        33.);
+  assignedDiscTrapezoidBoundsObject = DiscTrapezoidBoundsObject;
+  BOOST_CHECK_EQUAL(assignedDiscTrapezoidBoundsObject,
+                    DiscTrapezoidBoundsObject);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}  // namespace Test
+
+}  // namespace Acts
diff --git a/Tests/UnitTests/Core/Surfaces/DiscTrapezoidalBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/DiscTrapezoidalBoundsTests.cpp
index 9b3ab57e814b8f6423a656f40d009045bd98eb01..82deb4d79d28fb46213abda4777f73bcf82d69d8 100644
--- a/Tests/UnitTests/Core/Surfaces/DiscTrapezoidalBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/DiscTrapezoidalBoundsTests.cpp
@@ -12,7 +12,7 @@
 
 #include <limits>
 
-#include "Acts/Surfaces/DiscTrapezoidalBounds.hpp"
+#include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
@@ -20,108 +20,106 @@ namespace Acts {
 
 namespace Test {
 BOOST_AUTO_TEST_SUITE(Surfaces)
-/// Unit tests for DiscTrapezoidalBounds constrcuctors
-BOOST_AUTO_TEST_CASE(DiscTrapezoidalBoundsConstruction) {
+/// Unit tests for DiscTrapezoidBounds constrcuctors
+BOOST_AUTO_TEST_CASE(DiscTrapezoidBoundsConstruction) {
   double minHalfX(1.0), maxHalfX(5.0), rMin(2.0), rMax(6.0), averagePhi(0.0),
       stereo(0.1);
   // test default construction
-  // DiscTrapezoidalBounds defaultConstructedDiscTrapezoidalBounds;  should be
+  // DiscTrapezoidBounds defaultConstructedDiscTrapezoidBounds;  should be
   // deleted
   //
   /// Test construction with dimensions and default stereo
   BOOST_CHECK_EQUAL(
-      DiscTrapezoidalBounds(minHalfX, maxHalfX, rMin, rMax, averagePhi).type(),
+      DiscTrapezoidBounds(minHalfX, maxHalfX, rMin, rMax, averagePhi).type(),
       SurfaceBounds::DiscTrapezoidal);
   //
   /// Test construction with all dimensions
   BOOST_CHECK_EQUAL(
-      DiscTrapezoidalBounds(minHalfX, maxHalfX, rMin, rMax, averagePhi, stereo)
+      DiscTrapezoidBounds(minHalfX, maxHalfX, rMin, rMax, averagePhi, stereo)
           .type(),
       SurfaceBounds::DiscTrapezoidal);
   //
   /// Copy constructor
-  DiscTrapezoidalBounds original(minHalfX, maxHalfX, rMin, rMax, averagePhi);
-  DiscTrapezoidalBounds copied(original);
+  DiscTrapezoidBounds original(minHalfX, maxHalfX, rMin, rMax, averagePhi);
+  DiscTrapezoidBounds copied(original);
   BOOST_CHECK_EQUAL(copied.type(), SurfaceBounds::DiscTrapezoidal);
 }
 
-/// Unit tests for DiscTrapezoidalBounds properties
-BOOST_AUTO_TEST_CASE(DiscTrapezoidalBoundsProperties) {
+/// Unit tests for DiscTrapezoidBounds properties
+BOOST_AUTO_TEST_CASE(DiscTrapezoidBoundsProperties) {
   double minHalfX(1.0), maxHalfX(5.0), rMin(2.0), rMax(6.0),
       averagePhi(0.0) /*, stereo(0.1)*/;
   /// Test clone
-  DiscTrapezoidalBounds discTrapezoidalBoundsObject(minHalfX, maxHalfX, rMin,
-                                                    rMax, averagePhi);
-  auto pClonedDiscTrapezoidalBounds = discTrapezoidalBoundsObject.clone();
-  BOOST_CHECK_NE(pClonedDiscTrapezoidalBounds, nullptr);
-  delete pClonedDiscTrapezoidalBounds;
+  DiscTrapezoidBounds discTrapezoidBoundsObject(minHalfX, maxHalfX, rMin, rMax,
+                                                averagePhi);
+  auto pClonedDiscTrapezoidBounds = discTrapezoidBoundsObject.clone();
+  BOOST_CHECK_NE(pClonedDiscTrapezoidBounds, nullptr);
+  delete pClonedDiscTrapezoidBounds;
   //
   /// Test type() (redundant; already used in constructor confirmation)
-  BOOST_CHECK_EQUAL(discTrapezoidalBoundsObject.type(),
+  BOOST_CHECK_EQUAL(discTrapezoidBoundsObject.type(),
                     SurfaceBounds::DiscTrapezoidal);
   //
   /// Test distanceToBoundary
   Vector2D origin(0., 0.);
   Vector2D outside(30., 0.);
   Vector2D inSurface(2., 0.0);
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.distanceToBoundary(origin), 2.0,
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.distanceToBoundary(origin), 2.0,
                   1e-6);
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.distanceToBoundary(outside), 24.0,
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.distanceToBoundary(outside), 24.0,
                   1e-6);
   //
   /// Test dump
   boost::test_tools::output_test_stream dumpOuput;
-  discTrapezoidalBoundsObject.toStream(dumpOuput);
+  discTrapezoidBoundsObject.toStream(dumpOuput);
   BOOST_CHECK(dumpOuput.is_equal(
-      "Acts::DiscTrapezoidalBounds:  (innerRadius, outerRadius, hMinX, "
+      "Acts::DiscTrapezoidBounds:  (innerRadius, outerRadius, hMinX, "
       "hMaxX, hlengthY, hPhiSector, averagePhi, rCenter, stereo) = "
       "(2.0000000, 6.0000000, 1.0000000, 5.0000000, 0.7922870, 0.9851108, "
       "0.0000000, 2.5243378, 0.0000000)"));
   //
   /// Test inside
-  BOOST_CHECK(
-      discTrapezoidalBoundsObject.inside(inSurface, BoundaryCheck(true)));
-  BOOST_CHECK(
-      !discTrapezoidalBoundsObject.inside(outside, BoundaryCheck(true)));
+  BOOST_CHECK(discTrapezoidBoundsObject.inside(inSurface, BoundaryCheck(true)));
+  BOOST_CHECK(!discTrapezoidBoundsObject.inside(outside, BoundaryCheck(true)));
   //
   /// Test rMin
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.rMin(), rMin, 1e-6);
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.rMin(), rMin, 1e-6);
   //
   /// Test rMax
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.rMax(), rMax, 1e-6);
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.rMax(), rMax, 1e-6);
   //
   /// Test averagePhi
-  CHECK_SMALL(discTrapezoidalBoundsObject.averagePhi(), 1e-9);
+  CHECK_SMALL(discTrapezoidBoundsObject.averagePhi(), 1e-9);
   //
   /// Test rCenter (redundant; not configurable)
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.rCenter(), 2.524337798, 1e-6);
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.rCenter(), 2.524337798, 1e-6);
   //
   /// Test halfPhiSector (redundant; not configurable)
-  CHECK_SMALL(discTrapezoidalBoundsObject.stereo(), 1e-6);
+  CHECK_SMALL(discTrapezoidBoundsObject.stereo(), 1e-6);
   //
   /// Test minHalflengthX
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.minHalflengthX(), minHalfX, 1e-6);
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.minHalflengthX(), minHalfX, 1e-6);
   //
   /// Test maxHalflengthX
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.maxHalflengthX(), maxHalfX, 1e-6);
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.maxHalflengthX(), maxHalfX, 1e-6);
   //
   /// Test halflengthY
-  CHECK_CLOSE_REL(discTrapezoidalBoundsObject.halflengthY(), 0.792286991, 1e-6);
+  CHECK_CLOSE_REL(discTrapezoidBoundsObject.halflengthY(), 0.792286991, 1e-6);
 }
-/// Unit test for testing DiscTrapezoidalBounds assignment
-BOOST_AUTO_TEST_CASE(DiscTrapezoidalBoundsAssignment) {
+/// Unit test for testing DiscTrapezoidBounds assignment
+BOOST_AUTO_TEST_CASE(DiscTrapezoidBoundsAssignment) {
   double minHalfX(1.0), maxHalfX(5.0), rMin(2.0), rMax(6.0), averagePhi(0.0),
       stereo(0.1);
-  DiscTrapezoidalBounds discTrapezoidalBoundsObject(minHalfX, maxHalfX, rMin,
-                                                    rMax, averagePhi, stereo);
+  DiscTrapezoidBounds discTrapezoidBoundsObject(minHalfX, maxHalfX, rMin, rMax,
+                                                averagePhi, stereo);
   // operator == not implemented in this class
   //
   /// Test assignment
-  DiscTrapezoidalBounds assignedDiscTrapezoidalBoundsObject(2.1, 6.6, 3.4, 4.2,
-                                                            33.);
-  assignedDiscTrapezoidalBoundsObject = discTrapezoidalBoundsObject;
-  BOOST_CHECK_EQUAL(assignedDiscTrapezoidalBoundsObject,
-                    discTrapezoidalBoundsObject);
+  DiscTrapezoidBounds assignedDiscTrapezoidBoundsObject(2.1, 6.6, 3.4, 4.2,
+                                                        33.);
+  assignedDiscTrapezoidBoundsObject = discTrapezoidBoundsObject;
+  BOOST_CHECK_EQUAL(assignedDiscTrapezoidBoundsObject,
+                    discTrapezoidBoundsObject);
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Surfaces/EllipseBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/EllipseBoundsTests.cpp
index 87dafd0bedeb4173935267e9b6b75bd795ff039d..a6b384ceac230db0b25ba43eb8d6e4e536da0ecb 100644
--- a/Tests/UnitTests/Core/Surfaces/EllipseBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/EllipseBoundsTests.cpp
@@ -82,11 +82,12 @@ BOOST_AUTO_TEST_CASE(EllipseBoundsProperties) {
   BOOST_CHECK_EQUAL(ellipseBoundsObject.averagePhi(), averagePhi);
   //
   /// Test vertices
-  std::vector<Vector2D> expectedVertices{{15, 0}, {0, 20}, {-15, 0}, {0, -20}};
-  const auto& actualVertices = ellipseBoundsObject.vertices();
-  BOOST_CHECK_EQUAL_COLLECTIONS(actualVertices.cbegin(), actualVertices.cend(),
-                                expectedVertices.cbegin(),
-                                expectedVertices.cend());
+  // std::vector<Vector2D> expectedVertices{{15, 0}, {0, 20}, {-15, 0}, {0,
+  // -20}}; const auto& actualVertices = ellipseBoundsObject.vertices(4);
+  // BOOST_CHECK_EQUAL_COLLECTIONS(actualVertices.cbegin(),
+  // actualVertices.cend(),
+  //                              expectedVertices.cbegin(),
+  //                              expectedVertices.cend());
   //
   /// Test boundingBox
   BOOST_CHECK_EQUAL(ellipseBoundsObject.boundingBox(),
diff --git a/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp b/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c1bdcfff233a5ca8aba8cf32e9dff1e03117275
--- /dev/null
+++ b/Tests/UnitTests/Core/Surfaces/PolyhedronSurfacesTests.cpp
@@ -0,0 +1,587 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/data/test_case.hpp>
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+// Helper
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Tests/CommonHelpers/ObjTestWriter.hpp"
+
+// The class to test
+#include "Acts/Geometry/Extent.hpp"
+#include "Acts/Geometry/Polyhedron.hpp"
+
+// Cone surface
+#include "Acts/Surfaces/ConeBounds.hpp"
+#include "Acts/Surfaces/ConeSurface.hpp"
+
+// Cylinder surface
+#include "Acts/Surfaces/CylinderBounds.hpp"
+#include "Acts/Surfaces/CylinderSurface.hpp"
+
+// Disc Surface
+#include "Acts/Surfaces/AnnulusBounds.hpp"
+#include "Acts/Surfaces/DiscSurface.hpp"
+#include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
+#include "Acts/Surfaces/RadialBounds.hpp"
+
+// Plane Surface
+#include "Acts/Surfaces/ConvexPolygonBounds.hpp"
+#include "Acts/Surfaces/DiamondBounds.hpp"
+#include "Acts/Surfaces/EllipseBounds.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Surfaces/TrapezoidBounds.hpp"
+
+// Straw Surface
+#include "Acts/Surfaces/LineBounds.hpp"
+#include "Acts/Surfaces/StrawSurface.hpp"
+
+#include "Acts/Utilities/ObjHelper.hpp"
+#include "Acts/Utilities/Units.hpp"
+
+#include <fstream>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+namespace Acts {
+
+using namespace UnitLiterals;
+
+using IdentifiedPolyderon = std::tuple<std::string, bool, Polyhedron>;
+
+namespace Test {
+
+// Create a test context
+GeometryContext tgContext = GeometryContext();
+
+std::vector<std::tuple<std::string, bool, unsigned int>> testModes = {
+    {"", false, 72}, {"Triangulate", true, 72}, {"Extremas", false, 1}};
+
+auto transform = std::make_shared<Transform3D>(Transform3D::Identity());
+
+BOOST_AUTO_TEST_SUITE(Surfaces)
+
+/// Unit tests for Cone Surfaces
+BOOST_AUTO_TEST_CASE(ConeSurfacePolyhedrons) {
+  std::vector<IdentifiedPolyderon> testTypes;
+
+  double hzpmin = 10_mm;
+  double hzpos = 35_mm;
+  double hzneg = -20_mm;
+  double alpha = 0.234;
+  double phiSector = 0.358;
+  ObjTestWriter::writeSectorPlanesObj("ConeSectorPlanes", phiSector, 0., hzpos,
+                                      hzpos);
+
+  for (const auto& mode : testModes) {
+    unsigned int segments = std::get<unsigned int>(mode);
+    std::string modename = std::get<std::string>(mode);
+    bool modetrg = std::get<bool>(mode);
+
+    /// The full cone on one side
+    auto cone = std::make_shared<ConeBounds>(alpha, 0_mm, hzpos);
+    auto oneCone = Surface::makeShared<ConeSurface>(transform, cone);
+    auto oneConePh = oneCone->polyhedronRepresentation(tgContext, segments);
+    size_t expectedFaces = segments < 4 ? 4 : segments;
+    BOOST_CHECK_EQUAL(oneConePh.faces.size(), expectedFaces);
+    BOOST_CHECK_EQUAL(oneConePh.vertices.size(), expectedFaces + 1);
+    // Check the extent in space
+    double r = hzpos * std::tan(alpha);
+    auto extent = oneConePh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0_mm, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0_mm, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hzpos, 1e-6);
+    testTypes.push_back({"ConeOneFull" + modename, modetrg, oneConePh});
+
+    /// The full cone on one side
+    auto conePiece = std::make_shared<ConeBounds>(alpha, hzpmin, hzpos);
+    auto oneConePiece = Surface::makeShared<ConeSurface>(transform, conePiece);
+    auto oneConePiecePh =
+        oneConePiece->polyhedronRepresentation(tgContext, segments);
+    expectedFaces = segments < 4 ? 4 : segments;
+    // Check the extent in space
+    double rmin = hzpmin * std::tan(alpha);
+    extent = oneConePiecePh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, rmin, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, hzpmin, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hzpos, 1e-6);
+    testTypes.push_back(
+        {"ConeOnePieceFull" + modename, modetrg, oneConePiecePh});
+
+    // The full cone on both sides
+    auto coneBoth = std::make_shared<ConeBounds>(alpha, hzneg, hzpos);
+    auto twoCones = Surface::makeShared<ConeSurface>(transform, coneBoth);
+    auto twoConesPh = twoCones->polyhedronRepresentation(tgContext, segments);
+    expectedFaces = segments < 4 ? 8 : 2 * segments;
+    BOOST_CHECK_EQUAL(twoConesPh.faces.size(), expectedFaces);
+    BOOST_CHECK_EQUAL(twoConesPh.vertices.size(), expectedFaces + 1);
+    extent = twoConesPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0_mm, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, hzneg, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hzpos, 1e-6);
+    testTypes.push_back({"ConesTwoFull" + modename, modetrg, twoConesPh});
+
+    // A centered sectoral cone on both sides
+    auto sectoralBoth =
+        std::make_shared<ConeBounds>(alpha, hzneg, hzpos, phiSector, 0.);
+    auto sectoralCones =
+        Surface::makeShared<ConeSurface>(transform, sectoralBoth);
+    auto sectoralConesPh =
+        sectoralCones->polyhedronRepresentation(tgContext, segments);
+    extent = sectoralConesPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0_mm, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, hzneg, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hzpos, 1e-6);
+    testTypes.push_back({"ConesSectoral" + modename, modetrg, sectoralConesPh});
+  }
+  ObjTestWriter::writeObj(testTypes);
+}
+
+/// Unit tests for Cylinder Surfaces
+BOOST_AUTO_TEST_CASE(CylinderSurfacePolyhedrons) {
+  double r = 25_mm;
+  double hZ = 35_mm;
+
+  double phiSector = 0.458;
+  double averagePhi = -1.345;
+  ObjTestWriter::writeSectorPlanesObj("CylinderCentralSectorPlanes", phiSector,
+                                      0., 1.5 * r, 1.5 * hZ);
+  ObjTestWriter::writeSectorPlanesObj("CylinderShiftedSectorPlanes", phiSector,
+                                      averagePhi, 1.5 * r, 1.5 * hZ);
+
+  std::vector<IdentifiedPolyderon> testTypes;
+
+  for (const auto& mode : testModes) {
+    unsigned int segments = std::get<unsigned int>(mode);
+    std::string modename = std::get<std::string>(mode);
+    bool modetrg = std::get<bool>(mode);
+
+    size_t expectedFaces = segments < 4 ? 4 : segments;
+    size_t expectedVertices = segments < 4 ? 8 : 2 * segments;
+
+    /// The full cone on one side
+    auto cylinder = std::make_shared<CylinderBounds>(r, hZ);
+    auto fullCylinder =
+        Surface::makeShared<CylinderSurface>(transform, cylinder);
+    auto fullCylinderPh =
+        fullCylinder->polyhedronRepresentation(tgContext, segments);
+
+    BOOST_CHECK_EQUAL(fullCylinderPh.faces.size(), expectedFaces);
+    BOOST_CHECK_EQUAL(fullCylinderPh.vertices.size(), expectedVertices);
+    // Check the extent in space
+    auto extent = fullCylinderPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, -hZ, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hZ, 1e-6);
+    testTypes.push_back({"CylinderFull" + modename, modetrg, fullCylinderPh});
+
+    /// The full cone on one side
+    auto sectorCentered = std::make_shared<CylinderBounds>(r, phiSector, hZ);
+    auto centerSectoredCylinder =
+        Surface::makeShared<CylinderSurface>(transform, sectorCentered);
+    auto centerSectoredCylinderPh =
+        centerSectoredCylinder->polyhedronRepresentation(tgContext, segments);
+
+    // Check the extent in space
+    extent = centerSectoredCylinderPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, r * std::cos(phiSector), 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -r * std::sin(phiSector), 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, r * std::sin(phiSector), 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, -hZ, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hZ, 1e-6);
+    testTypes.push_back({"CylinderSectorCentered" + modename, modetrg,
+                         centerSectoredCylinderPh});
+
+    /// The full cone on one side
+    auto sectorShifted =
+        std::make_shared<CylinderBounds>(r, averagePhi, phiSector, hZ);
+    auto shiftedSectoredCylinder =
+        Surface::makeShared<CylinderSurface>(transform, sectorShifted);
+    auto shiftedSectoredCylinderPh =
+        shiftedSectoredCylinder->polyhedronRepresentation(tgContext, segments);
+
+    // Check the extent in space
+    extent = shiftedSectoredCylinderPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, r, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, -hZ, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, hZ, 1e-6);
+    testTypes.push_back({"CylinderSectorShifted" + modename, modetrg,
+                         shiftedSectoredCylinderPh});
+  }
+
+  ObjTestWriter::writeObj(testTypes);
+}
+
+/// Unit tests for Disc Surfaces
+BOOST_AUTO_TEST_CASE(DiscSurfacePolyhedrons) {
+  std::vector<IdentifiedPolyderon> testTypes;
+
+  double innerR = 10_mm;
+  double outerR = 25_mm;
+
+  double phiSector = 0.345;
+  double averagePhi = -1.0;
+
+  double cphi = std::cos(phiSector);
+  double sphi = std::sin(phiSector);
+
+  std::pair<Vector3D, Vector3D> lineA = {
+      Vector3D(0., 0., 0.), Vector3D(outerR * cphi, outerR * sphi, 0.)};
+  std::pair<Vector3D, Vector3D> lineB = {
+      Vector3D(0., 0., 0.), Vector3D(outerR * cphi, -outerR * sphi, 0.)};
+  ObjTestWriter::writeSectorLinesObj("DiscSectorLines", lineA, lineB);
+
+  double minPhi = averagePhi - phiSector;
+  double maxPhi = averagePhi + phiSector;
+  lineA = {Vector3D(0., 0., 0.),
+           Vector3D(outerR * std::cos(minPhi), outerR * std::sin(minPhi), 0.)};
+  lineB = {Vector3D(0., 0., 0.),
+           Vector3D(outerR * std::cos(maxPhi), outerR * std::sin(maxPhi), 0.)};
+  ObjTestWriter::writeSectorLinesObj("DiscSectorLinesShifted", lineA, lineB);
+
+  for (const auto& mode : testModes) {
+    unsigned int segments = std::get<unsigned int>(mode);
+    std::string modename = std::get<std::string>(mode);
+    bool modetrg = std::get<bool>(mode);
+
+    // Full disc
+    auto disc = std::make_shared<RadialBounds>(0_mm, outerR);
+    auto fullDisc = Surface::makeShared<DiscSurface>(transform, disc);
+    auto fullDiscPh = fullDisc->polyhedronRepresentation(tgContext, segments);
+
+    unsigned int expectedVertices = segments > 4 ? segments : 4;
+    unsigned int expectedFaces = 1;
+
+    BOOST_CHECK_EQUAL(fullDiscPh.faces.size(), expectedFaces);
+    BOOST_CHECK_EQUAL(fullDiscPh.vertices.size(), expectedVertices);
+
+    auto extent = fullDiscPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+
+    testTypes.push_back({"DiscFull" + modename, modetrg, fullDiscPh});
+
+    // Ring disc
+    auto radial = std::make_shared<RadialBounds>(innerR, outerR);
+    auto radialDisc = Surface::makeShared<DiscSurface>(transform, radial);
+    auto radialPh = radialDisc->polyhedronRepresentation(tgContext, segments);
+    extent = radialPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, innerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back({"DiscRing" + modename, modetrg, radialPh});
+
+    // Sectoral disc - around 0.
+    auto sector = std::make_shared<RadialBounds>(0., outerR, phiSector);
+    auto sectorDisc = Surface::makeShared<DiscSurface>(transform, sector);
+    auto sectorPh = sectorDisc->polyhedronRepresentation(tgContext, segments);
+    extent = sectorPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -outerR * std::sin(phiSector),
+                    1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, outerR * std::sin(phiSector),
+                    1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back({"DiscSectorCentered" + modename, modetrg, sectorPh});
+
+    // Sectoral ring - around 0.
+    auto sectorRing = std::make_shared<RadialBounds>(innerR, outerR, phiSector);
+    auto sectorRingDisc =
+        Surface::makeShared<DiscSurface>(transform, sectorRing);
+    auto sectorRingDiscPh =
+        sectorRingDisc->polyhedronRepresentation(tgContext, segments);
+    extent = sectorRingDiscPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, innerR * std::cos(phiSector),
+                    1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -outerR * std::sin(phiSector),
+                    1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, outerR * std::sin(phiSector),
+                    1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, innerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back(
+        {"DiscRingSectorCentered" + modename, modetrg, sectorRingDiscPh});
+
+    // Sectoral disc - shifted
+    auto sectorRingShifted =
+        std::make_shared<RadialBounds>(innerR, outerR, averagePhi, phiSector);
+    auto sectorRingDiscShifted =
+        Surface::makeShared<DiscSurface>(transform, sectorRingShifted);
+    auto sectorRingDiscShiftedPh =
+        sectorRingDiscShifted->polyhedronRepresentation(tgContext, segments);
+    extent = sectorRingDiscShiftedPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, innerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back(
+        {"DiscRingSectorShifted" + modename, modetrg, sectorRingDiscShiftedPh});
+
+    // Trapezoid for a disc
+    double halfXmin = 10_mm;
+    double halfXmax = 20_mm;
+    auto trapezoidDisc = std::make_shared<DiscTrapezoidBounds>(
+        halfXmin, halfXmax, innerR, outerR, 0.);
+    auto trapezoidDiscSf =
+        Surface::makeShared<DiscSurface>(transform, trapezoidDisc);
+    auto trapezoidDiscSfPh =
+        trapezoidDiscSf->polyhedronRepresentation(tgContext, segments);
+    extent = trapezoidDiscSfPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, innerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back(
+        {"DiscTrapezoidCentered" + modename, modetrg, trapezoidDiscSfPh});
+
+    auto trapezoidDiscShifted = std::make_shared<DiscTrapezoidBounds>(
+        halfXmin, halfXmax, innerR, outerR, averagePhi);
+    auto trapezoidDiscShiftedSf =
+        Surface::makeShared<DiscSurface>(transform, trapezoidDiscShifted);
+    auto trapezoidDiscShiftedSfPh =
+        trapezoidDiscShiftedSf->polyhedronRepresentation(tgContext, segments);
+    extent = trapezoidDiscShiftedSfPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, innerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, outerR, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back(
+        {"DiscTrapezoidShifted" + modename, modetrg, trapezoidDiscShiftedSfPh});
+
+    double minRadius = 7.;
+    double maxRadius = 12.;
+    double minPhiA = 0.75;
+    double maxPhiA = 1.4;
+
+    Vector2D offset(-2., 2.);
+
+    auto annulus = std::make_shared<AnnulusBounds>(minRadius, maxRadius,
+                                                   minPhiA, maxPhiA, offset);
+    auto annulusDisc = Surface::makeShared<DiscSurface>(transform, annulus);
+    auto annulusDiscPh =
+        annulusDisc->polyhedronRepresentation(tgContext, segments);
+
+    testTypes.push_back(
+        {"DiscAnnulus" + modename, modetrg, trapezoidDiscShiftedSfPh});
+  }
+
+  ObjTestWriter::writeObj(testTypes);
+}
+
+/// Unit tests for Plane Surfaces
+BOOST_AUTO_TEST_CASE(PlaneSurfacePolyhedrons) {
+  std::vector<IdentifiedPolyderon> testTypes;
+
+  double rhX = 10_mm;
+  double rhY = 25_mm;
+  double shiftY = 50_mm;
+  auto rectangular = std::make_shared<RectangleBounds>(rhX, rhY);
+
+  // Special test for shifted plane to check rMin/rMax
+  Vector3D shift(0., shiftY, 0.);
+  auto shiftedTransform =
+      std::make_shared<Transform3D>(Transform3D::Identity());
+  shiftedTransform->pretranslate(shift);
+  auto shiftedPlane =
+      Surface::makeShared<PlaneSurface>(shiftedTransform, rectangular);
+  auto shiftedPh = shiftedPlane->polyhedronRepresentation(tgContext, 1);
+  auto shiftedExtent = shiftedPh.extent();
+  // Let's check the extent
+  CHECK_CLOSE_ABS(shiftedExtent.ranges[binX].first, -rhX, 1e-6);
+  CHECK_CLOSE_ABS(shiftedExtent.ranges[binX].second, rhX, 1e-6);
+  CHECK_CLOSE_ABS(shiftedExtent.ranges[binY].first, -rhY + shiftY, 1e-6);
+  CHECK_CLOSE_ABS(shiftedExtent.ranges[binY].second, rhY + shiftY, 1e-6);
+
+  for (const auto& mode : testModes) {
+    unsigned int segments = std::get<unsigned int>(mode);
+    std::string modename = std::get<std::string>(mode);
+    bool modetrg = std::get<bool>(mode);
+
+    /// Rectangular Plane
+    auto rectangularPlane =
+        Surface::makeShared<PlaneSurface>(transform, rectangular);
+    auto rectangularPh =
+        rectangularPlane->polyhedronRepresentation(tgContext, segments);
+    auto extent = rectangularPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -rhX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, rhX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -rhY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, rhY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second,
+                    std::sqrt(rhX * rhX + rhY * rhY), 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    BOOST_CHECK(rectangularPh.vertices.size() == 4);
+    BOOST_CHECK(rectangularPh.faces.size() == 1);
+    std::vector<size_t> expectedRect = {0, 1, 2, 3};
+    BOOST_CHECK(rectangularPh.faces[0] == expectedRect);
+    testTypes.push_back({"PlaneRectangle" + modename, modetrg, rectangularPh});
+
+    /// Trapezoidal Plane
+    double thX1 = 10_mm;
+    double thX2 = 25_mm;
+    double thY = 35_mm;
+
+    auto trapezoid = std::make_shared<TrapezoidBounds>(thX1, thX2, thY);
+    auto trapezoidalPlane =
+        Surface::makeShared<PlaneSurface>(transform, trapezoid);
+    auto trapezoidalPh =
+        trapezoidalPlane->polyhedronRepresentation(tgContext, segments);
+    extent = trapezoidalPh.extent();
+
+    double thX = std::max(thX1, thX2);
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -thX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, thX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -thY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, thY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second,
+                    std::sqrt(thX * thX + thY * thY), 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    BOOST_CHECK(trapezoidalPh.vertices.size() == 4);
+    BOOST_CHECK(trapezoidalPh.faces.size() == 1);
+    std::vector<size_t> expectedTra = {0, 1, 2, 3};
+    BOOST_CHECK(trapezoidalPh.faces[0] == expectedTra);
+    testTypes.push_back({"PlaneTrapezoid" + modename, modetrg, trapezoidalPh});
+
+    /// Ring-like ellispoidal Plane
+    double rMaxX = 30_mm;
+    double rMaxY = 40_mm;
+    auto ellipse = std::make_shared<EllipseBounds>(0., 0., rMaxX, rMaxY);
+    auto ellipsoidPlane = Surface::makeShared<PlaneSurface>(transform, ellipse);
+    auto ellispoidPh =
+        ellipsoidPlane->polyhedronRepresentation(tgContext, segments);
+    extent = ellispoidPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -rMaxX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, rMaxX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -rMaxY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -rMaxY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, rMaxY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+
+    testTypes.push_back({"PlaneFullEllipse" + modename, modetrg, ellispoidPh});
+
+    double rMinX = 10_mm;
+    double rMinY = 20_mm;
+    auto ellipseRing =
+        std::make_shared<EllipseBounds>(rMinX, rMaxX, rMinY, rMaxY);
+    auto ellipsoidRingPlane =
+        Surface::makeShared<PlaneSurface>(transform, ellipseRing);
+    auto ellispoidRingPh =
+        ellipsoidRingPlane->polyhedronRepresentation(tgContext, segments);
+
+    extent = ellispoidPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -rMaxX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, rMaxX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -rMaxY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, rMinX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second, rMaxY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+
+    // BOOST_CHECK(ellispoidPh.vertices.size()==72);
+    testTypes.push_back(
+        {"PlaneRingEllipse" + modename, modetrg, ellispoidRingPh});
+
+    /// ConvextPolygonBounds test
+    std::vector<Vector2D> vtxs = {
+        Vector2D(-40_mm, -10_mm), Vector2D(-10_mm, -30_mm),
+        Vector2D(30_mm, -20_mm),  Vector2D(10_mm, 20_mm),
+        Vector2D(-20_mm, 50_mm),  Vector2D(-30_mm, 30_mm)};
+
+    auto sextagon = std::make_shared<ConvexPolygonBounds<6>>(vtxs);
+    auto sextagonPlane = Surface::makeShared<PlaneSurface>(transform, sextagon);
+    auto sextagonPlanePh =
+        sextagonPlane->polyhedronRepresentation(tgContext, segments);
+    testTypes.push_back({"PlaneSextagon" + modename, modetrg, sextagonPlanePh});
+
+    /// Diamond shaped plane
+    double hMinX = 10_mm;
+    double hMedX = 20_mm;
+    double hMaxX = 15_mm;
+    double hMinY = 40_mm;
+    double hMaxY = 50_mm;
+    auto diamond =
+        std::make_shared<DiamondBounds>(hMinX, hMedX, hMaxX, hMinY, hMaxY);
+    auto diamondPlane = Surface::makeShared<PlaneSurface>(transform, diamond);
+    auto diamondPh =
+        diamondPlane->polyhedronRepresentation(tgContext, segments);
+    BOOST_CHECK(diamondPh.vertices.size() == 6);
+    BOOST_CHECK(diamondPh.faces.size() == 1);
+    extent = diamondPh.extent();
+    CHECK_CLOSE_ABS(extent.ranges[binX].first, -hMedX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binX].second, hMedX, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].first, -hMinY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binY].second, hMaxY, 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binR].second,
+                    std::sqrt(hMaxX * hMaxX + hMaxY * hMaxY), 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].first, 0., 1e-6);
+    CHECK_CLOSE_ABS(extent.ranges[binZ].second, 0., 1e-6);
+    testTypes.push_back({"PlaneDiamond" + modename, modetrg, diamondPh});
+  }
+  ObjTestWriter::writeObj(testTypes);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}  // namespace Test
+
+}  // namespace Acts
\ No newline at end of file
diff --git a/Tests/UnitTests/Core/Surfaces/RectangleBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/RectangleBoundsTests.cpp
index 7ebf7760492f9d08156ebda3e479602c690d2759..e7d4e7c3d6056f5b353d18b8ac7be466d34fe89b 100644
--- a/Tests/UnitTests/Core/Surfaces/RectangleBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/RectangleBoundsTests.cpp
@@ -82,7 +82,7 @@ BOOST_AUTO_TEST_CASE(RectangleBoundsProperties) {
   CHECK_CLOSE_ABS(rect.max(), Vector2D(halfX, halfY), 1e-6);
 
   const std::vector<Vector2D> coords = {
-      {10., -5.}, {10., 5.}, {-10., 5.}, {-10., -5.}};
+      {-10., -5.}, {10., -5.}, {10., 5.}, {-10., 5.}};
   // equality, ensure ordering is ok
   const auto& rectVertices = rect.vertices();
   BOOST_CHECK_EQUAL_COLLECTIONS(coords.cbegin(), coords.cend(),
diff --git a/Tests/UnitTests/Core/Surfaces/SurfaceStub.hpp b/Tests/UnitTests/Core/Surfaces/SurfaceStub.hpp
index 9605ffaa5fd204fcf1dbd6f6b943b140307044b1..04171db5ca32ede03d45d5e5882b7502910521e7 100644
--- a/Tests/UnitTests/Core/Surfaces/SurfaceStub.hpp
+++ b/Tests/UnitTests/Core/Surfaces/SurfaceStub.hpp
@@ -100,6 +100,16 @@ class SurfaceStub : public Surface {
   /// Simply return true to check a method can be called on a constructed object
   bool constructedOk() const { return true; }
 
+  /// Return a Polyhedron for the surfaces
+  Polyhedron polyhedronRepresentation(const GeometryContext& /*gctx*/,
+                                      size_t /*lseg */) const final {
+    std::vector<Vector3D> vertices;
+    std::vector<std::vector<size_t>> faces;
+    std::vector<std::vector<size_t>> triangularMesh;
+
+    return Polyhedron(vertices, faces, triangularMesh);
+  }
+
  private:
   /// the bounds of this surface
   std::shared_ptr<const PlanarBounds> m_bounds;
diff --git a/Tests/UnitTests/Core/Utilities/VisualizationTests.cpp b/Tests/UnitTests/Core/Utilities/VisualizationTests.cpp
index 1d7319e33c221e12b04f88d9167b576a79c68e7b..7c68dae4409e1acf0056b2a9cd36c6911a77e9df 100644
--- a/Tests/UnitTests/Core/Utilities/VisualizationTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/VisualizationTests.cpp
@@ -168,10 +168,6 @@ BOOST_AUTO_TEST_CASE(obj_output_test) {
 
   BOOST_CHECK(output.is_equal(exp));
 
-  obj.clear();
-
-  BOOST_CHECK_THROW(obj.line({1, 0, 0}, {0, 1, 0}), std::runtime_error);
-
   obj.clear();
   obj.face({{1, 0, 0}, {1, 1, 0}, {0, 1, 0}});
   output << obj;