Opencv中的PCA实现

// EigenTest.cpp : Defines the entry point for the console application.
//
 
#include "stdafx.h"
 
void MyPCA( const Mat &_data, int dim, Mat &eigenvalues, Mat &eigenvectors );
 
void printMat( Mat _data )
{
    Mat data = cv::Mat_<double>(_data);
    for ( int i=0; i<data.rows; i++ )
    {
        for ( int j=0; j< data.cols; j++ )
        {
            cout << data.at<double>(i,j) << "    ";
        }
        cout << endl;
    }
}
 
int _tmain(int argc, _TCHAR* argv[])
{
    float A[ 60 ]=
    { 1.5 , 2.3 ,  1.5 , 2.3 ,  1.5 , 2.3 , 
    3.0 , 1.7 , 3.0 , 1.7 , 3.0 , 1.7 , 
    1.2 , 2.9 , 1.2 , 2.9 , 1.2 , 2.9 , 
    2.1 , 2.2 , 2.1 , 2.2 ,2.1 , 2.2 ,
    3.1 , 3.1 , 3.1 , 3.1 , 3.1 , 3.1 , 
    1.3 , 2.7 , 1.3 , 2.7 , 1.3 , 2.7 , 
    2.0 , 1.7 , 2.0 , 1.7 , 2.0 , 1.7 , 
    1.0 , 2.0 , 1.0 , 2.0 , 1.0 , 2.0 , 
    0.5 , 0.6 , 0.5 , 0.6 , 0.5 , 0.6 , 
    1.0 , 0.9 , 1.0 , 0.9 ,1.0 , 0.9 }; 
 
    Mat m = Mat::zeros( 10, 6, CV_32F );
 
    
    for ( int i=0; i<10; i++ )
    {
        for ( int j=0; j<6; j++ )
        {
            m.at<float>(i,j) = A[i*6+j];
        }
    }
 
    cout << "origin data: " << endl;
    printMat( m );
 
    /// OPENCV PCA
    PCA pca( m, noArray(), CV_PCA_DATA_AS_ROW );
 
    cout << "pca.eigenvalues:" << endl;
    printMat( pca.eigenvalues );
 
    cout << "pca.eigenvectors:" << endl;
    printMat( pca.eigenvectors );
 
 
    /// My PCA
    Mat eigenvalues;
    Mat eigenvectors;
 
    MyPCA( m, 3, eigenvalues, eigenvectors ); 
 
    cout << "eigenvalues:" << endl;
    printMat( eigenvalues );
 
    cout << "eigenvectors:" << endl;
    printMat( eigenvectors );
 
    return 0;
}
 
 
void MyPCA( const Mat &_data, int dim, Mat &eigenvalues, Mat &eigenvectors  )
{
    //#define ShowPcaMsg
 
    assert( dim>0 );
 
    Mat data =  cv::Mat_<double>(_data);
 
#ifdef ShowPcaMsg
    printMat( data );
#endif
 
    int N = data.rows;
    int D = data.cols;
 
    if ( dim>D )
    {
        dim = D;
    }
 
    /// mean
    Mat m = Mat::zeros( 1, D, data.type() );
 
    for ( int j=0; j<D; j++ )
    {
        for ( int i=0; i<N; i++ )
        {
            m.at<double>(0,j) += data.at<double>(i,j);
        }
    }
 
    m = m/N;
 
    /// 离散度矩阵计算
    Mat S =  Mat::zeros( N, D, data.type() );
    for ( int i=0; i<N; i++ )
    {
        for ( int j=0; j<D; j++ )
        {
            S.at<double>(i,j) = data.at<double>(i,j) - m.at<double>(0,j); ///< data-mean
        }
    }
 
    Mat C = S.t() * S /(N); ///< 除以N，或者N-1
 
 
#ifdef ShowPcaMsg
    printMat( C  );
#endif
 
    /// 求特征值特征向量
    eigen( C, eigenvalues, eigenvectors );
}