1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
|
---|
5 | //
|
---|
6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
9 |
|
---|
10 | #ifndef EIGEN_SPLINE_FITTING_H
|
---|
11 | #define EIGEN_SPLINE_FITTING_H
|
---|
12 |
|
---|
13 | #include <numeric>
|
---|
14 |
|
---|
15 | #include "SplineFwd.h"
|
---|
16 |
|
---|
17 | #include <Eigen/QR>
|
---|
18 |
|
---|
19 | namespace Eigen
|
---|
20 | {
|
---|
21 | /**
|
---|
22 | * \brief Computes knot averages.
|
---|
23 | * \ingroup Splines_Module
|
---|
24 | *
|
---|
25 | * The knots are computed as
|
---|
26 | * \f{align*}
|
---|
27 | * u_0 & = \hdots = u_p = 0 \\
|
---|
28 | * u_{m-p} & = \hdots = u_{m} = 1 \\
|
---|
29 | * u_{j+p} & = \frac{1}{p}\sum_{i=j}^{j+p-1}\bar{u}_i \quad\quad j=1,\hdots,n-p
|
---|
30 | * \f}
|
---|
31 | * where \f$p\f$ is the degree and \f$m+1\f$ the number knots
|
---|
32 | * of the desired interpolating spline.
|
---|
33 | *
|
---|
34 | * \param[in] parameters The input parameters. During interpolation one for each data point.
|
---|
35 | * \param[in] degree The spline degree which is used during the interpolation.
|
---|
36 | * \param[out] knots The output knot vector.
|
---|
37 | *
|
---|
38 | * \sa Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data
|
---|
39 | **/
|
---|
40 | template <typename KnotVectorType>
|
---|
41 | void KnotAveraging(const KnotVectorType& parameters, DenseIndex degree, KnotVectorType& knots)
|
---|
42 | {
|
---|
43 | knots.resize(parameters.size()+degree+1);
|
---|
44 |
|
---|
45 | for (DenseIndex j=1; j<parameters.size()-degree; ++j)
|
---|
46 | knots(j+degree) = parameters.segment(j,degree).mean();
|
---|
47 |
|
---|
48 | knots.segment(0,degree+1) = KnotVectorType::Zero(degree+1);
|
---|
49 | knots.segment(knots.size()-degree-1,degree+1) = KnotVectorType::Ones(degree+1);
|
---|
50 | }
|
---|
51 |
|
---|
52 | /**
|
---|
53 | * \brief Computes chord length parameters which are required for spline interpolation.
|
---|
54 | * \ingroup Splines_Module
|
---|
55 | *
|
---|
56 | * \param[in] pts The data points to which a spline should be fit.
|
---|
57 | * \param[out] chord_lengths The resulting chord lenggth vector.
|
---|
58 | *
|
---|
59 | * \sa Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data
|
---|
60 | **/
|
---|
61 | template <typename PointArrayType, typename KnotVectorType>
|
---|
62 | void ChordLengths(const PointArrayType& pts, KnotVectorType& chord_lengths)
|
---|
63 | {
|
---|
64 | typedef typename KnotVectorType::Scalar Scalar;
|
---|
65 |
|
---|
66 | const DenseIndex n = pts.cols();
|
---|
67 |
|
---|
68 | // 1. compute the column-wise norms
|
---|
69 | chord_lengths.resize(pts.cols());
|
---|
70 | chord_lengths[0] = 0;
|
---|
71 | chord_lengths.rightCols(n-1) = (pts.array().leftCols(n-1) - pts.array().rightCols(n-1)).matrix().colwise().norm();
|
---|
72 |
|
---|
73 | // 2. compute the partial sums
|
---|
74 | std::partial_sum(chord_lengths.data(), chord_lengths.data()+n, chord_lengths.data());
|
---|
75 |
|
---|
76 | // 3. normalize the data
|
---|
77 | chord_lengths /= chord_lengths(n-1);
|
---|
78 | chord_lengths(n-1) = Scalar(1);
|
---|
79 | }
|
---|
80 |
|
---|
81 | /**
|
---|
82 | * \brief Spline fitting methods.
|
---|
83 | * \ingroup Splines_Module
|
---|
84 | **/
|
---|
85 | template <typename SplineType>
|
---|
86 | struct SplineFitting
|
---|
87 | {
|
---|
88 | typedef typename SplineType::KnotVectorType KnotVectorType;
|
---|
89 |
|
---|
90 | /**
|
---|
91 | * \brief Fits an interpolating Spline to the given data points.
|
---|
92 | *
|
---|
93 | * \param pts The points for which an interpolating spline will be computed.
|
---|
94 | * \param degree The degree of the interpolating spline.
|
---|
95 | *
|
---|
96 | * \returns A spline interpolating the initially provided points.
|
---|
97 | **/
|
---|
98 | template <typename PointArrayType>
|
---|
99 | static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree);
|
---|
100 |
|
---|
101 | /**
|
---|
102 | * \brief Fits an interpolating Spline to the given data points.
|
---|
103 | *
|
---|
104 | * \param pts The points for which an interpolating spline will be computed.
|
---|
105 | * \param degree The degree of the interpolating spline.
|
---|
106 | * \param knot_parameters The knot parameters for the interpolation.
|
---|
107 | *
|
---|
108 | * \returns A spline interpolating the initially provided points.
|
---|
109 | **/
|
---|
110 | template <typename PointArrayType>
|
---|
111 | static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters);
|
---|
112 | };
|
---|
113 |
|
---|
114 | template <typename SplineType>
|
---|
115 | template <typename PointArrayType>
|
---|
116 | SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters)
|
---|
117 | {
|
---|
118 | typedef typename SplineType::KnotVectorType::Scalar Scalar;
|
---|
119 | typedef typename SplineType::ControlPointVectorType ControlPointVectorType;
|
---|
120 |
|
---|
121 | typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
|
---|
122 |
|
---|
123 | KnotVectorType knots;
|
---|
124 | KnotAveraging(knot_parameters, degree, knots);
|
---|
125 |
|
---|
126 | DenseIndex n = pts.cols();
|
---|
127 | MatrixType A = MatrixType::Zero(n,n);
|
---|
128 | for (DenseIndex i=1; i<n-1; ++i)
|
---|
129 | {
|
---|
130 | const DenseIndex span = SplineType::Span(knot_parameters[i], degree, knots);
|
---|
131 |
|
---|
132 | // The segment call should somehow be told the spline order at compile time.
|
---|
133 | A.row(i).segment(span-degree, degree+1) = SplineType::BasisFunctions(knot_parameters[i], degree, knots);
|
---|
134 | }
|
---|
135 | A(0,0) = 1.0;
|
---|
136 | A(n-1,n-1) = 1.0;
|
---|
137 |
|
---|
138 | HouseholderQR<MatrixType> qr(A);
|
---|
139 |
|
---|
140 | // Here, we are creating a temporary due to an Eigen issue.
|
---|
141 | ControlPointVectorType ctrls = qr.solve(MatrixType(pts.transpose())).transpose();
|
---|
142 |
|
---|
143 | return SplineType(knots, ctrls);
|
---|
144 | }
|
---|
145 |
|
---|
146 | template <typename SplineType>
|
---|
147 | template <typename PointArrayType>
|
---|
148 | SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree)
|
---|
149 | {
|
---|
150 | KnotVectorType chord_lengths; // knot parameters
|
---|
151 | ChordLengths(pts, chord_lengths);
|
---|
152 | return Interpolate(pts, degree, chord_lengths);
|
---|
153 | }
|
---|
154 | }
|
---|
155 |
|
---|
156 | #endif // EIGEN_SPLINE_FITTING_H
|
---|