SplineFitting.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_SPLINE_FITTING_H
00011 #define EIGEN_SPLINE_FITTING_H
00012 
00013 #include <numeric>
00014 
00015 #include "SplineFwd.h"
00016 
00017 #include <Eigen/QR>
00018 
00019 namespace Eigen
00020 {
00040   template <typename KnotVectorType>
00041   void KnotAveraging(const KnotVectorType& parameters, DenseIndex degree, KnotVectorType& knots)
00042   {
00043     typedef typename KnotVectorType::Scalar Scalar;
00044 
00045     knots.resize(parameters.size()+degree+1);      
00046 
00047     for (DenseIndex j=1; j<parameters.size()-degree; ++j)
00048       knots(j+degree) = parameters.segment(j,degree).mean();
00049 
00050     knots.segment(0,degree+1) = KnotVectorType::Zero(degree+1);
00051     knots.segment(knots.size()-degree-1,degree+1) = KnotVectorType::Ones(degree+1);
00052   }
00053 
00063   template <typename PointArrayType, typename KnotVectorType>
00064   void ChordLengths(const PointArrayType& pts, KnotVectorType& chord_lengths)
00065   {
00066     typedef typename KnotVectorType::Scalar Scalar;
00067 
00068     const DenseIndex n = pts.cols();
00069 
00070     // 1. compute the column-wise norms
00071     chord_lengths.resize(pts.cols());
00072     chord_lengths[0] = 0;
00073     chord_lengths.rightCols(n-1) = (pts.array().leftCols(n-1) - pts.array().rightCols(n-1)).matrix().colwise().norm();
00074 
00075     // 2. compute the partial sums
00076     std::partial_sum(chord_lengths.data(), chord_lengths.data()+n, chord_lengths.data());
00077 
00078     // 3. normalize the data
00079     chord_lengths /= chord_lengths(n-1);
00080     chord_lengths(n-1) = Scalar(1);
00081   }
00082 
00087   template <typename SplineType>
00088   struct SplineFitting
00089   {
00090     typedef typename SplineType::KnotVectorType KnotVectorType;
00091 
00100     template <typename PointArrayType>
00101     static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree);
00102 
00112     template <typename PointArrayType>
00113     static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters);
00114   };
00115 
00116   template <typename SplineType>
00117   template <typename PointArrayType>
00118   SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters)
00119   {
00120     typedef typename SplineType::KnotVectorType::Scalar Scalar;      
00121     typedef typename SplineType::BasisVectorType BasisVectorType;
00122     typedef typename SplineType::ControlPointVectorType ControlPointVectorType;      
00123 
00124     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
00125 
00126     KnotVectorType knots;
00127     KnotAveraging(knot_parameters, degree, knots);
00128 
00129     DenseIndex n = pts.cols();
00130     MatrixType A = MatrixType::Zero(n,n);
00131     for (DenseIndex i=1; i<n-1; ++i)
00132     {
00133       const DenseIndex span = SplineType::Span(knot_parameters[i], degree, knots);
00134 
00135       // The segment call should somehow be told the spline order at compile time.
00136       A.row(i).segment(span-degree, degree+1) = SplineType::BasisFunctions(knot_parameters[i], degree, knots);
00137     }
00138     A(0,0) = 1.0;
00139     A(n-1,n-1) = 1.0;
00140 
00141     HouseholderQR<MatrixType> qr(A);
00142 
00143     // Here, we are creating a temporary due to an Eigen issue.
00144     ControlPointVectorType ctrls = qr.solve(MatrixType(pts.transpose())).transpose();
00145 
00146     return SplineType(knots, ctrls);
00147   }
00148 
00149   template <typename SplineType>
00150   template <typename PointArrayType>
00151   SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree)
00152   {
00153     KnotVectorType chord_lengths; // knot parameters
00154     ChordLengths(pts, chord_lengths);
00155     return Interpolate(pts, degree, chord_lengths);
00156   }
00157 }
00158 
00159 #endif // EIGEN_SPLINE_FITTING_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:26:27