[Insight-users] Fwd: Registration Speed using Normalized Mutual Information

Xabier Artaechevarria Artieda xabiarta at unav.es
Mon Jan 22 06:43:26 EST 2007


Hi Alexia,
this interesting topic was already discussed in a previous thread:  
http://public.kitware.com/pipermail/insight-users/2006-May/017911.html

As Luis suggested, it is possible to modify the Mattes method to get  
the normalized mutual information.

Regards,
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



More information about the Insight-users mailing list