[Insight-users] registration 3D
Gael Pentang
gael.pentang at med.uni-duesseldorf.de
Mon Feb 21 09:54:15 EST 2011
//>/
/>/ 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110221/58514e27/attachment.htm>
More information about the Insight-users
mailing list