[Insight-users] Image Registration on Regions of Images

[sg13]

Hello, 

I am a biomedical engineering student at Rensselaer Polytechnic Institute
and I am using the ITK library in my research.  I am trying to perform image
registration on 20x20 regions of the fixed image.  I read in a region, then
move the start point of the region over 10 each time, followed by moving
down by 10.  (using nested for-loops)  I have previous code that is written
in Matlab (without using ITK) and was converting to c++ in hopes of speeding
up the run time.  After running the code I wrote however, it takes about the
same amount of time to run.  My code uses the example image registration
code as a basis, and then I modified it to my specific needs, thus I was
expecting results.  My c++ code is below and any help would be greatly
appreciated!

Thanks!
Samantha


#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkNormalizedCorrelationImageToImageMetric.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImage.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"

#include <windows.h>

class CommandIterationUpdate : public itk::Command 
{
public:
  typedef  CommandIterationUpdate   Self;
  typedef  itk::Command             Superclass;
  typedef itk::SmartPointer<Self>  Pointer;
  itkNewMacro(Self);

protected:
  CommandIterationUpdate() {};

public:

  typedef itk::RegularStepGradientDescentOptimizer     OptimizerType;
  typedef const OptimizerType                         *OptimizerPointer;

  void Execute(itk::Object *caller, const itk::EventObject & event)
  {
    Execute( (const itk::Object *)caller, event);
  }

  void Execute(const itk::Object * object, const itk::EventObject & event)
  {
    OptimizerPointer optimizer = 
                         dynamic_cast< OptimizerPointer >( object );

    if( ! itk::IterationEvent().CheckEvent( &event ) )
      {
      return;
      }

    std::cout << optimizer->GetCurrentIteration() << " = ";
    std::cout << optimizer->GetValue() << " : ";
    std::cout << optimizer->GetCurrentPosition() << std::endl;
  }
   
};


int main( int argc, char *argv[] )
{
  float initialTime = GetTickCount();
  //Check Number of Input Parameters
  if( argc != 5 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile ";
    std::cerr << "outputImagefile [differenceImageAfter]";
    std::cerr << "[differenceImageBefore]" << std::endl;
    return EXIT_FAILURE;
    }
  
  //Initializations
  const unsigned int Dimension = 2;
  typedef float PixelType;
 
  typedef itk::Image<PixelType, Dimension>  FixedImageType;
  typedef itk::Image<PixelType, Dimension>  MovingImageType;
 
  //Maps original image onto deformed image
  typedef itk::TranslationTransform<double, Dimension> TransformType;
  typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
  typedef itk::NormalizedCorrelationImageToImageMetric<FixedImageType,
MovingImageType> MetricType;
  //Performs BSpline Interpolation
  typedef itk::BSplineInterpolateImageFunction<MovingImageType, double>
InterpolatorType;
  typedef itk::ImageRegistrationMethod<FixedImageType, MovingImageType>
RegistrationType;
                                    
  MetricType::Pointer metric = MetricType::New();
  TransformType::Pointer transform = TransformType::New();
  OptimizerType::Pointer optimizer = OptimizerType::New();
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  RegistrationType::Pointer registration = RegistrationType::New();

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

  //Read in Original and Deformed Image
  typedef itk::ImageFileReader<FixedImageType> FixedImageReaderType;
  typedef itk::ImageFileReader<MovingImageType> MovingImageReaderType;
  FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();

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

  typedef itk::RegionOfInterestImageFilter <FixedImageType, FixedImageType >
RegionFilterType;
  RegionFilterType::Pointer getRegion = RegionFilterType::New();
  typedef RegistrationType::ParametersType ParametersType;
  FixedImageType::RegionType region;
  FixedImageType::SizeType size;
  FixedImageType::IndexType startPoint;
  size[0] = 20;
  size[1] = 20;
  region.SetSize(size);
  unsigned int xStart = 20;
  unsigned int yStart = 20;
  unsigned int xItr = 20;
  unsigned int yItr = 20;

  //Equivalent of GridStep in the Matlab code
  optimizer->SetMaximumStepLength(5);  
  std::cout << "Max Step Length = 10" << std::endl;
  optimizer->SetMinimumStepLength(0.0001);  
 
  //Maximum Number of Iterations in case optimizer doesn't reach above
precision
  optimizer->SetNumberOfIterations(1000);
		 
  //Set the Initial Parameters
  ParametersType initialParameters(transform->GetNumberOfParameters());
  initialParameters[0] = 6;  // Initial offset in mm along X
  initialParameters[1] = -7;  // Initial offset in mm along Y
 
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();

  ParametersType finalParameters;
  double xTrans, yTrans;
  unsigned int numItr;
  double bestVal;

  for(unsigned int x=0; x != xItr; x++)
  {
	  startPoint[0] = xStart + (x*10);
	  for(unsigned int y=0; y != yItr; y++)
	  {
		  startPoint[1] = yStart + (y*10);
		  region.SetIndex(startPoint);

		  getRegion->SetRegionOfInterest(region);
		  getRegion->SetInput(fixedImageReader->GetOutput());
		  getRegion->Update();

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

		  //Set Image Region as only the region currently loaded into memory
		  //In our case, the whole image for right now
		 
registration->SetFixedImageRegion(getRegion->GetOutput()->GetBufferedRegion());
		  
		  //Set the x and y transform parameters
		  registration->SetInitialTransformParameters(initialParameters);	 

		  optimizer->AddObserver(itk::IterationEvent(), observer);

		  try 
			{ 
			registration->Update(); 
			} 
		  catch( itk::ExceptionObject & err ) 
			{ 
			std::cerr << "ExceptionObject caught !" << std::endl; 
			std::cerr << err << std::endl; 
			return EXIT_FAILURE;
			} 
		  
		  finalParameters = registration->GetLastTransformParameters();
		  xTrans = finalParameters[0];
		  yTrans = finalParameters[1];
		  numItr = optimizer->GetCurrentIteration();
		  bestVal = optimizer->GetValue();

	  }
  }

  typedef itk::ResampleImageFilter<MovingImageType, FixedImageType>   
ResampleFilterType;
 
  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
  resampler->SetInput(movingImageReader->GetOutput());
 
  resampler->SetTransform(registration->GetOutput()->Get());
 
  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
  resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resampler->SetOutputOrigin(fixedImage->GetOrigin());
  resampler->SetOutputSpacing(fixedImage->GetSpacing());
  resampler->SetOutputDirection(fixedImage->GetDirection());
  resampler->SetDefaultPixelValue(100);
 
  typedef unsigned char OutputPixelType;
  typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
  typedef itk::CastImageFilter<FixedImageType, OutputImageType>
CastFilterType;
  typedef itk::ImageFileWriter<OutputImageType>  WriterType;
 
  WriterType::Pointer writer =  WriterType::New();
  CastFilterType::Pointer caster =  CastFilterType::New();

  writer->SetFileName(argv[3]);

  caster->SetInput(resampler->GetOutput());
  writer->SetInput(caster->GetOutput());
  writer->Update();
 
  typedef itk::SubtractImageFilter<FixedImageType, FixedImageType,
FixedImageType> DifferenceFilterType;

  DifferenceFilterType::Pointer difference = DifferenceFilterType::New();

  difference->SetInput1(fixedImageReader->GetOutput());
  difference->SetInput2(resampler->GetOutput());
 
  typedef itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>  
RescalerType;

  RescalerType::Pointer intensityRescaler = RescalerType::New();
  
  intensityRescaler->SetInput(difference->GetOutput());
  intensityRescaler->SetOutputMinimum(0);
  intensityRescaler->SetOutputMaximum(255);

  resampler->SetDefaultPixelValue(1);
 
  WriterType::Pointer writer2 = WriterType::New();
  writer2->SetInput(intensityRescaler->GetOutput());  

  if(argc > 4)
    {
    writer2->SetFileName(argv[4]);
    writer2->Update();
    }

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

  if(argc > 5)
    {
    writer2->SetFileName(argv[5]);
    writer2->Update();
    }
  
  float finalTime = GetTickCount();
  //calculate total run time in minutes
  float runTime = (finalTime-initialTime)*.001;
  std::cout << " Run Time: " << runTime << " seconds" << std::endl;
  std::getchar();

  return EXIT_SUCCESS;
}

-- 
View this message in context: http://old.nabble.com/Image-Registration-on-Regions-of-Images-tp31939711p31939711.html
Sent from the ITK - Users mailing list archive at Nabble.com.