freTransformInitializer.h

Go to the documentation of this file.

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

  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.1.0.0
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: freTransformInitializer.h,v $
  Language:  C++


  Date:      $Date: 2004/10/03 $
  Version:   $Revision: 1.0 $

  Copyright (c) 2007 Ralf o Floca (Department of Medical Informatics,
  Institute for Medical Biometry and Informatics, University of Heidelberg,
  Germany). All rights reserved.
  See FREECopyright.txt or http://www.mi.med.uni-hd.de/free/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 __freTransformInitializer_h
#define __freTransformInitializer_h

#include "freElementals.h"
#include "freImageTypes.h"
#include "freExceptions.h"

#include "itkImage.h"
#include "itkVector.h"
#include "itkPoint.h"
#include "itkImageMomentsCalculator.h"
#include "itkAffineTransform.h"

namespace FREE
{

enum RotationAxisType 
{
        RA_X = 0, 
        RA_Y = 1, 
        RA_Z = 2, 
};

// Define a postfix increment operator for RotationAxisType.
inline RotationAxisType operator++( RotationAxisType &ra, int )
{
   return ra = (RotationAxisType)(ra + 1);
}

// Define a postfix increment operator for RotationAxisType.
inline RotationAxisType operator--( RotationAxisType &ra, int )
{
   return ra = (RotationAxisType)(ra - 1);
}

template <typename TPixel, int iDimension>
class CenterInitializer
{
public: 
  typedef itk::Point<ScalarType,iDimension> PointType;
  typedef itk::Vector<ScalarType,iDimension> VectorType;
  typedef itk::Image<TPixel,iDimension> TImage; //is same like InternalImage of free

  void GeometryOn()
  { 
          if (m_bUseMoments)
          {
            m_bUseMoments = false;
            m_bInitialized = false;
          }
  };
  
  void MomentsOn()
  { 
          if (!m_bUseMoments)
          {
            m_bUseMoments = true;
            m_bInitialized = false;
          }
  };

  void SetImage(TImage* pImage)
  { 
          m_pImage = pImage;
    m_bInitialized = false;
  };

  void ComputeCenter()
  {
          if( m_bUseMoments )
    {
            // Here use the center of mass for each image.
            typedef itk::ImageMomentsCalculator< TImage >  ImageCalculatorType;

            typename ImageCalculatorType::Pointer calculator = ImageCalculatorType::New();

            try
            {
                    calculator->SetImage(  m_pImage );
                    calculator->Compute();
    
                    m_Center = calculator->GetCenterOfGravity();
            }
            catch(...)
            {
                    //itk throws error if the mass is zero, so the image is all black ->
                    //use geometrical center.
                    m_bUseMoments = false;
            }
    }
        
          if( !m_bUseMoments )
    {
            // Here use the geometrical center of each image.

            const typename TImage::SpacingType& rSpacing = m_pImage->GetSpacing();
            const typename TImage::PointType& rOrigin = m_pImage->GetOrigin();
   
            typename TImage::SizeType size = m_pImage->GetLargestPossibleRegion().GetSize();
    
            for( unsigned int iIndex=0; iIndex<iDimension; iIndex++)
      {
                    m_Center[iIndex] = rOrigin[iIndex] + rSpacing[iIndex] * (size[iIndex] /2.0);
      }
    }
          m_bInitialized = true;
  };

  void GetCenter(VectorType& rVector)
  {
          if (!m_bInitialized) ComputeCenter();

          rVector = m_Center;
  };

  void GetCenter(ParameterArrayType& rVector)
  {
          if (!m_bInitialized) ComputeCenter();

          rVector.SetSize(iDimension);
          for (int iIndex = 0; iIndex<iDimension; iIndex++)
          {
            rVector.SetElement(iIndex,m_Center[iIndex]);
          }
  };

  CenterInitializer()
  {
          m_bUseMoments = true;
          m_bInitialized = false;
          m_pImage = NULL;
  };

  ~CenterInitializer() {};

private:
  bool m_bUseMoments;
  bool m_bInitialized;

  VectorType m_Center;

  TImage* m_pImage;
};

template <int iDimension>
class RotationInitializer
{
public: 
  typedef itk::Vector<ScalarType,iDimension> VectorType;

  void ClearSchedule(bool bToDefault = true) 
  { 
    m_Schedule.clear();
          if(bToDefault)
          {
      m_Schedule.push_back(RA_Z);
            m_Schedule.push_back(RA_Y);
            m_Schedule.push_back(RA_X);
          }
          m_bInitialized = false;
  };
  
  void AddToSchedule(RotationAxisType axis)
  { 
          m_Schedule.push_back(axis);
          m_bInitialized = false;
  };

  void SetMovingCenter(VectorType center)
  { 
          m_MovingCenter.Fill(0.0);
          for (int iIndex=0; iIndex<iDimension; iIndex++) m_MovingCenter.SetNthComponent(iIndex,center[iIndex]);
          m_bInitialized = false;
  };

  void SetFixedCenter(VectorType center)
  { 
          m_FixedCenter.Fill(0);
          for (int iIndex=0; iIndex<iDimension; iIndex++) m_FixedCenter.SetNthComponent(iIndex,center[iIndex]);
          m_bInitialized = false;
  };

  void SetCenterPair(VectorType movingCenter,VectorType fixedCenter)
  { 
    SetFixedCenter(fixedCenter);
    SetMovingCenter(movingCenter);
  };

  void SetCenterPair(ParameterArrayType movingCenter,ParameterArrayType fixedCenter)
  { 
          VectorType movingCenterVector;
          VectorType fixedCenterVector;

          for (int iIndex = 0; iIndex<iDimension; iIndex++)
          {
            movingCenterVector.SetNthComponent(iIndex,movingCenter[iIndex]);
            fixedCenterVector.SetNthComponent(iIndex,fixedCenter[iIndex]);
          };

    SetFixedCenter(fixedCenterVector);
    SetMovingCenter(movingCenterVector);
  };

  void ClearReferences()
  {
          this->m_FixedReferences.RemoveAll();
          this->m_MovingReferences.RemoveAll();
          m_bInitialized = false;
  };

  void AddReferencePair(VectorType movingReference,VectorType fixedReference)
  {
          Vector3DType movingVector3D;
          Vector3DType fixedVector3D;

          movingVector3D.Fill(0);
          fixedVector3D.Fill(0);

          for (int iIndex=0; iIndex<iDimension; iIndex++)
          {
            movingVector3D.SetNthComponent(iIndex,movingReference[iIndex]);
            fixedVector3D.SetNthComponent(iIndex,fixedReference[iIndex]);
          }

          m_FixedReferences.push_back(fixedVector3D);
          m_MovingReferences.push_back(movingVector3D);
          m_bInitialized = false;
  };

  void AddReferencePair(ParameterArrayType movingReference,ParameterArrayType fixedReference)
  { 
          VectorType movingReferenceVector;
          VectorType fixedReferenceVector;

          for (int iIndex = 0; iIndex<iDimension; iIndex++)
          {
            movingReferenceVector.SetNthComponent(iIndex,movingReference[iIndex]);
            fixedReferenceVector.SetNthComponent(iIndex,fixedReference[iIndex]);
          };

    AddReferencePair(movingReferenceVector,fixedReferenceVector);
  };

  void ComputeRotations()
  {
    typedef std::vector<double> AnglesArrayType;
          AnglesArrayType xAngles;
          AnglesArrayType yAngles;
          AnglesArrayType zAngles;

          for (int iIndex=0; iIndex<m_FixedReferences.size(); iIndex++)
          {
            //calculate ref points vector relative to the center points
            Vector3DType movingRef = m_MovingReferences[iIndex] - m_MovingCenter;
            Vector3DType fixedRef = m_FixedReferences[iIndex] - m_FixedCenter;
          
            //Calculating vector length for smart composing
            Vector3DType::ValueType lengthFixed = fixedRef.GetNorm();
            Vector3DType smartRef= fixedRef;

            //normalize
            movingRef.Normalize();
            fixedRef.Normalize();

            for (int iScheduleIndex=0; iScheduleIndex<m_Schedule.size(); iScheduleIndex++)
            {
                  RotationAxisType axis = m_Schedule[iScheduleIndex];

                    Vector3DType movingTmpRef =movingRef;
                    Vector3DType fixedTmpRef = fixedRef;
                    movingTmpRef[axis]=0;
                    fixedTmpRef[axis]=0;

                    //Calculating vector length for smart composing
                    Vector3DType smartTmpRef = smartRef;
                    smartTmpRef[axis]=0;
                    Vector3DType::ValueType length2DFixed = smartTmpRef.GetNorm();
                    ScalarType rotationWeight = length2DFixed/lengthFixed;

                    //Preparing angel calculation
                    movingTmpRef.Normalize();
                    fixedTmpRef.Normalize();

                    Vector3DType rotationAxis;
                    rotationAxis.Fill(0);
                    rotationAxis[axis]=1;

                    //calculate angel between ref vectors.
                    ScalarType dAngel = CalculateRotationAngel(fixedTmpRef,movingTmpRef,rotationAxis);
                
                    /* Smart composing should prevent faulty extreme rotations, caused by a small distances between
                     * the reference point and the estimated center of gravity, so that small estimation errors
                     * could have a large impact in the angel. (e.g. if the 2D projection of the reference vector is
                     * only three pixels long, a missestimation of three Pixel could mean a rotation of 45. If
                     * the reference vector has a length of 30 pixel the impact of missestimation would be much smaller.
                     * To compensate this the Angles are weightened by ratio between the length of the reference vector
                     * and his 2D projection.*/
                    if (m_bUseRotationWeight) dAngel *= rotationWeight;//smart composing

                    //compute matrix for the rotation
                    typedef itk::AffineTransform<ScalarType,3> CalculateTransformType;
                    CalculateTransformType::Pointer oneRotation = CalculateTransformType::New();
                    oneRotation->Rotate3D(rotationAxis,dAngel);

                    //update fixedRef for next computation step
                    movingRef = oneRotation->TransformVector(movingRef);

                    //update original fixed referencefixedRef for smart computation in the next step
                    smartRef = oneRotation->TransformVector(smartRef);

                    switch (axis)
                    {
                      case RA_X: xAngles.push_back(dAngel); break;
                      case RA_Y: yAngles.push_back(dAngel); break;
                      case RA_Z: zAngles.push_back(dAngel); break;
                    };
            };
          };

          for (int iIndex=0; iIndex<xAngles.size(); iIndex++)
          {
            m_Angles[RA_X] += xAngles[iIndex];
          }

          for (int iIndex=0; iIndex<yAngles.size(); iIndex++)
          {
            m_Angles[RA_Y] += yAngles[iIndex];
          }

          for (int iIndex=0; iIndex<zAngles.size(); iIndex++)
          {
            m_Angles[RA_Z] += zAngles[iIndex];
          }

          if (xAngles.size()>0) m_Angles[RA_X] = m_Angles[RA_X]/ xAngles.size(); 
          if (yAngles.size()>0) m_Angles[RA_Y] = m_Angles[RA_Y]/ yAngles.size(); 
          if (zAngles.size()>0) m_Angles[RA_Z] = m_Angles[RA_Z]/ zAngles.size(); 

          m_bInitialized = true;
  };

  double GetRotation(RotationAxisType axis)
  {
          if (!m_bInitialized) ComputeRotations();

          return m_Angles[axis];
  };

  bool GetUseRotationWeight() {return m_bUseRotationWeight;};

  void SetUseRotationWeight(bool bUseWeight) 
  {
          m_bUseRotationWeight = bUseWeight;
        m_bInitialized = true;
  };

  RotationInitializer()
  {
          ClearSchedule();
          m_Angles[RA_X] = 0.0;
          m_Angles[RA_Y] = 0.0;
          m_Angles[RA_Z] = 0.0;

          m_bInitialized = false;
          m_bUseRotationWeight = false;
  };

  ~RotationInitializer() {};

private:
  double m_Angles[3];
  bool m_bInitialized;
  bool m_bUseRotationWeight;

  typedef std::vector<Vector3DType> VectorArrayType;
  VectorArrayType m_FixedReferences;
  VectorArrayType m_MovingReferences;

  typedef std::vector<RotationAxisType> ScheduleType;
  
  ScheduleType m_Schedule;

  Vector3DType m_FixedCenter;
  Vector3DType m_MovingCenter;

  double CalculateRotationAngel(Vector3DType vector1, Vector3DType vector2,
                                                                Vector3DType rotationAxis)
  {
          vector1.Normalize();
          vector2.Normalize();

          double dAngel = acos(vector1*vector2);

          //calculate the matrix
          typedef itk::AffineTransform<ScalarType,3> CalculateTransformType;
          CalculateTransformType::Pointer calcTransform = CalculateTransformType::New();

          calcTransform->Rotate3D(rotationAxis,dAngel);
          
          double dTestangel = acos((calcTransform->TransformVector(vector2))*vector1);

          if (dTestangel>dAngel)
          {
            //angel of scalar product wasn't the rotation from vector2 to vector1, but reverse
            dAngel = -1*dAngel;
          };
  
          return dAngel;
  };
};

} //end of namespace free
#endif

---