[ITK] Metric NaN and multi resolution registration
Matias Montroull
matimontg at gmail.com
Wed Dec 16 09:35:39 EST 2015
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/4aa6c5a4/attachment-0001.html>
More information about the Community
mailing list