[ITK] Metric NaN and multi resolution registration
Matias Montroull
matimontg at gmail.com
Mon Dec 7 10:06:35 EST 2015
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/20151207/a2f36de9/attachment-0001.html>
More information about the Community
mailing list