[Insight-users] Fwd: Registration Speed using Normalized Mutual
Information
Luis Ibanez
luis.ibanez at kitware.com
Wed Jan 17 16:42:28 EST 2007
Hi Alexia,
The NormalizedMutualInformationHistogramImageToImageMetric
is computed from the joint histogram. It does not use
subsampling. You may want to try the MutualInformationMetric
(which is the Viola Wells implementation), or the Mattes
Mutual Information metric.
One minute is a suspiciously long time for one iteration...
Please find attached a version of a 3D registration using
Mattes Mutual Information. Give it a try and let us know
how it works for you.
Regards,
Luis
-------------------------------
Alexia Rodríguez Ruano wrote:
>
>>
>> Hi
>> I'm modifying the example ImageRegistration8.cxx to reguister 3D
>> images using Normalized Mutual Information.
>> The time it takes to calculate each iteration is too slow (1 minute
>> each one), and I am working with release version.
>
>
> The images are identical but one of them is translated and rotated and
> with a size of 181x271x181
> The number of histogram bins is 32
>
>> How can I solve that??
>
>
> The NormalizedMutualInformationHistogramImageToImageMetric is
> implemented directly from the histogram or follows any specific method
> (Viola, Mattes, etc)
>
>
>> -----------------------------------------
>> ...
>> ....
>> const unsigned int Dimension = 3;
>> typedef unsigned char PixelType;
>>
>> typedef itk::Image< PixelType, Dimension > FixedImageType;
>> typedef itk::Image< PixelType, Dimension > MovingImageType;
>> typedef itk::VersorRigid3DTransform< double > TransformType;
>> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
>> typedef itk::NormalizedMutualInformationHistogramImageToImageMetric<
>> FixedImageType,
>> MovingImageType >
>> MetricType;
>> typedef itk::ImageRegistrationMethod<
>> FixedImageType,
>> MovingImageType >
>> 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 );
>>
>> unsigned int numberOfHistogramBins = *argv6;
>> MetricType::HistogramType::SizeType histogramSize;
>> histogramSize[0] = numberOfHistogramBins;
>> histogramSize[1] = numberOfHistogramBins;
>> metric->SetHistogramSize( histogramSize );
>>
>> TransformType::Pointer transform = TransformType::New();
>> registration->SetTransform( transform );
>>
>> const unsigned int numberOfParameters =
>> transform->GetNumberOfParameters();
>> typedef MetricType::ScalesType ScalesType;
>> ScalesType scales( numberOfParameters );
>> const double transscale = 1.0 / 1000.0;
>> scales.Fill( 1.0 );
>> scales[0] = 1.0;
>> scales[1] = 1.0;
>> scales[2] = 1.0;
>> scales[3] = transscale;
>> scales[4] = transscale;
>> scales[5] = transscale;
>> metric->SetDerivativeStepLengthScales(scales);
>>
>> typedef itk::ImportImageFilter< PixelType, 3 >
>> FixedImageImportFilterType;
>> typedef itk::ImportImageFilter< PixelType, 3 >
>> MovingImageImportFilterType;
>>
>> FixedImageImportFilterType::Pointer FixedImageImport =
>> FixedImageImportFilterType::New();
>> MovingImageImportFilterType::Pointer MovingImageImport =
>> MovingImageImportFilterType::New();
>>
>> FixedImageImportFilterType::SizeType size;
>> size[0] = argv4[0]; // size along X
>> size[1] = argv4[1]; // size along Y
>> size[2] = argv4[2]; // size along Z
>> FixedImageImportFilterType::IndexType start;
>> start.Fill( 0 );
>> FixedImageImportFilterType::RegionType region;
>> region.SetIndex( start );
>> region.SetSize( size );
>> FixedImageImport->SetRegion( region );
>> const unsigned int numberOfPixels = size[0]*size[1]*size[2];
>>
>> MovingImageImportFilterType::SizeType sizem;
>> sizem[0] = argv5[0]; // size along X
>> sizem[1] = argv5[1]; // size along Y
>> sizem[2] = argv5[2]; // size along Z
>> MovingImageImportFilterType::IndexType startm;
>> startm.Fill( 0 );
>> MovingImageImportFilterType::RegionType regionm;
>> regionm.SetIndex( startm );
>> regionm.SetSize( sizem );
>> MovingImageImport->SetRegion( regionm );
>> const unsigned int numberOfPixelsm = sizem[0]*sizem[1]*sizem[2];
>>
>> const bool importImageFilterWillOwnTheBuffer = false;
>> FixedImageImport->SetImportPointer( argv1,
>> numberOfPixels,importImageFilterWillOwnTheBuffer );
>> MovingImageImport->SetImportPointer( argv2,
>> numberOfPixelsm,importImageFilterWillOwnTheBuffer );
>>
>> registration->SetFixedImage( FixedImageImport->GetOutput() );
>> registration->SetMovingImage( MovingImageImport->GetOutput() );
>>
>> FixedImageImport->Update();
>> registration->SetFixedImageRegion(
>> FixedImageImport->GetOutput()->GetBufferedRegion() );
>>
>> typedef itk::CenteredTransformInitializer< TransformType,
>> FixedImageType,
>> MovingImageType
>> >
>> TransformInitializerType;
>>
>> TransformInitializerType::Pointer initializer =
>>
>> TransformInitializerType::New();
>>
>> initializer->SetTransform( transform );
>> initializer->SetFixedImage( FixedImageImport->GetOutput() );
>> initializer->SetMovingImage( MovingImageImport->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 );
>> std::cout<< "transform->SetRotation( rotation )" << std::endl;
>> 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.2 );
>> optimizer->SetMinimumStepLength( 0.001 );
>> optimizer->SetRelaxationFactor( 0.90 );//Alexia
>> optimizer->SetNumberOfIterations( 100 );
>> optimizer->MaximizeOn();
>>
>> CommandIterationUpdate::Pointer observer =
>> CommandIterationUpdate::New();
>>
>> optimizer->AddObserver( itk::IterationEvent(), observer );
>>
>> registration->ReleaseDataFlagOn();
>>
>> try
>> {
>> registration->StartRegistration();
>> }
>> catch( itk::ExceptionObject & err )
>> {
>> std::cerr << "ExceptionObject caught !" << std::endl;
>> std::cerr << err << std::endl;
>> return -1;
>> }
>> ...
>> ...
>> ...
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: MultiModalityRegistration1.cxx,v $
Language: C++
Date: $Date: 2006/12/20 22:15:40 $
Version: $Revision: 1.16 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageRegistrationMethod.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkOrientedImage.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 "itkTransformFileWriter.h"
#include "metaCommand.h"
#include "metaOutput.h"
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
//
#include "itkCommand.h"
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;
}
};
int usage_examples()
{
std::cout << std::endl << "Usage Examples:" << std::endl;
std::cout << " -h" << std::endl;
std::cout << " Print usage statement and examples" << std::endl;
std::cout << std::endl;
std::cout << " FixedImageFilename" << std::endl;
std::cout << " MovingImageFilename" << std::endl;
std::cout << " OutputImageFilename" << std::endl;
std::cout << " TransformFilename" << std::endl;
return 0;
}
int main( int argc, char *argv[] )
{
MetaCommand command;
command.AddField("fixedImage","FixedImageFilename",MetaCommand::STRING,true);
command.AddField("movingImage","MovingImageFilename",MetaCommand::STRING,true);
command.AddField("outputImage","OutputImageFilename",MetaCommand::STRING,true);
command.SetOption("outputTransform","t",false,"TransformFilename",
MetaCommand::STRING);
// This should be put before the parsing
MetaOutput output;
output.SetMetaCommand(&command);
if ( !command.Parse(argc,argv) )
{
return EXIT_FAILURE;
}
std::string fixedImageFileName = command.GetValueAsString("fixedImage");
std::string movingImageFileName = command.GetValueAsString("movingImage");
std::string outputImageFileName = command.GetValueAsString("outputImage");
const unsigned int Dimension = 3;
typedef signed short PixelType;
typedef itk::OrientedImage< PixelType, Dimension > FixedImageType;
typedef itk::OrientedImage< 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::ImageRegistrationMethod<
FixedImageType,
MovingImageType > 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 );
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName( fixedImageFileName );
movingImageReader->SetFileName( movingImageFileName );
try
{
movingImageReader->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
try
{
fixedImageReader->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
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->GeometryOn();
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() );
FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
FixedImageType::SpacingType spacing = fixedImage->GetSpacing();
FixedImageType::SizeType size = fixedImage->GetLargestPossibleRegion().GetSize();
const double translationScaleFactor = 10.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = 1.0 / ( translationScaleFactor * spacing[0] * size[0] );
optimizerScales[4] = 1.0 / ( translationScaleFactor * spacing[1] * size[1] );
optimizerScales[5] = 1.0 / ( translationScaleFactor * spacing[2] * size[2] );
optimizer->SetScales( optimizerScales );
optimizer->SetMaximumStepLength( 0.2000 );
optimizer->SetMinimumStepLength( 0.0001 );
optimizer->SetNumberOfIterations( 300 );
optimizer->MaximizeOff();
const unsigned long numberOfImagePixels =
fixedImage->GetLargestPossibleRegion().GetNumberOfPixels();
const unsigned long numberOfSpatialSamples =
static_cast<unsigned long>( numberOfImagePixels * 0.010 );
metric->SetNumberOfSpatialSamples( numberOfSpatialSamples );
metric->SetNumberOfHistogramBins( 64 );
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 -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();
// Print out results
//
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;
output.AddFloatField("MetricValue","Final Value of the Metric",bestValue);
output.AddIntField("Iterations","Number of iterations",numberOfIterations);
output.Write();
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;
if(command.GetOptionWasSet("outputTransform"))
{
std::string outputTransform = command.GetValueAsString("outputTransform");
itk::TransformFileWriter::Pointer transformWriter = itk::TransformFileWriter::New();
transformWriter->SetFileName(outputTransform.c_str());
transformWriter->SetInput(transform);
transformWriter->Update();
}
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() );
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resampler->SetOutputOrigin( fixedImage->GetOrigin() );
resampler->SetOutputSpacing( fixedImage->GetSpacing() );
resampler->SetOutputDirection( fixedImage->GetDirection() );
resampler->SetDefaultPixelValue( 0 );
typedef itk::ImageFileWriter< FixedImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( outputImageFileName );
writer->UseCompressionOn();
writer->SetInput( resampler->GetOutput() );
writer->Update();
return 0;
}
More information about the Insight-users
mailing list