Skip to content

Commit

Permalink
feat: move spline to helper in Visualization (#3950)
Browse files Browse the repository at this point in the history
This PR moves the spline-interpolation from `Examples/Plugins/Obj` to `Core/Visualization` because it can be used in other visualisation areas as well.

I've added also a UnitTest since this is from some `experimetnal/usupported` Eigen directory, but it works nicely!
  • Loading branch information
asalzburger authored Dec 6, 2024
1 parent 5c41085 commit a646e06
Show file tree
Hide file tree
Showing 6 changed files with 211 additions and 52 deletions.
96 changes: 96 additions & 0 deletions Core/include/Acts/Visualization/Interpolation3D.hpp
Original file line number Diff line number Diff line change
@@ -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
10 changes: 8 additions & 2 deletions Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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.
Expand Down
60 changes: 11 additions & 49 deletions Examples/Io/Obj/src/ObjSimHitWriter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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)
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 3 additions & 1 deletion Examples/Python/src/Output.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,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<>());
Expand Down
1 change: 1 addition & 0 deletions Tests/UnitTests/Core/Visualization/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
Expand Down
92 changes: 92 additions & 0 deletions Tests/UnitTests/Core/Visualization/Interpolation3DTests.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
// 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_CASE(SplineInterpolationErrors) {
std::vector<std::array<double, 3u>> inputs;

// Test with single point
inputs.push_back({0., 0., 0.});
auto result = Acts::Interpolation3D::spline(inputs, 10);
BOOST_CHECK_EQUAL(result.size(), 1);

// Test with two points
inputs.push_back({1., 1., 1.});
result = Acts::Interpolation3D::spline(inputs, 10);
BOOST_CHECK_EQUAL(result.size(), 2);
}

BOOST_AUTO_TEST_SUITE_END()

} // namespace Acts::Test

0 comments on commit a646e06

Please sign in to comment.