[Fwd: Re: [Insight-users] Fwd: Registration Speed using Normalized
Mutual Information]
Marius Staring
marius at isi.uu.nl
Mon Jan 22 09:16:53 EST 2007
Maybe useful for others too:
Hi Xabi,
a collegea of mine (Stefan Klein) reimplemented the Mattes MI, making it
a bit more general. As a result it was also quite easy to extend it, and
indeed the itkParzenWindowNormalizedMutualInformationImageToImageMetric
implements Normalized MI, similar to what is done in the Mattes MI version.
For this class
- It is possible to set the number of histogram bins. This can be done
separately for the fixed and the moving image.
- It is possible to set the spline order of the parzen windows. This can
again be done separately.
- It is also possible to set the number of spatial samples. But we do it
differently from the itk implementation. We have defined a Sampler base
class from which for example a RandomSampler and a FullSampler are
derived. This is then set into an AdvancedImageToImageMetric,
equivalently to setting the fixed image, moving image, etc. The
itkParzenWindowNormalizedMutualInformationImageToImageMetric inherits
from this AdvancedImageToImageMetric. So if you want to use a subset of
randomly chosen samples, you make a RandomSampler, set the number of
desired samples, and then set it into the Normalized MI metric.
- The resulting Normalised MI is just as fast as the Mattes MI.
- But if you want to use it you are either stuck with our modified
framework (Samplers, AdvancedImageToImageMetric, etc) or you have to
change the code. Then again, the modified framework is created fully in itk.
Hope that helps, regards,
Marius
Xabier Artaechevarria Artieda wrote:
> Hi Marius,
>
> I am a PhD student at the Centre for Applied Medical Research of the
> University of Navarra in Spain. As part of my working I am trying to
> segment microCT images using an atlas and I am trying registration
> methods for that.
> I have question which is related to the ITK thread and to your work
> with Elastix, which I have seen on the web:
>
> The itkParzenWindowNormalizedMutualInformatioImageToImageMetric, is it
> a modification of the Mattes method to calculate the normalized MI? If
> so, is it possible to set the number of spatial samples, apart from
> the number of histogram bins?
>
> My problem is that I want to compute the NMI in a fast way, and I
> think this class might be very helpful.
>
> Thanks for your attention,
> Xabi
>
> --Xabier Artaechevarria
> Cancer Imaging Laboratory
> Centre for Applied Medical Research
> www.cima.es
>
>
>
> Marius Staring <marius at isi.uu.nl> ha escrito:
>
>> 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
>>
>> _______________________________________________
>> Insight-users mailing list
>> Insight-users at itk.org
>> http://www.itk.org/mailman/listinfo/insight-users
>
>
>
> ----------------------------------------------------------------
> Este mensaje ha sido enviado desde https://webmail.unav.es
>
>
--
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
--
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