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:
  using Self = GaussianDistanceKernel;
  using Superclass = KernelFunctionBase<TRealValueType>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  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 override
    { return ( std::exp( u * m_OneOverMinusTwoSigmaSqr )); }

protected:
  GaussianDistanceKernel() {}
  ~GaussianDistanceKernel() override {}
  void PrintSelf(std::ostream& os, Indent indent) const override
    { 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:
  ITK_DISALLOW_COPY_AND_ASSIGN(VectorFieldPCA);

  using Self = VectorFieldPCA;
  using Superclass = Object;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkNewMacro(Self);

  itkTypeMacro(VectorFieldPCA, Object);

  using InputPointSetType = TPointSetType;

  using InputPointType = typename InputPointSetType::PointType;

  using PointsContainer = typename InputPointSetType::PointsContainer;
  using PointsContainerIterator = typename PointsContainer::Iterator;

  using InputPointSetPointer = typename InputPointSetType::Pointer;

  using InputPointSetConstPointer = typename InputPointSetType::ConstPointer;

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

  using VectorFieldType = vnl_matrix< TVectorFieldElementType >;
  using VectorFieldSetType = VectorContainer< unsigned int, VectorFieldType >;

  using VectorFieldSetTypePointer = typename VectorFieldSetType::Pointer;
  using VectorFieldSetTypeConstPointer = typename VectorFieldSetType::ConstPointer;

  using MatrixType = vnl_matrix< TPCType >;
  using VectorType = vnl_vector< TPCType >;

  using BasisSetType = VectorContainer< unsigned int, MatrixType >;
  using ResSetType = VectorContainer< unsigned int, VectorType >;

  using BasisSetTypePointer = typename BasisSetType::Pointer;
  using KernelFunctionPointer = typename KernelFunctionType::Pointer;

  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();
  ~VectorFieldPCA() override {};
  void PrintSelf(std::ostream& os, Indent indent) const override;

  void KernelPCA();

  void ComputeMomentumSCP();

private:
  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

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkVectorFieldPCA.hxx"
#endif

#endif