itkParticleShapeStatistics.h

/*=========================================================================
   Program:   ShapeWorks: Particle-based Shape Correspondence & Visualization
   Module:    $RCSfile: itkParticleShapeStatistics.h,v $
   Date:      $Date: 2011/03/24 01:17:41 $
   Version:   $Revision: 1.3 $
   Author:    $Author: wmartin $

   Copyright (c) 2009 Scientific Computing and Imaging Institute.
   See ShapeWorksLicense.txt for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.
   =========================================================================*/
#ifndef _ParticleShapeStatistics_h
#define _ParticleShapeStatistics_h

#include <iostream>
#include <vector>
#include "vnl/vnl_vector.h"
#include "vnl/algo/vnl_symmetric_eigensystem.h"
#include "vnl/vnl_matrix.h"
#include "itkParticlePositionReader.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/algo/vnl_matrix_inverse.h"
#include <iostream>
#include <fstream>
#include <string>
#include <cstdio>
#include "itkParticlePositionWriter.h"

template <unsigned int VDimension = 3>
class ITK_EXPORT ParticleShapeStatistics
{
public:
  ParticleShapeStatistics() {}
  ~ParticleShapeStatistics() {}

  itkStaticConstMacro( Dimension, unsigned int, VDimension );

  int ImportPoints( std::vector<vnl_vector<double> > points );

  int ReadPointFiles( const char* fname );

  int ReloadPointFiles();

  int WriteCSVFile( const char* s );
  int WriteCSVFile( const std::string &s )
  { return this->WriteCSVFile( s.c_str() ); }
  int WriteCSVFile2( const char* s );
  int WriteCSVFile2( const std::string &s )
  { return this->WriteCSVFile2( s.c_str() ); }

  int ComputeModes();

  int PrincipalComponentProjections();

  int FisherLinearDiscriminant( unsigned int );

  inline int SampleSize() const
  { return m_numSamples; }
  inline int Group1SampleSize() const
  { return m_numSamples1; }
  inline int Group2SampleSize() const
  { return m_numSamples2; }

  inline int NumberOfDimensions() const
  { return m_numDimensions; }

  inline int GroupID( unsigned int i ) const
  { return m_groupIDs[i]; }
  const std::vector<int> &GroupID() const
  { return m_groupIDs; }

  const vnl_matrix<double> &Eigenvectors() const
  { return m_eigenvectors; }
  const vnl_vector<double> &Eigenvalues() const
  { return m_eigenvalues; }

  const vnl_vector<double> &Mean() const
  { return m_mean; }
  const vnl_vector<double> &Group1Mean() const
  { return m_mean1; }
  const vnl_vector<double> &Group2Mean() const
  { return m_mean2; }

  // Returns group2 - group1 mean
  const vnl_vector<double> &NormalizedGroupDifference() const
  { return m_groupdiffnorm; }
  const vnl_vector<double> &GroupDifference() const
  { return m_groupdiff; }

  int ComputeMedianShape( const int );

  double L1Norm( unsigned int a, unsigned int b );

  const vnl_matrix<double> &PCALoadings() const
  { return m_principals; }

  const vnl_vector<double> &FishersLDA() const
  { return m_fishersLD; }

  const vnl_matrix<double> &ShapeMatrix() const
  {return m_shapes; }

  const vnl_vector<double> &Shape( unsigned int i ) const
  { return m_shapes.get_column( i ); }

  int SimpleLinearRegression( const std::vector<double> &y,
                              const std::vector<double> &x,
                              double &a, double &b ) const;

protected:
  unsigned int m_numSamples1;
  unsigned int m_numSamples2;
  unsigned int m_numSamples;
  unsigned int m_domainsPerShape;
  unsigned int m_numDimensions;
  std::vector<int> m_groupIDs;

  vnl_matrix<double> m_pooled_covariance;
  vnl_matrix<double> m_eigenvectors;
  vnl_vector<double> m_eigenvalues;
  vnl_vector<double> m_mean;
  vnl_vector<double> m_mean1;
  vnl_vector<double> m_mean2;
  vnl_matrix<double> m_pointsMinusMean;
  vnl_matrix<double> m_shapes;
  vnl_matrix<double> m_projectedPMM1;
  vnl_matrix<double> m_projectedPMM2;
  vnl_vector<double> m_projectedMean1;
  vnl_vector<double> m_projectedMean2;
  std::vector<double> m_fishersProjection;
  std::vector<double> m_percentVarByMode;
  vnl_vector<double> m_fishersLD;
  int m_top95;
  vnl_matrix<double> m_principals;

  vnl_vector<double> m_groupdiff;
  vnl_vector<double> m_groupdiffnorm;

  // used to keep the points' files that needs to be reloaded when new updates come in.
  std::vector< std::string > m_pointsfiles;
};

#if ITK_TEMPLATE_EXPLICIT
# include "Templates/itkParticleShapeStatistics+-.h"
#endif

#if ITK_TEMPLATE_TXX
# include "itkParticleShapeStatistics.txx"
#endif

#endif // ifndef _ParticleShapeStatistics_h