[Insight-users] Image Piece Overlap Registration
Luis Ibanez
luis.ibanez at kitware.com
Sat Jul 29 18:38:33 EDT 2006
Hi Kevin,
> Kevin H. Hobbs wrote:
> I'll consider an IJ submission.
Please do so !!
This is the kind of contribution that is a perfect fit for the IJ.
It is code that will certainly be useful to many other users.
Thanks
Luis
-----------------------------------------------------------
Kevin H. Hobbs wrote:
> On Mon, 2006-07-17 at 11:44 -0400, Luis Ibanez wrote:
>
>>We will be interested to hear about your findings. BTW that will be
>>a nice submission to the Insight Journal.
>
>
> The attached two programs seem to work for my project.
>
> "align.cxx" takes two images with origin (0,0,0). It sets the correct
> origin for the second image starting from a user supplied guess. It uses
> only the regions of the two images that overlap with the guess.
>
> "ImageMaxStitch.cxx" has a streaming VTK to ITK to VTK pipeline that
> combines two images into one. The streaming may not be working
> correctly.
>
> I'll consider an IJ submission.
>
>
> ------------------------------------------------------------------------
>
> #include "itkPoint.h"
> #include "itkImage.h"
> #include "vtkImageData.h"
>
> #include "vtkITKUtility.h"
>
> #include "vtkXMLImageDataReader.h"
> #include "vtkImageClip.h"
> #include "itkImageFileReader.h"
>
> #include "vtkImageExport.h"
> #include "itkVTKImageImport.h"
>
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkTranslationTransform.h"
> #include "itkResampleImageFilter.h"
> #include "itkMaximumImageFilter.h"
>
> #include "itkVTKImageExport.h"
> #include "vtkImageImport.h"
>
> #include "vtkExtentTranslator.h"
> #include "vtkIntArray.h"
> #include "vtkIdTypeArray.h"
> #include "vtkSortDataArray.h"
> #include "vtkTableExtentTranslator.h"
>
> #include "vtkXMLImageDataWriter.h"
> #include "itkImageFileWriter.h"
>
> int main( int argc, char *argv[] )
> {
>
> if( argc < 5 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0] << std::endl;
> std::cerr << " in_file1 " << std::endl;
> std::cerr << " in_file2 " << std::endl;
> std::cerr << " out_file " << std::endl;
> std::cerr << " pieces " << std::endl;
> return EXIT_FAILURE;
> }
>
> char * in_file1 = argv[1];
> char * in_file2 = argv[2];
> char * out_file = argv[3];
> int pieces = atoi( argv[4] );
>
> // ITK Image Type
> const unsigned int Dimension = 3;
> typedef unsigned char PixelType;
> typedef itk::Image< PixelType, Dimension > ImageType;
>
> // ITK Region Type
> typedef ImageType::RegionType RegionType;
>
> // ITK Import Type
> typedef itk::VTKImageImport < ImageType > ITKImportType;
>
> // ITK Transform Type
> typedef itk::TranslationTransform
> < double, Dimension > TransformType;
>
> // ITK Resample Type
> typedef itk::ResampleImageFilter
> < ImageType, ImageType > ResampleFilterType;
>
> // ITK Max Type
> typedef itk::MaximumImageFilter
> < ImageType, ImageType, ImageType > MaxFilterType;
>
> typedef itk::VTKImageExport< ImageType > ITKExportType;
>
> // VTK Reader 1
> vtkXMLImageDataReader * image1_reader
> = vtkXMLImageDataReader::New();
> image1_reader->ReleaseDataFlagOn();
> image1_reader->SetFileName( in_file1 );
>
> // VTK Clip 1 Just to update origin and spacing.
> vtkImageClip * clip1 = vtkImageClip::New();
> clip1->SetOutputWholeExtent(0,1,0,1,0,1);
> clip1->SetInputConnection( image1_reader->GetOutputPort() );
> clip1->Update();
>
> // VTK Export 1
> vtkImageExport * vtk_export1 = vtkImageExport::New();
> vtk_export1->SetInputConnection
> ( image1_reader->GetOutputPort() );
>
> // ITK Import 1
> ITKImportType::Pointer itk_import1 = ITKImportType::New();
> ConnectPipelines( vtk_export1, itk_import1 );
>
> // VTK Reader 2
> vtkXMLImageDataReader * image2_reader
> = vtkXMLImageDataReader::New();
> image2_reader->ReleaseDataFlagOn();
> image2_reader->SetFileName( in_file2 );
>
> // VTK Clip 2 Just to update origin and spacing.
> vtkImageClip * clip2 = vtkImageClip::New();
> clip2->SetOutputWholeExtent(0,1,0,1,0,1);
> clip2->SetInputConnection( image2_reader->GetOutputPort() );
> clip2->Update();
>
> // VTK Export 2
> vtkImageExport * vtk_export2 = vtkImageExport::New();
> vtk_export2->SetInputConnection
> ( image2_reader->GetOutputPort() );
>
> // ITK Import 2
> ITKImportType::Pointer itk_import2 = ITKImportType::New();
> ConnectPipelines( vtk_export2, itk_import2 );
>
> // Image 1 Information
> vtkImageData * image1 = image1_reader->GetOutput();
> //image1->Print(cout);
> double * image1_spacing = image1->GetSpacing();
> double * image1_origin = image1->GetOrigin();
> int * whole_extent1 = image1->GetWholeExtent();
> int image1_dim[3];
> image1_dim[0] = whole_extent1[1] + 1;
> image1_dim[1] = whole_extent1[3] + 1;
> image1_dim[2] = whole_extent1[5] + 1;
> ImageType::PointType image1_corner;
> image1_corner[0]
> = image1_origin[0] + image1_dim[0] * image1_spacing[0];
> image1_corner[1]
> = image1_origin[1] + image1_dim[1] * image1_spacing[1];
> image1_corner[2]
> = image1_origin[2] + image1_dim[2] * image1_spacing[2];
>
> // Image 2 Information
> vtkImageData * image2 = image2_reader->GetOutput();
> //image2->Print(cout);
> double * image2_spacing = image2->GetSpacing();
> double * image2_origin = image2->GetOrigin();
> int * whole_extent2 = image2->GetWholeExtent();
> int image2_dim[3];
> image2_dim[0] = whole_extent2[1] + 1;
> image2_dim[1] = whole_extent2[3] + 1;
> image2_dim[2] = whole_extent2[5] + 1;
> ImageType::PointType image2_corner;
> image2_corner[0]
> = image2_origin[0] + image2_dim[0] * image2_spacing[0];
> image2_corner[1]
> = image2_origin[1] + image2_dim[1] * image2_spacing[1];
> image2_corner[2]
> = image2_origin[2] + image2_dim[2] * image2_spacing[2];
>
> // Output Image Information
> ImageType::SpacingType output_spacing;
> output_spacing[0] = fmin( image1_spacing[0], image2_spacing[0] );
> output_spacing[1] = fmin( image1_spacing[1], image2_spacing[1] );
> output_spacing[2] = fmin( image1_spacing[2], image2_spacing[2] );
>
> ImageType::PointType output_origin;
> output_origin[0] = fmin( image1_origin[0], image2_origin[0] );
> output_origin[1] = fmin( image1_origin[1], image2_origin[1] );
> output_origin[2] = fmin( image1_origin[2], image2_origin[2] );
>
> ImageType::PointType output_corner;
> output_corner[0] = fmax( image1_corner[0], image2_corner[0] );
> output_corner[1] = fmax( image1_corner[1], image2_corner[1] );
> output_corner[2] = fmax( image1_corner[2], image2_corner[2] );
>
> ImageType::SizeType output_dim;
> output_dim[0] = (int)( ( output_corner[0] - output_origin[0] )
> / output_spacing[0] );
> output_dim[1] = (int)( ( output_corner[1] - output_origin[1] )
> / output_spacing[0] );
> output_dim[2] = (int)( ( output_corner[2] - output_origin[2] )
> / output_spacing[2] );
>
> // ITK Transform
> TransformType::Pointer transform = TransformType::New();
> transform->SetIdentity();
>
> // ITK Resample 1
> ResampleFilterType::Pointer resample1
> = ResampleFilterType::New();
> resample1->SetTransform( transform );
> resample1->SetSize( output_dim );
> resample1->SetOutputOrigin( output_origin );
> resample1->SetOutputSpacing( output_spacing );
> resample1->SetDefaultPixelValue( 0 );
> resample1->SetInput( itk_import1->GetOutput() );
>
> // ITK Resample 2
> ResampleFilterType::Pointer resample2
> = ResampleFilterType::New();
> resample2->SetTransform( transform );
> resample2->SetSize( output_dim );
> resample2->SetOutputOrigin( output_origin );
> resample2->SetOutputSpacing( output_spacing );
> resample2->SetDefaultPixelValue( 0 );
> resample2->SetInput( itk_import2->GetOutput() );
>
> // ITK Max
> MaxFilterType::Pointer max_filter = MaxFilterType::New();
> max_filter->SetInput1( resample1->GetOutput() );
> max_filter->SetInput2( resample2->GetOutput() );
>
> /*
> // ITK Writer
> typedef itk::ImageFileWriter< ImageType > itkWriterType;
> itkWriterType::Pointer itk_writer = itkWriterType::New();
> itk_writer->SetFileName( out_file );
> itk_writer->SetInput( max_filter->GetOutput() );
> itk_writer->Update();
> */
>
> // ITK Export
> ITKExportType::Pointer itk_export = ITKExportType::New();
> itk_export->SetInput( max_filter->GetOutput() );
>
> // VTK Import
> vtkImageImport * vtk_import = vtkImageImport::New();
> ConnectPipelines( itk_export, vtk_import );
>
> // Default Extent Translator
> vtkExtentTranslator * default_translator
> = vtkExtentTranslator::New();
>
> // Table Extent Translator
> vtkTableExtentTranslator * table_translator
> = vtkTableExtentTranslator::New();
> table_translator->SetNumberOfPieces( pieces );
> table_translator->SetNumberOfPiecesInTable( pieces );
>
> vtkIdType piece;
> int piece_extent[6];
> vtkIdType key;
> int whole_extent[6];
> whole_extent[0] = 0;
> whole_extent[1] = output_dim[0] - 1;
> whole_extent[2] = 0;
> whole_extent[3] = output_dim[1] - 1;
> whole_extent[4] = 0;
> whole_extent[5] = output_dim[2] - 1;
>
> vtkIdTypeArray * keys = vtkIdTypeArray::New();
> keys->SetNumberOfComponents( 1 );
> keys->SetNumberOfTuples( pieces );
>
> vtkIntArray * extents = vtkIntArray::New();
> extents->SetNumberOfComponents( 6 );
> extents->SetNumberOfTuples( pieces );
>
> for( piece = 0; piece < pieces; ++piece )
> {
> default_translator->PieceToExtentThreadSafe(
> piece, pieces,
> 0, whole_extent, piece_extent,
> 3, 1);
> key = piece_extent[4]
> * whole_extent[1] * whole_extent[3]
> + piece_extent[2] * whole_extent[1]
> + piece_extent[0];
> keys->SetTupleValue( piece, &key );
> extents->SetTupleValue( piece, piece_extent );
> }
>
> vtkSortDataArray::Sort( keys, extents );
>
> for( piece = 0; piece < pieces; ++piece )
> {
>
> extents->GetTupleValue( piece, piece_extent );
> table_translator->SetExtentForPiece(
> piece, piece_extent);
> }
>
> // VTK Image Writer
> vtkXMLImageDataWriter * writer
> = vtkXMLImageDataWriter::New();
> writer->SetFileName( out_file );
> writer->SetNumberOfPieces( pieces );
> writer->SetExtentTranslator( table_translator );
> writer->SetInputConnection( vtk_import->GetOutputPort() );
>
> writer->Update();
>
>
> return 0;
> }
>
>
> ------------------------------------------------------------------------
>
> #include "itkImage.h"
>
> #include "itkImageFileReader.h"
>
> #include "itkRegionOfInterestImageFilter.h"
>
> #include "itkTranslationTransform.h"
>
> #include "itkImageRegistrationMethod.h"
>
> #include "itkNormalizedCorrelationImageToImageMetric.h"
>
> #include "itkLinearInterpolateImageFunction.h"
>
> #include "itkAmoebaOptimizer.h"
>
> #include "itkImageFileWriter.h"
>
> #include "itkCommand.h"
>
>
> class CommandIterationUpdateAmoeba : public itk::Command
> {
> public:
> typedef CommandIterationUpdateAmoeba Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
> protected:
> CommandIterationUpdateAmoeba()
> {
> m_IterationNumber=0;
> }
> public:
> typedef itk::AmoebaOptimizer 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::FunctionEvaluationIterationEvent().CheckEvent( &event ) )
> {
> std::cout << m_IterationNumber++ << " ";
> std::cout << optimizer->GetCachedValue() << " ";
> std::cout << optimizer->GetCachedCurrentPosition() << std::endl;
> }
> }
> private:
> unsigned long m_IterationNumber;
> };
>
>
> int main( int argc, char *argv[] )
> {
>
> if( argc < 9 )
> {
> 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 init_step min_step" << 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 init_step = atof( argv[4] );
> double min_step = atof( argv[5] );
> double trans_x = atof( argv[6] );
> double trans_y = atof( argv[7] );
> double trans_z = atof( argv[8] );
>
> const unsigned int Dimension = 3;
> typedef unsigned char PixelType;
>
> typedef itk::Image< PixelType, Dimension > ImageType;
>
> typedef itk::TranslationTransform
> < double, Dimension> TransformType;
>
> typedef itk::AmoebaOptimizer OptimizerType;
>
> typedef itk::NormalizedCorrelationImageToImageMetric
> < 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();
>
> metric->ComputeGradientOff();
>
> 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;
>
> fixed_roi_filter = 0; // Delete
> fixedImageReader = 0; // 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;
>
> moving_roi_filter = 0; // Delete
> movingImageReader = 0; // Delete
>
>
> registration->SetFixedImageRegion
> ( fixed_roi->GetBufferedRegion() );
>
> registration->SetFixedImage( fixed_roi );
> registration->SetMovingImage( moving_roi );
>
> typedef RegistrationType::ParametersType ParametersType;
> ParametersType initialParameters
> ( transform->GetNumberOfParameters() );
> initialParameters[0] = -trans_x;
> initialParameters[1] = -trans_y;
> initialParameters[2] = -trans_z;
> registration->SetInitialTransformParameters( initialParameters );
>
> typedef OptimizerType::ParametersType AmoebaParametersType;
> ParametersType initialStep
> ( transform->GetNumberOfParameters() );
> initialParameters[0] = init_step;
> initialParameters[1] = init_step;
> initialParameters[2] = init_step;
> optimizer->SetInitialSimplexDelta( initialParameters );
>
> optimizer->SetParametersConvergenceTolerance( min_step );
> optimizer->AutomaticInitialSimplexOff();
> optimizer->SetMaximumNumberOfIterations( 2000 );
>
>
> CommandIterationUpdateAmoeba::Pointer observer =
> CommandIterationUpdateAmoeba::New();
> optimizer->AddObserver( itk::FunctionEvaluationIterationEvent(), 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 );
>
> std::cout << "Writing output image...";
>
> try
> {
> moving_writer->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return EXIT_FAILURE;
> }
>
> std::cout << "Done" << std::endl;
>
>
>
> moving_writer = 0; // Delete
> moving_image = 0; // Delete
> movingImageReReader = 0; // Delete
> observer = 0; // Delete
> moving_roi = 0; // Delete
> fixed_roi = 0; // Delete
> transform = 0; // Delete
> registration = 0; // Delete
> interpolator = 0; // Delete
> optimizer = 0; // Delete
> metric = 0; // Delete
>
>
> 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