Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | { |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 |