/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: $
  Language:  C++
  Date:      $Date: $
  Version:   $Revision: $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __itkMeshVectorFieldPCA_h
#define __itkMeshVectorFieldPCA_h

#include "itkObject.h"
#include "vnl/vnl_vector.h"
#include "vnl/vnl_matrix.h"

namespace itk
{

/** \class MeshVectorFieldPCA
 * \brief Decomposes a mesh into directions along basis components.
 * 
 * This calculator produces a set of basis functions composed of the
 * principal components of the vector values at the mesh vertices.
 *  
 * This class is templated over the input mesh type.
 *
 * \warning This method assumes that the input mesh consists of vector pixel
 * types.
 *
 * \author Michael Bowers
 *
 * \ingroup Operators
 */
template <
    class TInputMesh,
    typename TVectorFieldElementType,
    typename TPCType
  >
class ITK_EXPORT MeshVectorFieldPCA : public Object 
{
public:
  /** Standard class typedefs. */
  typedef MeshVectorFieldPCA              Self;
  typedef Object                          Superclass;
  typedef SmartPointer              Pointer;
  typedef SmartPointer        ConstPointer;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(MeshVectorFieldPCA, Object);

  /** Type definitions for the mesh. */
  typedef TInputMesh   InputMeshType;

  /** Definitions for points on the mesh. */
  typedef typename InputMeshType::PointType     InputPointType;
  
  /** Pointer types for the mesh. */
  typedef typename TInputMesh::Pointer  InputMeshPointer;
  
  /** Const Pointer type for the mesh. */
  typedef typename TInputMesh::ConstPointer InputMeshConstPointer;

  /** Input Mesh dimension */
  itkStaticConstMacro(InputMeshDimension, unsigned int,
                      TInputMesh::PointDimension); 

  /** type for the vector fields. */
  typedef vnl_matrix< TVectorFieldElementType >               VectorFieldType;
  typedef VectorContainer< unsigned int, VectorFieldType >    VectorFieldSetType;

  typedef typename VectorFieldSetType::Pointer       VectorFieldSetTypePointer;
  typedef typename VectorFieldSetType::ConstPointer  VectorFieldSetTypeConstPointer;

  /** types for the output. */
  typedef vnl_matrix< TPCType > MatrixType;
  typedef vnl_vector< TPCType > VectorType;

  typedef VectorContainer< unsigned int, MatrixType > BasisSetType;
  typedef VectorContainer< unsigned int, VectorType > ResSetType;

  typedef typename BasisSetType::Pointer       BasisSetTypePointer;

  /** Set and get the input mesh. */
//  itkSetConstObjectMacro(Mesh, InputMeshPointer);
//  itkGetConstObjectMacro(Mesh, InputMeshPointer);
  itkSetMacro(Mesh, InputMeshPointer);
  itkGetMacro(Mesh, InputMeshPointer);
  
  /** Set and get the vector fields for the analysis. */
  itkSetMacro(VectorFieldSet, VectorFieldSetTypePointer);
  itkGetMacro(VectorFieldSet, VectorFieldSetTypePointer);
  
  /** Set and get the PCA count. */
  itkSetMacro( ComponentCount, unsigned int );
  itkGetMacro( ComponentCount, unsigned int );

  /** Set and get the Kernel sigma. */
  itkSetMacro( KernelSigma, double );
  itkGetMacro( KernelSigma, double );
  
  /** Compute the PCA decomposition of the input mesh. */
  void Compute(void);

  /** Return the res. */
  itkGetConstReferenceMacro(Projection, MatrixType);
  itkGetConstReferenceMacro(AveVectorField, MatrixType);
  itkGetConstReferenceMacro(PCAEigenValues, VectorType);
  itkGetConstObjectMacro(BasisVectors, BasisSetType);
/*
  virtual const BasisSetType & GetBasisVectors() const
  {
    return this->m_BasisVectors;
  }
  itkGetConstReferenceMacro(BasisVectors, BasisSetType);
*/

protected:
  MeshVectorFieldPCA();
  virtual ~MeshVectorFieldPCA() {};
  void PrintSelf(std::ostream& os, Indent indent) const;

  void                      KernelPCA(void);
  void                      computeMomentumSCP(void);
  
private:
    
  MeshVectorFieldPCA(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
  
  VectorType                m_PCAEigenValues;

  BasisSetTypePointer       m_BasisVectors;
  VectorFieldSetTypePointer m_VectorFieldSet;
  InputMeshPointer          m_Mesh;
  double                    m_KernelSigma;

  // problem dimensions
  unsigned int              m_ComponentCount;
  unsigned int              m_SetSize;
  unsigned int              m_VertexCount;
  unsigned int              m_PointDim;

  MatrixType                m_Projection;
  MatrixType                m_V0;
  MatrixType                m_AveVectorField;
  MatrixType                m_K;

  bool                      m_PCACalculated;

};

} // end namespace itk


#ifndef ITK_MANUAL_INSTANTIATION
#include "itkMeshVectorFieldPCA.txx"
#endif

#endif