Developer Documentation

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
PCA.cc
1 /*===========================================================================*\
2 * *
3 * OpenFlipper *
4 * Copyright (c) 2001-2015, RWTH-Aachen University *
5 * Department of Computer Graphics and Multimedia *
6 * All rights reserved. *
7 * www.openflipper.org *
8 * *
9 *---------------------------------------------------------------------------*
10 * This file is part of OpenFlipper. *
11 *---------------------------------------------------------------------------*
12 * *
13 * Redistribution and use in source and binary forms, with or without *
14 * modification, are permitted provided that the following conditions *
15 * are met: *
16 * *
17 * 1. Redistributions of source code must retain the above copyright notice, *
18 * this list of conditions and the following disclaimer. *
19 * *
20 * 2. Redistributions in binary form must reproduce the above copyright *
21 * notice, this list of conditions and the following disclaimer in the *
22 * documentation and/or other materials provided with the distribution. *
23 * *
24 * 3. Neither the name of the copyright holder nor the names of its *
25 * contributors may be used to endorse or promote products derived from *
26 * this software without specific prior written permission. *
27 * *
28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 \*===========================================================================*/
41
42 /*===========================================================================*\
43 * *
44 * $Revision$ *
45 * $LastChangedBy$ *
46 * $Date$ *
47 * *
48 \*===========================================================================*/
49
50
51
52
53 //=============================================================================
54 //
55 // CLASS PCA - IMPLEMENTATION
56 //
57 //=============================================================================
58
59 #define PCA_C
60
61 //== INCLUDES =================================================================
62
63 #include "PCA.hh"
64
65 extern "C" {
66
67 typedef void (*U_fp)(...);
68
69 int dsyev_( char *jobz, char *uplo, long int *n, double *a,
70 long int *lda, double *w, double *work, long int *lwork,
71 long int *info);
72
73 }
74
75 // typedef gmm::dense_matrix< double > Matrix;
76 // typedef std::vector< double > Vector;
77
79 static double * p_work_double = 0;
80
82 static int m_workSize_double = 0;
83
84 //== NAMESPACES ===============================================================
85
86 namespace Pca {
87
88 //== IMPLEMENTATION ==========================================================
89
90 template < typename VectorT >
91 PCA< VectorT >::PCA( std::vector< VectorT >& _points , VectorT& _first , VectorT& _second , VectorT& _third)
92 {
93 pca(_points,_first,_second,_third);
94 }
95
96 template < typename VectorT >
97 PCA< VectorT >::~PCA()
98 {
99 if ( p_work_double != 0){
100 delete [] p_work_double;
101 p_work_double = 0;
102 m_workSize_double = 0;
103 }
104 }
105
106 template < typename VectorT >
107 inline
108 VectorT
109 PCA< VectorT >::
110 center_of_gravity(const std::vector< VectorT >& _points ) {
111 // compute cog of Points
112 VectorT cog(0.0, 0.0, 0.0);
113
114 for (uint i = 0 ; i < _points.size() ; ++i ) {
115 cog = cog + _points[i];
116 }
117
118 cog = cog / ( typename VectorT::value_type )_points.size();
119
120 return cog;
121 }
122
123 template < typename VectorT >
124 bool
125 PCA< VectorT >::SymRightEigenproblem( Matrix &_mat_A, Matrix & _mat_VR,
126 std::vector< double > & _vec_EV )
127 {
128 char jobz, uplo;
129 long int n, lda;
130 long int lwork, info;
131 double size;
132
133 jobz = 'V';
134 uplo = 'U';
135 n = gmm::mat_nrows( _mat_A );
136 lda = n;
137
138 info = 0;
139 // compute optimal size of working array
140 lwork = -1;
141 dsyev_( &jobz, &uplo, &n, &_mat_A( 0, 0 ), &lda, &_vec_EV[ 0 ],
142 &size, &lwork, &info );
143
144 if( info != 0 )
145 return false;
146
147
148 if( (long int) size > m_workSize_double )
149 {
150 if( p_work_double)
151 delete [] p_work_double;
152
153 m_workSize_double = (long int) size;
154
155 p_work_double= new double[ m_workSize_double ];
156 }
157 lwork = m_workSize_double;
158
159 // compute eigenvalues / eigenvectors
160 dsyev_( &jobz, &uplo, &n, &_mat_A( 0, 0 ), &lda, &_vec_EV[ 0 ],
161 p_work_double, &lwork, &info );
162
163 if( info != 0 )
164 return false;
165
166 // copy eigenvectors to matrix
167 gmm::copy( _mat_A, _mat_VR );
168
169 return true;
170 }
171
172
173 template < typename VectorT >
174 void
175 PCA< VectorT >::pca(std::vector< VectorT >& _points , VectorT& _first , VectorT& _second , VectorT& _third)
176 {
177 // compute Mean of Points
178 const VectorT cog = center_of_gravity(_points);
179
180 //build covariance Matrix
181 Matrix points(3 , _points.size() );
182
183 for ( uint i = 0 ; i < _points.size() ; ++i)
184 {
185 points(0 , i ) = ( _points[i] - cog) [0];
186 points(1 , i ) = ( _points[i] - cog) [1];
187 points(2 , i ) = ( _points[i] - cog) [2];
188 }
189
190 Matrix cov(3,3);
191 gmm::mult(points,gmm::transposed(points),cov );
192
193 Matrix EV(3,3);
194 std::vector< double > ew(3);
195
196
197 if ( !SymRightEigenproblem( cov, EV ,ew ) )
198 std::cerr << "Error: Could not compute Eigenvectors for PCA" << std::endl;
199
200 int maximum,middle,minimum;
201
202 if ( (ew[0] > ew[1]) && (ew[0] > ew[2]) ) {
203 maximum = 0;
204 } else {
205 if ( (ew[1] > ew[0]) && (ew[1] > ew[2]) )
206 maximum = 1;
207 else
208 maximum = 2;
209 }
210
211 if ( (ew[0] < ew[1]) && (ew[0] < ew[2]) )
212 minimum = 0;
213 else if ( (ew[1] < ew[0]) && (ew[1] < ew[2]) )
214 minimum = 1;
215 else
216 minimum = 2;
217
218 if ( (minimum != 0 ) && ( maximum != 0 ) )
219 middle = 0;
220 else
221 if ( (minimum != 1 ) && ( maximum != 1 ) )
222 middle = 1;
223 else
224 middle = 2;
225
226 _first = VectorT( EV(0,maximum) , EV(1,maximum) , EV(2,maximum) );
227 _second = VectorT( EV(0,middle) , EV(1,middle) , EV(2,middle ) );
228
229 _first = _first.normalize();
230 _second = _second.normalize();
231 _third = _first % _second;
232
233 //remember the eigenvalues
234 lastEigenValues_.clear();
235 lastEigenValues_.push_back( ew[maximum] );
236 lastEigenValues_.push_back( ew[middle] );
237 lastEigenValues_.push_back( ew[minimum] );
238 }
239 //-----------------------------------------------------------------------------
240
241 template < typename VectorT >
242 inline
243 std::vector<double>&
244 PCA< VectorT >::getLastEigenValues()
245 {
246 return lastEigenValues_;
247 }
248
249
250 //=============================================================================
251
252 } //namespace Pca
253
254 //=============================================================================
VectorT
Pca::PCA::pca
void pca(std::vector< VectorT > &_points, VectorT &_first, VectorT &_second, VectorT &_third)
Definition: PCA.cc:175
PCA.hh
Classes for doing Principal component analysis.
Pca::PCA::PCA
PCA()
Constructor.
Definition: PCA.hh:97
Pca::PCA::~PCA
~PCA()
Destructor.
Definition: PCA.cc:97
Pca::PCA::center_of_gravity
VectorT center_of_gravity(const std::vector< VectorT > &_points)
Definition: PCA.cc:110
Pca
Namespace for principal Component Analysis.
Definition: PCA.cc:86
Pca::PCA
Class for principal component Analysis.
Definition: PCA.hh:89