[Insight-users] Image Piece Overlap Registration
Luis Ibanez
luis.ibanez at kitware.com
Thu Jul 20 19:50:55 EDT 2006
Hi Kevin,
Then Normalized correlation metric uses a good deal of memory
for storing the gradient of the moving image. (As do most of the
other image metrics in ITK).
Depending on the optimizer that you use, this gradient may or may not
be necessary.
In your specific case, the RegularStepGradientDescentOptimizer, does
take advantage of this gradient, so you are actually using it during
the registration.
If memory is an issue, you may want to disable the computation of the
gradient and use a different optimizer. For example an optimizer like
the AmoebaOptimizer will not require the metric to compute its
derivatives.
Note also that in your case you should crop your images before proceeding
to the registration. This is because the overlapping region between the
images is the only portion that is useful for the purpose of the
registration.
Let's say that your images are 2048 x 248 each, and that when
you overlap them correctly they share a rectangular section of 20x2048.
That section of 20x2048 is the only one that you want to use for the
registration. The other non-overlapping regions are an unnecessary burden
for the metric.
So, to summarize, use the itkRegionOfInterestImageFilter in order to
extract the useful sections of the fixed and moving images, and then
use only those sections as input fixed and input moving images for the
registration framework.
Regards,
Luis
-----------------------
Kevin H. Hobbs wrote:
> On Mon, 2006-07-17 at 11:44 -0400, Luis Ibanez wrote:
>
>> Since your images are of the same
>>modality you probably should try the NormalizedCorrelation image Metric
>>first.
>
>
> Does this metric use a bunch of doubles? It seams to use a lot of RAM.
>
> Other than that I have something working. The code is attached.
>
>
> ------------------------------------------------------------------------
>
> #include "itkImage.h"
>
> #include "itkImageFileReader.h"
>
> #include "itkRegionOfInterestImageFilter.h"
>
> #include "itkTranslationTransform.h"
>
> #include "itkImageRegistrationMethod.h"
>
> #include "itkNormalizedCorrelationImageToImageMetric.h"
> #include "itkMeanSquaresImageToImageMetric.h"
>
> #include "itkLinearInterpolateImageFunction.h"
>
> #include "itkRegularStepGradientDescentOptimizer.h"
>
> #include "itkResampleImageFilter.h"
>
> #include "itkSubtractImageFilter.h"
>
> #include "itkImageFileWriter.h"
>
> #include "itkCommand.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::cout << optimizer->GetCurrentStepLength() << std::endl;
> }
> };
>
>
> int main( int argc, char *argv[] )
> {
>
> if( argc < 12 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0] << std::endl;
> std::cerr << " fixed_file moving_file " << std::endl;
> std::cerr << " out_file max_step min_step" << std::endl;
> std::cerr << " scale_x scale_y scale_z" << std::endl;
> std::cerr << " trans_x trans_y trans_z " << std::endl;
> return EXIT_FAILURE;
> }
>
> char * fixed_file = argv[1];
> char * moving_file = argv[2];
> char * out_file = argv[3];
> double max_step = atof( argv[4] );
> double min_step = atof( argv[5] );
> double scale_x = atof( argv[6] );
> double scale_y = atof( argv[7] );
> double scale_z = atof( argv[8] );
> double trans_x = atof( argv[9] );
> double trans_y = atof( argv[10] );
> double trans_z = atof( argv[11] );
>
> const unsigned int Dimension = 3;
> typedef unsigned char PixelType;
>
> typedef itk::Image< PixelType, Dimension > ImageType;
>
> typedef itk::TranslationTransform
> < double, Dimension> TransformType;
>
> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>
>
> typedef itk::NormalizedCorrelationImageToImageMetric
> < ImageType, ImageType > MetricType;
> /*
> typedef itk::MeanSquaresImageToImageMetric
> < ImageType, ImageType > MetricType;
> */
>
> typedef itk:: LinearInterpolateImageFunction
> < ImageType, double > InterpolatorType;
>
> typedef itk::ImageRegistrationMethod
> < ImageType, ImageType > RegistrationType;
>
> MetricType::Pointer metric = MetricType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
>
> TransformType::Pointer transform = TransformType::New();
> registration->SetTransform( transform );
>
> typedef itk::ImageFileReader< ImageType > ImageReaderType;
> ImageReaderType::Pointer fixedImageReader = ImageReaderType::New();
> ImageReaderType::Pointer movingImageReader = ImageReaderType::New();
>
> fixedImageReader->SetFileName( fixed_file );
> movingImageReader->SetFileName( moving_file );
>
> fixedImageReader->Update();
>
> ImageType::RegionType fixed_region
> = fixedImageReader->GetOutput()->GetLargestPossibleRegion();
> ImageType::SizeType fixed_size = fixed_region.GetSize();
> ImageType::IndexType fixed_index = fixed_region.GetIndex();
> ImageType::SpacingType fixed_spacing
> = fixedImageReader->GetOutput()->GetSpacing();
> ImageType::PointType fixed_origin
> = fixedImageReader->GetOutput()->GetOrigin();
>
> fixed_size[0] -= (long unsigned int)( fabs( trans_x ) / fixed_spacing[0] );
> fixed_size[1] -= (long unsigned int)( fabs( trans_y ) / fixed_spacing[1] );
> fixed_size[2] -= (long unsigned int)( fabs( trans_z ) / fixed_spacing[2] );
>
> std::cout << "Fixed x size = " << fixed_size[0] << std::endl;
> std::cout << "Fixed y size = " << fixed_size[1] << std::endl;
> std::cout << "Fixed z size = " << fixed_size[2] << std::endl;
>
> if ( trans_x > 0.0 )
> fixed_index[0] += (long int)( trans_x / fixed_spacing[0] );
> if ( trans_y > 0.0 )
> fixed_index[1] += (long int)( trans_y / fixed_spacing[1] );
> if ( trans_z > 0.0 )
> fixed_index[2] += (long int)( trans_z / fixed_spacing[2] );
>
> std::cout << "Fixed x start = " << fixed_index[0] << std::endl;
> std::cout << "Fixed y start = " << fixed_index[1] << std::endl;
> std::cout << "Fixed z start = " << fixed_index[2] << std::endl;
>
> fixed_region.SetSize( fixed_size );
> fixed_region.SetIndex( fixed_index );
>
> typedef itk::RegionOfInterestImageFilter
> < ImageType, ImageType > ROIFilterType;
> ROIFilterType::Pointer fixed_roi_filter = ROIFilterType::New();
> fixed_roi_filter->SetRegionOfInterest( fixed_region );
> fixed_roi_filter->SetInput( fixedImageReader->GetOutput() );
> fixed_roi_filter->Update();
> ImageType::Pointer fixed_roi = fixed_roi_filter->GetOutput();
>
> std::cout << "Deleting used fixed filters." << std::endl;
> fixedImageReader->Delete();
> fixed_roi_filter->Delete();
>
> movingImageReader->Update();
>
> ImageType::RegionType moving_region
> = movingImageReader->GetOutput()->GetLargestPossibleRegion();
> ImageType::SizeType moving_size = moving_region.GetSize();
> ImageType::IndexType moving_index = moving_region.GetIndex();
> ImageType::SpacingType moving_spacing
> = movingImageReader->GetOutput()->GetSpacing();
> ImageType::PointType moving_origin
> = movingImageReader->GetOutput()->GetOrigin();
>
> moving_size[0] -= (long unsigned int)( fabs( trans_x ) / moving_spacing[0] );
> moving_size[1] -= (long unsigned int)( fabs( trans_y ) / moving_spacing[1] );
> moving_size[2] -= (long unsigned int)( fabs( trans_z ) / moving_spacing[2] );
>
> std::cout << "Moving x size = " << moving_size[0] << std::endl;
> std::cout << "Moving y size = " << moving_size[1] << std::endl;
> std::cout << "Moving z size = " << moving_size[2] << std::endl;
>
> if ( trans_x < 0.0 )
> moving_index[0] -= (long int)( trans_x / moving_spacing[0] );
> if ( trans_y < 0.0 )
> moving_index[1] -= (long int)( trans_y / moving_spacing[1] );
> if ( trans_z < 0.0 )
> moving_index[2] -= (long int)( trans_z / moving_spacing[2] );
>
> std::cout << "Moving x start = " << moving_index[0] << std::endl;
> std::cout << "Moving y start = " << moving_index[1] << std::endl;
> std::cout << "Moving z start = " << moving_index[2] << std::endl;
>
> moving_region.SetSize( moving_size );
> moving_region.SetIndex( moving_index );
>
>
> ROIFilterType::Pointer moving_roi_filter = ROIFilterType::New();
> moving_roi_filter->SetRegionOfInterest( moving_region );
> moving_roi_filter->SetInput( movingImageReader->GetOutput() );
> moving_roi_filter->Update();
> ImageType::Pointer moving_roi = moving_roi_filter->GetOutput();
>
> std::cout << "Deleting used moving filters." << std::endl;
> movingImageReader->Delete();
> moving_roi_filter->Delete();
>
>
> std::cout << "Writing the ROIs." << std::endl;
> typedef itk::ImageFileWriter< ImageType > ROIWriterType;
> ROIWriterType::Pointer roi_f_writer = ROIWriterType::New();
> ROIWriterType::Pointer roi_m_writer = ROIWriterType::New();
> roi_f_writer->SetFileName( "f_roi.vtk" );
> roi_m_writer->SetFileName( "m_roi.vtk" );
> roi_f_writer->SetInput( fixed_roi );
> roi_m_writer->SetInput( moving_roi );
> roi_f_writer->Update();
> roi_m_writer->Update();
> roi_f_writer->Delete();
> roi_m_writer->Delete();
>
>
> registration->SetFixedImageRegion
> ( fixed_roi->GetBufferedRegion() );
>
> registration->SetFixedImage( fixed_roi );
> registration->SetMovingImage( moving_roi );
>
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales( Dimension );
> optimizerScales[0] = scale_x;
> optimizerScales[1] = scale_y;
> optimizerScales[2] = scale_z;
>
> optimizer->SetScales( optimizerScales );
>
> typedef RegistrationType::ParametersType ParametersType;
> ParametersType initialParameters( transform->GetNumberOfParameters() );
>
> initialParameters[0] = -trans_x;
> initialParameters[1] = -trans_y;
> initialParameters[2] = -trans_z;
>
> registration->SetInitialTransformParameters( initialParameters );
>
> optimizer->SetRelaxationFactor( 0.6 );
> optimizer->SetMaximumStepLength( max_step );
> optimizer->SetMinimumStepLength( min_step );
> optimizer->SetNumberOfIterations( 2000 );
>
>
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return EXIT_FAILURE;
> }
>
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
> double final_trans[Dimension];
> final_trans[0] = -finalParameters[0];
> final_trans[1] = -finalParameters[1];
> final_trans[2] = -finalParameters[2];
>
> std::cout << " Translation X = " << final_trans[0] << std::endl;
> std::cout << " Translation Y = " << final_trans[1] << std::endl;
> std::cout << " Translation Z = " << final_trans[2] << std::endl;
>
>
> ImageReaderType::Pointer movingImageReReader
> = ImageReaderType::New();
> movingImageReReader->SetFileName( moving_file );
> movingImageReReader->Update();
>
> ImageType::Pointer moving_image
> = movingImageReReader->GetOutput();
> moving_image->DisconnectPipeline();
> moving_image->SetOrigin( final_trans );
>
> typedef itk::ImageFileWriter< ImageType > MovingWriterType;
> MovingWriterType::Pointer moving_writer = MovingWriterType::New();
> moving_writer->SetFileName( out_file );
> moving_writer->SetInput( moving_image );
>
> moving_writer->Update();
>
>
> return 0;
> }
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list