[ITK] Metric NaN and multi resolution registration
Luigi Riba
ribaluigi at gmail.com
Wed Dec 16 09:44:26 EST 2015
Hello,
in fact I am surprised by this warning since I am using the
CenteredTransformInitializer. This is the piece of code that involves the
initialization of the transform:
typedef itk::CenteredTransformInitializer<TransformType, ImageType,
> ImageType> TransformInitializerType;
> TransformInitializerType::Pointer transformInitializer =
> TransformInitializerType::New();
> transformInitializer->SetTransform(transform);
> transformInitializer->SetFixedImage(fixedReader->GetOutput());
> transformInitializer->SetMovingImage(movingReader->GetOutput());
> transformInitializer->GeometryOn();
> transformInitializer->InitializeTransform();
...
> registration->SetInitialTransform(transform);
Thanks for the support,
Luigi
2015-12-16 15:35 GMT+01:00 Matias Montroull <matimontg at gmail.com>:
> The second warning looks like centers not aligned properly therefore the
> registration fails. Are you using CenteredTransformInitializer?
>
> El mié., 16 de dic. de 2015 a la(s) 06:24, Luigi Riba <ribaluigi at gmail.com>
> escribió:
>
>> Dear Matias and ITK community,
>>
>> thank you very much for your quick and detailed reply.
>>
>> I have tried a little bit of tinkering starting from your suggestion but
>> I didn't overcome the errors.
>>
>> I am registering two images using the v4 framework with the following
>> header
>> #include <itkImageRegistrationMethodv4.h>
>> #include <itkMeanSquaresImageToImageMetricv4.h>
>> #include <itkRegularStepGradientDescentOptimizerv4.h>
>>
>> I have initialized the transform with the transform initializer and
>> computed the optimizer parameters with
>> RegistratonParametersFromPhysicalShift and after a little of playing I left
>> the relaxation factor to .5.
>>
>> *The registration process runs correctly at full resolution. *
>>
>> If I choose a multi resolution approach, like
>>
>> const unsigned int numberOfLevels = 3;
>>> RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
>>> shrinkFactorsPerLevel.SetSize(numberOfLevels);
>>> shrinkFactorsPerLevel[0] = 4;
>>> shrinkFactorsPerLevel[1] = 2;
>>> shrinkFactorsPerLevel[2] = 1;
>>> RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
>>> smoothingSigmasPerLevel.SetSize(numberOfLevels);
>>> smoothingSigmasPerLevel[0] = 2;
>>> smoothingSigmasPerLevel[1] = 1;
>>> smoothingSigmasPerLevel[2] = 0;
>>
>>
>> the code stops working properly.
>>
>> In fact, I get the following warnings:
>>
>> - RegistrationParameterScalesFromPhysicalShift (00000000021F2630):
>> Variation in any parameter won't change a voxel position. The default
>> scales (1.0) are used to avoid division-by-zero.
>> - WARNING: In
>> f:\libs\insighttoolkit-4.8.2\modules\numerics\optimizersv4\include\itkObjectToObjectMetric.hxx,
>> line 529 MeanSquaresImageToImageMetricv4 (00000000650FCB40): No valid
>> points were found during metric evaluation. For image metrics, verify that
>> the images overlap appropriately. For instance, you can align the image
>> centers by translation. For point-set metrics, verify that the fixed
>> points, once transformed into the virtual domain space, actually lie within
>> the virtual domain.
>>
>> the metric is NaN and the registration process ends before its first step.
>>
>> Could you explain me the rationale behind this errors? I am probably
>> overlooking something trivial.
>>
>> Best,
>>
>> Luigi
>>
>>
>>
>> 2015-12-07 16:06 GMT+01:00 Matias Montroull <matimontg at gmail.com>:
>>
>>> Actually the min and max steplenght are taken into account, I use these:
>>> Min: 0.0001
>>> Max: 0.1
>>>
>>> El lun., 7 de dic. de 2015 a la(s) 12:04, Matias Montroull <
>>> matimontg at gmail.com> escribió:
>>>
>>>> Hi Luigi, it could be for various reasons, mailny the transform
>>>> initializer you use plus the number of samples and translationscale
>>>> Try this code I have which I've done for MultiResolution and works fine
>>>> for me:
>>>>
>>>> When you enter the parameters, try this:
>>>> Relaxation: 0.99
>>>> Number of Bins: 50
>>>> Samples 3 (this is actually translated to 3% in my code as you will see
>>>> below)
>>>> maximum and minimum steplenght is not taken into account so you can put
>>>> 0 and 0
>>>> translation scale 0.0001
>>>>
>>>> #include "itkAffineTransform.h"
>>>> #include "itkRegistrationParameterScalesFromPhysicalShift.h"
>>>> #include "itkCenteredTransformInitializer.h"
>>>> #include "itkLandmarkBasedTransformInitializer.h"
>>>> #include "itkMultiResolutionImageRegistrationMethod.h"
>>>> #include "itkMattesMutualInformationImageToImageMetric.h"
>>>> #include "itkMeanSquaresImageToImageMetric.h"
>>>> #include "itkRegularStepGradientDescentOptimizer.h"
>>>> #include "itkRecursiveMultiResolutionPyramidImageFilter.h"
>>>> #include "itkImage.h"
>>>> #include "itkImageFileReader.h"
>>>> #include "itkImageFileWriter.h"
>>>> #include "itkResampleImageFilter.h"
>>>> #include "itkCastImageFilter.h"
>>>> #include "itkCheckerBoardImageFilter.h"
>>>> #include "itkRescaleIntensityImageFilter.h"
>>>>
>>>> #include "itkNormalizeImageFilter.h"
>>>>
>>>> #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() : m_CumulativeIterationIndex(0) {};
>>>> public:
>>>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>>>> typedef const OptimizerType * OptimizerPointer;
>>>> void Execute(itk::Object *caller, const itk::EventObject & event)
>>>> ITK_OVERRIDE
>>>> {
>>>> Execute((const itk::Object *)caller, event);
>>>> }
>>>> void Execute(const itk::Object * object, const itk::EventObject &
>>>> event) ITK_OVERRIDE
>>>> {
>>>> OptimizerPointer optimizer = static_cast<OptimizerPointer>(object);
>>>> if (!(itk::IterationEvent().CheckEvent(&event)))
>>>> {
>>>> return;
>>>> }
>>>> std::cout << optimizer->GetCurrentIteration() << " ";
>>>> std::cout << optimizer->GetValue() << " ";
>>>> std::cout << optimizer->GetCurrentPosition() << " " <<
>>>> m_CumulativeIterationIndex++ << std::endl;
>>>> }
>>>> private:
>>>> unsigned int m_CumulativeIterationIndex;
>>>> };
>>>>
>>>> template <typename TRegistration>
>>>> class RegistrationInterfaceCommand : public itk::Command
>>>> {
>>>> public:
>>>> typedef RegistrationInterfaceCommand Self;
>>>> typedef itk::Command Superclass;
>>>> typedef itk::SmartPointer<Self> Pointer;
>>>> itkNewMacro(Self);
>>>> protected:
>>>> RegistrationInterfaceCommand() {};
>>>> public:
>>>> typedef TRegistration RegistrationType;
>>>> typedef RegistrationType *
>>>> RegistrationPointer;
>>>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>>>> typedef OptimizerType * OptimizerPointer;
>>>> void Execute(itk::Object * object, const itk::EventObject & event)
>>>> ITK_OVERRIDE
>>>> {
>>>> if (!(itk::IterationEvent().CheckEvent(&event)))
>>>> {
>>>> return;
>>>> }
>>>> RegistrationPointer registration =
>>>> static_cast<RegistrationPointer>(object);
>>>> OptimizerPointer optimizer =
>>>> static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());
>>>> std::cout << "-------------------------------------" << std::endl;
>>>> std::cout << "MultiResolution Level : "
>>>> << registration->GetCurrentLevel() << std::endl;
>>>> std::cout << std::endl;
>>>> if (registration->GetCurrentLevel() == 0)
>>>> {
>>>> optimizer->SetMaximumStepLength(1.0);
>>>> optimizer->SetMinimumStepLength(0.001);
>>>> }
>>>> else
>>>> {
>>>> optimizer->SetMaximumStepLength(optimizer->GetMaximumStepLength() /
>>>> 4.0);
>>>> optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() /
>>>> 10.0);
>>>> }
>>>> }
>>>> void Execute(const itk::Object *, const itk::EventObject &) ITK_OVERRIDE
>>>> { return; }
>>>> };
>>>> int main(int argc, char *argv[])
>>>> {
>>>> if (argc < 10)
>>>> {
>>>> std::cerr << "Faltan Parametros " << std::endl;
>>>> std::cerr << "Uso: " << argv[0];
>>>> std::cerr << " imagenfija imagenflotante ";
>>>> std::cerr << " salida";
>>>> std::cerr << " Iteraciones Relaxation ";
>>>> std::cerr << " NumeroBins Samples%(1/100) ";
>>>> std::cerr << " StepMax(No Habilitado) StepMin(No Habilitado)
>>>> TranslationFactor";
>>>> std::cerr << " " << std::endl;
>>>>
>>>> return EXIT_FAILURE;
>>>> }
>>>>
>>>> const unsigned int Dimension = 3;
>>>> typedef signed short PixelType;
>>>> typedef itk::Image< PixelType, Dimension > FixedImageType;
>>>> typedef itk::Image< PixelType, Dimension > MovingImageType;
>>>> typedef float InternalPixelType;
>>>> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
>>>>
>>>> typedef itk::AffineTransform< double, Dimension > TransformType;
>>>>
>>>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>>>> typedef itk::LinearInterpolateImageFunction<
>>>> InternalImageType,
>>>> double > InterpolatorType;
>>>> typedef itk::MattesMutualInformationImageToImageMetric<
>>>> InternalImageType,
>>>> InternalImageType > MetricType;
>>>> typedef OptimizerType::ScalesType OptimizerScalesType;
>>>> typedef itk::MultiResolutionImageRegistrationMethod<
>>>> InternalImageType,
>>>> InternalImageType > RegistrationType;
>>>> OptimizerType::Pointer optimizer = OptimizerType::New();
>>>> InterpolatorType::Pointer interpolator = InterpolatorType::New();
>>>> RegistrationType::Pointer registration = RegistrationType::New();
>>>> MetricType::Pointer metric = MetricType::New();
>>>>
>>>>
>>>> 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(argv[1]);
>>>> movingImageReader->SetFileName(argv[2]);
>>>>
>>>> fixedImageReader->Update();
>>>> movingImageReader->Update();
>>>>
>>>> typedef itk::CastImageFilter<
>>>> FixedImageType, InternalImageType > FixedCastFilterType;
>>>> typedef itk::CastImageFilter<
>>>> MovingImageType, InternalImageType > MovingCastFilterType;
>>>> FixedCastFilterType::Pointer fixedCaster = FixedCastFilterType::New();
>>>> MovingCastFilterType::Pointer movingCaster =
>>>> MovingCastFilterType::New();
>>>> fixedCaster->SetInput(fixedImageReader->GetOutput());
>>>> movingCaster->SetInput(movingImageReader->GetOutput());
>>>>
>>>> registration->SetFixedImage(fixedCaster->GetOutput());
>>>> registration->SetMovingImage(movingCaster->GetOutput());
>>>> fixedCaster->Update();
>>>>
>>>> registration->SetFixedImageRegion(fixedCaster->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::ParametersType ParametersType;
>>>> const unsigned int numberOfParameters =
>>>> transform->GetNumberOfParameters();
>>>> ParametersType parameters(numberOfParameters);
>>>> registration->SetInitialTransformParameters(transform->GetParameters());
>>>>
>>>> OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());
>>>>
>>>> double translationScale = atof(argv[10]);
>>>> optimizerScales[0] = 1.0;
>>>> optimizerScales[1] = 1.0;
>>>> optimizerScales[2] = 1.0;
>>>> optimizerScales[3] = 1.0;
>>>> optimizerScales[4] = 1.0;
>>>> optimizerScales[5] = 1.0;
>>>> optimizerScales[6] = 1.0;
>>>> optimizerScales[7] = 1.0;
>>>> optimizerScales[8] = 1.0;
>>>> optimizerScales[9] = translationScale;
>>>> optimizerScales[10] = translationScale;
>>>> optimizerScales[11] = translationScale;
>>>>
>>>> optimizer->SetScales(optimizerScales);
>>>>
>>>> FixedImageType::RegionType fixedImageRegion =
>>>> fixedCaster->GetOutput()->GetBufferedRegion();
>>>> const unsigned int numberOfPixels =
>>>> fixedImageRegion.GetNumberOfPixels();
>>>> double porcentage_samples = atof(argv[7]);
>>>> const unsigned int numberOfSamples = static_cast<unsigned
>>>> int>(numberOfPixels * porcentage_samples);
>>>>
>>>> optimizer->SetNumberOfIterations(atof(argv[4]));
>>>> optimizer->SetRelaxationFactor(atof(argv[5]));
>>>> metric->SetNumberOfHistogramBins(atof(argv[6]));
>>>> metric->SetNumberOfSpatialSamples(numberOfSamples);
>>>> optimizer->SetMaximumStepLength(atof(argv[8]));
>>>> optimizer->SetMinimumStepLength(atof(argv[9]));
>>>>
>>>>
>>>> optimizer->MinimizeOn();
>>>> CommandIterationUpdate::Pointer observer =
>>>> CommandIterationUpdate::New();
>>>> optimizer->AddObserver(itk::IterationEvent(), observer);
>>>>
>>>>
>>>> typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
>>>> CommandType::Pointer command = CommandType::New();
>>>> registration->AddObserver(itk::IterationEvent(), command);
>>>> registration->SetNumberOfLevels(1);
>>>> registration->SetOptimizer(optimizer);
>>>> registration->SetInterpolator(interpolator);
>>>> registration->SetMetric(metric);
>>>>
>>>> FixedImageType::Pointer imagenfija = fixedImageReader->GetOutput();
>>>> MovingImageType::Pointer imagenflotante =
>>>> movingImageReader->GetOutput();
>>>>
>>>>
>>>> typedef RegistrationType::ParametersType ParametersType;
>>>> ParametersType InitialParameters =
>>>> registration->GetInitialTransformParameters();
>>>> std::cout << "TranslationScale: " << translationScale << std::endl;
>>>> std::cout << "MaximumStep: " << optimizer->GetMaximumStepLength() <<
>>>> std::endl;
>>>> std::cout << "MinimunStep: " << optimizer->GetMinimumStepLength() <<
>>>> std::endl;
>>>> std::cout << "Number of Bins: " << metric->GetNumberOfHistogramBins()
>>>> << std::endl;
>>>> std::cout << "Number of Samples: " <<
>>>> metric->GetNumberOfSpatialSamples() << std::endl;
>>>> std::cout << "Iteraciones: " << optimizer->GetNumberOfIterations() <<
>>>> std::endl;
>>>> std::cout << "Relaxation: " << optimizer->GetRelaxationFactor() <<
>>>> std::endl;
>>>> std::cout << "Initial Parameters: " << InitialParameters << std::endl;
>>>> std::cout << "Origen Fija: " << imagenfija->GetOrigin() << std::endl;
>>>> std::cout << "Origen Flotante: " << imagenflotante->GetOrigin() <<
>>>> std::endl;
>>>> std::cout << "Samples: " << numberOfSamples << std::endl;
>>>>
>>>> try
>>>> {
>>>> registration->Update();
>>>> std::cout << "Optimizer stop condition: "
>>>> << registration->GetOptimizer()->GetStopConditionDescription()
>>>> << std::endl;
>>>> }
>>>> catch (itk::ExceptionObject & err)
>>>> {
>>>> std::cout << "ExceptionObject caught !" << std::endl;
>>>> std::cout << err << std::endl;
>>>> return EXIT_FAILURE;
>>>> }
>>>> std::cout << "Optimizer Stopping Condition = "
>>>> << optimizer->GetStopCondition() << std::endl;
>>>>
>>>> ParametersType finalParameters =
>>>> registration->GetLastTransformParameters();
>>>>
>>>> double TranslationAlongX = finalParameters[9];
>>>> double TranslationAlongY = finalParameters[10];
>>>> double TranslationAlongZ = finalParameters[11];
>>>> unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>>>> double bestValue = optimizer->GetValue();
>>>>
>>>> std::cout << "Resultado = " << std::endl;
>>>> std::cout << " Traslacion X = " << TranslationAlongX << std::endl;
>>>> std::cout << " Traslacion Y = " << TranslationAlongY << std::endl;
>>>> std::cout << " Traslacion Z = " << TranslationAlongZ << std::endl;
>>>> std::cout << " Iteraciones = " << numberOfIterations << std::endl;
>>>> std::cout << " Metrica = " << bestValue << std::endl;
>>>>
>>>>
>>>> typedef itk::ResampleImageFilter<
>>>> MovingImageType,
>>>> FixedImageType > ResampleFilterType;
>>>> TransformType::Pointer finalTransform = TransformType::New();
>>>> finalTransform->SetParameters(finalParameters);
>>>> finalTransform->SetFixedParameters(transform->GetFixedParameters());
>>>> ResampleFilterType::Pointer resample = ResampleFilterType::New();
>>>> resample->SetTransform(finalTransform);
>>>> resample->SetInput(movingImageReader->GetOutput());
>>>> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
>>>> PixelType backgroundGrayLevel = 0; //estaba en 100
>>>>
>>>> resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
>>>> resample->SetOutputOrigin(fixedImage->GetOrigin());
>>>> resample->SetOutputSpacing(fixedImage->GetSpacing());
>>>> resample->SetOutputDirection(fixedImage->GetDirection());
>>>> resample->SetDefaultPixelValue(backgroundGrayLevel);
>>>> typedef signed short OutputPixelType;
>>>> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
>>>> typedef itk::CastImageFilter<
>>>> FixedImageType,
>>>> OutputImageType > CastFilterType;
>>>> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>>>> WriterType::Pointer writer = WriterType::New();
>>>> CastFilterType::Pointer caster = CastFilterType::New();
>>>> writer->SetFileName(argv[3]);
>>>> caster->SetInput(resample->GetOutput());
>>>> writer->SetInput(caster->GetOutput());
>>>> writer->Update();
>>>> Sleep(5000);
>>>> /*
>>>> typedef itk::CheckerBoardImageFilter< FixedImageType >
>>>> CheckerBoardFilterType;
>>>> CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
>>>> checker->SetInput1(fixedImage);
>>>> checker->SetInput2(resample->GetOutput());
>>>> caster->SetInput(checker->GetOutput());
>>>> writer->SetInput(caster->GetOutput());
>>>> resample->SetDefaultPixelValue(0);
>>>>
>>>> resample->SetTransform(finalTransform);
>>>> if (argc > 10)
>>>> {
>>>> writer->SetFileName(argv[11]);
>>>> writer->Update();
>>>> }*/
>>>> return EXIT_SUCCESS;
>>>> }
>>>>
>>>>
>>>> El lun., 7 de dic. de 2015 a la(s) 11:43, Luigi Riba <
>>>> ribaluigi at gmail.com> escribió:
>>>>
>>>>> Dear ITK community,
>>>>>
>>>>> today I was trying to implement a multi-resolution registration
>>>>> process for 3d images.
>>>>>
>>>>> I started from the example in the tutorial. Unfortunately, I've
>>>>> received the following warning:
>>>>>
>>>>>> WARNING: In
>>>>>> f:\libs\insighttoolkit-4.8.2\modules\numerics\optimizersv4\include\itkObjectToObjectMetric.hxx,
>>>>>> line 529
>>>>>> MeanSquaresImageToImageMetricv4 (000000006506C020): No valid points
>>>>>> were found during metric evaluation. For image metrics, verify that the
>>>>>> images overlap appropriately. For instance, you can align the image centers
>>>>>> by translation. For point-set metrics, verify that the fixed points, once
>>>>>> transformed into the virtual domain space, actually lie within the virtual
>>>>>> domain.
>>>>>
>>>>>
>>>>> These are the parameters I have used for the multi resolution process:
>>>>>> const unsigned int numberOfLevels = 3;
>>>>>> RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
>>>>>> shrinkFactorsPerLevel.SetSize(numberOfLevels);
>>>>>> shrinkFactorsPerLevel[0] = 4;
>>>>>> shrinkFactorsPerLevel[1] = 2;
>>>>>> shrinkFactorsPerLevel[2] = 1;
>>>>>> RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
>>>>>> smoothingSigmasPerLevel.SetSize(numberOfLevels);
>>>>>> smoothingSigmasPerLevel[0] = 2;
>>>>>> smoothingSigmasPerLevel[1] = 1;
>>>>>> smoothingSigmasPerLevel[2] = 0;
>>>>>> registration->SetNumberOfLevels(numberOfLevels);
>>>>>> registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
>>>>>> registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
>>>>>
>>>>>
>>>>> Could you please me explain what is the meaning of this warning?
>>>>>
>>>>> By the way, I am working with affine transformations and I have used
>>>>> the transform initializer.
>>>>>
>>>>> Best,
>>>>>
>>>>> Luigi
>>>>> _______________________________________________
>>>>> Community mailing list
>>>>> Community at itk.org
>>>>> http://public.kitware.com/mailman/listinfo/community
>>>>>
>>>> --
>>>> Matias
>>>>
>>> --
>>> Matias
>>>
>>
>> --
> Matias
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20151216/1664d067/attachment-0001.html>
More information about the Community
mailing list