LeastSquares.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN2_LEASTSQUARES_H
00026 #define EIGEN2_LEASTSQUARES_H
00027 
00097 template<typename VectorType>
00098 void linearRegression(int numPoints,
00099                       VectorType **points,
00100                       VectorType *result,
00101                       int funcOfOthers )
00102 {
00103   typedef typename VectorType::Scalar Scalar;
00104   typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
00105   const int size = points[0]->size();
00106   result->resize(size);
00107   HyperplaneType h(size);
00108   fitHyperplane(numPoints, points, &h);
00109   for(int i = 0; i < funcOfOthers; i++)
00110     result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
00111   for(int i = funcOfOthers; i < size; i++)
00112     result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
00113 }
00114 
00142 template<typename VectorType, typename HyperplaneType>
00143 void fitHyperplane(int numPoints,
00144                    VectorType **points,
00145                    HyperplaneType *result,
00146                    typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
00147 {
00148   typedef typename VectorType::Scalar Scalar;
00149   typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
00150   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
00151   ei_assert(numPoints >= 1);
00152   int size = points[0]->size();
00153   ei_assert(size+1 == result->coeffs().size());
00154 
00155   // compute the mean of the data
00156   VectorType mean = VectorType::Zero(size);
00157   for(int i = 0; i < numPoints; ++i)
00158     mean += *(points[i]);
00159   mean /= numPoints;
00160 
00161   // compute the covariance matrix
00162   CovMatrixType covMat = CovMatrixType::Zero(size, size);
00163   VectorType remean = VectorType::Zero(size);
00164   for(int i = 0; i < numPoints; ++i)
00165   {
00166     VectorType diff = (*(points[i]) - mean).conjugate();
00167     covMat += diff * diff.adjoint();
00168   }
00169 
00170   // now we just have to pick the eigen vector with smallest eigen value
00171   SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
00172   result->normal() = eig.eigenvectors().col(0);
00173   if (soundness)
00174     *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);
00175 
00176   // let's compute the constant coefficient such that the
00177   // plane pass trough the mean point:
00178   result->offset() = - (result->normal().cwise()* mean).sum();
00179 }
00180 
00181 
00182 #endif // EIGEN2_LEASTSQUARES_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:56