[ITK] Metric NaN and multi resolution registration
Luigi Riba
ribaluigi at gmail.com
Wed Dec 16 04:24:18 EST 2015
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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20151216/f922761a/attachment-0001.html>
More information about the Community
mailing list