itkKernelTransform2.h

Go to the documentation of this file.

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

Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkKernelTransform2.h,v $
Language:  C++
Date:      $Date: 2006-11-28 14:22:18 $
Version:   $Revision: 1.1 $

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 __itkKernelTransform2_h
#define __itkKernelTransform2_h

#include "itkAdvancedTransform.h"
#include "itkPoint.h"
#include "itkVector.h"
#include "itkMatrix.h"
#include "itkPointSet.h"
#include <deque>
#include <math.h>
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/algo/vnl_svd.h"
#include "vnl/vnl_sample.h"

namespace itk
{

template <class TScalarType, // probably only float and double make sense here
  unsigned int NDimensions>   // Number of dimensions
class KernelTransform2
  : public AdvancedTransform<TScalarType, NDimensions,NDimensions>
{
public:

  typedef KernelTransform2                  Self;
  typedef AdvancedTransform<
    TScalarType, NDimensions, NDimensions > Superclass;
  typedef SmartPointer<Self>                Pointer;
  typedef SmartPointer<const Self>          ConstPointer;

  itkTypeMacro( KernelTransform2, AdvancedTransform );

  itkNewMacro( Self );

  itkStaticConstMacro(SpaceDimension, unsigned int, NDimensions);

  typedef typename Superclass::ScalarType       ScalarType;
  typedef typename Superclass::ParametersType   ParametersType;
  typedef typename Superclass::JacobianType     JacobianType;
  typedef typename Superclass::InputPointType   InputPointType;
  typedef typename Superclass::OutputPointType  OutputPointType;
  typedef typename Superclass::InputVectorType  InputVectorType;
  typedef typename Superclass::OutputVectorType OutputVectorType;

  typedef typename Superclass
    ::NonZeroJacobianIndicesType                    NonZeroJacobianIndicesType;
  typedef typename Superclass::SpatialJacobianType  SpatialJacobianType;
  typedef typename Superclass
    ::JacobianOfSpatialJacobianType                 JacobianOfSpatialJacobianType;
  typedef typename Superclass::SpatialHessianType   SpatialHessianType;
  typedef typename Superclass
    ::JacobianOfSpatialHessianType                  JacobianOfSpatialHessianType;
  typedef typename Superclass::InternalMatrixType   InternalMatrixType;

  typedef DefaultStaticMeshTraits< TScalarType,
    NDimensions, NDimensions, TScalarType, TScalarType>       PointSetTraitsType;
  typedef PointSet<InputPointType, NDimensions,
    PointSetTraitsType>                                       PointSetType;
  typedef typename PointSetType::Pointer                      PointSetPointer;
  typedef typename PointSetType::PointsContainer              PointsContainer;
  typedef typename PointSetType::PointsContainerIterator      PointsIterator;
  typedef typename PointSetType::PointsContainerConstIterator PointsConstIterator;

  typedef VectorContainer<unsigned long,InputVectorType>      VectorSetType;
  typedef typename VectorSetType::Pointer                     VectorSetPointer;

  typedef vnl_matrix_fixed<TScalarType, NDimensions, NDimensions> IMatrixType;

  virtual bool GetHasNonZeroSpatialHessian( void ) const
  {
    return true;
  }
  virtual bool HasNonZeroJacobianOfSpatialHessian( void ) const
  {
    return true;
  }

  virtual unsigned int GetNumberOfParameters(void) const
  {
    return ( this->m_SourceLandmarks->GetNumberOfPoints() * SpaceDimension );
  }

  itkGetObjectMacro( SourceLandmarks, PointSetType );

  virtual void SetSourceLandmarks( PointSetType * );

  itkGetObjectMacro( TargetLandmarks, PointSetType );

  virtual void SetTargetLandmarks( PointSetType * );

  itkGetObjectMacro( Displacements, VectorSetType );

  void ComputeWMatrix( void );

  void ComputeLInverse( void );

  virtual OutputPointType TransformPoint( const InputPointType & thisPoint ) const;

  virtual const JacobianType & GetJacobian( const InputPointType & point ) const;

  virtual void GetJacobian(
    const InputPointType &,
    JacobianType &,
    NonZeroJacobianIndicesType & ) const;

  virtual void SetIdentity( void );

  virtual void SetParameters( const ParametersType & );

  virtual void SetFixedParameters( const ParametersType & );

  virtual void UpdateParameters( void );

  virtual const ParametersType & GetParameters( void ) const;

  virtual const ParametersType & GetFixedParameters( void ) const;

  virtual void SetStiffness( double stiffness )
  {
    this->m_Stiffness = stiffness > 0 ? stiffness : 0.0;
    this->m_LMatrixComputed = false;
    this->m_LInverseComputed = false;
    this->m_WMatrixComputed = false;
  }
  itkGetMacro( Stiffness, double );

  virtual void SetAlpha( TScalarType Alpha ) {};
  virtual TScalarType GetAlpha( void ) const { return -1.0; }

  itkSetMacro( PoissonRatio, TScalarType );
  virtual const TScalarType GetPoissonRatio( void ) const
  {
    return this->m_PoissonRatio;
  };

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

public:
  typedef vnl_matrix_fixed<TScalarType, NDimensions, NDimensions> GMatrixType;

  typedef vnl_matrix<TScalarType> LMatrixType;

  typedef vnl_matrix<TScalarType> KMatrixType;

  typedef vnl_matrix<TScalarType> PMatrixType;

  typedef vnl_matrix<TScalarType> YMatrixType;

  typedef vnl_matrix<TScalarType> WMatrixType;

  typedef vnl_matrix<TScalarType> DMatrixType;

  typedef vnl_matrix_fixed<TScalarType,NDimensions,NDimensions> AMatrixType;

  typedef vnl_vector_fixed<TScalarType,NDimensions> BMatrixType;

  typedef vnl_matrix_fixed<TScalarType, 1, NDimensions> RowMatrixType;

  typedef vnl_matrix_fixed<TScalarType, NDimensions, 1> ColumnMatrixType;

  PointSetPointer m_SourceLandmarks;

  PointSetPointer m_TargetLandmarks;

protected:

  virtual void ComputeG( const InputVectorType & landmarkVector,
    GMatrixType & GMatrix ) const;

  virtual void ComputeReflexiveG( PointsIterator, GMatrixType & GMatrix ) const;

  virtual void ComputeDeformationContribution(
    const InputPointType & inputPoint,
    OutputPointType & result ) const;

  void ComputeK( void );

  void ComputeL( void );

  void ComputeP( void );

  void ComputeY( void );

  void ComputeD( void );

  void ReorganizeW( void );

  double m_Stiffness;

  VectorSetPointer m_Displacements;

  LMatrixType m_LMatrix;

  LMatrixType m_LMatrixInverse;

  KMatrixType m_KMatrix;

  PMatrixType m_PMatrix;

  YMatrixType m_YMatrix;

  WMatrixType m_WMatrix;

  DMatrixType m_DMatrix;

  AMatrixType m_AMatrix;

  BMatrixType m_BVector;

  //GMatrixType m_GMatrix;

  bool m_WMatrixComputed;
  bool m_LMatrixComputed;
  bool m_LInverseComputed;

  IMatrixType m_I;

  NonZeroJacobianIndicesType m_NonZeroJacobianIndices;

  mutable NonZeroJacobianIndicesType m_NonZeroJacobianIndicesTemp;

private:
  KernelTransform2(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  TScalarType m_PoissonRatio;

};

} // end namespace itk

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

#endif // __itkKernelTransform2_h