67 typedef void (*U_fp)(...);
69 int dsyev_(
char *jobz,
char *uplo,
long int *n,
double *a,
70 long int *lda,
double *w,
double *work,
long int *lwork,
79 static double * p_work_double = 0;
82 static int m_workSize_double = 0;
90 template <
typename VectorT >
93 pca(_points,_first,_second,_third);
96 template <
typename VectorT >
99 if ( p_work_double != 0){
100 delete [] p_work_double;
102 m_workSize_double = 0;
106 template <
typename VectorT >
114 for (uint i = 0 ; i < _points.size() ; ++i ) {
115 cog = cog + _points[i];
118 cog = cog / (
typename VectorT::value_type )_points.size();
123 template <
typename VectorT >
126 std::vector< double > & _vec_EV )
130 long int lwork, info;
135 n = gmm::mat_nrows( _mat_A );
141 dsyev_( &jobz, &uplo, &n, &_mat_A( 0, 0 ), &lda, &_vec_EV[ 0 ],
142 &size, &lwork, &info );
148 if( (
long int) size > m_workSize_double )
151 delete [] p_work_double;
153 m_workSize_double = (
long int) size;
155 p_work_double=
new double[ m_workSize_double ];
157 lwork = m_workSize_double;
160 dsyev_( &jobz, &uplo, &n, &_mat_A( 0, 0 ), &lda, &_vec_EV[ 0 ],
161 p_work_double, &lwork, &info );
167 gmm::copy( _mat_A, _mat_VR );
173 template <
typename VectorT >
178 const VectorT cog = center_of_gravity(_points);
181 Matrix points(3 , _points.size() );
183 for ( uint i = 0 ; i < _points.size() ; ++i)
185 points(0 , i ) = ( _points[i] - cog) [0];
186 points(1 , i ) = ( _points[i] - cog) [1];
187 points(2 , i ) = ( _points[i] - cog) [2];
191 gmm::mult(points,gmm::transposed(points),cov );
194 std::vector< double > ew(3);
197 if ( !SymRightEigenproblem( cov, EV ,ew ) )
198 std::cerr <<
"Error: Could not compute Eigenvectors for PCA" << std::endl;
200 int maximum,middle,minimum;
202 if ( (ew[0] > ew[1]) && (ew[0] > ew[2]) ) {
205 if ( (ew[1] > ew[0]) && (ew[1] > ew[2]) )
211 if ( (ew[0] < ew[1]) && (ew[0] < ew[2]) )
213 else if ( (ew[1] < ew[0]) && (ew[1] < ew[2]) )
218 if ( (minimum != 0 ) && ( maximum != 0 ) )
221 if ( (minimum != 1 ) && ( maximum != 1 ) )
226 _first =
VectorT( EV(0,maximum) , EV(1,maximum) , EV(2,maximum) );
227 _second =
VectorT( EV(0,middle) , EV(1,middle) , EV(2,middle ) );
229 _first = _first.normalize();
230 _second = _second.normalize();
231 _third = _first % _second;
234 lastEigenValues_.clear();
235 lastEigenValues_.push_back( ew[maximum] );
236 lastEigenValues_.push_back( ew[middle] );
237 lastEigenValues_.push_back( ew[minimum] );
241 template <
typename VectorT >
246 return lastEigenValues_;
void pca(std::vector< VectorT > &_points, VectorT &_first, VectorT &_second, VectorT &_third)
Classes for doing Principal component analysis.
VectorT center_of_gravity(const std::vector< VectorT > &_points)
Namespace for principal Component Analysis.
Class for principal component Analysis.