[Insight-users] Fwd: Registration Speed using Normalized
Mutual Information
Alexia Rodríguez Ruano
arodriguez at mce.hggm.es
Fri Jan 19 10:12:33 EST 2007
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;
>}
More information about the Insight-users
mailing list