feat: move spline to helper in Visualization

Core/include/Acts/Visualization/Interpolation3D.hpp

@@ -0,0 +1,96 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+
+#include <unsupported/Eigen/Splines>
+
+namespace Acts::Interpolation3D {
+
+/// @brief Helper function to interpolate points using a spline
+/// from Eigen
+///
+/// The only requirement is that the input trajectory type has
+/// a method empty() and size() and that the elements can be
+/// accessed with operator[] and have themselves a operator[] to
+/// access the coordinates.
+///
+/// @tparam input_trajectory_type input trajectory type
+///
+/// @param inputsRaw input vector points
+/// @param nPoints number of interpolation points
+/// @param keepOriginalHits keep the original hits in the trajectory
+///
+/// @return std::vector<Acts::Vector3> interpolated points
+template <typename trajectory_type>
+trajectory_type spline(const trajectory_type& inputsRaw, std::size_t nPoints,
+                       bool keepOriginalHits = false) {
+  trajectory_type output;
+  if (inputsRaw.empty()) {
+    return output;
+  }
+
+  using InputVectorType = typename trajectory_type::value_type;
+
+  std::vector<Vector3> inputs;
+  // If input type is a vector of Vector3 we can use it directly
+  if constexpr (std::is_same_v<trajectory_type, std::vector<Vector3>>) {
+    inputs = inputsRaw;
+  } else {
+    inputs.reserve(inputsRaw.size());
+    for (const auto& input : inputsRaw) {
+      inputs.push_back(Vector3(input[0], input[1], input[2]));
+    }
+  }
+
+  // Don't do anything if we have less than 3 points or less interpolation
+  // points than input points
+  if (inputsRaw.size() < 3 || nPoints <= inputsRaw.size()) {
+    return inputsRaw;
+  } else {
+    Eigen::MatrixXd points(3, inputs.size());
+    for (std::size_t i = 0; i < inputs.size(); ++i) {
+      points.col(i) = inputs[i].transpose();
+    }
+    Eigen::Spline<double, 3> spline3D =
+        Eigen::SplineFitting<Eigen::Spline<double, 3>>::Interpolate(points, 2);
+
+    double step = 1. / (nPoints - 1);
+    for (std::size_t i = 0; i < nPoints; ++i) {
+      double t = i * step;
+      InputVectorType point;
+      point[0] = spline3D(t)[0];
+      point[1] = spline3D(t)[1];
+      point[2] = spline3D(t)[2];
+      output.push_back(point);
+    }
+  }
+  // If we want to keep the original hits, we add them to the output
+  // (first and last are there anyway)
+  if (keepOriginalHits) {
+    output.insert(output.begin(), inputsRaw.begin() + 1, inputsRaw.end() - 1);
+    // We need to sort the output in distance to first
+    std::sort(output.begin(), output.end(),
+              [&inputs](const auto& a, const auto& b) {
+                const auto ifront = inputs.front();
+                double da2 = (a[0] - ifront[0]) * (a[0] - ifront[0]) +
+                             (a[1] - ifront[1]) * (a[1] - ifront[1]) +
+                             (a[2] - ifront[2]) * (a[2] - ifront[2]);
+                double db2 = (b[0] - ifront[0]) * (b[0] - ifront[0]) +
+                             (b[1] - ifront[1]) * (b[1] - ifront[1]) +
+                             (b[2] - ifront[2]) * (b[2] - ifront[2]);
+                return da2 < db2;
+              });
+  }
+
+  return output;
+}
+
+}  // namespace Acts::Interpolation3D

Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp

@@ -32,6 +32,9 @@ struct AlgorithmContext;
 ///     event000000002-<stem>.obj
 ///     event000000002-<stem>_trajectory.obj
 ///
+///
+/// The trajectory can be smoothed using a spline interpolation, where
+/// nInterpolatedPoints points are added between each hit.
 class ObjSimHitWriter : public WriterT<SimHitContainer> {
  public:
   struct Config {
@@ -49,8 +52,11 @@ class ObjSimHitWriter : public WriterT<SimHitContainer> {
     double momentumThreshold = 0.05 * Acts::UnitConstants::GeV;
     /// Momentum threshold for trajectories
     double momentumThresholdTraj = 0.05 * Acts::UnitConstants::GeV;
-    /// Number of points to interpolate between hits
-    std::size_t nInterpolatedPoints = 10;
+    /// Number of points to interpolated between hits to smooth the
+    /// trajectory view in the obj file.
+    std::size_t nInterpolatedPoints = 4;
+    /// Keep the original hits in the trajectory file
+    bool keepOriginalHits = false;
   };
 
   /// Construct the particle writer.

Examples/Io/Obj/src/ObjSimHitWriter.cpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Common.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Visualization/Interpolation3D.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsExamples/Utilities/Paths.hpp"
@@ -22,47 +23,6 @@
 #include <unordered_map>
 #include <vector>
 
-#include <unsupported/Eigen/Splines>
-
-namespace {
-
-/// @brief Helper function to interpolate points
-///
-/// @tparam input_vector_type
-/// @param inputs input vector points
-/// @param nPoints number of interpolation points
-///
-/// @return std::vector<Acts::Vector3> interpolated points
-template <typename input_vector_type>
-std::vector<Acts::Vector3> interpolatedPoints(
-    const std::vector<input_vector_type>& inputs, std::size_t nPoints) {
-  std::vector<Acts::Vector3> output;
-
-  if (nPoints < 2) {
-    // No interpolation done return simply the output vector
-    for (const auto& input : inputs) {
-      output.push_back(input.template head<3>());
-    }
-
-  } else {
-    Eigen::MatrixXd points(3, inputs.size());
-    for (std::size_t i = 0; i < inputs.size(); ++i) {
-      points.col(i) = inputs[i].template head<3>().transpose();
-    }
-    Eigen::Spline<double, 3> spline3D =
-        Eigen::SplineFitting<Eigen::Spline<double, 3>>::Interpolate(points, 2);
-
-    double step = 1. / (nPoints - 1);
-    for (std::size_t i = 0; i < nPoints; ++i) {
-      double t = i * step;
-      output.push_back(spline3D(t));
-    }
-  }
-  return output;
-}
-
-}  // namespace
-
 ActsExamples::ObjSimHitWriter::ObjSimHitWriter(
     const ActsExamples::ObjSimHitWriter::Config& config,
     Acts::Logging::Level level)
@@ -152,15 +112,17 @@ ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::writeT(
       }
       osHits << '\n';
 
-      // Interpolate the points
-      std::vector<Acts::Vector3> trajectory;
-      if (pHits.size() < 3) {
-        for (const auto& hit : pHits) {
-          trajectory.push_back(hit.template head<3>());
-        }
+      // Interpolate the points, a minimum number of 3 hits is necessary for
+      // that
+      std::vector<Acts::Vector4> trajectory;
+      if (pHits.size() < 3 || m_cfg.nInterpolatedPoints == 0) {
+        trajectory = pHits;
       } else {
-        trajectory =
-            interpolatedPoints(pHits, pHits.size() * m_cfg.nInterpolatedPoints);
+        // The total number of points is the number of hits times the number of
+        // interpolated points plus the number of hits
+        trajectory = Acts::Interpolation3D::spline(
+            pHits, pHits.size() * (m_cfg.nInterpolatedPoints + 1) - 1,
+            m_cfg.keepOriginalHits);
       }
 
       osTrajectory << "o particle_trajectory_" << pId << std::endl;

Examples/Python/src/Output.cpp

@@ -114,7 +114,9 @@ void addOutput(Context& ctx) {
 
   ACTS_PYTHON_DECLARE_WRITER(ActsExamples::ObjSimHitWriter, mex,
                              "ObjSimHitWriter", inputSimHits, outputDir,
-                             outputStem, outputPrecision, drawConnections);
+                             outputStem, outputPrecision, drawConnections,
+                             momentumThreshold, momentumThresholdTraj,
+                             nInterpolatedPoints, keepOriginalHits);
 
   {
     auto c = py::class_<ViewConfig>(m, "ViewConfig").def(py::init<>());

Tests/UnitTests/Core/Visualization/CMakeLists.txt

@@ -1,5 +1,6 @@
 add_unittest(Visualization3D Visualization3DTests.cpp)
 add_unittest(EventDataView3D EventDataView3DTests.cpp)
+add_unittest(Interpolation3D Interpolation3DTests.cpp)
 add_unittest(PrimitivesView3D PrimitivesView3DTests.cpp)
 add_unittest(SurfaceView3D SurfaceView3DTests.cpp)
 add_unittest(TrackingGeometryView3D TrackingGeometryView3DTests.cpp)

Tests/UnitTests/Core/Visualization/Interpolation3DTests.cpp

@@ -0,0 +1,78 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Visualization/Interpolation3D.hpp"
+
+#include <numbers>
+
+namespace Acts::Test {
+
+BOOST_AUTO_TEST_SUITE(Visualization)
+
+BOOST_AUTO_TEST_CASE(SplineInterpolationEigen) {
+  /// Define the input vector
+  double R = 10.;
+  std::vector<Acts::Vector3> inputs;
+
+  // Interpolate the points options
+  std::vector<Acts::Vector3> trajectory;
+
+  // Empty in empty out check
+  trajectory = Acts::Interpolation3D::spline(inputs, 10);
+  BOOST_CHECK(trajectory.empty());
+
+  for (double phi = 0; phi < 2 * std::numbers::pi;
+       phi += std::numbers::pi / 4) {
+    inputs.push_back(Acts::Vector3(R * cos(phi), R * sin(phi), 0.));
+  }
+
+  // (0) - nothing happens
+  trajectory = Acts::Interpolation3D::spline(inputs, 1);
+  // Check input and output size are the same
+  BOOST_CHECK_EQUAL(trajectory.size(), inputs.size());
+
+  // (1) - interpolate between the points with 12 points in total
+  trajectory = Acts::Interpolation3D::spline(inputs, 12);
+  // Check the output size is correct
+  BOOST_CHECK_EQUAL(trajectory.size(), 12);
+
+  for (const auto& point : trajectory) {
+    // Check the interpolated points are on the circle
+    // with a tolerance of course
+    CHECK_CLOSE_ABS(point.norm(), R, 0.1);
+    // Verify points remain in the XY plane
+    CHECK_CLOSE_ABS(point.z(), 0., 0.1);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(SplineInterpolationArray) {
+  /// Define the input vector
+  std::vector<std::array<double, 3u>> inputs;
+
+  for (double x = 0; x < 10; x += 1) {
+    inputs.push_back({x, x * x, 0.});
+  }
+
+  // This time we keep the original hits
+  auto trajectory = Acts::Interpolation3D::spline(inputs, 100, true);
+
+  // Check the outpu type is correct
+  constexpr bool isOutput =
+      std::is_same_v<decltype(trajectory), decltype(inputs)>;
+  BOOST_CHECK(isOutput);
+
+  // Check the output size is correct
+  BOOST_CHECK_EQUAL(trajectory.size(), 108);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}  // namespace Acts::Test