1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@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 EIGEN2_LEASTSQUARES_H
|
---|
11 | #define EIGEN2_LEASTSQUARES_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | /** \ingroup LeastSquares_Module
|
---|
16 | *
|
---|
17 | * \leastsquares_module
|
---|
18 | *
|
---|
19 | * For a set of points, this function tries to express
|
---|
20 | * one of the coords as a linear (affine) function of the other coords.
|
---|
21 | *
|
---|
22 | * This is best explained by an example. This function works in full
|
---|
23 | * generality, for points in a space of arbitrary dimension, and also over
|
---|
24 | * the complex numbers, but for this example we will work in dimension 3
|
---|
25 | * over the real numbers (doubles).
|
---|
26 | *
|
---|
27 | * So let us work with the following set of 5 points given by their
|
---|
28 | * \f$(x,y,z)\f$ coordinates:
|
---|
29 | * @code
|
---|
30 | Vector3d points[5];
|
---|
31 | points[0] = Vector3d( 3.02, 6.89, -4.32 );
|
---|
32 | points[1] = Vector3d( 2.01, 5.39, -3.79 );
|
---|
33 | points[2] = Vector3d( 2.41, 6.01, -4.01 );
|
---|
34 | points[3] = Vector3d( 2.09, 5.55, -3.86 );
|
---|
35 | points[4] = Vector3d( 2.58, 6.32, -4.10 );
|
---|
36 | * @endcode
|
---|
37 | * Suppose that we want to express the second coordinate (\f$y\f$) as a linear
|
---|
38 | * expression in \f$x\f$ and \f$z\f$, that is,
|
---|
39 | * \f[ y=ax+bz+c \f]
|
---|
40 | * for some constants \f$a,b,c\f$. Thus, we want to find the best possible
|
---|
41 | * constants \f$a,b,c\f$ so that the plane of equation \f$y=ax+bz+c\f$ fits
|
---|
42 | * best the five above points. To do that, call this function as follows:
|
---|
43 | * @code
|
---|
44 | Vector3d coeffs; // will store the coefficients a, b, c
|
---|
45 | linearRegression(
|
---|
46 | 5,
|
---|
47 | &points,
|
---|
48 | &coeffs,
|
---|
49 | 1 // the coord to express as a function of
|
---|
50 | // the other ones. 0 means x, 1 means y, 2 means z.
|
---|
51 | );
|
---|
52 | * @endcode
|
---|
53 | * Now the vector \a coeffs is approximately
|
---|
54 | * \f$( 0.495 , -1.927 , -2.906 )\f$.
|
---|
55 | * Thus, we get \f$a=0.495, b = -1.927, c = -2.906\f$. Let us check for
|
---|
56 | * instance how near points[0] is from the plane of equation \f$y=ax+bz+c\f$.
|
---|
57 | * Looking at the coords of points[0], we see that:
|
---|
58 | * \f[ax+bz+c = 0.495 * 3.02 + (-1.927) * (-4.32) + (-2.906) = 6.91.\f]
|
---|
59 | * On the other hand, we have \f$y=6.89\f$. We see that the values
|
---|
60 | * \f$6.91\f$ and \f$6.89\f$
|
---|
61 | * are near, so points[0] is very near the plane of equation \f$y=ax+bz+c\f$.
|
---|
62 | *
|
---|
63 | * Let's now describe precisely the parameters:
|
---|
64 | * @param numPoints the number of points
|
---|
65 | * @param points the array of pointers to the points on which to perform the linear regression
|
---|
66 | * @param result pointer to the vector in which to store the result.
|
---|
67 | This vector must be of the same type and size as the
|
---|
68 | data points. The meaning of its coords is as follows.
|
---|
69 | For brevity, let \f$n=Size\f$,
|
---|
70 | \f$r_i=result[i]\f$,
|
---|
71 | and \f$f=funcOfOthers\f$. Denote by
|
---|
72 | \f$x_0,\ldots,x_{n-1}\f$
|
---|
73 | the n coordinates in the n-dimensional space.
|
---|
74 | Then the resulting equation is:
|
---|
75 | \f[ x_f = r_0 x_0 + \cdots + r_{f-1}x_{f-1}
|
---|
76 | + r_{f+1}x_{f+1} + \cdots + r_{n-1}x_{n-1} + r_n. \f]
|
---|
77 | * @param funcOfOthers Determines which coord to express as a function of the
|
---|
78 | others. Coords are numbered starting from 0, so that a
|
---|
79 | value of 0 means \f$x\f$, 1 means \f$y\f$,
|
---|
80 | 2 means \f$z\f$, ...
|
---|
81 | *
|
---|
82 | * \sa fitHyperplane()
|
---|
83 | */
|
---|
84 | template<typename VectorType>
|
---|
85 | void linearRegression(int numPoints,
|
---|
86 | VectorType **points,
|
---|
87 | VectorType *result,
|
---|
88 | int funcOfOthers )
|
---|
89 | {
|
---|
90 | typedef typename VectorType::Scalar Scalar;
|
---|
91 | typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
|
---|
92 | const int size = points[0]->size();
|
---|
93 | result->resize(size);
|
---|
94 | HyperplaneType h(size);
|
---|
95 | fitHyperplane(numPoints, points, &h);
|
---|
96 | for(int i = 0; i < funcOfOthers; i++)
|
---|
97 | result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
|
---|
98 | for(int i = funcOfOthers; i < size; i++)
|
---|
99 | result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
|
---|
100 | }
|
---|
101 |
|
---|
102 | /** \ingroup LeastSquares_Module
|
---|
103 | *
|
---|
104 | * \leastsquares_module
|
---|
105 | *
|
---|
106 | * This function is quite similar to linearRegression(), so we refer to the
|
---|
107 | * documentation of this function and only list here the differences.
|
---|
108 | *
|
---|
109 | * The main difference from linearRegression() is that this function doesn't
|
---|
110 | * take a \a funcOfOthers argument. Instead, it finds a general equation
|
---|
111 | * of the form
|
---|
112 | * \f[ r_0 x_0 + \cdots + r_{n-1}x_{n-1} + r_n = 0, \f]
|
---|
113 | * where \f$n=Size\f$, \f$r_i=retCoefficients[i]\f$, and we denote by
|
---|
114 | * \f$x_0,\ldots,x_{n-1}\f$ the n coordinates in the n-dimensional space.
|
---|
115 | *
|
---|
116 | * Thus, the vector \a retCoefficients has size \f$n+1\f$, which is another
|
---|
117 | * difference from linearRegression().
|
---|
118 | *
|
---|
119 | * In practice, this function performs an hyper-plane fit in a total least square sense
|
---|
120 | * via the following steps:
|
---|
121 | * 1 - center the data to the mean
|
---|
122 | * 2 - compute the covariance matrix
|
---|
123 | * 3 - pick the eigenvector corresponding to the smallest eigenvalue of the covariance matrix
|
---|
124 | * The ratio of the smallest eigenvalue and the second one gives us a hint about the relevance
|
---|
125 | * of the solution. This value is optionally returned in \a soundness.
|
---|
126 | *
|
---|
127 | * \sa linearRegression()
|
---|
128 | */
|
---|
129 | template<typename VectorType, typename HyperplaneType>
|
---|
130 | void fitHyperplane(int numPoints,
|
---|
131 | VectorType **points,
|
---|
132 | HyperplaneType *result,
|
---|
133 | typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
|
---|
134 | {
|
---|
135 | typedef typename VectorType::Scalar Scalar;
|
---|
136 | typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
|
---|
137 | EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
|
---|
138 | ei_assert(numPoints >= 1);
|
---|
139 | int size = points[0]->size();
|
---|
140 | ei_assert(size+1 == result->coeffs().size());
|
---|
141 |
|
---|
142 | // compute the mean of the data
|
---|
143 | VectorType mean = VectorType::Zero(size);
|
---|
144 | for(int i = 0; i < numPoints; ++i)
|
---|
145 | mean += *(points[i]);
|
---|
146 | mean /= numPoints;
|
---|
147 |
|
---|
148 | // compute the covariance matrix
|
---|
149 | CovMatrixType covMat = CovMatrixType::Zero(size, size);
|
---|
150 | for(int i = 0; i < numPoints; ++i)
|
---|
151 | {
|
---|
152 | VectorType diff = (*(points[i]) - mean).conjugate();
|
---|
153 | covMat += diff * diff.adjoint();
|
---|
154 | }
|
---|
155 |
|
---|
156 | // now we just have to pick the eigen vector with smallest eigen value
|
---|
157 | SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
|
---|
158 | result->normal() = eig.eigenvectors().col(0);
|
---|
159 | if (soundness)
|
---|
160 | *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);
|
---|
161 |
|
---|
162 | // let's compute the constant coefficient such that the
|
---|
163 | // plane pass trough the mean point:
|
---|
164 | result->offset() = - (result->normal().cwise()* mean).sum();
|
---|
165 | }
|
---|
166 |
|
---|
167 | } // end namespace Eigen
|
---|
168 |
|
---|
169 | #endif // EIGEN2_LEASTSQUARES_H
|
---|