[Insight-users] registration 3D
Gael Pentang
gael.pentang at med.uni-duesseldorf.de
Tue Feb 22 04:02:41 EST 2011
Hi Luis,
thanks for ur reply.
Actually I am getting an exception abort at the first iteration.
I have already verified that the imported images are OK.
In a first step, I did the registration without using the
normalizeImagerFilter and the DiscreteGaussianImageFilter.
Meaning that I run the registration on the original images that I
imported without filtering them and the program works fine.
But when I try to use the Filters before registration, I am getting
that exception.
To precise my question:
I am using the _versor3DRigidTransform_ and it needs the
_CenteredTransformInitializer_ to initialize the transform.
initializer->SetFixedImage( importFilter1->GetOutput() );
initializer->SetMovingImage( importFilter2->GetOutput() );
Now I am confused, what should be the input of the initializer? Should I
use the filtered data as input to initialize the transform or should I
use the original data without any processing to initialize the transform.
Actually in my code, I am using the original data (and it works fine),
but when I try to use the filtered data, the registration abort with an
exception at the first iteration.
Thanks in advance for your reply.
Gael.
Am 22/02/2011 00:51, schrieb Luis Ibanez:
> 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
>>
>>
--
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110222/f5560612/attachment.htm>
More information about the Insight-users
mailing list