PCA.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <MathAlgo/PCA.hh>
#include <MathAlgo/SVD.hh>
#include <Base/Common/BIASpragma.hh>

using namespace BIAS;
using namespace std; 


void PCA::ComputeReductionMatrix(const vector<vector<Vector<PCAType> > >& vec, 
                                 BIAS::Matrix<PCAType> &matrix, 
                                 int reductionSize){
  int counter = 0;
  for (unsigned int i=0; i<vec.size(); i++) {
    counter += vec[i].size();
  }
  vector<Vector<PCAType> > newvec;
  newvec.reserve(counter);
  for (unsigned int i=0; i<vec.size(); i++) {
    const vector<Vector<PCAType> > &partvec = vec[i];
    for (unsigned int j=0; j<partvec.size(); j++) {
      newvec.push_back(partvec[j]);
    }
  }
  ComputeReductionMatrix(newvec, matrix, reductionSize);

}

void PCA::ComputeReductionMatrix(const vector<BIAS::Vector<PCAType> > &vec, 
                                 BIAS::Matrix<PCAType> &matrix, 
                                 int reductionSize){
  if (reductionSize==-1) reductionSize = reductionSize_;
  BIASASSERT(reductionSize>0);

  const int dimensions=vec[0].size();
     
  
  //covariance-matrix and matrix U,V,S for SVD
  Matrix<PCAType> cov(dimensions, dimensions, 0.0f);
      
  //compute mean
  ComputeMean(vec, mean_);
  ComputeScatter(vec, mean_, cov);
  ComputeReductionMatrix(cov, matrix, false, reductionSize);
}

 

void PCA::ComputeReductionMatrix(BIAS::Matrix<PCAType> &cov,
                                 BIAS::Matrix<PCAType> &matrix, 
                                 bool normalize,
                                 int reductionSize){
  if (reductionSize==-1) reductionSize = reductionSize_;
  BIASASSERT(reductionSize>0);
  BIASASSERT(cov.num_rows() == cov.num_cols());

  Matrix<double> covd(cov.num_rows(), cov.num_cols(), 0.0);   
  for (int i=0; i<cov.num_cols(); i++){
    for (int j=0; j<=i; j++){
      // make symmetric (scale by 2 does not matter)
      covd[i][j] = covd[j][i] = double( cov[j][i] ) + double( cov[i][j] ) ;
    }
  }

  Matrix<double> U,V;
  //make SVD
  SVD svd(covd);
  S_ = svd.GetS();
  U  = svd.GetVT();

  // if rows are normalize we should use mahalanobis distance in transformed
  // space instead of norm2, so do not normalize here ...
  if (normalize) U.NormalizeRows();

  matrix.newsize(cov.num_rows(), reductionSize);
  for (int i=0;i<matrix.num_cols();i++){
    for (int j=0;j<matrix.num_rows();j++){
      matrix[j][i]= PCAType( U[i][j] );//<1.0e-13?0.0f:U[j][i];
    }
  }
 
}

void PCA::ComputeScatter(const vector<BIAS::Vector<PCAType> > &vec,
                         const BIAS::Vector<PCAType> &mean,
                         BIAS::Matrix<PCAType> &cov) {
  // compute covariance matrix
  for (int k=0; k<(int)vec.size(); k++){
    const BIAS::Vector<PCAType>& veck = vec[k];
    for (int i=0; i<veck.size(); i++){
      const PCAType &vecki = veck[i];
      for (int j=0; j<=i; j++){
        cov[i][j] += (vecki-mean[i]) * (veck[j]-mean[j]);
      }
    }
  }
  for (int i=0; i<cov.num_rows(); i++){
    for (int j=i+1; j<cov.num_cols(); j++){
      cov[i][j] = cov[j][i];
    }
  }
}

void PCA::ComputeMean(const vector<BIAS::Vector<PCAType> > &vec, 
                      BIAS::Vector<PCAType> &mean){
  BIASASSERT(!vec.empty());
  // sum up
  mean.newsize(vec[0].size());
  mean.SetZero();
  for (int k=0; k<(int)vec.size(); k++){
    const BIAS::Vector<PCAType>& veck = vec[k];
    for (int i=0; i<veck.size(); i++){
      mean[i] += veck[i];
    }
  }
  // divide
  const PCAType vecsize = 1.0f / PCAType(vec.size());
  for (int i=0; i<mean.size(); i++){
    mean[i] *= vecsize;
  }
}

void PCA::SetReductionSize(int size){
  reductionSize_ = size;
}

void PCA::GetMean(BIAS::Vector<PCAType> &mean){
  mean = mean_;
}