itkSwitchableAffineTransform.h

Go to the documentation of this file.

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

  NifTK: A software platform for medical image computing.

  Copyright (c) University College London (UCL). All rights reserved.

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.

  See LICENSE.txt in the top level directory for details.

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

#ifndef itkSwitchableAffineTransform_h
#define itkSwitchableAffineTransform_h

#include <iostream>
#include <itkMatrix.h>
#include "itkUCLBaseTransform.h"
#include <itkExceptionObject.h>
#include <itkMacro.h>
#include <itkAffineTransform.h>
#include "itkMatrixLinearCombinationFunctions.h"

namespace itk
{

template <
  class TScalarType=double,         // Data type for scalars 
  unsigned int NInputDimensions=3,  // Number of dimensions in the input space
  unsigned int NOutputDimensions=3> // Number of dimensions in the output space
class ITK_EXPORT SwitchableAffineTransform 
  : public UCLBaseTransform< TScalarType, NInputDimensions, NOutputDimensions >
{
public:
  typedef SwitchableAffineTransform<TScalarType, 
                     NInputDimensions, 
                     NInputDimensions>          Self;
  typedef UCLBaseTransform< TScalarType,
                     NInputDimensions,
                     NOutputDimensions >        Superclass;
  typedef SmartPointer<Self>                    Pointer;
  typedef SmartPointer<const Self>              ConstPointer;
  
  itkTypeMacro( SwitchableAffineTransform, UCLBaseTransform );
  
  itkStaticConstMacro(InputSpaceDimension, unsigned int, NInputDimensions);
  itkStaticConstMacro(OutputSpaceDimension, unsigned int, NOutputDimensions);
  itkStaticConstMacro(ParametersDimension, unsigned int, (NInputDimensions+1 * NOutputDimensions+1)-1);
  
  typedef typename Superclass::ParametersType                  ParametersType;

  typedef typename Superclass::JacobianType                    JacobianType;

  typedef typename Superclass::ScalarType                      ScalarType;

  typedef Vector<TScalarType,
                 itkGetStaticConstMacro(InputSpaceDimension)>  InputVectorType;
  typedef Vector<TScalarType,
                 itkGetStaticConstMacro(OutputSpaceDimension)> OutputVectorType;
  
  typedef CovariantVector<TScalarType,
                          itkGetStaticConstMacro(InputSpaceDimension)>  
                                                    InputCovariantVectorType;
  typedef CovariantVector<TScalarType,
                          itkGetStaticConstMacro(OutputSpaceDimension)>  
                                                    OutputCovariantVectorType;
  
  typedef vnl_vector_fixed<TScalarType,
                           itkGetStaticConstMacro(InputSpaceDimension)> 
                                                    InputVnlVectorType;
  typedef vnl_vector_fixed<TScalarType,
                           itkGetStaticConstMacro(OutputSpaceDimension)> 
                                                    OutputVnlVectorType;
  
  typedef Point<TScalarType,
                itkGetStaticConstMacro(InputSpaceDimension)>   
                                                    InputPointType;
  typedef Point<TScalarType,
                itkGetStaticConstMacro(OutputSpaceDimension)>  
                                                    OutputPointType;
  
  typedef Matrix<TScalarType, itkGetStaticConstMacro(OutputSpaceDimension),
                 itkGetStaticConstMacro(InputSpaceDimension)>  
                                                    MatrixType;

  typedef Matrix<TScalarType, itkGetStaticConstMacro(InputSpaceDimension),
                 itkGetStaticConstMacro(OutputSpaceDimension)> 
                                                    InverseMatrixType;
                 
  typedef Matrix<TScalarType, NOutputDimensions+1, NInputDimensions+1>          FullAffineMatrixType; 
  typedef AffineTransform<TScalarType, NInputDimensions>                        FullAffineTransformType;
  typedef typename FullAffineTransformType::Pointer                             FullAffineTransformPointer;
  
  typedef InputPointType                            CenterType;

  typedef Array<double>                             TranslationType;

  typedef Array<double>                             RotationType;
  
  typedef Array<double>                             ScaleType;
                                                      
  typedef Array<double>                             SkewType;                     

  typedef Array<double>                             RelativeParameterWeightingType;

  void SetIdentity( void );

  const MatrixType & GetMatrix() const { return m_Matrix; }

  const OutputVectorType & GetOffset(void) const { return m_Offset; }
  
  FullAffineTransformType* GetFullAffineTransform()
    {
      ParametersType parameters;
      parameters.SetSize((NInputDimensions+1)*(NInputDimensions+1));
      
      FullAffineMatrixType matrix = GetFullAffineMatrix();
      unsigned int counter = 0;
      
      for (unsigned int i = 0; i < NOutputDimensions+1; i++)
      {
        for (unsigned int j = 0; j < NInputDimensions+1; j++)
        {
          parameters.SetElement(counter, matrix(i,j));
          counter++;
          //std::cout << "Matt: matrix[" << i << "," << j << "]" << "=" << matrix(i,j) << ", parameter=" << parameters.GetElement(counter) << std::endl;
        }
      }
      
      m_AffineMatrixTransform->SetParameters(parameters);
      return m_AffineMatrixTransform.GetPointer();
    }
  
  FullAffineMatrixType GetFullAffineMatrix() const
  {
    FullAffineMatrixType fullAffineMatrix; 
    
    fullAffineMatrix.SetIdentity(); 
    for (unsigned int i = 0; i < NOutputDimensions; i++)
    {
      for (unsigned int j = 0; j < NInputDimensions; j++)
      {
        fullAffineMatrix(i,j) = this->m_Matrix(i,j); 
      }
    }
    for (unsigned int i = 0; i < NOutputDimensions; i++)
    {
      fullAffineMatrix(i,NInputDimensions) = this->m_Offset[i]; 
    }
    return fullAffineMatrix; 
  }
  
  void SetFullAffineMatrix(const FullAffineMatrixType& fullAffineMatrix)
  {
    this->m_Matrix.SetIdentity(); 
    for (unsigned int i = 0; i < NOutputDimensions; i++)
    {
      for (unsigned int j = 0; j < NInputDimensions; j++)
      {
        this->m_Matrix(i,j) = fullAffineMatrix(i,j); 
      }
    }
    for (unsigned int i = 0; i < NOutputDimensions; i++)
    {
      this->m_Offset[i] = fullAffineMatrix(i,NInputDimensions); 
    }
    this->m_MatrixMTime.Modified(); 
  }
  
  void HalfTransformationMatrix()
  {
    FullAffineMatrixType fullAffineMatrix = GetFullAffineMatrix(); 
    
    FullAffineMatrixType sqaureRoot(MatrixLinearCombinationFunctions<typename FullAffineMatrixType::InternalMatrixType>::ComputeMatrixSquareRoot(fullAffineMatrix.GetVnlMatrix(), 0.001)); 
    
    SetFullAffineMatrix(sqaureRoot); 
  }

  void SetCenter(const InputPointType & center)
      { m_Center = center;
        this->ComputeMatrixAndOffset();
        this->Modified(); 
        return; 
      }

  const InputPointType & GetCenter() const
      { return m_Center; }

  void SetTranslation(const TranslationType & translation)
      { m_Translation = translation; 
        this->ComputeMatrixAndOffset();
        this->Modified(); 
        return; }

  const TranslationType & GetTranslation(void) const
      { return m_Translation; }

  void SetRotation(const RotationType & rotation)
      { m_Rotation = rotation;
        this->ComputeMatrixAndOffset();
        this->Modified(); 
        return; }

  const RotationType & GetRotation(void) const
      { return m_Rotation; }

  void SetScale(const ScaleType & scale)
      { m_Scale = scale;
        this->ComputeMatrixAndOffset();
        this->Modified(); 
        return; }

  const ScaleType & GetScale(void) const
      { return m_Scale; }

  void SetSkew(const SkewType & skew)
      { m_Skew = skew;
        this->ComputeMatrixAndOffset();
        this->Modified(); 
        return; }

  const SkewType & GetSkew(void) const
      { return m_Skew; }

  void SetTranslationRelativeWeighting(const RelativeParameterWeightingType &weighting)
  {
    m_TranslationRelativeWeighting = weighting; 
  }
  void SetTranslationRelativeWeighting(const double &weighting)
  {
    m_TranslationRelativeWeighting.Fill( weighting ); 
  }

  void SetRotationRelativeWeighting(const RelativeParameterWeightingType &weighting)
  {
    m_RotationRelativeWeighting = weighting; 
  }
  void SetRotationRelativeWeighting(const double &weighting)
  {
    m_RotationRelativeWeighting.Fill( weighting ); 
  }
  void SetScaleRelativeWeighting(const RelativeParameterWeightingType &weighting)
  {
    m_ScaleRelativeWeighting = weighting; 
  }
  void SetScaleRelativeWeighting(const double &weighting)
  {
    m_ScaleRelativeWeighting.Fill( weighting ); 
  }
  void SetSkewRelativeWeighting(const RelativeParameterWeightingType &weighting)
  {
    m_SkewRelativeWeighting = weighting; 
  }
  void SetSkewRelativeWeighting(const double &weighting)
  {
    m_SkewRelativeWeighting.Fill( weighting ); 
  }



  void SetParameters( const ParametersType & parameters );

  const ParametersType& GetParameters(void) const;
  
  void SetFixedParameters(const ParametersType& parameters)
  {
    if (parameters.GetSize() != 1+NInputDimensions)
      itkExceptionMacro("SwitchableAffineTransform: number of expected fixed parameters does not match.");
    SetNumberOfDOF(static_cast<int>(parameters.GetElement(0))); 
    for (unsigned int d = 0; d < NInputDimensions; d++)
    {
      this->m_Center[d] = parameters.GetElement(1+d); 
    }
  }
  
  const ParametersType& GetFixedParameters() const
  {
    // Store the dof and centre of the transformation. 
    this->m_FixedParameters.SetSize(1+NInputDimensions); 
    
    this->m_FixedParameters.SetElement(0, this->GetNumberOfDOF()); 
    for (unsigned int d = 0; d < NInputDimensions; d++)
    {
      this->m_FixedParameters.SetElement(d+1, this->m_Center[d]); 
    }
    return this->m_FixedParameters; 
  }

  const RelativeParameterWeightingType& GetRelativeParameterWeightingFactors();

  OutputPointType     TransformPoint(const InputPointType & point) const;
  OutputVectorType    TransformVector(const InputVectorType & vector) const;
  OutputVnlVectorType TransformVector(const InputVnlVectorType & vector) const;
  OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &vector) const;
  
  itkSetMacro( OptimiseRotation, bool );
  itkGetConstMacro( OptimiseRotation, bool );
  itkBooleanMacro( OptimiseRotation );

  itkSetMacro( OptimiseTranslation, bool );
  itkGetConstMacro( OptimiseTranslation, bool );
  itkBooleanMacro( OptimiseTranslation );

  itkSetMacro( OptimiseScale, bool );
  itkGetConstMacro( OptimiseScale, bool );
  itkBooleanMacro( OptimiseScale );

  itkSetMacro( OptimiseSkew, bool );
  itkGetConstMacro( OptimiseSkew, bool );
  itkBooleanMacro( OptimiseSkew );

  virtual itk::TransformBase::NumberOfParametersType GetNumberOfParameters() const 
    { 
      return this->GetNumberOfDOF();
    }
  
  unsigned int GetNumberOfDOF() const;
  
  void SetNumberOfDOF(int number);

  void SetRigid()
    {
      OptimiseRotationOn();
      OptimiseTranslationOn();
      OptimiseScaleOff();
      OptimiseSkewOff();
      this->Modified();
    } 
  
  void SetRigidPlusScale()
    {
      OptimiseRotationOn();
      OptimiseTranslationOn();
      OptimiseScaleOn();
      OptimiseSkewOff();
      this->Modified();
    }   
  
  void SetFullAffine()
    {
      OptimiseRotationOn();
      OptimiseTranslationOn();
      OptimiseScaleOn();
      OptimiseSkewOn();
      this->Modified();
    } 

  void SetJustScale()
    {
      OptimiseRotationOff();
      OptimiseTranslationOff();
      OptimiseScaleOn();
      OptimiseSkewOff();
      this->Modified();
    }

  void SetJustRotation()
    {
      OptimiseRotationOn();
      OptimiseTranslationOff();
      OptimiseScaleOff();
      OptimiseSkewOff();
      this->Modified();
    }

  void SetJustTranslation()
    {
      OptimiseRotationOff();
      OptimiseTranslationOn();
      OptimiseScaleOff();
      OptimiseSkewOff();
      this->Modified();
    }

  virtual const JacobianType GetJacobian(const InputPointType & point ) const = 0;
  
  bool GetInverse(Self* inverse) const; 

  virtual void TransformPoint(const InputPointType  &input, OutputPointType &output ) const;
  
protected:

  SwitchableAffineTransform(unsigned int outputDims,
                            unsigned int paramDims);
  SwitchableAffineTransform();      
  
  virtual ~SwitchableAffineTransform() = 0;

  void PrintSelf(std::ostream &s, Indent indent) const;

  virtual void ComputeMatrixAndOffset(void) {};
  
  virtual void ComputeParametersFromMatrixAndOffset(void) {};
  
  void SetDefaultRelativeParameterWeightings( void );

  const InverseMatrixType & GetInverseMatrix( void ) const;

  const InverseMatrixType & GetVarInverseMatrix( void ) const
    { return m_InverseMatrix; };
  
  void SetVarInverseMatrix(const InverseMatrixType & matrix) const
    { m_InverseMatrix = matrix; m_InverseMatrixMTime.Modified(); };
  
  bool InverseMatrixIsOld(void) const
    { if(m_MatrixMTime != m_InverseMatrixMTime)
        { return true; } else { return false; } };

  MatrixType                   m_Matrix;        // Matrix of the transformation
  OutputVectorType             m_Offset;        // Offset of the transformation  
  
private:

  SwitchableAffineTransform(const Self & other); // Purposely not implemented
  const Self & operator=( const Self & );        // Purposely not implemented

  InputPointType              m_Center;
  TranslationType             m_Translation;
  RotationType                m_Rotation;
  ScaleType                   m_Scale;
  SkewType                    m_Skew;
  bool                        m_OptimiseRotation;
  bool                        m_OptimiseTranslation;
  bool                        m_OptimiseScale;
  bool                        m_OptimiseSkew;

  RelativeParameterWeightingType m_TranslationRelativeWeighting;
  RelativeParameterWeightingType m_RotationRelativeWeighting;
  RelativeParameterWeightingType m_ScaleRelativeWeighting;
  RelativeParameterWeightingType m_SkewRelativeWeighting;
  RelativeParameterWeightingType m_AllRelativeWeightings;

  mutable InverseMatrixType   m_InverseMatrix;  // Inverse of the matrix
  mutable bool                m_Singular;       // Is m_Inverse singular?
  
  FullAffineTransformPointer  m_AffineMatrixTransform;
  
  TimeStamp                   m_MatrixMTime;
  mutable TimeStamp           m_InverseMatrixMTime;

}; // class SwitchableAffineTransform

}  // namespace itk

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

#endif /* __itkSwitchableAffineTransform_h */