preCICE v3.1.1
Loading...
Searching...
No Matches
Bspline.cpp
Go to the documentation of this file.
1#include "math/Bspline.hpp"
3#include "utils/assertion.hpp"
4
5#include <Eigen/Core>
6#include <algorithm>
7#include <cstdlib>
8#include <unsupported/Eigen/Splines>
10#include "utils/assertion.hpp"
11
12namespace precice::math {
13
14Bspline::Bspline(Eigen::VectorXd ts, const Eigen::MatrixXd &xs, int splineDegree)
15{
16
17 PRECICE_ASSERT(ts.size() >= 2, "Interpolation requires at least 2 samples");
18#if EIGEN_VERSION_AT_LEAST(3, 4, 0)
19 PRECICE_ASSERT(std::is_sorted(ts.begin(), ts.end()), "Timestamps must be sorted");
20#else
21 PRECICE_ASSERT(std::is_sorted(ts.data(), ts.data() + ts.size()), "Timestamps must be sorted");
22#endif
23
24 // organize data in columns. Each column represents one sample in time.
25 PRECICE_ASSERT(xs.cols() == ts.size());
26 _ndofs = xs.rows(); // number of dofs. Each dof needs its own interpolant.
27 _tsMin = ts(0);
28 _tsMax = ts(ts.size() - 1);
29 auto relativeTime = [tsMin = _tsMin, tsMax = _tsMax](double t) -> double { return (t - tsMin) / (tsMax - tsMin); };
30 ts = ts.unaryExpr(relativeTime);
31
32 //The code for computing the knots and the control points is copied from Eigens bspline interpolation with minor modifications, https://gitlab.com/libeigen/eigen/-/blob/master/unsupported/Eigen/src/Splines/SplineFitting.h
33
34 Eigen::KnotAveraging(ts, splineDegree, _knots);
35 Eigen::DenseIndex n = xs.cols();
36 Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n, n);
37
38 for (Eigen::DenseIndex i = 1; i < n - 1; ++i) {
39 //Attempt at hack... the spline dimension is not used here explicitly
40 const Eigen::DenseIndex span = Eigen::Spline<double, 1>::Span(ts[i], splineDegree, _knots);
41
42 // The segment call should somehow be told the spline order at compile time.
43 A.row(i).segment(span - splineDegree, splineDegree + 1) = Eigen::Spline<double, 1>::BasisFunctions(ts[i], splineDegree, _knots);
44 }
45 A(0, 0) = 1.0;
46 A(n - 1, n - 1) = 1.0;
47
48 auto qr = A.householderQr();
49
50 Eigen::MatrixXd controls = Eigen::MatrixXd::Zero(n, _ndofs);
51
52 for (int i = 0; i < _ndofs; i++) {
53 controls.col(i) = qr.solve(xs.row(i).transpose());
54 }
55 _ctrls = std::move(controls);
56}
57
58Eigen::VectorXd Bspline::interpolateAt(double t) const
59{
60 // transform t to the relative interval [0; 1]
61 const double tRelative = std::clamp((t - _tsMin) / (_tsMax - _tsMin), 0.0, 1.0);
62
63 Eigen::VectorXd interpolated(_ndofs);
64 constexpr int splineDimension = 1;
65
66 for (int i = 0; i < _ndofs; i++) {
67 interpolated[i] = Eigen::Spline<double, splineDimension>(_knots, _ctrls.col(i))(tRelative)[0];
68 }
69
70 return interpolated;
71}
72} // namespace precice::math
Eigen::Vector2d ts
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T clamp(T... args)
Eigen::MatrixXd _ctrls
Definition Bspline.hpp:31
Eigen::VectorXd _knots
Definition Bspline.hpp:30
Eigen::VectorXd interpolateAt(double t) const
Samples the B-Spline interpolation.
Definition Bspline.cpp:58
Bspline(Eigen::VectorXd ts, const Eigen::MatrixXd &xs, int splineDegree)
Initialises the B-Spline interpolation with the given data (x0,t0), (x1,t1), ..., (xn,...
Definition Bspline.cpp:14
A C++ 11 implementation of the non-owning C++20 std::span type.
Definition span.hpp:284
T is_sorted(T... args)
provides general mathematical constants and functions.
Definition barycenter.cpp:9