[Insight-users] Fine-registration gone wrong

[vanderLaak J.]

Hi all,
Currently I'm trying to register a set of 50 images which were taken by
transmitted light microscopy from serial histological sections
(immunohistochemically stained for bloodvessels, without any
counterstain). Sections were manually 'optically' (using the microscope
stage) aligned as much as possible (including rotation). The result of
this is a set of RGB images showing only the vascular wall in brown,
which are already pretty well registered. I would like to have ITK do a
fine-registration, preferably an affine-transform, to account for small
stretches in the tissue. I tried various registration metrics,
transforms and interpolators, and came up with the routine below, which
produces good results for a number of the images. However, in some cases
unpredictable (to me, at least ;-) behaviour occurs: although images are
pretty close to good alignment the routine doesn't converge and produces
results that are way off (I could send examples of these, but don't want
to burden this list too much). Below is the code-fragment, and part of
the output of the registration which produces incorrect results. My
questions: is it possible to put a limit to the registration paramaters
(e.g. x translation no more than 20 pixels). Does anyone have a clue
what's wrong with my code (or images)? How does the registration deal
with RGB images? Can I detect the procedure going wrong, in order to
(worst case) keep the manual registration and continue registration with
the next images? In the output I often see that the last metric value is
not the smallest, why doesn't the registration process take the set of
parameters with lowest (highest) metric????
Regards,
Jeroen

Code:
  const unsigned int Dimension =3D 2;
  typedef float PixelType;

  typedef itk::Image< PixelType, Dimension >  FixedImageType;
  typedef itk::Image< PixelType, Dimension >  MovingImageType;
  typedef itk::CenteredAffineTransform< double, Dimension >
TransformType;
  typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
  typedef itk::NormalizedCorrelationImageToImageMetric< FixedImageType,
MovingImageType > MetricType;
  typedef itk:: BSplineInterpolateImageFunction<MovingImageType,double>
InterpolatorType;
  typedef itk::ImageRegistrationMethod<FixedImageType,MovingImageType >
RegistrationType;

  MetricType::Pointer metric =3D MetricType::New();
  OptimizerType::Pointer optimizer =3D OptimizerType::New();
  InterpolatorType::Pointer interpolator =3D InterpolatorType::New();
  RegistrationType::Pointer registration =3D RegistrationType::New();

  registration->SetMetric( metric );
  registration->SetOptimizer( optimizer );
  registration->SetInterpolator( interpolator );

  TransformType::Pointer transform =3D TransformType::New();
  registration->SetTransform(transform);

  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
  FixedImageReaderType::Pointer  fixedImageReader  =3D
FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader =3D
MovingImageReaderType::New();

  fixedImageReader->SetFileName(  argv[1] );
  movingImageReader->SetFileName( argv[2] );

  registration->SetFixedImage(fixedImageReader->GetOutput());
  registration->SetMovingImage(movingImageReader->GetOutput());
  fixedImageReader->Update();

  registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );

  typedef itk::CenteredTransformInitializer<TransformType,
FixedImageType,MovingImageType > TransformInitializerType;
  TransformInitializerType::Pointer initializer =3D
TransformInitializerType::New();

  initializer->SetTransform(transform);
  initializer->SetFixedImage(fixedImageReader->GetOutput());
  initializer->SetMovingImage(movingImageReader->GetOutput());
  initializer->MomentsOn();
  initializer->InitializeTransform();

  std::cout << "Initial Transform Parameters " << std::endl;
  std::cout << transform->GetParameters() << std::endl;

  registration->SetInitialTransformParameters(
transform->GetParameters() );

  typedef OptimizerType::ScalesType OptimizerScalesType;
  OptimizerScalesType optimizerScales(transform->GetNumberOfParameters()
);
  const double translationScale =3D 1.0 / 1000.0;

  optimizerScales[0] =3D 1.0;
  optimizerScales[1] =3D  1.0;
  optimizerScales[2] =3D  1.0;
  optimizerScales[3] =3D  1.0;
  optimizerScales[4] =3D  translationScale;
  optimizerScales[5] =3D  translationScale;
  optimizerScales[6] =3D  translationScale;
  optimizerScales[7] =3D  translationScale;

  optimizer->SetScales( optimizerScales );

  optimizer->SetMaximumStepLength( .1 );
  optimizer->SetMinimumStepLength( 0.001 );=20
  optimizer->SetNumberOfIterations( 400 );
  optimizer->MinimizeOn();

  CommandIterationUpdate::Pointer observer =3D
CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer );

  try=20
    { registration->StartRegistration();  }=20
  catch( itk::ExceptionObject & err )=20
    { std::cerr << "ExceptionObject caught !" << std::endl;=20
      std::cerr << err << std::endl;=20
      return -1; }=20

  OptimizerType::ParametersType finalParameters =3D
registration->GetLastTransformParameters();

  const double finalRotationCenterX =3D finalParameters[4];
  const double finalRotationCenterY =3D finalParameters[5];
  const double finalTranslationX    =3D finalParameters[6];
  const double finalTranslationY    =3D finalParameters[7];

  const unsigned int numberOfIterations =3D
optimizer->GetCurrentIteration();
  const double bestValue =3D optimizer->GetValue();

  std::cout << "Result =3D " << std::endl;
  std::cout << " Center X      =3D " << finalRotationCenterX  <<
std::endl;
  std::cout << " Center Y      =3D " << finalRotationCenterY  <<
std::endl;
  std::cout << " Translation X =3D " << finalTranslationX  << std::endl;
  std::cout << " Translation Y =3D " << finalTranslationY  << std::endl;
  std::cout << " Iterations    =3D " << numberOfIterations << std::endl;
  std::cout << " Metric value  =3D " << bestValue          << std::endl;


Ouput:
Initial Transform Parameters=20
[1, 0, 0, 1, 385.213, 286.714, 0.907753, -0.0524303]
0   -0.981916   [1.03616, -0.0124404, 0.00123626, 1.02807, 385.213,
286.714, 0.927529, 0.0333425]
1   -0.977545   [1.00046, 0.0363107, 0.00624682, 1.0067, 385.211,
286.713, 0.988304, 0.0799002]
2   -0.977815   [0.999092, 0.00879042, -0.00218918, 1.011, 385.21,
286.714, 0.953614, 0.10103]
3   -0.980874   [0.991963, -0.0173192, -0.00861846, 1.00402, 385.21,
286.715, 0.913292, 0.108212]
4   -0.980221   [0.995119, -0.0216787, -0.0217679, 0.994541, 385.21,
286.714, 0.877521, 0.0777477]
5   -0.980084   [1.00374, -0.0234345, -0.0352056, 0.987303, 385.209,
286.714, 0.868013, 0.0319443]
6   -0.979511   [1.00972, -0.0188327, -0.0443085, 0.979458, 385.207,
286.713, 0.865454, -0.015898]
7   -0.978891   [1.00964, -0.00596139, -0.048326, 0.971347, 385.205,
286.712, 0.871117, -0.0629617]
...
203   -0.969197   [1.4177, 0.383298, -0.0274694, 0.987543, 384.906,
286.44, 2.18153, -1.08256]
204   -0.969208   [1.41863, 0.383454, -0.0273697, 0.987456, 384.907,
286.441, 2.18077, -1.08345]
205   -0.969217   [1.4199, 0.383883, -0.0268136, 0.987596, 384.907,
286.441, 2.18039, -1.0831]
206   -0.969232   [1.42031, 0.384117, -0.0273033, 0.987269, 384.907,
286.441, 2.18001, -1.08439]
Result =3D=20
 Center X      =3D 384.907
 Center Y      =3D 286.441
 Translation X =3D 2.18001
 Translation Y =3D -1.08439
 Iterations    =3D 208
 Metric value  =3D -0.969233