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