[Insight-users] Intensity-based 2D-3D registration problems

[Jian Wu]

Hi,
Is there anyone still working on or has been worked on intensity-based 2D-3D 
registration? It seems itk code intensityBased2D3DRegistration.cxx is full 
of bugs and has not been updated since 2005/01/02. I'm wondering who is in 
charge of maintaining this piece of code right now. I've tried to contact 
the authors a few times but got not response. I fixed a few bugs but there 
are more. I wish to contact the people who are still working on the project 
and those who are in charge of code maintaining.

My testing environment:
ITK version: 2.6
OS: Window XP Professional
Compiler: Microsoft VC++ 2005 Standard Edition

Here are the corrections I suggest.

=======================================================================
File Name: IntensityBased2D3DRegistration.cxx
Version: Revision: 1.4

Line Number: 222
Original:
  std::cerr << "       <-step float float>      Maximum and minimum step 
sizes [default: 4 and 0.01]\n";
Replaced by:
  std::cerr << "       <-step float float>      Maximum and minimum step 
sizes [default: 4 and 1]\n";

Line Number: 244
Original:
  char *fileOutput = NULL;
Replaced by:
  char *fileOutput = "DRR_Registered.img";

Line Number: 723
Original:
  imageReader2D->GetOutput()->SetOrigin( origin2D );
Replaced by:
  caster2D->GetOutput()->SetOrigin( origin2D );

Line Number: 753-759
Original:
  std::vector<double> sampledRes2D;
  std::vector<ImageRegionType2D> imageRegionPyramid2D;

  imagePyramid2D->SetNumberOfLevels( nScales );

  sampledRes2D.reserve( nScales );
  imageRegionPyramid2D.reserve( nScales );
Replaced by:
  std::vector<double> sampledRes2D( nScales);
  std::vector<ImageRegionType2D> imageRegionPyramid2D( nScales);

  imagePyramid2D->SetNumberOfLevels( nScales );

Line Number: 842-846
Original:
  std::vector<ImageRegionType3D> imageRegionPyramid3D;

  imagePyramid3D->SetNumberOfLevels( nScales );

  imageRegionPyramid3D.reserve( nScales );
Replaced by:
  std::vector<ImageRegionType3D> imageRegionPyramid3D( nScales);

  imagePyramid3D->SetNumberOfLevels( nScales );

Line Number: 994
Original:
  itk::Optimizer::ScalesType weightings( 6 );
Replaced by:
  itk::Optimizer::ScalesType weightings( transform->GetNumberOfParameters() 
);

Line Number: 1001
Original:
  weightings[5] = 1.;
Replaced by:
  weightings[5] = 1.;
  weightings[6] = 1.;
  weightings[7] = 1.;
  weightings[8] = 1.;

Line Number: 1074-1076
Original:
  const double TranslationAlongX = finalParameters[0];
  const double TranslationAlongY = finalParameters[1];
  const double TranslationAlongZ = finalParameters[2];
Replaced by:
  const double RotationAlongX = finalParameters[0];
  const double RotationAlongY = finalParameters[1];
  const double RotationAlongZ = finalParameters[2];
  const double RotationCenterX = finalParameters[3];
  const double RotationCenterY = finalParameters[4];
  const double RotationCenterZ = finalParameters[5];
  const double TranslationAlongX = finalParameters[6];
  const double TranslationAlongY = finalParameters[7];
  const double TranslationAlongZ = finalParameters[8];

Line Number: 1083-1085
Original:
  std::cout << " Translation X = " << TranslationAlongX  << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
  std::cout << " Translation Z = " << TranslationAlongZ  << std::endl;
Replaced by:
  std::cout << " Rotation Along X = " << RotationAlongX  << std::endl;
  std::cout << " Rotation Along Y = " << RotationAlongY  << std::endl;
  std::cout << " Rotation Along Z = " << RotationAlongZ  << std::endl;
  std::cout << " Rotation Center X = " << RotationCenterX  << std::endl;
  std::cout << " Rotation Center Y = " << RotationCenterY  << std::endl;
  std::cout << " Rotation Center Z = " << RotationCenterZ  << std::endl;
  std::cout << " Translation X = " << TranslationAlongX  << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
  std::cout << " Translation Z = " << TranslationAlongZ  << std::endl;
=======================================================================

=======================================================================
File Name: itkRayCastInterpolateImageFunction.txx
Revision: 1.18
Line Number: 731-738
Original:
  m_CurrentRayPositionInMM[0] =
    RayPosn[0] + 0.5*m_VoxelDimensionInX*(double)m_NumberOfVoxelsInX;

  m_CurrentRayPositionInMM[1] =
    RayPosn[1] + 0.5*m_VoxelDimensionInY*(double)m_NumberOfVoxelsInY;

  m_CurrentRayPositionInMM[2] =
RayPosn[2] + 0.5*m_VoxelDimensionInZ*(double)m_NumberOfVoxelsInZ;
Replaced by:
  m_CurrentRayPositionInMM[0] = RayPosn[0];

  m_CurrentRayPositionInMM[1] = RayPosn[1];

  m_CurrentRayPositionInMM[2] = RayPosn[2];
=======================================================================

The following are the problems that still troubles me.

1)
Even though I made the change of the code in line 723 as below
  // imageReader2D->GetOutput()->SetOrigin( origin2D );
  caster2D->GetOutput()->SetOrigin( origin2D );
The customized origin of the fixed image was reset to (0, 0, 0) when 
preparing the multiresolution pyramid, i.e., after the execution of
  this->PreparePyramids(); //(line 230 in 
itkMultiResolutionImageRegistrationMethod.txx)
Please enlighten me how to prevent the origin being reset. Thanks!

2)
There is a problem with the metric--GradientDifferenceImageToImageMetric. 
The member function GetDerivative() always returned zero. (See lines 461-484 
in itkGradientDifferenceImageToImageMetric.txx). It appears the variables 
valuep1 and valuep0 always have identical values no matter how the 
transformation parameters were changed. I guess the outputs from the 
TransformMovingImageFilter and the MovedSobelFilters had not been properly 
updated. I could not figure out how to correct it. This problem caused the 
registration stopped at the first iteration and went nowhere.

Jian Wu
Nuclear and Radiological Engineering
University of Florida