[Insight-users] initialization transform parameters not taken into account

[Renaud Isabelle]

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;  
 
Isabelle Renaud
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120119/641dd766/attachment.htm>