Developer Documentation
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 >
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
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>&
245 {
246  return lastEigenValues_;
247 }
248 
249 
250 //=============================================================================
251 
252 } //namespace Pca
253 
254 //=============================================================================
VectorT center_of_gravity(const std::vector< VectorT > &_points)
Definition: PCA.cc:110
Class for principal component Analysis.
Definition: PCA.hh:89
Namespace for principal Component Analysis.
Definition: PCA.cc:86
~PCA()
Destructor.
Definition: PCA.cc:97
void pca(std::vector< VectorT > &_points, VectorT &_first, VectorT &_second, VectorT &_third)
Definition: PCA.cc:175
PCA()
Constructor.
Definition: PCA.hh:97
Classes for doing Principal component analysis.