[Insight-users] rectangular image registration issue

[Sabine Husse]

Hi all,

I am trying to register 2D images (spect modality) witk itk Mattes mutual
information. Visually, offsets between both images are very small
(translation of ~ 5 mms). I am initializing transform parameters to 0 for
the 3 parameters of euler 2D transform. However, after registration I got a
rotation around 90 degres. Plus, when I am observing optimizer iterations,
first iteration is already 1.6 radian for rotation parameter, even if my
initial rotation parameter was 0.

Here is the observer:

MultiResolution Level : 0
0 = -0.372419 : [1.7683, -0.0955618, -0.929501]
1 = -0.0383217 : [1.92583, -2.08874, -0.978926]
2 = -0.0362538 : [1.69325, -1.10031, -2.21542]
3 = -0.0377709 : [1.56321, -1.28384, -3.79953]
4 = -0.0419196 : [1.66387, -2.53705, -3.55923]
141 = -0.042436 : [1.51798, -17.6144, -6.53557]
142 = -0.0423536 : [1.52173, -17.635, -6.52587]
143 = -0.042424 : [1.51635, -17.6507, -6.54188]

Iterations    = 145
Metric value  = -0.0424182
Stop Condition  = RegularStepGradientDescentOptimizer: Step too small after
144
iterations. Current step (0.0184467) is less than minimum step (0.02).
Result =
Translation X = 17.6507 mm
Translation Y = 6.54188 mm
Rotation X = -86.8805 degres

Here is part of the print out of registration:
  CurrentLevel: 1
  InitialTransformParameters: [0, 0, 0]
  InitialTransformParametersOfNextLevel: [0, 0, 0]
  LastTransformParameters: [1.51635, -17.6507, -6.54188]
  FixedImageRegion: ImageRegion (103F17EC)
  Dimension: 2
  Index: [0, 0]
  Size: [256, 1024]

  FixedImageRegion at level 0: ImageRegion (04DA4C98)
  Dimension: 2
  Index: [0, 0]
  Size: [128, 512]

  FixedImageRegion at level 1: ImageRegion (04DA4CAC)
  Dimension: 2
  Index: [0, 0]
  Size: [256, 1024]


It's clear that I am missing something. Please help. Here is my code:

  ImageType::RegionType fixedImageRegion = fixedImage->GetBufferedRegion();
  ImageType::RegionType movingImageRegion =
movingImage->GetBufferedRegion();  typedef itk::Euler2DTransform< 
double > TransformType;
  TransformType::Pointer transform = TransformType::New();
  registration_->SetTransform(transform);

  // Initialize transform parameters
  int transformParametersTotal =
registration_->GetTransform()->GetNumberOfParameters();
  typedef RegistrationType::ParametersType ParametersType;
  ParametersType initialParameters(transformParametersTotal);
  initialParameters.Fill(0.0);
  registration_->SetInitialTransformParameters(initialParameters);

  metric_->SetNumberOfHistogramBins(50);
  metric_->UseAllPixelsOn();

  registration_->SetFixedImage(fixedImage);
  registration_->SetMovingImage(movingImage);    
registration_->SetFixedImageRegion(fixedImageRegion);

  typedef OptimizerType::ScalesType ScalesType;
  int transformParametersTotal =
registration_->GetTransform()->GetNumberOfParameters();
  ScalesType parametersScales(transformParametersTotal);
  parametersScales.Fill( 1.0 );

  for (int j = 1; j < transformParametersTotal; j++ ){
    parametersScales[j] = 1.0 / 1000;
  }        optimizer_->SetScales(parametersScales);

  // The initial step length is defined with SetMaximumStepLength()
  optimizer_->SetMaximumStepLength(2.0);
  optimizer_->SetMinimumStepLength(0.02);

  // Define number of iterations per multi-resolution level
  optimizer_->SetNumberOfIterations( 200);

  optimizer_->SetRelaxationFactor(0.8);    
OptimizerIterationUpdate::Pointer observer =
OptimizerIterationUpdate::New();
  optimizer_->AddObserver( itk::IterationEvent(), observer );

   typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
   CommandType::Pointer command = CommandType::New();
   registration_->AddObserver( itk::IterationEvent(), command );

    registration_->Update();
    // Output parameters
    ParametersType finalParameters =
registration_->GetLastTransformParameters();
     result[0] = ;
      result[1] =
      result[2] = ;                    // Print out results
    std::cout << "Result = " << std::endl;
    std::cout << " Translation X = " << - finalParameters[1]  << std::endl;
    std::cout << " Translation Y = " <<- finalParameters[2]; << std::endl;
     std::cout << " Rotation X = " << - finalParameters[0] * RAD_TO_DEG <<
std::endl;

Sabine Husse