[Insight-users] Multiresolution registration and
MattesMutualInformation
Christian Kollee
sichkoll at stud.uni-erlangen.de
Mon Feb 13 12:45:20 EST 2006
Hi,
I try to do a 3D rigid registration between two volumes. For testing
purposes I constructed two 256x256x256 volumes and extended the
ImageRegistration8 example. Each volume includes a smaller black cube
(100x100x100). One cube is shifted (25x25x25) compared to the other one.
If I use a single resolution approach everything is fine and the
registration works. But if I switch to a multiresolution approach, I
always get the following error message:
itk::ExceptionObject (0x81c71e8)
Location: "Unknown"
File: /home/cip/adm/sichkoll/itk-2.4.1/Code/Algorithms/itkMattesMutualInformationImageToImageMetric.txx
Line: 986
Description: itk::ERROR: MattesMutualInformationImageToImageMetric(0x81c6fc8): Too many samples map outside moving image buffer
I checked the transformation parameters as suggested a few days ago and
I think they are correct. Maybe someone here is able to tell me how to
do the multiresolution right.
Regards,
Christian
-------------- next part --------------
#include <iostream>
#include <ctime>
#include "itkImageRegistrationMethod.h"
#include "itkMultiResolutionImageRegistrationMethod.h"
#include "itkMultiResolutionPyramidImageFilter.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImage.h"
#include "itkVersorRigid3DTransform.h"
#include "itkCenteredTransformInitializer.h"
#include "itkVersorRigid3DTransformOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkCommand.h"
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
typedef RegistrationInterfaceCommand Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
RegistrationInterfaceCommand() {};
public:
typedef TRegistration RegistrationType;
typedef RegistrationType * RegistrationPointer;
typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
typedef OptimizerType * OptimizerPointer;
void Execute(itk::Object * object, const itk::EventObject & event)
{
if( !(itk::IterationEvent().CheckEvent( &event )) )
{
return;
}
RegistrationPointer registration =
dynamic_cast<RegistrationPointer>( object );
OptimizerPointer optimizer = dynamic_cast< OptimizerPointer >(
registration->GetOptimizer() );
if ( registration->GetCurrentLevel() == 0 )
{
optimizer->SetMaximumStepLength( optimizer->GetCurrentStepLength() );
optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() );
}
else
{
optimizer->SetMaximumStepLength(
optimizer->GetCurrentStepLength() );
optimizer->SetMinimumStepLength(
optimizer->GetMinimumStepLength() / 10.0 );
}
}
void Execute(const itk::Object * , const itk::EventObject & )
{ return; }
};
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)
{
std::cerr << "Optimizer step: ";
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;
}
};
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile " << std::endl;
return 1;
}
const unsigned int Dimension = 3;
typedef float PixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
typedef itk::VersorRigid3DTransform< double > TransformType;
typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
typedef itk::MattesMutualInformationImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::MultiResolutionImageRegistrationMethod<
// typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
typedef itk::MultiResolutionPyramidImageFilter<
FixedImageType,
FixedImageType > FixedImagePyramidType;
typedef itk::MultiResolutionPyramidImageFilter<
MovingImageType,
MovingImageType > MovingImagePyramidType;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
FixedImagePyramidType::Pointer fixedImagePyramid =
FixedImagePyramidType::New();
MovingImagePyramidType::Pointer movingImagePyramid =
MovingImagePyramidType::New();
registration->SetFixedImagePyramid( fixedImagePyramid );
registration->SetMovingImagePyramid( movingImagePyramid );
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
TransformType::Pointer transform = TransformType::New();
registration->SetTransform( transform );
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName( argv[1] );
movingImageReader->SetFileName( argv[2] );
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
fixedImageReader->Update();
registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );
typedef itk::CenteredTransformInitializer< TransformType,
FixedImageType,
MovingImageType
> TransformInitializerType;
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform( transform );
initializer->SetFixedImage( fixedImageReader->GetOutput() );
initializer->SetMovingImage( movingImageReader->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;
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 );
optimizer->MinimizeOn();
// metric->UseAllPixelsOn();
metric->SetNumberOfSpatialSamples( 100000 );
metric->SetNumberOfHistogramBins( 10 );
metric->ReinitializeSeed( 39234823 );
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
CommandType::Pointer command = CommandType::New();
registration->AddObserver( itk::IterationEvent(), command );
registration->SetNumberOfLevels( 3 );
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
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();
std::cout << std::endl << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " versor X = " << versorX << std::endl;
std::cout << " versor Y = " << versorY << std::endl;
std::cout << " versor Z = " << versorZ << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Translation Z = " << finalTranslationZ << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
transform->SetParameters( finalParameters );
TransformType::MatrixType matrix = transform->GetRotationMatrix();
TransformType::OffsetType offset = transform->GetOffset();
std::cout << "Matrix = " << std::endl << matrix << std::endl;
std::cout << "Offset = " << std::endl << offset << std::endl;
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( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resampler->SetOutputOrigin( fixedImage->GetOrigin() );
resampler->SetOutputSpacing( fixedImage->GetSpacing() );
resampler->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] );
caster->SetInput( resampler->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
return EXIT_SUCCESS;
}
More information about the Insight-users
mailing list