Developer Documentation

PCAT_impl.hh
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
45
46
47 //=============================================================================
48 //
49 // CLASS PCA - IMPLEMENTATION
50 //
51 //=============================================================================
52
53 #define PCA_C
54
55 //== INCLUDES =================================================================
56
57 #include "PCA.hh"
58
59 extern "C" {
60
61 typedef void (*U_fp)(...);
62
63 int dsyev_( char *jobz, char *uplo, long int *n, double *a,
64 long int *lda, double *w, double *work, long int *lwork,
65 long int *info);
66
67 }
68
69 // typedef gmm::dense_matrix< double > Matrix;
70 // typedef std::vector< double > Vector;
71
73 static double * p_work_double = 0;
74
76 static int m_workSize_double = 0;
77
78 //== NAMESPACES ===============================================================
79
80 namespace Pca {
81
82 //== IMPLEMENTATION ==========================================================
83
84 template < typename VectorT >
85 PCA< VectorT >::PCA( std::vector< VectorT >& _points , VectorT& _first , VectorT& _second , VectorT& _third)
86 {
87 pca(_points,_first,_second,_third);
88 }
89
90 template < typename VectorT >
91 PCA< VectorT >::~PCA()
92 {
93 if ( p_work_double != 0){
94 delete [] p_work_double;
95 p_work_double = 0;
96 m_workSize_double = 0;
97 }
98 }
99
100 template < typename VectorT >
101 inline
102 VectorT
103 PCA< VectorT >::
104 center_of_gravity(const std::vector< VectorT >& _points ) {
105 // compute cog of Points
106 VectorT cog(0.0, 0.0, 0.0);
107
108 for (uint i = 0 ; i < _points.size() ; ++i ) {
109 cog = cog + _points[i];
110 }
111
112 cog = cog / ( typename VectorT::value_type )_points.size();
113
114 return cog;
115 }
116
117 template < typename VectorT >
118 bool
119 PCA< VectorT >::SymRightEigenproblem( Matrix &_mat_A, Matrix & _mat_VR,
120 std::vector< double > & _vec_EV )
121 {
122 char jobz, uplo;
123 long int n, lda;
124 long int lwork, info;
125 double size;
126
127 jobz = 'V';
128 uplo = 'U';
129 n = gmm::mat_nrows( _mat_A );
130 lda = n;
131
132 info = 0;
133 // compute optimal size of working array
134 lwork = -1;
135 dsyev_( &jobz, &uplo, &n, &_mat_A( 0, 0 ), &lda, &_vec_EV[ 0 ],
136 &size, &lwork, &info );
137
138 if( info != 0 )
139 return false;
140
141
142 if( (long int) size > m_workSize_double )
143 {
144 if( p_work_double)
145 delete [] p_work_double;
146
147 m_workSize_double = (long int) size;
148
149 p_work_double= new double[ m_workSize_double ];
150 }
151 lwork = m_workSize_double;
152
153 // compute eigenvalues / eigenvectors
154 dsyev_( &jobz, &uplo, &n, &_mat_A( 0, 0 ), &lda, &_vec_EV[ 0 ],
155 p_work_double, &lwork, &info );
156
157 if( info != 0 )
158 return false;
159
160 // copy eigenvectors to matrix
161 gmm::copy( _mat_A, _mat_VR );
162
163 return true;
164 }
165
166
167 template < typename VectorT >
168 void
169 PCA< VectorT >::pca(std::vector< VectorT >& _points , VectorT& _first , VectorT& _second , VectorT& _third)
170 {
171 // compute Mean of Points
172 const VectorT cog = center_of_gravity(_points);
173
174 //build covariance Matrix
175 Matrix points(3 , _points.size() );
176
177 for ( uint i = 0 ; i < _points.size() ; ++i)
178 {
179 points(0 , i ) = ( _points[i] - cog) [0];
180 points(1 , i ) = ( _points[i] - cog) [1];
181 points(2 , i ) = ( _points[i] - cog) [2];
182 }
183
184 Matrix cov(3,3);
185 gmm::mult(points,gmm::transposed(points),cov );
186
187 Matrix EV(3,3);
188 std::vector< double > ew(3);
189
190
191 if ( !SymRightEigenproblem( cov, EV ,ew ) )
192 std::cerr << "Error: Could not compute Eigenvectors for PCA" << std::endl;
193
194 int maximum,middle,minimum;
195
196 if ( (ew[0] > ew[1]) && (ew[0] > ew[2]) ) {
197 maximum = 0;
198 } else {
199 if ( (ew[1] > ew[0]) && (ew[1] > ew[2]) )
200 maximum = 1;
201 else
202 maximum = 2;
203 }
204
205 if ( (ew[0] < ew[1]) && (ew[0] < ew[2]) )
206 minimum = 0;
207 else if ( (ew[1] < ew[0]) && (ew[1] < ew[2]) )
208 minimum = 1;
209 else
210 minimum = 2;
211
212 if ( (minimum != 0 ) && ( maximum != 0 ) )
213 middle = 0;
214 else
215 if ( (minimum != 1 ) && ( maximum != 1 ) )
216 middle = 1;
217 else
218 middle = 2;
219
220 _first = VectorT( EV(0,maximum) , EV(1,maximum) , EV(2,maximum) );
221 _second = VectorT( EV(0,middle) , EV(1,middle) , EV(2,middle ) );
222
223 _first = _first.normalize();
224 _second = _second.normalize();
225 _third = _first % _second;
226
227 //remember the eigenvalues
228 lastEigenValues_.clear();
229 lastEigenValues_.push_back( ew[maximum] );
230 lastEigenValues_.push_back( ew[middle] );
231 lastEigenValues_.push_back( ew[minimum] );
232 }
233 //-----------------------------------------------------------------------------
234
235 template < typename VectorT >
236 inline
237 std::vector<double>&
238 PCA< VectorT >::getLastEigenValues()
239 {
240 return lastEigenValues_;
241 }
242
243
244 //=============================================================================
245
246 } //namespace Pca
247
248 //=============================================================================
Pca::PCA
Class for principal component Analysis.
Definition: PCA.hh:83
Pca::PCA::PCA
PCA()
Constructor.
Definition: PCA.hh:91
Pca::PCA::~PCA
~PCA()
Destructor.
Definition: PCAT_impl.hh:91
Pca::PCA::center_of_gravity
VectorT center_of_gravity(const std::vector< VectorT > &_points)
Definition: PCAT_impl.hh:104
VectorT
Pca
Namespace for principal Component Analysis.
Definition: PCA.hh:72
Pca::PCA::pca
void pca(std::vector< VectorT > &_points, VectorT &_first, VectorT &_second, VectorT &_third)
Definition: PCAT_impl.hh:169
PCA.hh
Classes for doing Principal component analysis.