diff --git a/Core/include/Acts/Utilities/Axis.hpp b/Core/include/Acts/Utilities/Axis.hpp
index 211938b0df7..a8601e7eff9 100644
--- a/Core/include/Acts/Utilities/Axis.hpp
+++ b/Core/include/Acts/Utilities/Axis.hpp
@@ -206,7 +206,7 @@ class Axis<AxisType::Equidistant, bdt> : public IAxis {
   ///       bounds, that is \f$l <= x < u\f$ if the value @c x lies within a
   ///       bin with lower bound @c l and upper bound @c u.
   /// @note Bin indices start at @c 1. The underflow bin has the index @c 0
-  ///       while the index <tt>nBins + 1</tt> indicates the overflow bin .
+  ///       while the index <tt>nBins + 1</tt> indicates the overflow bin.
   std::size_t getBin(double x) const final {
     return wrapBin(
         static_cast<int>(std::floor((x - getMin()) / getBinWidth()) + 1));
diff --git a/Core/include/Acts/Utilities/AxisDefinitions.hpp b/Core/include/Acts/Utilities/AxisDefinitions.hpp
index 03938f04ae9..13ae1272185 100644
--- a/Core/include/Acts/Utilities/AxisDefinitions.hpp
+++ b/Core/include/Acts/Utilities/AxisDefinitions.hpp
@@ -73,17 +73,15 @@ std::string axesDirectionName(const std::vector<AxisDirection>& manyDir);
 std::ostream& operator<<(std::ostream& os,
                          const std::vector<AxisDirection>& manyDir);
 
-/// Enum which determines how the axis handle its outer boundaries
-/// possible values values
+/// Enum which determines how the axis handle its outer boundaries possible
+/// values values
 enum class AxisBoundaryType {
-  /// Default behaviour: out of bounds
-  /// positions are filled into the over or underflow bins
+  /// Default behaviour: out of bounds positions are filled into the over or
+  /// underflow bins
   Open,
-  /// Out-of-bounds positions resolve to first/last bin
-  /// respectively
+  /// Out-of-bounds positions resolve to first/last bin respectively
   Bound,
-  /// Out-of-bounds positions resolve to the outermost
-  /// bin on the opposite side
+  /// Out-of-bounds positions resolve to the outermost bin on the opposite side
   Closed,
 };
 
diff --git a/Core/include/Acts/Utilities/ProtoAxis.hpp b/Core/include/Acts/Utilities/ProtoAxis.hpp
index 6a72ea9747a..0101fc68cfb 100644
--- a/Core/include/Acts/Utilities/ProtoAxis.hpp
+++ b/Core/include/Acts/Utilities/ProtoAxis.hpp
@@ -8,7 +8,6 @@
 
 #pragma once
 
-#include "Acts/Utilities/Axis.hpp"
 #include "Acts/Utilities/AxisDefinitions.hpp"
 #include "Acts/Utilities/Grid.hpp"
 #include "Acts/Utilities/IAxis.hpp"
@@ -17,6 +16,7 @@
 #include <vector>
 
 namespace Acts {
+
 /// This class is a pure run-time placeholder for the axis definition
 ///
 /// The IAxis allows via the visitor pattern to access the actual axis type
@@ -31,7 +31,7 @@ class ProtoAxis {
   ///
   /// @param abType the axis boundary type
   /// @param edges the bin edges (variable binning)
-  ProtoAxis(Acts::AxisBoundaryType abType, const std::vector<double>& edges);
+  ProtoAxis(AxisBoundaryType abType, const std::vector<double>& edges);
 
   /// Convenience constructors - for equidistant binning
   ///
@@ -65,8 +65,6 @@ class ProtoAxis {
   /// @return Reference to this ProtoAxis for assignment chaining
   ProtoAxis& operator=(ProtoAxis&&) = default;
 
-  ~ProtoAxis() = default;
-
   /// @brief Return the IAxis representation
   ///
   /// @return @c AxisType of this axis
@@ -96,6 +94,46 @@ class ProtoAxis {
   /// @return the string representation
   std::string toString() const;
 
+  /// Get the minimum edge of the axis
+  /// @return the minimum edge of the axis
+  double min() const;
+
+  /// Get the maximum edge of the axis
+  /// @return the maximum edge of the axis
+  double max() const;
+
+  /// Get the number of bins
+  /// @return the number of bins
+  std::size_t nBins() const;
+
+  /// Compute the center of a given bin
+  /// @param bin the bin number for which to compute the center
+  /// @return the center of the given bin
+  double binCenter(std::size_t bin) const;
+
+  /// Compute the width of a given bin
+  /// @param bin the bin number for which to compute the width
+  /// @return the width of the given bin
+  double binWidth(std::size_t bin) const;
+
+  /// Compute the bin number for a given value
+  /// @param value the value for which to compute the bin number
+  /// @return the bin number for the given value
+  std::size_t bin(double value) const;
+
+  /// Get the bin edges
+  /// @return the bin edges
+  std::vector<double> binEdges() const;
+
+  /// Check if two axes are equal
+  /// @param lhs first axis
+  /// @param rhs second axis
+  /// @return true if the axes are equal
+  friend bool operator==(const ProtoAxis& lhs, const ProtoAxis& rhs) {
+    return lhs.getAxis() == rhs.getAxis() &&
+           lhs.isAutorange() == rhs.isAutorange();
+  }
+
  private:
   /// @brief Hidden friend ostream operator
   /// @param os the output stream
@@ -167,7 +205,7 @@ std::unique_ptr<IGrid> makeGrid(const ProtoAxis& a, const ProtoAxis& b) {
 }
 
 /// A Directed proto axis
-struct DirectedProtoAxis : public ProtoAxis {
+struct DirectedProtoAxis final : public ProtoAxis {
  public:
   /// Convenience constructors - for variable binning
   ///
@@ -197,8 +235,6 @@ struct DirectedProtoAxis : public ProtoAxis {
   DirectedProtoAxis(AxisDirection axisDir, AxisBoundaryType abType,
                     std::size_t nbins);
 
-  virtual ~DirectedProtoAxis() = default;
-
   /// Access to the AxisDirection
   /// @return the axis direction
   AxisDirection getAxisDirection() const;
@@ -207,6 +243,17 @@ struct DirectedProtoAxis : public ProtoAxis {
   /// @return the string representation
   std::string toString() const;
 
+  /// Check if two axes are equal
+  /// @param lhs first axis
+  /// @param rhs second axis
+  /// @return true if the axes are equal
+  friend bool operator==(const DirectedProtoAxis& lhs,
+                         const DirectedProtoAxis& rhs) {
+    return operator==(static_cast<const ProtoAxis&>(lhs),
+                      static_cast<const ProtoAxis&>(rhs)) &&
+           lhs.getAxisDirection() == rhs.getAxisDirection();
+  }
+
  private:
   /// @brief Hidden friend ostream operator
   /// @param os the output stream
diff --git a/Core/src/Utilities/ProtoAxis.cpp b/Core/src/Utilities/ProtoAxis.cpp
index 38be4d5acd7..bf836ff02c7 100644
--- a/Core/src/Utilities/ProtoAxis.cpp
+++ b/Core/src/Utilities/ProtoAxis.cpp
@@ -96,6 +96,35 @@ void ProtoAxis::toStream(std::ostream& os) const {
   os << toString();
 }
 
+double ProtoAxis::min() const {
+  return getAxis().getMin();
+}
+
+double ProtoAxis::max() const {
+  return getAxis().getMax();
+}
+
+std::size_t ProtoAxis::nBins() const {
+  return getAxis().getNBins();
+}
+
+double ProtoAxis::binCenter(const std::size_t bin) const {
+  return 0.5 * (getAxis().getBinEdges().at(bin) +
+                getAxis().getBinEdges().at(bin + 1));
+}
+
+double ProtoAxis::binWidth(const std::size_t bin) const {
+  return getAxis().getBinEdges().at(bin + 1) - getAxis().getBinEdges().at(bin);
+}
+
+std::size_t ProtoAxis::bin(const double value) const {
+  return getAxis().getBin(value);
+}
+
+std::vector<double> ProtoAxis::binEdges() const {
+  return getAxis().getBinEdges();
+}
+
 std::string ProtoAxis::toString() const {
   std::stringstream ss;
   const auto& axis = getAxis();
diff --git a/Examples/Algorithms/Digitization/CMakeLists.txt b/Examples/Algorithms/Digitization/CMakeLists.txt
index 79772cfdce6..2f0ddf3a3fd 100644
--- a/Examples/Algorithms/Digitization/CMakeLists.txt
+++ b/Examples/Algorithms/Digitization/CMakeLists.txt
@@ -1,7 +1,7 @@
 acts_add_library(
     ExamplesDigitization
     src/DigitizationAlgorithm.cpp
-    src/DigitizationConfig.cpp
+    src/GeometricConfig.cpp
     src/DigitizationCoordinatesConverter.cpp
     src/MeasurementCreation.cpp
     src/DigitizationConfigurator.cpp
diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/GeometricConfig.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/GeometricConfig.hpp
index c084e8ff383..0a38fd478d3 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/GeometricConfig.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/GeometricConfig.hpp
@@ -10,7 +10,7 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "Acts/Utilities/Result.hpp"
 #include "ActsExamples/Digitization/Smearers.hpp"
 #include "ActsExamples/Framework/RandomNumbers.hpp"
@@ -28,18 +28,16 @@ namespace ActsExamples {
 
 /// Configuration struct for geometric digitization
 ///
-/// If this is defined, then the geometric digitization
-/// will create clusters with cells.
-/// The BinUtility defines the segmentation and which parameters
-/// are defined by this.
-///
+/// If this is defined, then the geometric digitization will create clusters
+/// with cells. The DirectedProtoAxis defines the segmentation and which
+/// parameters are defined by this.
 struct GeometricConfig {
   // The dimensions of the measurement
-  std::vector<Acts::BoundIndices> indices = {};
+  std::vector<Acts::BoundIndices> indices;
 
   /// The (multidimensional) binning definition for the segmentation of the
   /// sensor
-  Acts::BinUtility segmentation;
+  std::vector<Acts::DirectedProtoAxis> segmentation;
 
   /// The thickness of the sensor
   double thickness = 0.;
diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/ModuleClusters.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/ModuleClusters.hpp
index 443a223d9c0..30a6e32f15a 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/ModuleClusters.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/ModuleClusters.hpp
@@ -10,7 +10,6 @@
 
 #include "Acts/Clusterization/Clusterization.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
@@ -35,7 +34,7 @@ struct ModuleValue {
 
 class ModuleClusters {
  public:
-  ModuleClusters(Acts::BinUtility segmentation,
+  ModuleClusters(std::vector<Acts::DirectedProtoAxis> segmentation,
                  std::vector<Acts::BoundIndices> geoIndices, bool merge,
                  double nsigma, bool commonCorner)
       : m_segmentation(std::move(segmentation)),
@@ -49,7 +48,7 @@ class ModuleClusters {
   digitizedParameters();
 
  private:
-  Acts::BinUtility m_segmentation;
+  std::vector<Acts::DirectedProtoAxis> m_segmentation;
   std::vector<Acts::BoundIndices> m_geoIndices;
   std::vector<ModuleValue> m_moduleValues;
   bool m_merge;
diff --git a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
index fea230fcb36..ed8a66d91a1 100644
--- a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
+++ b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
@@ -11,7 +11,6 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
 #include "ActsExamples/Digitization/ModuleClusters.hpp"
 #include "ActsExamples/EventData/GeometryContainers.hpp"
 #include "ActsExamples/EventData/Index.hpp"
@@ -359,8 +358,6 @@ DigitizedParameters DigitizationAlgorithm::localParameters(
     RandomEngine& rng) const {
   DigitizedParameters dParameters;
 
-  const auto& binningData = geoCfg.segmentation.binningData();
-
   // For digital readout, the weight needs to be split in x and y
   std::array<double, 2u> pos = {0., 0.};
   std::array<double, 2u> totalWeight = {0., 0.};
@@ -383,12 +380,12 @@ DigitizedParameters DigitizationAlgorithm::localParameters(
           // only fill component of this row/column if not yet filled
           if (!componentChannels[ib].contains(bin[ib])) {
             totalWeight[ib] += weight;
-            pos[ib] += weight * binningData[ib].center(bin[ib]);
+            pos[ib] += weight * geoCfg.segmentation.at(ib).binCenter(bin[ib]);
             componentChannels[ib].insert(bin[ib]);
           }
         } else {
           totalWeight[ib] += weight;
-          pos[ib] += weight * binningData[ib].center(bin[ib]);
+          pos[ib] += weight * geoCfg.segmentation.at(ib).binCenter(bin[ib]);
         }
         // min max channels
         bmin[ib] = std::min(bmin[ib], static_cast<std::size_t>(bin[ib]));
diff --git a/Examples/Algorithms/Digitization/src/DigitizationConfigurator.cpp b/Examples/Algorithms/Digitization/src/DigitizationConfigurator.cpp
index 54e030f9796..3a352e284ac 100644
--- a/Examples/Algorithms/Digitization/src/DigitizationConfigurator.cpp
+++ b/Examples/Algorithms/Digitization/src/DigitizationConfigurator.cpp
@@ -15,8 +15,6 @@
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Surfaces/SurfaceBounds.hpp"
 #include "Acts/Surfaces/TrapezoidBounds.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Zip.hpp"
 #include "ActsExamples/Digitization/SmearingConfig.hpp"
 #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"
@@ -53,234 +51,245 @@ bool digiConfigMaybeEqual(DigiComponentsConfig &a, DigiComponentsConfig &b) {
 }  // namespace
 
 void DigitizationConfigurator::operator()(const Acts::Surface *surface) {
-  if (surface->isSensitive()) {
-    Acts::GeometryIdentifier geoId = surface->geometryId();
-    auto dInputConfig = inputDigiComponents.find((geoId));
-    if (dInputConfig != inputDigiComponents.end()) {
-      // The output config, copy over the smearing part
-      DigiComponentsConfig dOutputConfig;
-      dOutputConfig.smearingDigiConfig = dInputConfig->smearingDigiConfig;
-
-      if (!dInputConfig->geometricDigiConfig.indices.empty()) {
-        // Copy over what can be done
-        dOutputConfig.geometricDigiConfig.indices =
-            dInputConfig->geometricDigiConfig.indices;
-        dOutputConfig.geometricDigiConfig.thickness =
-            dInputConfig->geometricDigiConfig.thickness;
-        dOutputConfig.geometricDigiConfig.chargeSmearer =
-            dInputConfig->geometricDigiConfig.chargeSmearer;
-        dOutputConfig.geometricDigiConfig.digital =
-            dInputConfig->geometricDigiConfig.digital;
-
-        const Acts::SurfaceBounds &sBounds = surface->bounds();
-        auto boundValues = sBounds.values();
-
-        const auto &inputSegmentation =
-            dInputConfig->geometricDigiConfig.segmentation;
-        Acts::BinUtility outputSegmentation;
-
-        switch (sBounds.type()) {
-          // The module is a rectangle module
-          case Acts::SurfaceBounds::eRectangle: {
-            if (inputSegmentation.binningData()[0].binvalue ==
-                Acts::AxisDirection::AxisX) {
-              double minX = boundValues[Acts::RectangleBounds::eMinX];
-              double maxX = boundValues[Acts::RectangleBounds::eMaxX];
-              unsigned int nBins = static_cast<unsigned int>(std::round(
-                  (maxX - minX) / inputSegmentation.binningData()[0].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minX), static_cast<float>(maxX),
-                  Acts::open, Acts::AxisDirection::AxisX);
-            }
-            if (inputSegmentation.binningData()[0].binvalue ==
-                    Acts::AxisDirection::AxisY ||
-                inputSegmentation.dimensions() == 2) {
-              unsigned int accessBin =
-                  inputSegmentation.dimensions() == 2 ? 1 : 0;
-              double minY = boundValues[Acts::RectangleBounds::eMinY];
-              double maxY = boundValues[Acts::RectangleBounds::eMaxY];
-              unsigned int nBins = static_cast<unsigned int>(
-                  std::round((maxY - minY) /
-                             inputSegmentation.binningData()[accessBin].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minY), static_cast<float>(maxY),
-                  Acts::open, Acts::AxisDirection::AxisY);
-            }
-          } break;
-
-          // The module is a trapezoid module
-          case Acts::SurfaceBounds::eTrapezoid: {
-            if (inputSegmentation.binningData()[0].binvalue ==
-                Acts::AxisDirection::AxisX) {
-              double maxX = std::max(
-                  boundValues[Acts::TrapezoidBounds::eHalfLengthXnegY],
-                  boundValues[Acts::TrapezoidBounds::eHalfLengthXposY]);
-              unsigned int nBins = static_cast<unsigned int>(std::round(
-                  2 * maxX / inputSegmentation.binningData()[0].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, -static_cast<float>(maxX), static_cast<float>(maxX),
-                  Acts::open, Acts::AxisDirection::AxisX);
-            }
-            if (inputSegmentation.binningData()[0].binvalue ==
-                    Acts::AxisDirection::AxisY ||
-                inputSegmentation.dimensions() == 2) {
-              unsigned int accessBin =
-                  inputSegmentation.dimensions() == 2 ? 1 : 0;
-              double maxY = boundValues[Acts::TrapezoidBounds::eHalfLengthY];
-              unsigned int nBins = static_cast<unsigned int>(
-                  std::round((2 * maxY) /
-                             inputSegmentation.binningData()[accessBin].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, -static_cast<float>(maxY), static_cast<float>(maxY),
-                  Acts::open, Acts::AxisDirection::AxisY);
-            }
-          } break;
-
-          // The module is an annulus module
-          case Acts::SurfaceBounds::eAnnulus: {
-            if (inputSegmentation.binningData()[0].binvalue ==
-                Acts::AxisDirection::AxisR) {
-              double minR = boundValues[Acts::AnnulusBounds::eMinR];
-              double maxR = boundValues[Acts::AnnulusBounds::eMaxR];
-              unsigned int nBins = static_cast<unsigned int>(std::round(
-                  (maxR - minR) / inputSegmentation.binningData()[0].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minR), static_cast<float>(maxR),
-                  Acts::open, Acts::AxisDirection::AxisR);
-            }
-            if (inputSegmentation.binningData()[0].binvalue ==
-                    Acts::AxisDirection::AxisPhi ||
-                inputSegmentation.dimensions() == 2) {
-              unsigned int accessBin =
-                  inputSegmentation.dimensions() == 2 ? 1 : 0;
-              double averagePhi = boundValues[Acts::AnnulusBounds::eAveragePhi];
-              double minPhi =
-                  averagePhi - boundValues[Acts::AnnulusBounds::eMinPhiRel];
-              double maxPhi =
-                  averagePhi + boundValues[Acts::AnnulusBounds::eMaxPhiRel];
-              unsigned int nBins = static_cast<unsigned int>(
-                  std::round((maxPhi - minPhi) /
-                             inputSegmentation.binningData()[accessBin].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minPhi), static_cast<float>(maxPhi),
-                  Acts::open, Acts::AxisDirection::AxisPhi);
-            }
-
-          } break;
-
-          // The module is a Disc Trapezoid
-          case Acts::SurfaceBounds::eDiscTrapezoid: {
-            double minR = boundValues[Acts::DiscTrapezoidBounds::eMinR];
-            double maxR = boundValues[Acts::DiscTrapezoidBounds::eMaxR];
-
-            if (inputSegmentation.binningData()[0].binvalue ==
-                Acts::AxisDirection::AxisR) {
-              unsigned int nBins = static_cast<unsigned int>(std::round(
-                  (maxR - minR) / inputSegmentation.binningData()[0].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minR), static_cast<float>(maxR),
-                  Acts::open, Acts::AxisDirection::AxisR);
-            }
-            if (inputSegmentation.binningData()[0].binvalue ==
-                    Acts::AxisDirection::AxisPhi ||
-                inputSegmentation.dimensions() == 2) {
-              unsigned int accessBin =
-                  inputSegmentation.dimensions() == 2 ? 1 : 0;
-              double hxMinR =
-                  boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXminR];
-              double hxMaxR =
-                  boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXmaxR];
-
-              double averagePhi =
-                  boundValues[Acts::DiscTrapezoidBounds::eAveragePhi];
-              double alphaMinR = std::atan2(minR, hxMinR);
-              double alphaMaxR = std::atan2(maxR, hxMaxR);
-              double alpha = std::max(alphaMinR, alphaMaxR);
-              unsigned int nBins = static_cast<unsigned int>(std::round(
-                  2 * alpha / inputSegmentation.binningData()[accessBin].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(averagePhi - alpha),
-                  static_cast<float>(averagePhi + alpha), Acts::open,
-                  Acts::AxisDirection::AxisPhi);
-            }
-
-          } break;
-
-          case Acts::SurfaceBounds::eDisc: {
-            if (inputSegmentation.binningData()[0].binvalue ==
-                Acts::AxisDirection::AxisR) {
-              double minR = boundValues[Acts::RadialBounds::eMinR];
-              double maxR = boundValues[Acts::RadialBounds::eMaxR];
-              unsigned int nBins = static_cast<unsigned int>(std::round(
-                  (maxR - minR) / inputSegmentation.binningData()[0].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minR), static_cast<float>(maxR),
-                  Acts::open, Acts::AxisDirection::AxisR);
-            }
-            if (inputSegmentation.binningData()[0].binvalue ==
-                    Acts::AxisDirection::AxisPhi ||
-                inputSegmentation.dimensions() == 2) {
-              unsigned int accessBin =
-                  inputSegmentation.dimensions() == 2 ? 1 : 0;
-
-              double averagePhi = boundValues[Acts::RadialBounds::eAveragePhi];
-              double halfPhiSector =
-                  boundValues[Acts::RadialBounds::eHalfPhiSector];
-              double minPhi = averagePhi - halfPhiSector;
-              double maxPhi = averagePhi + halfPhiSector;
-
-              unsigned int nBins = static_cast<unsigned int>(
-                  std::round((maxPhi - minPhi) /
-                             inputSegmentation.binningData()[accessBin].step));
-              outputSegmentation += Acts::BinUtility(
-                  nBins, static_cast<float>(minPhi), static_cast<float>(maxPhi),
-                  Acts::open, Acts::AxisDirection::AxisPhi);
-            }
-
-          } break;
-
-          default:
-            break;
+  if (!surface->isSensitive()) {
+    return;
+  }
+
+  Acts::GeometryIdentifier geoId = surface->geometryId();
+  const auto dInputConfig = inputDigiComponents.find(geoId);
+  if (dInputConfig == inputDigiComponents.end()) {
+    return;
+  }
+
+  // The output config, copy over the smearing part
+  DigiComponentsConfig dOutputConfig;
+  dOutputConfig.smearingDigiConfig = dInputConfig->smearingDigiConfig;
+
+  if (!dInputConfig->geometricDigiConfig.indices.empty()) {
+    // Copy over what can be done
+    dOutputConfig.geometricDigiConfig.indices =
+        dInputConfig->geometricDigiConfig.indices;
+    dOutputConfig.geometricDigiConfig.thickness =
+        dInputConfig->geometricDigiConfig.thickness;
+    dOutputConfig.geometricDigiConfig.chargeSmearer =
+        dInputConfig->geometricDigiConfig.chargeSmearer;
+    dOutputConfig.geometricDigiConfig.digital =
+        dInputConfig->geometricDigiConfig.digital;
+
+    const Acts::SurfaceBounds &sBounds = surface->bounds();
+    const std::vector<double> boundValues = sBounds.values();
+
+    const auto &inputSegmentation =
+        dInputConfig->geometricDigiConfig.segmentation;
+    std::vector<Acts::DirectedProtoAxis> outputSegmentation;
+
+    switch (sBounds.type()) {
+      // The module is a rectangle module
+      case Acts::SurfaceBounds::eRectangle: {
+        if (inputSegmentation.at(0).getAxisDirection() ==
+            Acts::AxisDirection::AxisX) {
+          const double minX = boundValues[Acts::RectangleBounds::eMinX];
+          const double maxX = boundValues[Acts::RectangleBounds::eMaxX];
+
+          const unsigned int nBins =
+              inputSegmentation.at(0).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisX,
+                                          Acts::AxisBoundaryType::Open, minX,
+                                          maxX, nBins);
+        }
+        if (inputSegmentation.at(0).getAxisDirection() ==
+                Acts::AxisDirection::AxisY ||
+            inputSegmentation.size() == 2) {
+          const unsigned int accessBin = inputSegmentation.size() == 2 ? 1 : 0;
+
+          const double minY = boundValues[Acts::RectangleBounds::eMinY];
+          const double maxY = boundValues[Acts::RectangleBounds::eMaxY];
+
+          const unsigned int nBins =
+              inputSegmentation.at(accessBin).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisY,
+                                          Acts::AxisBoundaryType::Open, minY,
+                                          maxY, nBins);
         }
-        // Set the adapted segmentation class
-        dOutputConfig.geometricDigiConfig.segmentation = outputSegmentation;
-      }
-
-      // Compactify the output map where possible
-      if (compactify) {
-        // Check for a representing volume configuration, insert if not
-        // present
-        Acts::GeometryIdentifier volGeoId =
-            Acts::GeometryIdentifier().withVolume(geoId.volume());
-
-        auto volRep = volumeLayerComponents.find(volGeoId);
-        if (volRep != volumeLayerComponents.end() &&
-            digiConfigMaybeEqual(dOutputConfig, volRep->second)) {
-          // return if the volume representation already covers this one
-          return;
-        } else {
-          volumeLayerComponents[volGeoId] = dOutputConfig;
-          outputDigiComponents.push_back({volGeoId, dOutputConfig});
+      } break;
+
+      // The module is a trapezoid module
+      case Acts::SurfaceBounds::eTrapezoid: {
+        if (inputSegmentation.at(0).getAxisDirection() ==
+            Acts::AxisDirection::AxisX) {
+          const double maxX =
+              std::max(boundValues[Acts::TrapezoidBounds::eHalfLengthXnegY],
+                       boundValues[Acts::TrapezoidBounds::eHalfLengthXposY]);
+
+          const unsigned int nBins =
+              inputSegmentation.at(0).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisX,
+                                          Acts::AxisBoundaryType::Open, -maxX,
+                                          maxX, nBins);
         }
+        if (inputSegmentation.at(0).getAxisDirection() ==
+                Acts::AxisDirection::AxisY ||
+            inputSegmentation.size() == 2) {
+          const unsigned int accessBin = inputSegmentation.size() == 2 ? 1 : 0;
+
+          const double maxY = boundValues[Acts::TrapezoidBounds::eHalfLengthY];
+
+          const unsigned int nBins =
+              inputSegmentation.at(accessBin).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisY,
+                                          Acts::AxisBoundaryType::Open, -maxY,
+                                          maxY, nBins);
+        }
+      } break;
+
+      // The module is an annulus module
+      case Acts::SurfaceBounds::eAnnulus: {
+        if (inputSegmentation.at(0).getAxisDirection() ==
+            Acts::AxisDirection::AxisR) {
+          const double minR = boundValues[Acts::AnnulusBounds::eMinR];
+          const double maxR = boundValues[Acts::AnnulusBounds::eMaxR];
+
+          const unsigned int nBins =
+              inputSegmentation.at(0).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisR,
+                                          Acts::AxisBoundaryType::Open, minR,
+                                          maxR, nBins);
+        }
+        if (inputSegmentation.at(0).getAxisDirection() ==
+                Acts::AxisDirection::AxisPhi ||
+            inputSegmentation.size() == 2) {
+          const double averagePhi =
+              boundValues[Acts::AnnulusBounds::eAveragePhi];
+          const double minPhi =
+              averagePhi - boundValues[Acts::AnnulusBounds::eMinPhiRel];
+          const double maxPhi =
+              averagePhi + boundValues[Acts::AnnulusBounds::eMaxPhiRel];
+
+          const unsigned int nBins =
+              inputSegmentation.at(0).getAxis().getNBins();
 
-        // Check for a representing layer configuration, insert if not present
-        Acts::GeometryIdentifier volLayGeoId =
-            Acts::GeometryIdentifier(volGeoId).withLayer(geoId.layer());
-        auto volLayRep = volumeLayerComponents.find(volLayGeoId);
-
-        if (volLayRep != volumeLayerComponents.end() &&
-            digiConfigMaybeEqual(dOutputConfig, volLayRep->second)) {
-          return;
-        } else {
-          volumeLayerComponents[volLayGeoId] = dOutputConfig;
-          outputDigiComponents.push_back({volLayGeoId, dOutputConfig});
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisPhi,
+                                          Acts::AxisBoundaryType::Open, minPhi,
+                                          maxPhi, nBins);
         }
-      }
+      } break;
 
-      // Insert into the output list
-      outputDigiComponents.push_back({geoId, dOutputConfig});
+      // The module is a Disc Trapezoid
+      case Acts::SurfaceBounds::eDiscTrapezoid: {
+        const double minR = boundValues[Acts::DiscTrapezoidBounds::eMinR];
+        const double maxR = boundValues[Acts::DiscTrapezoidBounds::eMaxR];
+
+        if (inputSegmentation.at(0).getAxisDirection() ==
+            Acts::AxisDirection::AxisR) {
+          const unsigned int nBins =
+              inputSegmentation.at(0).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisR,
+                                          Acts::AxisBoundaryType::Open, minR,
+                                          maxR, nBins);
+        }
+        if (inputSegmentation.at(0).getAxisDirection() ==
+                Acts::AxisDirection::AxisPhi ||
+            inputSegmentation.size() == 2) {
+          const unsigned int accessBin = inputSegmentation.size() == 2 ? 1 : 0;
+
+          const double hxMinR =
+              boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXminR];
+          const double hxMaxR =
+              boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXmaxR];
+
+          const double averagePhi =
+              boundValues[Acts::DiscTrapezoidBounds::eAveragePhi];
+          const double alphaMinR = std::atan2(minR, hxMinR);
+          const double alphaMaxR = std::atan2(maxR, hxMaxR);
+          const double alpha = std::max(alphaMinR, alphaMaxR);
+
+          const unsigned int nBins =
+              inputSegmentation.at(accessBin).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(
+              Acts::AxisDirection::AxisPhi, Acts::AxisBoundaryType::Open,
+              averagePhi - alpha, averagePhi + alpha, nBins);
+        }
+      } break;
+
+      case Acts::SurfaceBounds::eDisc: {
+        if (inputSegmentation.at(0).getAxisDirection() ==
+            Acts::AxisDirection::AxisR) {
+          const double minR = boundValues[Acts::RadialBounds::eMinR];
+          const double maxR = boundValues[Acts::RadialBounds::eMaxR];
+
+          const unsigned int nBins =
+              inputSegmentation.at(0).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisR,
+                                          Acts::AxisBoundaryType::Open, minR,
+                                          maxR, nBins);
+        }
+        if (inputSegmentation.at(0).getAxisDirection() ==
+                Acts::AxisDirection::AxisPhi ||
+            inputSegmentation.size() == 2) {
+          const unsigned int accessBin = inputSegmentation.size() == 2 ? 1 : 0;
+
+          const double averagePhi =
+              boundValues[Acts::RadialBounds::eAveragePhi];
+          const double halfPhiSector =
+              boundValues[Acts::RadialBounds::eHalfPhiSector];
+          const double minPhi = averagePhi - halfPhiSector;
+          const double maxPhi = averagePhi + halfPhiSector;
+
+          const unsigned int nBins =
+              inputSegmentation.at(accessBin).getAxis().getNBins();
+
+          outputSegmentation.emplace_back(Acts::AxisDirection::AxisPhi,
+                                          Acts::AxisBoundaryType::Open, minPhi,
+                                          maxPhi, nBins);
+        }
+      } break;
+
+      default:
+        break;
+    }
+
+    // Set the adapted segmentation class
+    dOutputConfig.geometricDigiConfig.segmentation = outputSegmentation;
+  }
+
+  // Compactify the output map where possible
+  if (compactify) {
+    // Check for a representing volume configuration, insert if not
+    // present
+    Acts::GeometryIdentifier volGeoId =
+        Acts::GeometryIdentifier().withVolume(geoId.volume());
+
+    auto volRep = volumeLayerComponents.find(volGeoId);
+    if (volRep != volumeLayerComponents.end() &&
+        digiConfigMaybeEqual(dOutputConfig, volRep->second)) {
+      // return if the volume representation already covers this one
+      return;
     }
+    volumeLayerComponents[volGeoId] = dOutputConfig;
+    outputDigiComponents.push_back({volGeoId, dOutputConfig});
+
+    // Check for a representing layer configuration, insert if not present
+    Acts::GeometryIdentifier volLayGeoId =
+        Acts::GeometryIdentifier(volGeoId).withLayer(geoId.layer());
+    const auto volLayRep = volumeLayerComponents.find(volLayGeoId);
+    if (volLayRep != volumeLayerComponents.end() &&
+        digiConfigMaybeEqual(dOutputConfig, volLayRep->second)) {
+      return;
+    }
+    volumeLayerComponents[volLayGeoId] = dOutputConfig;
+    outputDigiComponents.push_back({volLayGeoId, dOutputConfig});
   }
+
+  // Insert into the output list
+  outputDigiComponents.push_back({geoId, dOutputConfig});
 }
 
 }  // namespace ActsExamples
diff --git a/Examples/Algorithms/Digitization/src/DigitizationConfig.cpp b/Examples/Algorithms/Digitization/src/GeometricConfig.cpp
similarity index 74%
rename from Examples/Algorithms/Digitization/src/DigitizationConfig.cpp
rename to Examples/Algorithms/Digitization/src/GeometricConfig.cpp
index 232ae598cc7..fc1e48cebb2 100644
--- a/Examples/Algorithms/Digitization/src/DigitizationConfig.cpp
+++ b/Examples/Algorithms/Digitization/src/GeometricConfig.cpp
@@ -6,7 +6,7 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-#include "ActsExamples/Digitization/DigitizationConfig.hpp"
+#include "ActsExamples/Digitization/GeometricConfig.hpp"
 
 namespace ActsExamples {
 
@@ -14,17 +14,18 @@ std::vector<double> GeometricConfig::variances(
     const std::array<std::size_t, 2u>& csizes,
     const std::array<std::size_t, 2u>& cmins) const {
   std::vector<double> rVariances;
-  for (const auto& bIndex : indices) {
+  for (std::size_t i = 0; i < indices.size(); ++i) {
+    const auto& bIndex = indices[i];
     double var = 0.;
     if (varianceMap.contains(bIndex)) {
       // Try to find the variance for this cluster size
-      std::size_t lsize =
+      const std::size_t lsize =
           std::min(csizes[bIndex], varianceMap.at(bIndex).size());
       var = varianceMap.at(bIndex).at(lsize - 1);
     } else {
       // Pitch size ofer / sqrt(12) as error instead
-      std::size_t ictr = cmins[bIndex] + csizes[bIndex] / 2;
-      var = std::pow(segmentation.binningData()[bIndex].width(ictr), 2) / 12.0;
+      const std::size_t ictr = cmins[i] + csizes[i] / 2;
+      var = std::pow(segmentation.at(bIndex).binWidth(ictr), 2) / 12.0;
     }
     rVariances.push_back(var);
   }
diff --git a/Examples/Algorithms/Digitization/src/ModuleClusters.cpp b/Examples/Algorithms/Digitization/src/ModuleClusters.cpp
index 879913db9b5..d82919754eb 100644
--- a/Examples/Algorithms/Digitization/src/ModuleClusters.cpp
+++ b/Examples/Algorithms/Digitization/src/ModuleClusters.cpp
@@ -275,7 +275,6 @@ ModuleValue ModuleClusters::squash(std::vector<ModuleValue>& values) {
   // Now do the geometric indices
   Cluster clus;
 
-  const auto& binningData = m_segmentation.binningData();
   Acts::Vector2 pos(0., 0.);
   Acts::Vector2 var(0., 0.);
 
@@ -301,10 +300,10 @@ ModuleValue ModuleClusters::squash(std::vector<ModuleValue>& values) {
     b1min = std::min(b1min, b1);
     b1max = std::max(b1max, b1);
 
-    float p0 = binningData[0].center(b0);
-    float w0 = binningData[0].width(b0);
-    float p1 = binningData[1].center(b1);
-    float w1 = binningData[1].width(b1);
+    const double p0 = m_segmentation[0].binCenter(b0);
+    const double w0 = m_segmentation[0].binWidth(b0);
+    const double p1 = m_segmentation[1].binCenter(b1);
+    const double w1 = m_segmentation[1].binWidth(b1);
 
     pos += Acts::Vector2(weights.at(i) * p0, weights.at(i) * p1);
     // Assume uniform distribution to compute error
diff --git a/Examples/Io/Json/src/JsonDigitizationConfig.cpp b/Examples/Io/Json/src/JsonDigitizationConfig.cpp
index 435dba3da9d..c1825d2eec5 100644
--- a/Examples/Io/Json/src/JsonDigitizationConfig.cpp
+++ b/Examples/Io/Json/src/JsonDigitizationConfig.cpp
@@ -9,7 +9,10 @@
 #include "ActsExamples/Io/Json/JsonDigitizationConfig.hpp"
 
 #include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Utilities/AxisDefinitions.hpp"
+#include "Acts/Utilities/BinUtility.hpp"
 #include "Acts/Utilities/BinningData.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "ActsExamples/Digitization/Smearers.hpp"
 #include "ActsExamples/Framework/RandomNumbers.hpp"
 #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"
@@ -131,7 +134,11 @@ void ActsExamples::to_json(nlohmann::json& j, const GeometricConfig& gdc) {
     indices.push_back(static_cast<std::size_t>(idx));
   }
   j["indices"] = indices;
-  j["segmentation"] = nlohmann::json(gdc.segmentation);
+  Acts::BinUtility segmentation;
+  for (const Acts::DirectedProtoAxis& axis : gdc.segmentation) {
+    segmentation += Acts::BinUtility(axis);
+  }
+  j["segmentation"] = segmentation;
   j["thickness"] = gdc.thickness;
   j["threshold"] = gdc.threshold;
   j["digital"] = gdc.digital;
@@ -144,7 +151,25 @@ void ActsExamples::from_json(const nlohmann::json& j, GeometricConfig& gdc) {
   for (const auto& jidx : j["indices"]) {
     gdc.indices.push_back(static_cast<Acts::BoundIndices>(jidx));
   }
-  from_json(j["segmentation"], gdc.segmentation);
+  Acts::BinUtility segmentation;
+  from_json(j["segmentation"], segmentation);
+  for (const Acts::BinningData& binData : segmentation.binningData()) {
+    const Acts::AxisDirection axisDir = binData.binvalue;
+    const Acts::AxisType axisType = binData.type == Acts::equidistant
+                                        ? Acts::AxisType::Equidistant
+                                        : Acts::AxisType::Variable;
+    const Acts::AxisBoundaryType abType = binData.option == Acts::open
+                                              ? Acts::AxisBoundaryType::Open
+                                              : Acts::AxisBoundaryType::Closed;
+    if (axisType == Acts::AxisType::Equidistant) {
+      gdc.segmentation.emplace_back(axisDir, abType, binData.min, binData.max,
+                                    binData.bins());
+    } else {
+      const std::vector<double> edges(binData.boundaries().begin(),
+                                      binData.boundaries().end());
+      gdc.segmentation.emplace_back(axisDir, abType, edges);
+    }
+  }
   gdc.thickness = j["thickness"];
   gdc.threshold = j["threshold"];
   gdc.digital = j["digital"];
diff --git a/Fatras/include/ActsFatras/Digitization/Channelizer.hpp b/Fatras/include/ActsFatras/Digitization/Channelizer.hpp
index 51d21386895..64b17d5eb40 100644
--- a/Fatras/include/ActsFatras/Digitization/Channelizer.hpp
+++ b/Fatras/include/ActsFatras/Digitization/Channelizer.hpp
@@ -10,14 +10,12 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "ActsFatras/Digitization/Segmentizer.hpp"
 #include "ActsFatras/Digitization/SurfaceDrift.hpp"
 #include "ActsFatras/Digitization/SurfaceMask.hpp"
 #include "ActsFatras/EventData/Hit.hpp"
 
-#include <numeric>
-
 namespace Acts {
 class Surface;
 }
@@ -51,8 +49,8 @@ class Channelizer {
   Acts::Result<std::vector<Segmentizer::ChannelSegment>> channelize(
       const Hit& hit, const Acts::Surface& surface,
       const Acts::GeometryContext& gctx, const Acts::Vector3& driftDir,
-      const Acts::BinUtility& segmentation, double thickness,
-      double minRelPerpDrift = 0.001) const;
+      const std::vector<Acts::DirectedProtoAxis>& segmentation,
+      double thickness, double minRelPerpDrift = 0.001) const;
 };
 
 }  // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp b/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp
index 79859335771..f0e3892baaa 100644
--- a/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp
+++ b/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp
@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 
 #include <array>
 #include <utility>
@@ -96,10 +97,10 @@ struct Segmentizer {
   /// @param segment The surface segment (cartesian coordinates)
   ///
   /// @return a vector of ChannelSegment objects
-  std::vector<ChannelSegment> segments(const Acts::GeometryContext& geoCtx,
-                                       const Acts::Surface& surface,
-                                       const Acts::BinUtility& segmentation,
-                                       const Segment2D& segment) const;
+  std::vector<ChannelSegment> segments(
+      const Acts::GeometryContext& geoCtx, const Acts::Surface& surface,
+      const std::vector<Acts::DirectedProtoAxis>& segmentation,
+      const Segment2D& segment) const;
 };
 
 }  // namespace ActsFatras
diff --git a/Fatras/src/Digitization/Channelizer.cpp b/Fatras/src/Digitization/Channelizer.cpp
index 1d180226b68..588fbdcf045 100644
--- a/Fatras/src/Digitization/Channelizer.cpp
+++ b/Fatras/src/Digitization/Channelizer.cpp
@@ -13,7 +13,7 @@ namespace ActsFatras {
 Acts::Result<std::vector<Segmentizer::ChannelSegment>> Channelizer::channelize(
     const Hit& hit, const Acts::Surface& surface,
     const Acts::GeometryContext& gctx, const Acts::Vector3& driftDir,
-    const Acts::BinUtility& segmentation, double thickness,
+    const std::vector<Acts::DirectedProtoAxis>& segmentation, double thickness,
     double minRelPerpDrift) const {
   // Drifted surface and scalor 2D to 3D segment
   // SurfaceDrift handles the surface-type-specific local frame internally
diff --git a/Fatras/src/Digitization/Segmentizer.cpp b/Fatras/src/Digitization/Segmentizer.cpp
index 37f9efbfe00..824f21c5d3b 100644
--- a/Fatras/src/Digitization/Segmentizer.cpp
+++ b/Fatras/src/Digitization/Segmentizer.cpp
@@ -10,27 +10,28 @@
 
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Surfaces/detail/IntersectionHelper2D.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 
 #include <algorithm>
 #include <cmath>
 #include <span>
+#include <vector>
 
 namespace ActsFatras {
 
 std::vector<Segmentizer::ChannelSegment> Segmentizer::segments(
     const Acts::GeometryContext& geoCtx, const Acts::Surface& surface,
-    const Acts::BinUtility& segmentation, const Segment2D& segment) const {
+    const std::vector<Acts::DirectedProtoAxis>& segmentation,
+    const Segment2D& segment) const {
   // Return if the segmentation is not two-dimensional
   // (strips need to have one bin along the strip)
-  if (segmentation.dimensions() != 2) {
+  if (segmentation.size() != 2) {
     return {};
   }
 
   // Start and end point
-  const auto& start = segment[0];
-  const auto& end = segment[1];
+  const Acts::Vector2& start = segment[0];
+  const Acts::Vector2& end = segment[1];
 
   // Full path length - the full channel
   auto segment2d = (end - start);
@@ -44,70 +45,70 @@ std::vector<Segmentizer::ChannelSegment> Segmentizer::segments(
     // unrolled readout frame (rPhi, z). Either way the cell boundaries are
     // axis-aligned straight lines and the stepping algorithm is identical.
     // Get the segmentation and convert it to lines & arcs
-    bstart = {static_cast<unsigned int>(segmentation.bin(start, 0)),
-              static_cast<unsigned int>(segmentation.bin(start, 1))};
-    bend = {static_cast<unsigned int>(segmentation.bin(end, 0)),
-            static_cast<unsigned int>(segmentation.bin(end, 1))};
+    bstart = {static_cast<unsigned int>(segmentation.at(0).bin(start[0]) - 1),
+              static_cast<unsigned int>(segmentation.at(1).bin(start[1]) - 1)};
+    bend = {static_cast<unsigned int>(segmentation.at(0).bin(end[0]) - 1),
+            static_cast<unsigned int>(segmentation.at(1).bin(end[1]) - 1)};
     // Fast single channel exit
     if (bstart == bend) {
       return {ChannelSegment(bstart, {start, end}, segment2d.norm())};
     }
     // The lines channel segment lines along x
     if (bstart[0] != bend[0]) {
-      double k = segment2d.y() / segment2d.x();
-      double d = start.y() - k * start.x();
+      const double k = segment2d.y() / segment2d.x();
+      const double d = start.y() - k * start.x();
 
-      const auto& xboundaries = segmentation.binningData()[0].boundaries();
-      std::span<const float> xbbounds(
+      const std::vector<double> xboundaries = segmentation.at(0).binEdges();
+      const std::span<const double> xbbounds(
           xboundaries.begin() + std::min(bstart[0], bend[0]) + 1,
           xboundaries.begin() + std::max(bstart[0], bend[0]) + 1);
-      for (const auto x : xbbounds) {
+      for (const double x : xbbounds) {
         cSteps.push_back(ChannelStep{
             {(bstart[0] < bend[0] ? 1 : -1), 0}, {x, k * x + d}, start});
       }
     }
     // The lines channel segment lines along y
     if (bstart[1] != bend[1]) {
-      double k = segment2d.x() / segment2d.y();
-      double d = start.x() - k * start.y();
-      const auto& yboundaries = segmentation.binningData()[1].boundaries();
-      std::span<const float> ybbounds(
+      const double k = segment2d.x() / segment2d.y();
+      const double d = start.x() - k * start.y();
+      const std::vector<double> yboundaries = segmentation.at(1).binEdges();
+      const std::span<const double> ybbounds(
           yboundaries.begin() + std::min(bstart[1], bend[1]) + 1,
           yboundaries.begin() + std::max(bstart[1], bend[1]) + 1);
-      for (const auto y : ybbounds) {
+      for (const double y : ybbounds) {
         cSteps.push_back(ChannelStep{
             {0, (bstart[1] < bend[1] ? 1 : -1)}, {k * y + d, y}, start});
       }
     }
 
   } else if (surface.type() == Acts::Surface::SurfaceType::Disc) {
-    Acts::Vector2 pstart(Acts::VectorHelpers::perp(start),
-                         Acts::VectorHelpers::phi(start));
-    Acts::Vector2 pend(Acts::VectorHelpers::perp(end),
-                       Acts::VectorHelpers::phi(end));
+    const Acts::Vector2 pstart(Acts::VectorHelpers::perp(start),
+                               Acts::VectorHelpers::phi(start));
+    const Acts::Vector2 pend(Acts::VectorHelpers::perp(end),
+                             Acts::VectorHelpers::phi(end));
 
     // Get the segmentation and convert it to lines & arcs
-    bstart = {static_cast<unsigned int>(segmentation.bin(pstart, 0)),
-              static_cast<unsigned int>(segmentation.bin(pstart, 1))};
-    bend = {static_cast<unsigned int>(segmentation.bin(pend, 0)),
-            static_cast<unsigned int>(segmentation.bin(pend, 1))};
+    bstart = {static_cast<unsigned int>(segmentation.at(0).bin(pstart[0]) - 1),
+              static_cast<unsigned int>(segmentation.at(1).bin(pstart[1]) - 1)};
+    bend = {static_cast<unsigned int>(segmentation.at(0).bin(pend[0]) - 1),
+            static_cast<unsigned int>(segmentation.at(1).bin(pend[1]) - 1)};
 
     // Fast single channel exit
     if (bstart == bend) {
       return {ChannelSegment(bstart, {start, end}, segment2d.norm())};
     }
 
-    double phistart = pstart[1];
-    double phiend = pend[1];
+    const double phistart = pstart[1];
+    const double phiend = pend[1];
 
     // The radial boundaries
     if (bstart[0] != bend[0]) {
-      const auto& rboundaries = segmentation.binningData()[0].boundaries();
-      std::span<const float> rbbounds(
+      const std::vector<double> rboundaries = segmentation.at(0).binEdges();
+      const std::span<const double> rbbounds(
           rboundaries.begin() + std::min(bstart[0], bend[0]) + 1,
           rboundaries.begin() + std::max(bstart[0], bend[0]) + 1);
-      for (const auto r : rbbounds) {
-        auto radIntersection =
+      for (const double r : rbbounds) {
+        const auto radIntersection =
             Acts::detail::IntersectionHelper2D::intersectCircleSegment(
                 r, std::min(phistart, phiend), std::max(phistart, phiend),
                 start, (end - start).normalized());
@@ -118,18 +119,18 @@ std::vector<Segmentizer::ChannelSegment> Segmentizer::segments(
     }
     // The phi boundaries
     if (bstart[1] != bend[1]) {
-      double referenceR =
+      const double referenceR =
           surface.referencePositionValue(geoCtx, Acts::AxisDirection::AxisR);
-      Acts::Vector2 origin = {0., 0.};
-      const auto& phiboundaries = segmentation.binningData()[1].boundaries();
-      std::span<const float> phibbounds(
+      const Acts::Vector2 origin = {0., 0.};
+      const std::vector<double> phiboundaries = segmentation.at(1).binEdges();
+      const std::span<const double> phibbounds(
           phiboundaries.begin() + std::min(bstart[1], bend[1]) + 1,
           phiboundaries.begin() + std::max(bstart[1], bend[1]) + 1);
 
-      for (const auto phi : phibbounds) {
+      for (const double phi : phibbounds) {
         Acts::Vector2 philine(referenceR * std::cos(phi),
                               referenceR * std::sin(phi));
-        auto phiIntersection =
+        const auto phiIntersection =
             Acts::detail::IntersectionHelper2D::intersectSegment(
                 origin, philine, start, (end - start).normalized());
         cSteps.push_back(ChannelStep{{0, (bstart[1] < bend[1] ? 1 : -1)},
diff --git a/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp b/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp
index cafc9f35073..c2714d00eaf 100644
--- a/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp
+++ b/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp
@@ -8,8 +8,8 @@
 
 #include <boost/test/unit_test.hpp>
 
-#include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/BinningData.hpp"
+#include "Acts/Utilities/AxisDefinitions.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "ActsExamples/Digitization/ModuleClusters.hpp"
 #include "ActsFatras/Digitization/Segmentizer.hpp"
 
@@ -19,16 +19,17 @@ using namespace ActsExamples;
 
 namespace ActsTests {
 
-DigitizedParameters makeDigitizationParameters(const Vector2 &position,
-                                               const Vector2 &variance,
-                                               const BinUtility &binUtility) {
-  auto [binX, binY, _] =
-      binUtility.binTriple((Vector3() << position, 0).finished());
-  Segmentizer::Bin2D bin = {static_cast<Segmentizer::Bin2D::value_type>(binX),
-                            static_cast<Segmentizer::Bin2D::value_type>(binY)};
-  Segmentizer::Segment2D segment = {position, position};
-  double activation = 1;
-  Cluster::Cell cell = {bin, segment, activation};
+DigitizedParameters makeDigitizationParameters(
+    const Vector2 &position, const Vector2 &variance,
+    const std::vector<DirectedProtoAxis> &segmentation) {
+  const std::size_t binX = segmentation[0].bin(position.x());
+  const std::size_t binY = segmentation[1].bin(position.y());
+  const Segmentizer::Bin2D bin = {
+      static_cast<Segmentizer::Bin2D::value_type>(binX),
+      static_cast<Segmentizer::Bin2D::value_type>(binY)};
+  const Segmentizer::Segment2D segment = {position, position};
+  const double activation = 1;
+  const Cluster::Cell cell = {bin, segment, activation};
 
   Cluster cluster;
   cluster.sizeLoc0 = 1;
@@ -46,22 +47,22 @@ DigitizedParameters makeDigitizationParameters(const Vector2 &position,
 
 auto testDigitizedParametersWithTwoClusters(bool merge, const Vector2 &firstHit,
                                             const Vector2 &secondHit) {
-  BinUtility binUtility;
-  binUtility += BinUtility(BinningData(
-      BinningOption::open, AxisDirection::AxisX, 20, -10.0f, 10.0f));
-  binUtility += BinUtility(BinningData(
-      BinningOption::open, AxisDirection::AxisY, 20, -10.0f, 10.0f));
+  std::vector<DirectedProtoAxis> segmentation;
+  segmentation.emplace_back(AxisDirection::AxisX, AxisBoundaryType::Open, -10.0,
+                            10.0, 20);
+  segmentation.emplace_back(AxisDirection::AxisY, AxisBoundaryType::Open, -10.0,
+                            10.0, 20);
   std::vector<Acts::BoundIndices> boundIndices = {eBoundLoc0, eBoundLoc1};
   double nsigma = 1;
   bool commonCorner = true;
 
-  ModuleClusters moduleClusters(binUtility, boundIndices, merge, nsigma,
+  ModuleClusters moduleClusters(segmentation, boundIndices, merge, nsigma,
                                 commonCorner);
 
-  moduleClusters.add(makeDigitizationParameters(firstHit, {1, 1}, binUtility),
+  moduleClusters.add(makeDigitizationParameters(firstHit, {1, 1}, segmentation),
                      0);
-  moduleClusters.add(makeDigitizationParameters(secondHit, {1, 1}, binUtility),
-                     1);
+  moduleClusters.add(
+      makeDigitizationParameters(secondHit, {1, 1}, segmentation), 1);
 
   return moduleClusters.digitizedParameters();
 }
diff --git a/Tests/UnitTests/Examples/Io/Json/JsonDigitizationConfigTests.cpp b/Tests/UnitTests/Examples/Io/Json/JsonDigitizationConfigTests.cpp
index 0c8a0ccbe9c..52185083849 100644
--- a/Tests/UnitTests/Examples/Io/Json/JsonDigitizationConfigTests.cpp
+++ b/Tests/UnitTests/Examples/Io/Json/JsonDigitizationConfigTests.cpp
@@ -13,8 +13,7 @@
 #include "Acts/EventData/detail/GenerateParameters.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Surfaces/StrawSurface.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "ActsExamples/Digitization/DigitizationConfig.hpp"
 #include "ActsExamples/Digitization/Smearers.hpp"
 #include "ActsExamples/Io/Json/JsonDigitizationConfig.hpp"
@@ -157,9 +156,11 @@ BOOST_AUTO_TEST_CASE(DigitizationConfigRoundTrip) {
 
   GeometricConfig gdc;
 
-  BinUtility segmentation;
-  segmentation += BinUtility(336, -8.4, 8.4, open, AxisDirection::AxisX);
-  segmentation += BinUtility(1280, -36, 36, open, AxisDirection::AxisY);
+  std::vector<DirectedProtoAxis> segmentation;
+  segmentation.emplace_back(AxisDirection::AxisX, AxisBoundaryType::Open, -8.4,
+                            8.4, 336);
+  segmentation.emplace_back(AxisDirection::AxisY, AxisBoundaryType::Open, -36,
+                            36, 1280);
 
   gdc.segmentation = segmentation;
   gdc.threshold = 0.01;
@@ -185,8 +186,8 @@ BOOST_AUTO_TEST_CASE(DigitizationConfigRoundTrip) {
   DigiComponentsConfig dcTest(dcJsonIn);
   BOOST_CHECK(dcTest.geometricDigiConfig.indices ==
               dcRef.geometricDigiConfig.indices);
-  BOOST_CHECK_EQUAL(dcTest.geometricDigiConfig.segmentation.dimensions(),
-                    dcRef.geometricDigiConfig.segmentation.dimensions());
+  BOOST_CHECK_EQUAL(dcTest.geometricDigiConfig.segmentation.size(),
+                    dcRef.geometricDigiConfig.segmentation.size());
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp b/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp
index fcc7a8cd7c2..7097a4337d1 100644
--- a/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp
+++ b/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp
@@ -10,10 +10,9 @@
 
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/Surfaces/CurvilinearSurface.hpp"
-#include "Acts/Surfaces/PlaneSurface.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "ActsFatras/Digitization/Channelizer.hpp"
-#include "ActsFatras/Digitization/SurfaceDrift.hpp"
 #include "ActsFatras/Digitization/SurfaceMask.hpp"
 
 #include <numeric>
@@ -24,7 +23,7 @@ using namespace ActsFatras;
 
 struct Helper {
   std::shared_ptr<Surface> surface;
-  BinUtility segmentation;
+  std::vector<DirectedProtoAxis> segmentation;
 
   GeometryContext gctx = GeometryContext::dangerouslyDefaultConstruct();
   double thickness = 125_um;
@@ -33,17 +32,17 @@ struct Helper {
   ActsFatras::Channelizer channelizer;
 
   Helper() {
-    surface = CurvilinearSurface(Vector3::Zero(), Vector3{0.0, 0.0, 1.0})
-                  .planeSurface();
+    surface =
+        CurvilinearSurface(Vector3::Zero(), Vector3{0.0, 0.0, 1.0}).surface();
 
     float pitchSize = 50_um;
     float min = -200_um;
     float max = 200_um;
     int bins = static_cast<int>((max - min) / pitchSize);
-    segmentation =
-        BinUtility(bins, min, max, BinningOption::open, AxisDirection::AxisX);
-    segmentation +=
-        BinUtility(bins, min, max, BinningOption::open, AxisDirection::AxisY);
+    segmentation.emplace_back(AxisDirection::AxisX, AxisBoundaryType::Open, min,
+                              max, bins);
+    segmentation.emplace_back(AxisDirection::AxisY, AxisBoundaryType::Open, min,
+                              max, bins);
   }
 
   auto channelize(const Vector3 &pos3, const Vector3 &dir3) const {
diff --git a/Tests/UnitTests/Fatras/Digitization/PlanarSurfaceTestBeds.hpp b/Tests/UnitTests/Fatras/Digitization/PlanarSurfaceTestBeds.hpp
index 8e22c46870d..5e9d3df014c 100644
--- a/Tests/UnitTests/Fatras/Digitization/PlanarSurfaceTestBeds.hpp
+++ b/Tests/UnitTests/Fatras/Digitization/PlanarSurfaceTestBeds.hpp
@@ -16,9 +16,9 @@
 #include "Acts/Surfaces/RadialBounds.hpp"
 #include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/TrapezoidBounds.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/AxisDefinitions.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 
 #include <numbers>
 #include <tuple>
@@ -32,7 +32,7 @@ using Randomizer = std::function<Acts::Vector2(double, double)>;
 
 using PlanarTestBed =
     std::tuple<std::string, std::shared_ptr<const Acts::Surface>,
-               Acts::BinUtility, Randomizer>;
+               std::vector<Acts::DirectedProtoAxis>, Randomizer>;
 
 /// Helper struct to create a testbed for Digitization steps
 struct PlanarSurfaceTestBeds {
@@ -51,10 +51,11 @@ struct PlanarSurfaceTestBeds {
     auto rectangle = std::make_shared<Acts::RectangleBounds>(xhalf, yhalf);
     auto rSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
         Acts::Transform3::Identity(), rectangle);
-    Acts::BinUtility pixelated(15, -xhalf, xhalf, Acts::open,
-                               Acts::AxisDirection::AxisX);
-    pixelated += Acts::BinUtility(26, -yhalf, yhalf, Acts::open,
-                                  Acts::AxisDirection::AxisY);
+    std::vector<Acts::DirectedProtoAxis> pixelated;
+    pixelated.emplace_back(Acts::AxisDirection::AxisX,
+                           Acts::AxisBoundaryType::Open, -xhalf, xhalf, 15);
+    pixelated.emplace_back(Acts::AxisDirection::AxisY,
+                           Acts::AxisBoundaryType::Open, -yhalf, yhalf, 26);
     RectangleRandom rRandom(xhalf * rScale, yhalf * rScale);
 
     // Cartesian strip test in Trapezoid
@@ -65,10 +66,12 @@ struct PlanarSurfaceTestBeds {
         std::make_shared<Acts::TrapezoidBounds>(xhalfminy, xhalfmaxy, yhalf);
     auto tSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
         Acts::Transform3::Identity(), trapezoid);
-    Acts::BinUtility stripsX(35, -xhalfmaxy, xhalfmaxy, Acts::open,
-                             Acts::AxisDirection::AxisX);
-    stripsX += Acts::BinUtility(1, -yhalf, yhalf, Acts::open,
-                                Acts::AxisDirection::AxisY);
+    std::vector<Acts::DirectedProtoAxis> stripsX;
+    stripsX.emplace_back(Acts::AxisDirection::AxisX,
+                         Acts::AxisBoundaryType::Open, -xhalfmaxy, xhalfmaxy,
+                         35);
+    stripsX.emplace_back(Acts::AxisDirection::AxisY,
+                         Acts::AxisBoundaryType::Open, -yhalf, yhalf, 1);
     TrapezoidRandom tRandom(xhalfminy * rScale, xhalfmaxy * rScale,
                             yhalf * rScale);
 
@@ -84,11 +87,12 @@ struct PlanarSurfaceTestBeds {
         std::make_shared<Acts::DiscTrapezoidBounds>(xmin, xmax, rmin, rmax);
     auto dtSurface = Acts::Surface::makeShared<Acts::DiscSurface>(
         Acts::Transform3::Identity(), discTrapezoid);
-    Acts::BinUtility stripsPhi(1, rmin, rmax, Acts::open,
-                               Acts::AxisDirection::AxisR);
-    stripsPhi += Acts::BinUtility(25, std::numbers::pi / 2. - alpha,
-                                  std::numbers::pi / 2. + alpha, Acts::open,
-                                  Acts::AxisDirection::AxisPhi);
+    std::vector<Acts::DirectedProtoAxis> stripsPhi;
+    stripsPhi.emplace_back(Acts::AxisDirection::AxisR,
+                           Acts::AxisBoundaryType::Open, rmin, rmax, 1);
+    stripsPhi.emplace_back(
+        Acts::AxisDirection::AxisPhi, Acts::AxisBoundaryType::Open,
+        std::numbers::pi / 2. - alpha, std::numbers::pi / 2. + alpha, 25);
     TrapezoidRandom dtRandom(xmin * rScale, xmax * rScale, rmin * irScale,
                              ymax * rScale);
 
@@ -97,12 +101,13 @@ struct PlanarSurfaceTestBeds {
         rmin, rmax, std::numbers::pi / 4., std::numbers::pi / 2.);
     auto dSurface = Acts::Surface::makeShared<Acts::DiscSurface>(
         Acts::Transform3::Identity(), discRadial);
-    Acts::BinUtility rphiseg(10, rmin, rmax, Acts::open,
-                             Acts::AxisDirection::AxisR);
-    rphiseg +=
-        Acts::BinUtility(20, (std::numbers::pi / 2. - std::numbers::pi / 4.),
-                         (std::numbers::pi / 2. + std::numbers::pi / 4.),
-                         Acts::open, Acts::AxisDirection::AxisPhi);
+    std::vector<Acts::DirectedProtoAxis> rphiseg;
+    rphiseg.emplace_back(Acts::AxisDirection::AxisR,
+                         Acts::AxisBoundaryType::Open, rmin, rmax, 10);
+    rphiseg.emplace_back(Acts::AxisDirection::AxisPhi,
+                         Acts::AxisBoundaryType::Open,
+                         (std::numbers::pi / 2. - std::numbers::pi / 4.),
+                         (std::numbers::pi / 2. + std::numbers::pi / 4.), 20);
 
     DiscRandom dRandom(
         rmin * irScale, rmax * rScale,
@@ -129,10 +134,11 @@ struct PlanarSurfaceTestBeds {
       rmax = std::max(rmax, r);
     });
 
-    Acts::BinUtility stripsPhiA(1, rmin, rmax, Acts::open,
-                                Acts::AxisDirection::AxisR);
-    stripsPhiA += Acts::BinUtility(12, phimin, phimax, Acts::open,
-                                   Acts::AxisDirection::AxisPhi);
+    std::vector<Acts::DirectedProtoAxis> stripsPhiA;
+    stripsPhiA.emplace_back(Acts::AxisDirection::AxisR,
+                            Acts::AxisBoundaryType::Open, rmin, rmax, 1);
+    stripsPhiA.emplace_back(Acts::AxisDirection::AxisPhi,
+                            Acts::AxisBoundaryType::Open, phimin, phimax, 12);
     AnnulusRandom aRandom(rmin * irScale, rmax * rScale, phimin * rScale,
                           phimax * rScale, aorigin.x(), aorigin.y());
 
diff --git a/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp b/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp
index fbab73704f6..325d817542d 100644
--- a/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp
+++ b/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp
@@ -18,8 +18,7 @@
 #include "Acts/Surfaces/RadialBounds.hpp"
 #include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
-#include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
 #include "ActsFatras/Digitization/Segmentizer.hpp"
 
 #include <cmath>
@@ -50,15 +49,18 @@ BOOST_AUTO_TEST_CASE(SegmentizerCartesian) {
                                                         rectangleBounds);
 
   // The segmentation
-  BinUtility pixelated(20, -1., 1., open, AxisDirection::AxisX);
-  pixelated += BinUtility(20, -1., 1., open, AxisDirection::AxisY);
+  std::vector<DirectedProtoAxis> segmentation;
+  segmentation.emplace_back(AxisDirection::AxisX, AxisBoundaryType::Open, -1.,
+                            1., 20);
+  segmentation.emplace_back(AxisDirection::AxisY, AxisBoundaryType::Open, -1.,
+                            1., 20);
 
   Segmentizer cl;
 
   // Test: Normal hit into the surface
   Vector2 nPosition(0.37, 0.76);
   auto nSegments =
-      cl.segments(geoCtx, *planeSurface, pixelated, {nPosition, nPosition});
+      cl.segments(geoCtx, *planeSurface, segmentation, {nPosition, nPosition});
   BOOST_CHECK_EQUAL(nSegments.size(), 1);
   BOOST_CHECK_EQUAL(nSegments[0].bin[0], 13);
   BOOST_CHECK_EQUAL(nSegments[0].bin[1], 17);
@@ -66,21 +68,21 @@ BOOST_AUTO_TEST_CASE(SegmentizerCartesian) {
   // Test: Inclined hit into the surface - negative x direction
   Vector2 ixPositionS(0.37, 0.76);
   Vector2 ixPositionE(0.02, 0.73);
-  auto ixSegments =
-      cl.segments(geoCtx, *planeSurface, pixelated, {ixPositionS, ixPositionE});
+  auto ixSegments = cl.segments(geoCtx, *planeSurface, segmentation,
+                                {ixPositionS, ixPositionE});
   BOOST_CHECK_EQUAL(ixSegments.size(), 4);
 
   // Test: Inclined hit into the surface - positive y direction
   Vector2 iyPositionS(0.37, 0.76);
   Vector2 iyPositionE(0.39, 0.91);
-  auto iySegments =
-      cl.segments(geoCtx, *planeSurface, pixelated, {iyPositionS, iyPositionE});
+  auto iySegments = cl.segments(geoCtx, *planeSurface, segmentation,
+                                {iyPositionS, iyPositionE});
   BOOST_CHECK_EQUAL(iySegments.size(), 3);
 
   // Test: Inclined hit into the surface - x/y direction
   Vector2 ixyPositionS(-0.27, 0.76);
   Vector2 ixyPositionE(-0.02, -0.73);
-  auto ixySegments = cl.segments(geoCtx, *planeSurface, pixelated,
+  auto ixySegments = cl.segments(geoCtx, *planeSurface, segmentation,
                                  {ixyPositionS, ixyPositionE});
   BOOST_CHECK_EQUAL(ixySegments.size(), 18);
 }
@@ -93,15 +95,18 @@ BOOST_AUTO_TEST_CASE(SegmentizerPolarRadial) {
       Surface::makeShared<DiscSurface>(Transform3::Identity(), radialBounds);
 
   // The segmentation
-  BinUtility strips(2, 5., 10., open, AxisDirection::AxisR);
-  strips += BinUtility(250, -0.25, 0.25, open, AxisDirection::AxisPhi);
+  std::vector<DirectedProtoAxis> segmentation;
+  segmentation.emplace_back(AxisDirection::AxisR, AxisBoundaryType::Open, 5.,
+                            10., 2);
+  segmentation.emplace_back(AxisDirection::AxisPhi, AxisBoundaryType::Open,
+                            -0.25, 0.25, 250);
 
   Segmentizer cl;
 
   // Test: Normal hit into the surface
   Vector2 nPosition(6.76, 0.5);
   auto nSegments =
-      cl.segments(geoCtx, *radialDisc, strips, {nPosition, nPosition});
+      cl.segments(geoCtx, *radialDisc, segmentation, {nPosition, nPosition});
   BOOST_CHECK_EQUAL(nSegments.size(), 1);
   BOOST_CHECK_EQUAL(nSegments[0].bin[0], 0);
   BOOST_CHECK_EQUAL(nSegments[0].bin[1], 161);
@@ -110,13 +115,14 @@ BOOST_AUTO_TEST_CASE(SegmentizerPolarRadial) {
   Vector2 sPositionS(6.76, 0.5);
   Vector2 sPositionE(7.03, -0.3);
   auto sSegment =
-      cl.segments(geoCtx, *radialDisc, strips, {sPositionS, sPositionE});
+      cl.segments(geoCtx, *radialDisc, segmentation, {sPositionS, sPositionE});
   BOOST_CHECK_EQUAL(sSegment.size(), 59);
 
   // Test: jump over R boundary, but stay in phi bin
   sPositionS = Vector2(6.76, 0.);
   sPositionE = Vector2(7.83, 0.);
-  sSegment = cl.segments(geoCtx, *radialDisc, strips, {sPositionS, sPositionE});
+  sSegment =
+      cl.segments(geoCtx, *radialDisc, segmentation, {sPositionS, sPositionE});
   BOOST_CHECK_EQUAL(sSegment.size(), 2);
 }
 
@@ -169,34 +175,32 @@ BOOST_DATA_TEST_CASE(
       }
       // 1 - write the grid
       grid.open("Segmentizer" + name + "Grid.csv");
-      if (segmentation.binningData()[0].binvalue == AxisDirection::AxisX &&
-          segmentation.binningData()[1].binvalue == AxisDirection::AxisY) {
-        double bxmin = segmentation.binningData()[0].min;
-        double bxmax = segmentation.binningData()[0].max;
-        double bymin = segmentation.binningData()[1].min;
-        double bymax = segmentation.binningData()[1].max;
-        const auto& xboundaries = segmentation.binningData()[0].boundaries();
-        const auto& yboundaries = segmentation.binningData()[1].boundaries();
-        for (const auto xval : xboundaries) {
+      if (segmentation[0].getAxisDirection() == AxisDirection::AxisX &&
+          segmentation[1].getAxisDirection() == AxisDirection::AxisY) {
+        double bxmin = segmentation[0].min();
+        double bxmax = segmentation[0].max();
+        double bymin = segmentation[1].min();
+        double bymax = segmentation[1].max();
+        const std::vector<double> xboundaries = segmentation[0].binEdges();
+        const std::vector<double> yboundaries = segmentation[1].binEdges();
+        for (const double xval : xboundaries) {
           csvHelper.writeLine(grid, {xval, bymin}, {xval, bymax});
         }
-        for (const auto yval : yboundaries) {
+        for (const double yval : yboundaries) {
           csvHelper.writeLine(grid, {bxmin, yval}, {bxmax, yval});
         }
-      } else if (segmentation.binningData()[0].binvalue ==
-                     AxisDirection::AxisR &&
-                 segmentation.binningData()[1].binvalue ==
-                     AxisDirection::AxisPhi) {
-        double brmin = segmentation.binningData()[0].min;
-        double brmax = segmentation.binningData()[0].max;
-        double bphimin = segmentation.binningData()[1].min;
-        double bphimax = segmentation.binningData()[1].max;
-        const auto& rboundaries = segmentation.binningData()[0].boundaries();
-        const auto& phiboundaries = segmentation.binningData()[1].boundaries();
-        for (const auto r : rboundaries) {
+      } else if (segmentation[0].getAxisDirection() == AxisDirection::AxisR &&
+                 segmentation[1].getAxisDirection() == AxisDirection::AxisPhi) {
+        double brmin = segmentation[0].min();
+        double brmax = segmentation[0].max();
+        double bphimin = segmentation[1].min();
+        double bphimax = segmentation[1].max();
+        const std::vector<double> rboundaries = segmentation[0].binEdges();
+        const std::vector<double> phiboundaries = segmentation[1].binEdges();
+        for (const double r : rboundaries) {
           csvHelper.writeArc(grid, r, bphimin, bphimax);
         }
-        for (const auto phi : phiboundaries) {
+        for (const double phi : phiboundaries) {
           double cphi = std::cos(phi);
           double sphi = std::sin(phi);
           csvHelper.writeLine(grid, {brmin * cphi, brmin * sphi},