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 // 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 EIGEN2_LEASTSQUARES_H
00011 #define EIGEN2_LEASTSQUARES_H
00012 
00013 namespace Eigen { 
00014 
00084 template<typename VectorType>
00085 void linearRegression(int numPoints,
00086                       VectorType **points,
00087                       VectorType *result,
00088                       int funcOfOthers )
00089 {
00090   typedef typename VectorType::Scalar Scalar;
00091   typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
00092   const int size = points[0]->size();
00093   result->resize(size);
00094   HyperplaneType h(size);
00095   fitHyperplane(numPoints, points, &h);
00096   for(int i = 0; i < funcOfOthers; i++)
00097     result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
00098   for(int i = funcOfOthers; i < size; i++)
00099     result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
00100 }
00101 
00129 template<typename VectorType, typename HyperplaneType>
00130 void fitHyperplane(int numPoints,
00131                    VectorType **points,
00132                    HyperplaneType *result,
00133                    typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
00134 {
00135   typedef typename VectorType::Scalar Scalar;
00136   typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
00137   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
00138   ei_assert(numPoints >= 1);
00139   int size = points[0]->size();
00140   ei_assert(size+1 == result->coeffs().size());
00141 
00142   // compute the mean of the data
00143   VectorType mean = VectorType::Zero(size);
00144   for(int i = 0; i < numPoints; ++i)
00145     mean += *(points[i]);
00146   mean /= numPoints;
00147 
00148   // compute the covariance matrix
00149   CovMatrixType covMat = CovMatrixType::Zero(size, size);
00150   VectorType remean = VectorType::Zero(size);
00151   for(int i = 0; i < numPoints; ++i)
00152   {
00153     VectorType diff = (*(points[i]) - mean).conjugate();
00154     covMat += diff * diff.adjoint();
00155   }
00156 
00157   // now we just have to pick the eigen vector with smallest eigen value
00158   SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
00159   result->normal() = eig.eigenvectors().col(0);
00160   if (soundness)
00161     *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);
00162 
00163   // let's compute the constant coefficient such that the
00164   // plane pass trough the mean point:
00165   result->offset() = - (result->normal().cwise()* mean).sum();
00166 }
00167 
00168 } // end namespace Eigen
00169 
00170 #endif // EIGEN2_LEASTSQUARES_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Wed Sep 16 2015 07:11:07