[ITK] Metric NaN and multi resolution registration
Matias Montroull
matimontg at gmail.com
Mon Dec 7 10:04:31 EST 2015
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20151207/c2c55f18/attachment-0001.html>
More information about the Community
mailing list