itkVectorFieldPCA.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
=========================================================================*/

#ifndef itkVectorFieldPCA_h
#define itkVectorFieldPCA_h

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

namespace itk
{

template< typename TRealValueType = double >
class ITK_EXPORT GaussianDistanceKernel : public KernelFunctionBase<TRealValueType>
{
public:
  typedef GaussianDistanceKernel                  Self;
  typedef KernelFunctionBase<TRealValueType>      Superclass;
  typedef SmartPointer<Self>                      Pointer;
  typedef SmartPointer<const Self>                ConstPointer;

  itkTypeMacro(GaussianDistanceKernel, KernelFunction);

  itkNewMacro(Self);

  void SetKernelSigma( double s )
    { m_KernelSigma = s;
      m_OneOverMinusTwoSigmaSqr = -1.0 / (2.0 * s * s); }
  itkGetMacro( KernelSigma, double );

  inline double Evaluate (const double& u) const
    { return ( vcl_exp( u * m_OneOverMinusTwoSigmaSqr )); }

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

private:
  double                    m_KernelSigma;
  double                    m_OneOverMinusTwoSigmaSqr;

};

template <
    typename TVectorFieldElementType,
    typename TPCType,
    typename TPointSetPixelType = float,
    typename TPointSetCoordRepType = float,
    typename KernelFunctionType = KernelFunctionBase<TPointSetCoordRepType>,
    class    TPointSetType = PointSet<TPointSetPixelType, 3, DefaultStaticMeshTraits< TPointSetPixelType, 3, 3, TPointSetCoordRepType > >
  >
class ITK_EXPORT VectorFieldPCA : public Object
{
public:
  typedef VectorFieldPCA                  Self;
  typedef Object                          Superclass;
  typedef SmartPointer<Self>              Pointer;
  typedef SmartPointer<const Self>        ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(VectorFieldPCA, Object);

  typedef TPointSetType   InputPointSetType;

  typedef typename InputPointSetType::PointType     InputPointType;

  typedef typename InputPointSetType::PointsContainer    PointsContainer;
  typedef typename PointsContainer::Iterator             PointsContainerIterator;

  typedef typename InputPointSetType::Pointer  InputPointSetPointer;

  typedef typename InputPointSetType::ConstPointer InputPointSetConstPointer;

  itkStaticConstMacro(InputMeshDimension, unsigned int,
                      TPointSetType::PointDimension);

  typedef vnl_matrix< TVectorFieldElementType >               VectorFieldType;
  typedef VectorContainer< unsigned int, VectorFieldType >    VectorFieldSetType;

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

  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;
  typedef typename KernelFunctionType::Pointer   KernelFunctionPointer;

  itkSetMacro(PointSet, InputPointSetPointer);
  itkGetMacro(PointSet, InputPointSetPointer);

  itkSetMacro(VectorFieldSet, VectorFieldSetTypePointer);
  itkGetMacro(VectorFieldSet, VectorFieldSetTypePointer);

  itkSetMacro( ComponentCount, unsigned int );
  itkGetMacro( ComponentCount, unsigned int );

  itkSetMacro( KernelFunction, KernelFunctionPointer );

  void Compute(void);

  itkGetConstReferenceMacro(AveVectorField, MatrixType);
  itkGetConstReferenceMacro(PCAEigenValues, VectorType);
  itkGetConstObjectMacro(BasisVectors, BasisSetType);

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

  void                      KernelPCA(void);
  void                      computeMomentumSCP(void);

private:

  VectorFieldPCA(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  VectorType                m_PCAEigenValues;

  BasisSetTypePointer       m_BasisVectors;
  VectorFieldSetTypePointer m_VectorFieldSet;
  InputPointSetPointer      m_PointSet;
  KernelFunctionPointer     m_KernelFunction;

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

  MatrixType                m_V0;
  MatrixType                m_AveVectorField;
  MatrixType                m_K;

  bool                      m_PCACalculated;

};

} // end namespace itk

#include "itkVectorFieldPCA.hxx"

#endif