LeastSquares.h

Go to the documentation of this file.

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN2_LEASTSQUARES_H
#define EIGEN2_LEASTSQUARES_H

namespace Eigen { 

template<typename VectorType>
void linearRegression(int numPoints,
                      VectorType **points,
                      VectorType *result,
                      int funcOfOthers )
{
  typedef typename VectorType::Scalar Scalar;
  typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
  const int size = points[0]->size();
  result->resize(size);
  HyperplaneType h(size);
  fitHyperplane(numPoints, points, &h);
  for(int i = 0; i < funcOfOthers; i++)
    result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
  for(int i = funcOfOthers; i < size; i++)
    result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
}

template<typename VectorType, typename HyperplaneType>
void fitHyperplane(int numPoints,
                   VectorType **points,
                   HyperplaneType *result,
                   typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
{
  typedef typename VectorType::Scalar Scalar;
  typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
  ei_assert(numPoints >= 1);
  int size = points[0]->size();
  ei_assert(size+1 == result->coeffs().size());

  // compute the mean of the data
  VectorType mean = VectorType::Zero(size);
  for(int i = 0; i < numPoints; ++i)
    mean += *(points[i]);
  mean /= numPoints;

  // compute the covariance matrix
  CovMatrixType covMat = CovMatrixType::Zero(size, size);
  VectorType remean = VectorType::Zero(size);
  for(int i = 0; i < numPoints; ++i)
  {
    VectorType diff = (*(points[i]) - mean).conjugate();
    covMat += diff * diff.adjoint();
  }

  // now we just have to pick the eigen vector with smallest eigen value
  SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
  result->normal() = eig.eigenvectors().col(0);
  if (soundness)
    *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);

  // let's compute the constant coefficient such that the
  // plane pass trough the mean point:
  result->offset() = - (result->normal().cwise()* mean).sum();
}

} // end namespace Eigen

#endif // EIGEN2_LEASTSQUARES_H