[Insight-users] registration question

[Ken Urish]

I have the mattes multimodality registration working great. The
checkerboard shows great alignment.

However, when I take the output images from the registration (ie the
moving images that are now aligned with teh fixed images) and do a
checkerboard with the 'fixed' instead of getting identical output as
the registration program its shifted by about 10 pixels and not
aligned???? The code is heavily based on the mattes registration
example.


In the registration, the checkerboard is pulled from the same pipeline
as the output images (immediately after resampling). In the 2nd
program that overlays the output folder and the 'fxed' folder I simply
put both straight into a checker board. I cant figure out why they are
not aligning up?

I attached code below.


This is for the checkerboard::

//----------------------------------------------------------------------------------
// DICOM SERIES READER
//----------------------------------------------------------------------------------
  typedef itk::ImageSeriesReader< FixedImageType > FixedReaderType;
  typedef itk::ImageSeriesReader< MovingImageType > MovingReaderType;
  typedef itk::GDCMSeriesFileNames NamesGeneratorType;
  typedef itk::GDCMImageIO  ImageIOType;

  ImageIOType::Pointer gdcmIO = ImageIOType::New();
  FixedReaderType::Pointer fixedReader  = FixedReaderType::New();
  MovingReaderType::Pointer movingReader  = MovingReaderType::New();
  NamesGeneratorType::Pointer fixedNamesGenerator = NamesGeneratorType::New();
  NamesGeneratorType::Pointer movingNamesGenerator = NamesGeneratorType::New();

  //reading fixed image files
  fixedNamesGenerator->SetInputDirectory( fixedDirectory );
  const FixedReaderType::FileNamesContainer & fixedFilenames =
fixedNamesGenerator->GetInputFileNames();
  fixedReader->SetImageIO( gdcmIO );
  fixedReader->SetFileNames( fixedFilenames );
  unsigned int numFixedFilenames =  fixedFilenames.size();
  std::cout << numFixedFilenames << std::endl;

  //reading moving image files
  movingNamesGenerator->SetInputDirectory( movingDirectory );
  const MovingReaderType::FileNamesContainer & movingFilenames =
movingNamesGenerator->GetInputFileNames();
  movingReader->SetImageIO( gdcmIO );
  movingReader->SetFileNames( movingFilenames );
  unsigned int numMovingFilenames =  movingFilenames.size();
  std::cout << numMovingFilenames << std::endl;

  //need update so register method can get region size on call
  try { fixedReader->Update(); }
  catch( itk::ExceptionObject & err )
    {
	std::cerr << "ITK says: What? You cant even read an image!" << std::endl;
    std::cerr << err << std::endl;
	getch();
    return EXIT_FAILURE;
    }


  //This is the important section because it is what is adjustng the contrast.
  //The whole reason for the fix.

  //Need to rescale because vtk image reader can only handle unsigned
short greysacle

  //typedef itk::CastImageFilter < FixedImageType, OutputImageType >
InputCastFilterType;

  typedef itk::RescaleIntensityImageFilter< FixedImageType,
OutputImageType >   InputCastFilterType;
  InputCastFilterType::Pointer inputCast = InputCastFilterType::New();
  inputCast->SetInput(fixedReader->GetOutput());
  inputCast->SetOutputMaximum (4000);
  inputCast->SetOutputMinimum (10);
  inputCast->Update();
  //std::cout << "getch: " << std::endl;
  //getch();

  typedef itk::RescaleIntensityImageFilter< FixedImageType,
OutputImageType >   MovingCastFilterType;
  MovingCastFilterType::Pointer movingCast = MovingCastFilterType::New();
  movingCast->SetOutputMaximum (1500);
  movingCast->SetOutputMinimum (10);
  movingCast->SetInput(movingReader->GetOutput());
  movingCast->Update();


//----------------------------------------------------------------------------------
// WRITERS
//----------------------------------------------------------------------------------

//First set up the checker board
typedef itk::FixedArray< unsigned char, 3 > ArrayType;
ArrayType checkerNum;
checkerNum[0] = xCheckerNum;
checkerNum[1] = yCheckerNum;
checkerNum[2] = zCheckerNum;

typedef itk::CheckerBoardImageFilter< FixedImageType > CheckerBoardFilterType;
  CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
  checker->SetInput1( movingCast->GetOutput());
  checker->SetInput2( inputCast->GetOutput() );
  checker->SetCheckerPattern (checkerNum);
  checker-> Update ();


  //Series writer declaration for After and Before shots
  typedef itk::ImageSeriesWriter< OutputImageType, OutputSliceType >
SeriesWriterType;

  SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
  fixedNamesGenerator->SetOutputDirectory(outputDirectory.c_str());
  seriesWriter->SetInput(checker->GetOutput() );
  seriesWriter->SetImageIO( gdcmIO );
  seriesWriter->SetFileNames( fixedNamesGenerator->GetOutputFileNames() );
  seriesWriter->SetMetaDataDictionaryArray(fixedReader->GetMetaDataDictionaryArray()
);

  try { seriesWriter->Update(); }
  catch( itk::ExceptionObject & err )
    {
	std::cerr << "ITK says: You wrong!" << std::endl;
    std::cerr << err << std::endl;
	getch();
    return EXIT_FAILURE;
    }



This is the last part of the registration:


  //Check region dimensions
  outputFile.open(txtFile,std::ios::out | std::ios::app);
  outputFile <<"Region registered: "<<regionRegister<<std::endl;
  std::cout <<"Region registered: "<<regionRegister<<std::endl;
  outputFile.close();

  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetInterpolator(  interpolator  );
  registration->SetFixedImage(    fixedReader->GetOutput()    );
  registration->SetMovingImage(   movingReader->GetOutput()   );
  registration->SetTransform(     transform );
  registration->SetFixedImageRegion(regionRegister);

  // Metric Type
  // One mechanism for bringing the Metric to its limit is to disable the
  // sampling and use all the pixels present in the FixedImageRegion.
  unsigned int numberOfSamples = 10000;
  const unsigned long numberOfImagePixels =
					fixedImagePtr->GetLargestPossibleRegion().GetNumberOfPixels();
  const unsigned long numberOfSpatialSamples =
					static_cast<unsigned long>( numberOfImagePixels * metricPercent );
  metric->SetNumberOfSpatialSamples( numberOfSpatialSamples );
  metric->SetNumberOfHistogramBins( numberOfBins );

  //Transform Initializer
  //Not necessary to call Update() on reader since the
  //CenteredTransformInitializer will do it as part of its computations.
  //If want geometrical centers call: GeometryOn()}
  //If want mass centers call: MomentsOn().
  //Computation of the center and translation is triggered by
InitializeTransform()
  typedef itk::CenteredTransformInitializer< TransformType, FixedImageType,
	                                         MovingImageType>
TransformInitializerType;
  TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
  initializer->SetTransform(   transform );
  initializer->SetFixedImage(  fixedReader->GetOutput() );
  initializer->SetMovingImage( movingReader->GetOutput() );
  initializer->MomentsOn();
  initializer->InitializeTransform();

  //  The rotation part of the transform is initialized using a
versor, which is simply a unit quaternion.
  //  The versor itself defines the type of the vector used to
indicate the rotation axis.
  //  This creates a versor object and initialize its parameters by passing a
  //  rotation axis and an angle.

  typedef TransformType::VersorType  VersorType;
  typedef VersorType::VectorType     VectorType;

  VersorType     rotation;
  VectorType     axis;
  axis[0] = 0.0;
  axis[1] = 0.0;
  axis[2] = 1.0;

  const double angle = 0;
  rotation.Set(  axis, angle  );
  transform->SetRotation( rotation );

  //  We now pass the parameters of the current transform as the initial
  //  parameters to be used when the registration process starts.
  registration->SetInitialTransformParameters( transform->GetParameters() );

  //  Another significant difference in the metric is that it computes the
  //  negative mutual information and hence we need to minimize the cost
  //  function in this case.
  typedef OptimizerType::ScalesType       OptimizerScalesType;
  OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );

  optimizerScales[0] = 1.0;
  optimizerScales[1] = 1.0;
  optimizerScales[2] = 1.0;
  optimizerScales[3] = 1.0 / ( translationScaleFactor*spacing[0]*size[0] );
  optimizerScales[4] = 1.0 / ( translationScaleFactor*spacing[1]*size[1] );
  optimizerScales[5] = 1.0 / ( translationScaleFactor*spacing[2]*size[2] );

  optimizer->SetScales( optimizerScales );
  optimizer->SetMaximumStepLength(  maxStepLength);
  optimizer->SetMinimumStepLength(  minStepLength);
  optimizer->SetNumberOfIterations( iterations);
  optimizer->MaximizeOff();

  // Create the Command observer and register it with the optimizer.
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer);

  //Lets see if we can register
  try { registration->StartRegistration(); }
  catch( itk::ExceptionObject & err )
    {
	std::cerr << "ITK says: You wrong - you will respect my authority!"
<< std::endl;
    std::cerr << err << std::endl;
	getch();
    return EXIT_FAILURE;
    }

  //Opening Observer File
  outputFile.open(txtFile,std::ios::out | std::ios::app);
  outputFile <<"Just finished registeration"<<std::endl;
  outputFile.close();

//----------------------------------------------------------------------------------
// NUMERICAL OUTPUT
//----------------------------------------------------------------------------------
  OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
  const double versorX              = finalParameters[0];
  const double versorY              = finalParameters[1];
  const double versorZ              = finalParameters[2];
  const double finalTranslationX    = finalParameters[3];
  const double finalTranslationY    = finalParameters[4];
  const double finalTranslationZ    = finalParameters[5];

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

  transform->SetParameters( finalParameters );
  TransformType::MatrixType matrix = transform->GetRotationMatrix();
  TransformType::OffsetType offset = transform->GetOffset();


  // Print out results
  std::cout << std::endl << std::endl;
  std::cout << "Result = " << std::endl;
  std::cout << " versor X      = " << versorX  << std::endl;
  std::cout << " versor Y      = " << versorY  << std::endl;
  std::cout << " versor Z      = " << versorZ  << std::endl;
  std::cout << " Translation X = " << finalTranslationX  << std::endl;
  std::cout << " Translation Y = " << finalTranslationY  << std::endl;
  std::cout << " Translation Z = " << finalTranslationZ  << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Best Metric   = " << bestValue          <<
std::endl<< std::endl;
  std::cout << " Matrix = " << std::endl << matrix << std::endl;
  std::cout << " Offset = " << std::endl << offset << std::endl;

  //Opening Observer File
  outputFile.open(txtFile,std::ios::out | std::ios::app);
  outputFile << "Result = " << std::endl;
  outputFile << " versor X      = " << versorX  << std::endl;
  outputFile << " versor Y      = " << versorY  << std::endl;
  outputFile << " versor Z      = " << versorZ  << std::endl;
  outputFile << " Translation X = " << finalTranslationX  << std::endl;
  outputFile << " Translation Y = " << finalTranslationY  << std::endl;
  outputFile << " Translation Z = " << finalTranslationZ  << std::endl;
  outputFile << " Iterations    = " << numberOfIterations << std::endl;
  outputFile << " Best Metric   = " << bestValue          <<
std::endl<< std::endl;
  outputFile << " Matrix = " << std::endl << matrix << std::endl;
  outputFile << " Offset = " << std::endl << offset << std::endl;
  outputFile.close();


//----------------------------------------------------------------------------------
// TRANSFORMS AND RESAMPLING
//----------------------------------------------------------------------------------
  TransformType::Pointer finalTransform = TransformType::New();
  finalTransform->SetCenter( transform->GetCenter() );
  finalTransform->SetParameters( finalParameters );

  typedef itk::ResampleImageFilter<MovingImageType,FixedImageType >
ResampleFilterType;
  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
  FixedImageType::Pointer fixedImage = fixedImagePtr;
  resampler->SetTransform( finalTransform );
  resampler->SetInput( movingImagePtr );
  resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
  resampler->SetOutputOrigin( fixedImage->GetOrigin() );
  resampler->SetOutputSpacing( fixedImage->GetSpacing() );
  resampler->SetDefaultPixelValue( 100 );


//----------------------------------------------------------------------------------
// WRITERS
//----------------------------------------------------------------------------------

  //Series writer declaration for After and Before shots
  typedef itk::CastImageFilter<FixedImageType, OutputImageType > CastFilterType;
  typedef itk::ImageSeriesWriter< OutputImageType, OutputSliceType >
SeriesWriterType;

  CastFilterType::Pointer  caster =  CastFilterType::New();
  caster->SetInput( resampler->GetOutput() );

  SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
  fixedNamesGenerator->SetOutputDirectory(outputDirectory.c_str());
  seriesWriter->SetInput(caster->GetOutput() );
  seriesWriter->SetImageIO( gdcmIO );
  seriesWriter->SetFileNames( fixedNamesGenerator->GetOutputFileNames() );
  seriesWriter->SetMetaDataDictionaryArray(fixedReader->GetMetaDataDictionaryArray()
);

  try { seriesWriter->Update(); }
  catch( itk::ExceptionObject & err )
    {
	std::cerr << "ITK says: You wrong!" << std::endl;
    std::cerr << err << std::endl;
	getch();
    return EXIT_FAILURE;
    }

  // Generate checkerboards before and after registration

  typedef itk::FixedArray< unsigned char, 3 > ArrayType;
  ArrayType checkerNum;
  checkerNum[0] = xCheckerNum;
  checkerNum[1] = yCheckerNum;
  checkerNum[2] = zCheckerNum;

  typedef itk::CheckerBoardImageFilter< FixedImageType > CheckerBoardFilterType;
  CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
  checker->SetInput1( fixedImage );
  checker->SetInput2( caster->GetOutput() );

  //This was the old code.
  //checker->SetInput2( resampler->GetOutput() );
  checker->SetCheckerPattern (checkerNum);

  caster->SetInput( checker->GetOutput() );
  resampler->SetDefaultPixelValue( 0 );

  // Before registration
  TransformType::Pointer identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  resampler->SetTransform( identityTransform );

  fixedNamesGenerator->SetOutputDirectory(beforeDirectory.c_str());
  seriesWriter->SetFileNames( fixedNamesGenerator->GetOutputFileNames() );
  try { seriesWriter->Update(); }
  catch( itk::ExceptionObject & err )
    {
	std::cerr << "ITK says: You wrong!" << std::endl;
    std::cerr << err << std::endl;
	getch();
    return EXIT_FAILURE;
    }

  // After registration
  resampler->SetTransform( finalTransform );
  fixedNamesGenerator->SetOutputDirectory(afterDirectory.c_str());
  seriesWriter->SetFileNames( fixedNamesGenerator->GetOutputFileNames() );

  try { seriesWriter->Update(); }
  catch( itk::ExceptionObject & err )
    {
	std::cerr << "ITK says: You wrong!" << std::endl;
    std::cerr << err << std::endl;
	getch();
    return EXIT_FAILURE;
    }
 	
  //Opening Observer File
  outputFile.open(txtFile,std::ios::out | std::ios::app);
  outputFile <<"Just wrote 'after' image series"<<std::endl;
  outputFile.close();