freCenteredAffineTransformController.cxx

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: freCenteredAffineTransformController.cxx,v $
  Language:  C++



  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.

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

#include "freCenteredAffineTransformController.h"

#include "freIntensityImageMediaController.h"
#include "freTransformSetupAdaptor.h"

#include "itkPoint.h"
#include "itkCenteredTransformInitializer.h"

namespace FREE
{



CenteredAffine2DTransformController::
CenteredAffine2DTransformController()
{
  this->UpdateControllerID(ControllerID::CenteredAffine2DTransformController);
  this->m_Description = "A affine transformation (2 Dimensions), composed of rotation (around a specified center), scaling, shearing and translation.";
};

void
CenteredAffine2DTransformController::
GenerateProfile(CtrlProfile::ControllerProfile& profile, 
                const SessionComponentCache* pComponentCache,
                bool bRegardOldSetup) const
{
  Superclass::GenerateProfile(profile,pComponentCache,bRegardOldSetup);

  //Parameters
  profile.Parameters().AddParameter(cParam_InitialTransformValues,CtrlProfile::Parameter::PVTDouble,"1..4: elements of the 2x2 matrix\n5..6: rotation center (x/y)\n7..8: elements of the translation vector.",8,"0");
  profile.Parameters().AddParameter(cParam_CurrentTransformValues,CtrlProfile::Parameter::PVTDouble,"1..4: elements of the 2x2 matrix\n5..6: rotation center (x/y)\n7..8: elements of the translation vector.",8,"0");
  profile.Parameters().AddParameter(cParam_TransformScale,CtrlProfile::Parameter::PVTDouble,"1..4: elements of the 2x2 matrix\n5..6: rotation center (x/y)\n7..8: elements of the translation vector.",8,"1",-1,true);

  profile.MediaMap().AddMedia("movingInitialImage",ControllerID::IntensityImage2DMediaController,DASet,2);
  profile.MediaMap().AddMedia("fixedInitialImage",ControllerID::IntensityImage2DMediaController,DASet,2);
};

bool
CenteredAffine2DTransformController::
SetInitialTransformValues( ComponentType* pTransformComponent,
                           SessionComponentCache* pComponentCache,
                           SessionInfo* pSessionInfo,
                           const int& iInitializeByITV) const
{
  int iCenterType;
  
  if (!pComponentCache->SetupIsAssigned()) throwCtrlExceptionMacro("","Cannot initialize transform values, transform setup is not assigned in cache.");
  ComponentSetup* pTransformSetup = pComponentCache->Setup();

  if (!pTransformSetup->Parameters().GetParameterValue("CenterType",iCenterType)) throwExceptionMacro("Missing parameter: " << "CenterType");
  
  TransformSetupAdaptor adaptor(pTransformSetup);
  ReferencePointsAccessor<2> pointsAccessor(pComponentCache);

        //2D version
        //Calculate the centers
        ParameterArrayType fixedCenter(2);
        ParameterArrayType movingCenter(2);

        if (2==(iInitializeByITV&2))
        { //Center is stored in the ITV and should be used
          // center point x/y
          ParameterArrayType itv = adaptor.GetInitialTransformParameters();
          for (int iIndex=0; iIndex<2; iIndex++)
          {
                movingCenter[iIndex] = itv[4+iIndex];
                fixedCenter[iIndex] = itv[4+iIndex] - itv[6+iIndex];
          };
        }
        else
        { 
          if (iCenterType==0)
          { //First paire of reference points are the center points
      if (!pointsAccessor.MovingPointsExist()) throwCtrlExceptionMacro("","Cannot set moving center via reference points. No moving reference point set defined. Please check media linkage.");
      if (!pointsAccessor.FixedPointsExist()) throwCtrlExceptionMacro("","Cannot set fixed center via reference points. No fixed reference point set defined. Please check media linkage.");
                fixedCenter = pointsAccessor.GetFixedPointArray(0);
                movingCenter = pointsAccessor.GetMovingPointArray(0);
          }
          else
          {
                //Center should be calculated
                CenterInitializer<PixelType,2> centerInitializer;
          
                if (iCenterType==1) centerInitializer.GeometryOn();
                else centerInitializer.MomentsOn();

    typedef ImageTypes<2>::InternalImageType InternalImageType;
    
    SessionAccessor::GenericMediaPointer smpGenericInitialFixedImage = SessionAccessor::GetLinkedMedia("fixedInitialImage",pComponentCache,pSessionInfo);
    if (smpGenericInitialFixedImage.IsNull()) throwCtrlExceptionMacro("","Cannot set fixed center via initial image. No initial fixed image linked. Please check media linkage.");
    SessionAccessor::GenericMediaPointer smpGenericInitialMovingImage = SessionAccessor::GetLinkedMedia("movingInitialImage",pComponentCache,pSessionInfo);
    if (smpGenericInitialMovingImage.IsNull()) throwCtrlExceptionMacro("","Cannot set moving center via initial image. No initial fixed image linked. Please check media linkage.");
    
    InternalImageType* pInitialFixedImage = 0;
    InternalImageType* pInitialMovingImage = 0;

    try
    {
      pInitialFixedImage = dynamic_cast<InternalImageType*>(smpGenericInitialFixedImage.GetPointer());
    }
    catchAllNPassMacro("Error while casting initial fixed image");
    if (!pInitialFixedImage) throwCtrlExceptionMacro("","Error while casting initial fixed image.");

    try
    {
      pInitialMovingImage = dynamic_cast<InternalImageType*>(smpGenericInitialMovingImage.GetPointer());
    }
    catchAllNPassMacro("Error while casting initial fixed image");
    if (!pInitialMovingImage) throwCtrlExceptionMacro("","Error while casting initial fixed image.");

    centerInitializer.SetImage(pInitialFixedImage);
                centerInitializer.GetCenter(fixedCenter);

                centerInitializer.SetImage(pInitialMovingImage);
                centerInitializer.GetCenter(movingCenter);
          }
        }

        //Calculate the translation
        ParameterArrayType translation(2);
        for (int iIndex=0; iIndex<2; iIndex++)
        {
          translation[iIndex] = movingCenter[iIndex]-fixedCenter[iIndex];
        };

    //add the calculated initial translation and center to the setup - to be present
    //for the Registrator - when using set initial position.
        for (int iIndex=0; iIndex<2; iIndex++)
        {
          pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,fixedCenter[iIndex],4+iIndex);
          pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,translation[iIndex],6+iIndex);
        };

        if (!(1==(iInitializeByITV&1)))
        { //Rotation should be calculated
    if (!pointsAccessor.MovingPointsExist()) throwCtrlExceptionMacro("","Cannot estimate rotation via reference points. No moving reference point set defined. Please check media linkage.");
    if (!pointsAccessor.FixedPointsExist()) throwCtrlExceptionMacro("","Cannot estimate rotation via via reference points. No fixed reference point set defined. Please check media linkage.");

          RotationInitializer<2> rotationInitializer;
          
          rotationInitializer.SetCenterPair(movingCenter,fixedCenter);

          //Add all references
          int iIndex = 0;
          if (iCenterType==0) iIndex = 1; //first reference is the center point;
          for (iIndex; iIndex<pointsAccessor.GetFixedPointsCount(); iIndex++)
          {
                rotationInitializer.AddReferencePair(pointsAccessor.GetMovingPointArray(iIndex),pointsAccessor.GetFixedPointArray(iIndex));
          };

          //set schedule
          rotationInitializer.ClearSchedule(false);
          rotationInitializer.AddToSchedule(RA_Z);

          rotationInitializer.ComputeRotations();

          //add angel to the parameters. Angel must be inverted, because itk uses back transform
          //(fixed to moving) and the initial calculation uses forward transform to allow an intuitive
          //parameter setting
          double dAngel = -1*rotationInitializer.GetRotation(RA_Z);
          //Set rotation matrix
          pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,cos(dAngel),0);
          pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,-sin(dAngel),1);
          pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,sin(dAngel),2);
          pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,cos(dAngel),3);
        }

  return true;
};



CenteredAffine3DTransformController::
CenteredAffine3DTransformController()
{
  this->UpdateControllerID(ControllerID::CenteredAffine3DTransformController);
  this->m_Description = "A affine transformation (3 Dimensions), composed of rotation (around a specified center), scaling, shearing and translation.";
};

void
CenteredAffine3DTransformController::
GenerateProfile(CtrlProfile::ControllerProfile& profile, 
                const SessionComponentCache* pComponentCache,
                bool bRegardOldSetup) const
{
  Superclass::GenerateProfile(profile,pComponentCache,bRegardOldSetup);

  //Parameters
  profile.Parameters().AddParameter("InitialRotationAxis",Parameter::PVTInteger,"Indicates the type of rotation used for initialisation, if the reference initialisation is used.\n0: The rotation axis is calculated as cross product of the reverence - centerpoint vectors.\n1: x axis\n2: y axis\n3: z axis\n4:Composing (in reading order) a z rotation, y rotation, x rotation\n5:Smart composing (like composing, but rotation is weightened).",1,"0");
  profile.Parameters().AddParameter(cParam_InitialTransformValues,CtrlProfile::Parameter::PVTDouble,"1..9: elements of the 3x3 matrix\n10..12: rotation center (x/y/z)\n13..15: elements of the translation vector.",15,"0");
  profile.Parameters().AddParameter(cParam_CurrentTransformValues,CtrlProfile::Parameter::PVTDouble,"1..9: elements of the 3x3 matrix\n10..12: rotation center (x/y/z)\n13..15: elements of the translation vector.",15,"0");
  profile.Parameters().AddParameter(cParam_TransformScale,CtrlProfile::Parameter::PVTDouble,"1..9: elements of the 3x3 matrix\n10..12: rotation center (x/y/z)\n13..15: elements of the translation vector.",15,"1",-1,true);

  profile.MediaMap().AddMedia("movingInitialImage",ControllerID::IntensityImage3DMediaController,DASet,3);
  profile.MediaMap().AddMedia("fixedInitialImage",ControllerID::IntensityImage3DMediaController,DASet,3);
};

bool
CenteredAffine3DTransformController::
SetInitialTransformValues( ComponentType* pTransformComponent,
                           SessionComponentCache* pComponentCache,
                           SessionInfo* pSessionInfo,
                           const int& iInitializeByITV) const
{
  int iCenterType;

  if (!pComponentCache->SetupIsAssigned()) throwCtrlExceptionMacro("","Cannot initialize transform values, transform setup is not assigned in cache.");
  ComponentSetup* pTransformSetup = pComponentCache->Setup();

  if (!pTransformSetup->Parameters().GetParameterValue("CenterType",iCenterType)) throwExceptionMacro("Missing parameter: " << "CenterType");

  //3D version
  TransformSetupAdaptor adaptor(pTransformSetup);
  ReferencePointsAccessor<3> pointsAccessor(pComponentCache);

  //Calculate the Centers
  ParameterArrayType fixedCenter(3);
  ParameterArrayType movingCenter(3);

  if (2==(iInitializeByITV&2))
  { //Center is stored in the ITV and should be used
    // center point x/y/z
    ParameterArrayType itv = adaptor.GetInitialTransformParameters();
    for (int iIndex=0; iIndex<3; iIndex++)
    {
      movingCenter[iIndex] = itv[9+iIndex];
      fixedCenter[iIndex] = itv[9+iIndex] - itv[12+iIndex];
    };
  }
  else
  { 
    if (iCenterType==0)
    { //First paire of reference points are the center points
      if (!pointsAccessor.MovingPointsExist()) throwCtrlExceptionMacro("","Cannot set moving center via reference points. No moving reference point set defined. Please check media linkage.");
      if (!pointsAccessor.FixedPointsExist()) throwCtrlExceptionMacro("","Cannot set fixed center via reference points. No fixed reference point set defined. Please check media linkage.");
      fixedCenter = pointsAccessor.GetFixedPointArray(0);
      movingCenter = pointsAccessor.GetMovingPointArray(0);
    }
    else
    {
      //Center should be calculated
      CenterInitializer<PixelType,3> centerInitializer;

      if (iCenterType==1) centerInitializer.GeometryOn();
      else centerInitializer.MomentsOn();

      typedef ImageTypes<3>::InternalImageType InternalImageType;
      
      SessionAccessor::GenericMediaPointer smpGenericInitialFixedImage = SessionAccessor::GetLinkedMedia("fixedInitialImage",pComponentCache,pSessionInfo);
      if (smpGenericInitialFixedImage.IsNull()) throwCtrlExceptionMacro("","Cannot set fixed center via initial image. No initial fixed image linked. Please check media linkage.");
      SessionAccessor::GenericMediaPointer smpGenericInitialMovingImage = SessionAccessor::GetLinkedMedia("movingInitialImage",pComponentCache,pSessionInfo);
      if (smpGenericInitialMovingImage.IsNull()) throwCtrlExceptionMacro("","Cannot set moving center via initial image. No initial fixed image linked. Please check media linkage.");
      
      InternalImageType* pInitialFixedImage = 0;
      InternalImageType* pInitialMovingImage = 0;

      try
      {
        pInitialFixedImage = dynamic_cast<InternalImageType*>(smpGenericInitialFixedImage.GetPointer());
      }
      catchAllNPassMacro("Error while casting initial fixed image");
      if (!pInitialFixedImage) throwCtrlExceptionMacro("","Error while casting initial fixed image.");

      try
      {
        pInitialMovingImage = dynamic_cast<InternalImageType*>(smpGenericInitialMovingImage.GetPointer());
      }
      catchAllNPassMacro("Error while casting initial fixed image");
      if (!pInitialMovingImage) throwCtrlExceptionMacro("","Error while casting initial fixed image.");

      centerInitializer.SetImage(pInitialFixedImage);
      centerInitializer.GetCenter(fixedCenter);

      centerInitializer.SetImage(pInitialMovingImage);
      centerInitializer.GetCenter(movingCenter);
    }
  }

  //Calculate the translation
  ParameterArrayType translation(3);
  for (int iIndex=0; iIndex<3; iIndex++)
  {
    translation[iIndex] = movingCenter[iIndex]-fixedCenter[iIndex];
  };

    //add the calculated initial translation and center to the setup - to be present
    //for the Registrator - when using set initial position.
  for (int iIndex=0; iIndex<3; iIndex++)
  {
    pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,fixedCenter[iIndex],9+iIndex);
    pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,translation[iIndex],12+iIndex);
  };

  if (!(1==(iInitializeByITV&1)))
  { //Rotation should be calculated
    if (!pointsAccessor.MovingPointsExist()) throwCtrlExceptionMacro("","Cannot estimate rotation via reference points. No moving reference point set defined. Please check media linkage.");
    if (!pointsAccessor.FixedPointsExist()) throwCtrlExceptionMacro("","Cannot estimate rotation via via reference points. No fixed reference point set defined. Please check media linkage.");

    RotationInitializer<3> rotationInitializer;

    rotationInitializer.SetCenterPair(movingCenter,fixedCenter);

    //Add all references
    int iIndex = 0;
    if (iCenterType==0) iIndex = 1; //first reference is the center point;
    for (iIndex; iIndex<pointsAccessor.GetFixedPointsCount(); iIndex++)
    {
      rotationInitializer.AddReferencePair(pointsAccessor.GetMovingPointArray(iIndex),pointsAccessor.GetFixedPointArray(iIndex));
    };

    //set schedule
    int iInitialRotationAxis;
    if (!pTransformSetup->Parameters().GetParameterValue("InitialRotationAxis",iInitialRotationAxis)) throwExceptionMacro("Missing parameter: " << "InitialRotationAxis");

    rotationInitializer.ClearSchedule(false);
    if ((iInitialRotationAxis<1) | (4 == (iInitialRotationAxis & 4))) rotationInitializer.AddToSchedule(RA_Z);
    if ((iInitialRotationAxis<1) | (2 == (iInitialRotationAxis & 2))) rotationInitializer.AddToSchedule(RA_Y);
    if ((iInitialRotationAxis<1) | (1 == (iInitialRotationAxis & 1))) rotationInitializer.AddToSchedule(RA_X);

    rotationInitializer.SetUseRotationWeight(iInitialRotationAxis==-1);

    rotationInitializer.ComputeRotations();

    typedef itk::AffineTransform<ScalarType,3> CalculateTransformType;
    CalculateTransformType::Pointer rotationTransform = CalculateTransformType::New();

    for ( RotationAxisType axis = RA_Z; axis>=RA_X; axis--)
    {
      Vector3DType rotationAxis;
      rotationAxis.Fill(0);
      rotationAxis[axis]=1;

      //add angel to the parameters. Angel must be inverted, because itk uses back transform
      //(fixed to moving) and the initial calculation uses forward transform to allow an intuitive
      double dAngel = -1*rotationInitializer.GetRotation(axis);

      //compute matrix for the rotation
      CalculateTransformType::Pointer oneRotation = CalculateTransformType::New();
      oneRotation->Rotate3D(rotationAxis,dAngel);

      //add the rotation of this step to the whole transformation
      rotationTransform->Compose(oneRotation,false);
    };

    //Set rotation matrix
    CalculateTransformType::MatrixType rotation = rotationTransform->GetMatrix();

    for (int iIndex=0; iIndex<9; iIndex++)
    {
      pTransformSetup->Parameters().SetParameterValue(cParam_CurrentTransformValues,rotation[(iIndex/3)][(iIndex%3)],iIndex);
    }
  };
  return true;
};



} //end of namespace free

---