[Insight-users] registration 3D
Luis Ibanez
luis.ibanez at kitware.com
Mon Feb 21 18:51:24 EST 2011
Hi Gael,
What sort of errors do you get ?
a) Segmentation fault ?
b) Exceptions abort ?
c) Warning messages ?
Do the errors happen in the first iteration ?
or do you get to see several iterations before the errors happen ?
You are using data importers in your code:
have you verified that the imported images make it to the
ITK level, before attempting to use them as input to the
registration ?
An easy way to do that is to write them to disk using the
ImageFileWriter.
Luis
-------------------------
On Mon, Feb 21, 2011 at 9:54 AM, Gael Pentang
<gael.pentang at med.uni-duesseldorf.de> wrote:
>>
>> Hi there,
>> I am trying to register two 3-D images,
>> using the MattesMutualInformationMetric with the versor3Drigid transform
>> and the
> VersorRigid3DTransformOptimizer optimizer .
> Actually after reading the images, I am using the normalizeImageFilter and
> the DiscreteGaussianImageFilter
> to preprocess the images before registration.
>
> But when I am calling the registration routine, I am getting errors.
> I guess there may be errors when I am calling the
> CenteredTransformInitializer, Could someone
> check my code and help me?
>
> Here is my code:
>
> // imreg.cpp : Defines the entry point for the DLL application.
> //
>
> // imreg.cpp : Defines the entry point for the DLL application.
> //
>
> #include "stdafx.h"
>
> // Header files for ITK image registration
> #include "itkImageRegistrationMethod.h"
> #include "itkMeanSquaresImageToImageMetric.h"
> #include "itkMattesMutualInformationImageToImageMetric.h"
> #include "itkMutualInformationImageToImageMetric.h"
> #include "itkMeanReciprocalSquareDifferenceImageToImageMetric.h"
> #include "itkNormalizedCorrelationImageToImageMetric.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkImage.h"
>
>
> #include "itkTimeProbesCollectorBase.h"
>
> #ifdef ITK_USE_REVIEW
> #include "itkMemoryProbesCollectorBase.h"
> #define itkProbesCreate() \
> itk::TimeProbesCollectorBase chronometer; \
> itk::MemoryProbesCollectorBase memorymeter
> #define itkProbesStart( text ) memorymeter.Start( text ); chronometer.Start(
> text )
> #define itkProbesStop( text ) chronometer.Stop( text ); memorymeter.Stop(
> text )
> #define itkProbesReport( stream ) chronometer.Report( stream );
> memorymeter.Report( stream )
> #else
> #define itkProbesCreate() \
> itk::TimeProbesCollectorBase chronometer
> #define itkProbesStart( text ) chronometer.Start( text )
> #define itkProbesStop( text ) chronometer.Stop( text )
> #define itkProbesReport( stream ) chronometer.Report( stream )
> #endif
>
> #include "itkBSplineDeformableTransform.h"
> #include "itkLBFGSBOptimizer.h"
> #include "itkVersorRigid3DTransform.h"
> #include "itkCenteredTransformInitializer.h"
> #include "itkVersorRigid3DTransformOptimizer.h"
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkSquaredDifferenceImageFilter.h"
>
> #include "itkImportImageFilter.h"
> #include "itkMultiThreader.h"
>
> // The following section of code implements a Command observer
> // that will monitor the evolution of the registration process.
> //
> #include "itkCommand.h"
>
> #include "itkNormalizeImageFilter.h"
> #include "itkDiscreteGaussianImageFilter.h"
>
> #define WM_MESSAGE_STATUS WM_USER +13
>
> 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::VersorRigid3DTransformOptimizer 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;
> }
> };
>
> BOOL APIENTRY DllMain( HANDLE hModule,
> DWORD ul_reason_for_call,
> LPVOID lpReserved
> )
> {
> return TRUE;
> }
>
>
> __declspec(dllexport) int RegisterImage2to1(short *pixdata1, short
> *pixdata2, int xdim, int ydim, int zdim){
> // define image dimension and data type
> const unsigned int ImageDimension = 3;
> typedef float PixelType;
> // define fixed and moving image
> typedef itk::Image< PixelType, ImageDimension > FixedImageType;
> typedef itk::Image< PixelType, ImageDimension > MovingImageType;
>
> typedef float InternalPixelType;
> typedef itk::Image< InternalPixelType, ImageDimension > InternalImageType;
>
> typedef itk::VersorRigid3DTransform< double > TransformType;
>
> const unsigned int SpaceDimension = ImageDimension;
> const unsigned int SplineOrder = 3;
> typedef double CoordinateRepType;
>
> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
>
> /*typedef itk::NormalizedCorrelationImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;*/
>
> /*typedef itk::MeanReciprocalSquareDifferenceImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;*/
>
> /*typedef itk::MeanSquaresImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
>
> typedef itk::MutualInformationImageToImageMetric<
> InternalImageType,
> InternalImageType >
> MetricType;*/
>
>
> typedef itk::MattesMutualInformationImageToImageMetric<
> InternalImageType,
> InternalImageType > MetricType;
>
> typedef itk:: LinearInterpolateImageFunction<
> InternalImageType,
> double > InterpolatorType;
>
> typedef itk::ImageRegistrationMethod<
> InternalImageType,
> InternalImageType > 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 );
>
> /*metric->SetFixedImageStandardDeviation( 0.4 );
> metric->SetMovingImageStandardDeviation( 0.4 );*/
>
> // define import filter #1 and #2 for data import to itk
> typedef itk::ImportImageFilter< PixelType, ImageDimension>
> ImportFilterType;
> ImportFilterType::Pointer importFilter1 = ImportFilterType::New();
> ImportFilterType::Pointer importFilter2 = ImportFilterType::New();
> ImportFilterType::SizeType size;
> size[0] = xdim;
> size[1] = ydim;
> size[2] = zdim;
> ImportFilterType::IndexType start;
> start.Fill( 0 );
> ImportFilterType::RegionType region;
> region.SetIndex( start );
> region.SetSize( size );
> importFilter1->SetRegion( region );
> importFilter2->SetRegion( region );
> double origin[ ImageDimension ];
> origin[0] = 0.0;
> origin[1] = 0.0;
> origin[2] = 0.0;
> importFilter1->SetOrigin( origin );
> importFilter2->SetOrigin( origin );
> double spacing[ ImageDimension ];
> spacing[0] = 1.0;
> spacing[1] = 1.0;
> spacing[2] = 1.0;
> importFilter1->SetSpacing( spacing);
> importFilter2->SetSpacing( spacing);
> // create buffers for images
> const unsigned int numberOfPixels = size[0]*size[1]*size[2];
> PixelType * localBuffer1 = new PixelType[ numberOfPixels ];
> PixelType * it1 = localBuffer1;
> PixelType * localBuffer2 = new PixelType[ numberOfPixels ];
> PixelType * it2 = localBuffer2;
> // set pixels into image #1: fixed image
> for(unsigned int i=0;i<size[0]*size[1]*size[2];i++){
> *it1++ = (float) pixdata1[i];
> // uncomment following lines for thresholding before registration
> /*if(pixdata1[i] > 16)
> *it1++ = (float) pixdata1[i];
> else
> *it1++ = 0;*/
> }
> importFilter1->SetImportPointer( localBuffer1, numberOfPixels, true);
> //registration->SetFixedImage( importFilter1->GetOutput() );
> importFilter1->Update();
> typedef FixedImageType::SpacingType SpacingType;
> typedef FixedImageType::PointType OriginType;
> typedef FixedImageType::RegionType RegionType;
> typedef FixedImageType::SizeType SizeType;
> FixedImageType::Pointer fixedImage = importFilter1->GetOutput();
> const SpacingType fixedSpacing = fixedImage->GetSpacing();
> const OriginType fixedOrigin = fixedImage->GetOrigin();
> const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
> const SizeType fixedSize = fixedRegion.GetSize();
> TransformType::InputPointType centerFixed;
> centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
> centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
> centerFixed[2] = fixedOrigin[2] + fixedSpacing[2] * fixedSize[2] / 2.0;
> // set pixels into image #2 moving image
> for(unsigned int i=0;i<size[0]*size[1]*size[2];i++){
> *it2++ = (float) pixdata2[i];
> // uncomment following lines for thresholding before registration
> /*if(pixdata2[i] > 16)
> *it2++ = (float) pixdata2[i];
> else
> *it2++ = 0;*/
> }
> importFilter2->SetImportPointer( localBuffer2, numberOfPixels, true);
> //registration->SetMovingImage( importFilter2->GetOutput() );
> importFilter2->Update();
> MovingImageType::Pointer movingImage = importFilter2->GetOutput();
> // set center of image
> const SpacingType movingSpacing = movingImage->GetSpacing();
> const OriginType movingOrigin = movingImage->GetOrigin();
> const RegionType movingRegion = movingImage->GetLargestPossibleRegion();
> const SizeType movingSize = movingRegion.GetSize();
> TransformType::InputPointType centerMoving;
> centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0;
> centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] / 2.0;
> centerMoving[2] = movingOrigin[2] + movingSpacing[2] * movingSize[2] / 2.0;
> //FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
> registration->SetFixedImageRegion( fixedRegion );
>
> /*typedef TransformType::RegionType RegionType;
> RegionType bsplineRegion;
> RegionType::SizeType gridSizeOnImage;
> RegionType::SizeType gridBorderSize;
> RegionType::SizeType totalGridSize;
>
> gridSizeOnImage.Fill( 5 ); // 12 ?
> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2
> upper )
> totalGridSize = gridSizeOnImage + gridBorderSize;
>
> bsplineRegion.SetSize( totalGridSize );
>
> //typedef TransformType::SpacingType SpacingType;
> //SpacingType spacing = fixedImage->GetSpacing();
> //typedef TransformType::OriginType OriginType;
> //OriginType origin = fixedImage->GetOrigin();;
>
> FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
>
> for(unsigned int r=0; r<ImageDimension; r++)
> {
> spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
> static_cast<double>(gridSizeOnImage[r] - 1);
> origin[r] -= spacing[r];
> }
>
> transform->SetGridSpacing( spacing );
> transform->SetGridOrigin( origin );
> transform->SetGridRegion( bsplineRegion );
> transform->SetGridDirection( fixedImage->GetDirection() );*/
>
> typedef itk::NormalizeImageFilter<
> FixedImageType,
> InternalImageType
> > FixedNormalizeFilterType;
> typedef itk::NormalizeImageFilter<
> MovingImageType,
> InternalImageType
> > MovingNormalizeFilterType;
>
> FixedNormalizeFilterType::Pointer fixedNormalizer =
> FixedNormalizeFilterType::New();
>
> MovingNormalizeFilterType::Pointer movingNormalizer =
> MovingNormalizeFilterType::New();
>
>
> typedef itk::DiscreteGaussianImageFilter<
> InternalImageType,
> InternalImageType
> > GaussianFilterType;
>
> GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
> GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
>
> fixedSmoother->SetVariance( 2.0 );
> movingSmoother->SetVariance( 2.0 );
>
> fixedNormalizer->SetInput( importFilter1->GetOutput() );
> movingNormalizer->SetInput( importFilter2->GetOutput() );
>
> fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
> movingSmoother->SetInput( movingNormalizer->GetOutput() );
>
> registration->SetFixedImage( fixedSmoother->GetOutput() );
> registration->SetMovingImage( movingSmoother->GetOutput() );
>
> fixedNormalizer->Update();
> FixedImageType::RegionType fixedImageRegion =
> fixedNormalizer->GetOutput()->GetBufferedRegion();
> registration->SetFixedImageRegion( fixedImageRegion );
>
> typedef itk::CenteredTransformInitializer< TransformType,
> FixedImageType,
> MovingImageType
> >
> TransformInitializerType;
>
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
>
> initializer->SetTransform( transform );
> initializer->SetFixedImage( importFilter1->GetOutput() );
> initializer->SetMovingImage( importFilter2->GetOutput() );
> initializer->MomentsOn();
> initializer->InitializeTransform();
> 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.0;
>
> rotation.Set( axis, angle );
>
> transform->SetRotation( rotation );
> registration->SetInitialTransformParameters( transform->GetParameters() );
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
> const double translationScale = 1.0 / 1000.0;
>
> optimizerScales[0] = 1.0;
> optimizerScales[1] = 1.0;
> optimizerScales[2] = 1.0;
> optimizerScales[3] = translationScale;
> optimizerScales[4] = translationScale;
> optimizerScales[5] = translationScale;
>
> optimizer->SetScales( optimizerScales );
>
> optimizer->SetMaximumStepLength( 0.2000 );
> optimizer->SetMinimumStepLength( 0.0001 );
>
> optimizer->SetNumberOfIterations( 200 );
>
> 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();
>
> 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();
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> TransformType::Pointer finalTransform = TransformType::New();
>
> finalTransform->SetCenter( transform->GetCenter() );
>
> finalTransform->SetParameters( finalParameters );
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>
> resampler->SetTransform( finalTransform );
> resampler->SetInput( movingImage );
>
>
>
> resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
> resampler->SetOutputOrigin( fixedImage->GetOrigin() );
> resampler->SetOutputSpacing( fixedImage->GetSpacing() );
> resampler->SetOutputDirection( fixedImage->GetDirection() );
> resampler->SetDefaultPixelValue( 0 );
> //resampler->SetDefaultPixelValue( 100 );
>
> movingImage = resampler->GetOutput();
> resampler->Update();
>
> // give back the registered image
> typedef itk::ImageRegionConstIterator< MovingImageType > IteratorType;
> IteratorType it( movingImage, region );
> it.GoToBegin();
> int k=0;
> while(! it.IsAtEnd()){
> pixdata2[k] = (short) it.Get();
> k++;
> ++it;
> }
>
> return EXIT_SUCCESS;
> return 0;
> };
>
>
>> Thanks in advance for your help
> Gael
>>
>
>
> --
> M.Sc. Comp. Eng. Gael Pentang
> Research Assistent
> University Dusseldorf
> Medical Faculty
> Department of Diagnostic and Interventional Radiology
> Moorenstrasse 5
> D-40225 Dusseldorf
> Germany
> Tel.: +49 211 8117430
> Mail: gael.pentang at med.uni-duesseldorf.de
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>
More information about the Insight-users
mailing list