[Insight-users] Fwd: Registration Speed using Normalized Mutual
Information
Marius Staring
marius at isi.uu.nl
Fri Jan 19 10:59:13 EST 2007
Hi Alexia,
I can see a couple of things:
- the NormalizedMutualInformationHistogramImageToImageMetric uses all
voxels, Mattes uses a subset. Registration time is aproximately linearly
dependent on the number of samples that you use.
- the NormalizedMutualInformation calls the GetValueAndDerivative() of
the HistogramImageToImageMetric, which simply calls GetValue() and then
GetDerivative(). As a result the function ComputeHistogram() is called
twice, which is once too many. So the implementation is not efficient.
If you would have used the BSplineTransform, the difference would even
be larger, since the MattesMutualInformation takes advantage of the
compact support of the B-splines, where NormalizedMutualInformation does
not.
The above is probably not an exhaustive list.
regards,
Marius
Alexia Rodríguez Ruano wrote:
> Hi,
> I have tried the version of 3D registration you attached, and it
> results incredibly fast compared with the
> NormalizedMutualInformationHistogramImageToImageMetric. I can
> understand that computing the joint histogram can be more slow, but i
> continue thinking that one minute is not a tolerable time for one
> iteration.
> If you have any idea of how that can be solved or why it happens?
> Regards,
> Alexia
>
>
> At 10:42 PM 17-01-07, Luis Ibanez wrote:
>
>> 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
>>
>>
>> /*=========================================================================
>>
>>
>> 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;
>> }
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
--
Marius Staring
Image Sciences Institute
University Medical Centre Utrecht
Heidelberglaan 100, 3584 CX Utrecht, The Netherlands
phone: +31 (0)30 250 3186, fax: +31 (0)30 251 3399
marius at isi.uu.nl, http://www.isi.uu.nl/People/Marius
More information about the Insight-users
mailing list