I want to combine two similarity measures in ITK. What should I do?
Evaluate both separately, then combine them using your formula. E.g.
auto metric1 = ...
auto metric2 = ...
double combinedMetric = metric1*0.2 + metric2*0.8;
I tried unsuccessfully. Do I have to write mymetric Mymetric.h and Mymetric.hxx integrates the two similarity measures?
I don’t think you have to write a new metric class. There is ObjectToObjectMultiMetricv4 which is worth trying if a simple combination didn’t work.
Can you share your attempt?
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
/*using OptimizerType = itk::PowellOptimizer;*/
using OptimizerType = itk::GradientDescentOptimizer;
using OptimizerPointer = const OptimizerType*;
void
Execute(itk::Object* caller, const itk::EventObject& event) override
{
Execute((const itk::Object*)caller, event);
}
void
Execute(const itk::Object* object, const itk::EventObject& event) override
{
auto optimizer = dynamic_cast<OptimizerPointer>(object);
auto optimizer1 = dynamic_cast<OptimizerPointer>(object);
if (!itk::IterationEvent().CheckEvent(&event))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue()*0.5 + optimizer1->GetValue() * 0.5 << " ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
}
};
int
main(int argc, char* argv[])
{
const char* fixedImageFile ="BrainProtonDensitySliceShifted13x17y.png";
const char* movingImageFile = "BrainT1SliceBorder20.png";
const char* outputImageFile = "output.png";
const char* checkerBoardBefore = "before.png";
const char* checkerBoardAfter = "after.png";
constexpr unsigned int Dimension = 2;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
// It is convenient to work with an internal image type because mutual
// information will perform better on images with a normalized statistical
// distribution. The fixed and moving images will be normalized and
// converted to this internal type.
using InternalPixelType = float;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
using TransformType = itk::TranslationTransform<double, Dimension>;
using OptimizerType = itk::GradientDescentOptimizer;
/* using OptimizerType = itk::PowellOptimizer;*/
using InterpolatorType = itk::LinearInterpolateImageFunction<InternalImageType, double>;
using RegistrationType = itk::ImageRegistrationMethod<InternalImageType, InternalImageType>;
using RegistrationType1 = itk::ImageRegistrationMethod<InternalImageType, InternalImageType>;
using MetricType = itk::MutualInformationImageToImageMetric<InternalImageType, InternalImageType>;
using MetricType1 = itk::GradientDifferenceImageToImageMetric<InternalImageType, InternalImageType>;
/* using MetricType = itk::NormalizedMutualInformationImageToImageMetric<InternalImageType, InternalImageType>;*/
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
OptimizerType::Pointer optimizer1 = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
RegistrationType1::Pointer registration1 = RegistrationType::New();
MetricType::Pointer metric = MetricType::New();
MetricType1::Pointer metric1 = MetricType1::New();
registration->SetOptimizer(optimizer);
registration->SetTransform(transform);
registration->SetInterpolator(interpolator);
registration->SetMetric(metric);
registration1->SetOptimizer(optimizer);
registration1->SetTransform(transform);
registration1->SetInterpolator(interpolator);
metric1->SetDerivativeDelta(0.5);
registration1->SetMetric(metric1);
// The metric requires a number of parameters to be selected, including
// the standard deviation of the Gaussian kernel for the fixed image
// density estimate, the standard deviation of the kernel for the moving
// image density and the number of samples use to compute the densities
// and entropy values. Experience has
// shown that a kernel standard deviation of 0.4 works well for images
// which have been normalized to a mean of zero and unit variance. We
// will follow this empirical rule in this example.
metric->SetFixedImageStandardDeviation(0.4);
metric->SetMovingImageStandardDeviation(0.4);
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName(movingImageFile);
movingImageReader->SetFileName(movingImageFile);
registration1->SetFixedImage(fixedImageReader->GetOutput());
registration1->SetMovingImage(movingImageReader->GetOutput());
registration1->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion());
using FixedNormalizeFilterType = itk::NormalizeImageFilter<FixedImageType, InternalImageType>;
using MovingNormalizeFilterType = itk::NormalizeImageFilter<MovingImageType, InternalImageType>;
FixedNormalizeFilterType::Pointer fixedNormalizer = FixedNormalizeFilterType::New();
MovingNormalizeFilterType::Pointer movingNormalizer = MovingNormalizeFilterType::New();
using GaussianFilterType = itk::DiscreteGaussianImageFilter<InternalImageType, InternalImageType>;
GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
fixedSmoother->SetVariance(2.0);
movingSmoother->SetVariance(2.0);
fixedNormalizer->SetInput(fixedImageReader->GetOutput());
movingNormalizer->SetInput(movingImageReader->GetOutput());
fixedSmoother->SetInput(fixedNormalizer->GetOutput());
movingSmoother->SetInput(movingNormalizer->GetOutput());
registration->SetFixedImage(fixedSmoother->GetOutput());
registration->SetMovingImage(movingSmoother->GetOutput());
fixedNormalizer->Update();
FixedImageType::RegionType fixedImageRegion = fixedNormalizer->GetOutput()->GetBufferedRegion();
registration->SetFixedImageRegion(fixedImageRegion);
using ParametersType = RegistrationType::ParametersType;
ParametersType initialParameters(transform->GetNumberOfParameters());
initialParameters[0] = 1.0; // Initial offset in mm along X
initialParameters[1] = 1.0; // Initial offset in mm along Y
registration->SetInitialTransformParameters(initialParameters);
registration1->SetInitialTransformParameters(initialParameters);
// We should now define the number of spatial samples to be considered in
// the metric computation. Note that we were forced to postpone this setting
// until we had done the preprocessing of the images because the number of
// samples is usually defined as a fraction of the total number of pixels in
// the fixed image.
//
// The number of spatial samples can usually be as low as $1\%$ of the total
// number of pixels in the fixed image. Increasing the number of samples
// improves the smoothness of the metric from one iteration to another and
// therefore helps when this metric is used in conjunction with optimizers
// that rely of the continuity of the metric values. The trade-off, of
// course, is that a larger number of samples result in longer computation
// times per every evaluation of the metric.
//
// It has been demonstrated empirically that the number of samples is not a
// critical parameter for the registration process. When you start fine
// tuning your own registration process, you should start using high values
// of number of samples, for example in the range of 20% to 50% of the
// number of pixels in the fixed image. Once you have succeeded to register
// your images you can then reduce the number of samples progressively until
// you find a good compromise on the time it takes to compute one evaluation
// of the Metric. Note that it is not useful to have very fast evaluations
// of the Metric if the noise in their values results in more iterations
// being required by the optimizer to converge.
// behavior of the metric values as the iterations progress.
const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
const auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.01);
metric->SetNumberOfSpatialSamples(numberOfSamples);
// For consistent results when regression testing.
metric->ReinitializeSeed(121212);
// Since larger values of mutual information indicate better matches than
// smaller values, we need to maximize the cost function in this example.
// By default the GradientDescentOptimizer class is set to minimize the
// value of the cost-function. It is therefore necessary to modify its
// default behavior by invoking the MaximizeOn() method.
// Additionally, we need to define the optimizer's step size using the
// SetLearningRate() method.
optimizer->SetLearningRate(15.0);
optimizer->SetNumberOfIterations(200);
optimizer->MaximizeOn();
optimizer1->SetLearningRate(15.0);
optimizer1->SetNumberOfIterations(200);
optimizer1->MaximizeOn();
// Note that large values of the learning rate will make the optimizer
// unstable. Small values, on the other hand, may result in the optimizer
// needing too many iterations in order to walk to the extrema of the cost
// function. The easy way of fine tuning this parameter is to start with
// small values, probably in the range of {5.0, 10.0}. Once the other
// registration parameters have been tuned for producing convergence, you
// may want to revisit the learning rate and start increasing its value until
// you observe that the optimization becomes unstable. The ideal value for
// this parameter is the one that results in a minimum number of iterations
// while still keeping a stable path on the parametric space of the
// optimization. Keep in mind that this parameter is a multiplicative factor
// applied on the gradient of the Metric. Therefore, its effect on the
// optimizer step length is proportional to the Metric values themselves.
// Metrics with large values will require you to use smaller values for the
// learning rate in order to maintain a similar optimizer behavior.
/* optimizer->SetLearningRate(15.0);*/
/* optimizer->SetMaximize(true);*/ // for NCC
/* optimizer->MaximizeOn();*/
//optimizer->SetMaximumIteration(100);
//optimizer->SetMaximumLineIteration(4); // for Powell's method
//optimizer->SetStepLength(0.5);
//optimizer->SetStepTolerance(0.02);
//optimizer->SetValueTolerance(0.001);
/* itk::Optimizer::ScalesType weightings(transform->GetNumberOfParameters());
weightings[0] = 1.;
weightings[1] = 1.;
optimizer->SetScales(weightings);*/
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
optimizer1->AddObserver(itk::IterationEvent(), observer);
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;
}
ParametersType finalParameters = registration->GetLastTransformParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = (optimizer->GetValue()*0.5+optimizer1->GetValue()*0.5);
// Print out results
std::cout << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
std::cout << " Numb. Samples = " << numberOfSamples << std::endl;
using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
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();
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType = itk::CastImageFilter<FixedImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName(outputImageFile);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
// Generate checkerboards before and after registration
using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
// Before registration
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
/* if (argc > 4)*/
writer->SetFileName(checkerBoardBefore);
writer->Update();
// After registration
resample->SetTransform(finalTransform);
/*if (argc > 5)*/
writer->SetFileName(checkerBoardAfter);
try
{
writer->Update();
}
catch (itk::ExceptionObject& error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
system("pause");
return EXIT_SUCCESS;
}
I did this, and the result seems wrong
Well, you failed to mention that you want to do registration using a combination of 2 metrics. For that, you indeed need a single metric, even if it is internally computed using two metrics. Using ObjectToObjectMultiMetricv4 is probably meant for this case.
I used it, but the following error occurred
My code:
template <typename TFilter>
class itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate : public itk::Command
{
public:
using Self = itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate() = default;
public:
void
Execute(itk::Object* caller, const itk::EventObject& event) override
{
Execute((const itk::Object*)caller, event);
}
void
Execute(const itk::Object* object, const itk::EventObject& event) override
{
if (typeid(event) != typeid(itk::IterationEvent))
{
return;
}
const auto* optimizer = dynamic_cast<const TFilter*>(object);
if (!optimizer)
{
itkGenericExceptionMacro("Error dynamic_cast failed");
}
std::cout << "It- " << optimizer->GetCurrentIteration() << " gradient: " << optimizer->GetGradient()
<< " metric value: " << optimizer->GetCurrentMetricValue()
<< " Params: " << const_cast<TFilter*>(optimizer)->GetCurrentPosition() << std::endl;
}
};
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
int
main()
{
constexpr int Dimension = 2;
using ImageType = itk::Image<float, Dimension>;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName("BrainProtonDensitySliceShifted13x17y.png");
movingImageReader->SetFileName("BrainT1SliceBorder20.png");
auto fixedImage = fixedImageReader->GetOutput();
auto movingImage = movingImageReader->GetOutput();
int numberOfIterations = 50;
/* if (argc > 1)
{
numberOfIterations = std::stoi(argv[1]);
}*/
// create an affine transform
using TranslationTransformType = itk::TranslationTransform<double, Dimension>;
auto translationTransform = TranslationTransformType::New();
translationTransform->SetIdentity();
using OptimizerType = itk::GradientDescentOptimizerv4;
auto optimizer = OptimizerType::New();
using RegistrationType = itk::ImageRegistrationMethodv4<ImageType, ImageType>;
auto registration = RegistrationType::New();
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
using MultiMetricType = itk::ObjectToObjectMultiMetricv4<Dimension, Dimension, ImageType>;
using MetricType = itk::CorrelationImageToImageMetricv4<ImageType, ImageType>;
auto correlationMetric = MetricType::New();
correlationMetric->SetFixedImage(fixedImage);
correlationMetric->SetMovingImage(movingImage);
correlationMetric->SetMovingTransform(translationTransform);
correlationMetric->Initialize();
/*translationTransform->SetIdentity();*/
// Test with two different metric types
/* using MeanSquaresMetricType = itk::MeanSquaresImageToImageMetricv4<ImageType, ImageType>;
auto meanSquaresMetric = MeanSquaresMetricType::New();
meanSquaresMetric->SetFixedImage(fixedImage);
meanSquaresMetric->SetMovingImage(movingImage);
meanSquaresMetric->SetMovingTransform(translationTransform);*/
using MattesMutualInformationMetricType = itk::MattesMutualInformationImageToImageMetricv4<ImageType,ImageType>;
auto MattesMutualInformationMetric = MattesMutualInformationMetricType::New();
MattesMutualInformationMetric->SetFixedImage(fixedImage);
MattesMutualInformationMetric->SetMovingImage(movingImage);
MattesMutualInformationMetric->SetMovingTransform(translationTransform);
MattesMutualInformationMetric->Initialize();
auto multiMetric2 = MultiMetricType::New();
multiMetric2->AddMetric(correlationMetric);
multiMetric2->AddMetric(MattesMutualInformationMetric);
multiMetric2->Initialize();
translationTransform->SetIdentity();
MultiMetricType::WeightsArrayType oriMetricWeights(2);
oriMetricWeights[0] = 0.5;
oriMetricWeights[1] = 0.5;
multiMetric2->SetMetricWeights(oriMetricWeights);
std::cout << "*** Multi-metric with different metric types: " << std::endl;
using ParametersType = TranslationTransformType::ParametersType;
ParametersType initialParameters(translationTransform->GetNumberOfParameters());
initialParameters[0] = 1.0; // Initial offset in mm along X
initialParameters[1] = 1.0; // Initial offset in mm along Y
translationTransform->SetParameters(initialParameters);
registration->SetInitialTransform(translationTransform);
registration->SetOptimizer(optimizer);
registration->SetMetric(multiMetric2);
optimizer->SetMetric(multiMetric2);
optimizer->SetNumberOfIterations(numberOfIterations);
optimizer->SetLearningRate(0.01);
optimizer->SetMaximumStepSizeInPhysicalUnits(1.0);
using CommandType = itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate<OptimizerType>;
auto observer = CommandType::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
optimizer->StartOptimization();
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;
system("pause");
return EXIT_FAILURE;
}
using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
OptimizerType::ParametersType finalParameters = translationTransform->GetParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
unsigned int numberOfIteration = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();
std::cout << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIteration << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
TranslationTransformType::Pointer finalTransform = TranslationTransformType::New();
finalTransform->SetParameters(finalParameters);
finalTransform->SetFixedParameters(translationTransform->GetFixedParameters());
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(finalTransform);
resample->SetInput(movingImageReader->GetOutput());
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType = itk::CastImageFilter<FixedImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName("out.png");
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
system("pause");
return EXIT_SUCCESS;
}
You came very close.
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
needed to be changed to
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
registration->SetFixedImage(1, fixedImage);
registration->SetMovingImage(1, movingImage);
There needs to be a separate pair of inputs for each metric in the multi-metric, as the different metrics don’t necessarily take the same inputs (e.g. first metric takes images, the second metric takes point sets). This should probably be noted in the documentation. Do you mind contributing that?
For posterity, here is a complete example which runs on my computer:
#include "itkImageRegistrationMethodv4.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCorrelationImageToImageMetricv4.h"
template <typename TFilter>
class itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate : public itk::Command
{
public:
using Self = itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate() = default;
public:
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
if (typeid(event) != typeid(itk::IterationEvent))
{
return;
}
const auto * optimizer = dynamic_cast<const TFilter *>(object);
if (!optimizer)
{
itkGenericExceptionMacro("Error dynamic_cast failed");
}
std::cout << "It- " << optimizer->GetCurrentIteration() << " gradient: " << optimizer->GetGradient()
<< " metric value: " << optimizer->GetCurrentMetricValue()
<< " Params: " << const_cast<TFilter *>(optimizer)->GetCurrentPosition() << std::endl;
}
};
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
int
main()
{
constexpr int Dimension = 2;
using ImageType = itk::Image<float, Dimension>;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName("C:/Dev/ITK-git/Examples/Data/BrainProtonDensitySliceShifted13x17y.png");
movingImageReader->SetFileName("C:/Dev/ITK-git/Examples/Data/BrainT1SliceBorder20.png");
auto fixedImage = fixedImageReader->GetOutput();
auto movingImage = movingImageReader->GetOutput();
int numberOfIterations = 50;
/* if (argc > 1)
{
numberOfIterations = std::stoi(argv[1]);
}*/
// create an affine transform
using TranslationTransformType = itk::TranslationTransform<double, Dimension>;
auto translationTransform = TranslationTransformType::New();
translationTransform->SetIdentity();
using OptimizerType = itk::GradientDescentOptimizerv4;
auto optimizer = OptimizerType::New();
using RegistrationType = itk::ImageRegistrationMethodv4<ImageType, ImageType>;
auto registration = RegistrationType::New();
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
registration->SetFixedImage(1, fixedImage);
registration->SetMovingImage(1, movingImage);
using MultiMetricType = itk::ObjectToObjectMultiMetricv4<Dimension, Dimension, ImageType>;
using MetricType = itk::CorrelationImageToImageMetricv4<ImageType, ImageType>;
auto correlationMetric = MetricType::New();
correlationMetric->SetFixedImage(fixedImage);
correlationMetric->SetMovingImage(movingImage);
correlationMetric->SetMovingTransform(translationTransform);
correlationMetric->Initialize();
/*translationTransform->SetIdentity();*/
// Test with two different metric types
/* using MeanSquaresMetricType = itk::MeanSquaresImageToImageMetricv4<ImageType, ImageType>;
auto meanSquaresMetric = MeanSquaresMetricType::New();
meanSquaresMetric->SetFixedImage(fixedImage);
meanSquaresMetric->SetMovingImage(movingImage);
meanSquaresMetric->SetMovingTransform(translationTransform);*/
using MattesMutualInformationMetricType = itk::MattesMutualInformationImageToImageMetricv4<ImageType, ImageType>;
auto MattesMutualInformationMetric = MattesMutualInformationMetricType::New();
MattesMutualInformationMetric->SetFixedImage(fixedImage);
MattesMutualInformationMetric->SetMovingImage(movingImage);
MattesMutualInformationMetric->SetMovingTransform(translationTransform);
MattesMutualInformationMetric->Initialize();
auto multiMetric2 = MultiMetricType::New();
multiMetric2->AddMetric(correlationMetric);
multiMetric2->AddMetric(MattesMutualInformationMetric);
multiMetric2->Initialize();
translationTransform->SetIdentity();
MultiMetricType::WeightsArrayType oriMetricWeights(2);
oriMetricWeights[0] = 0.5;
oriMetricWeights[1] = 0.5;
multiMetric2->SetMetricWeights(oriMetricWeights);
std::cout << "*** Multi-metric with different metric types: " << std::endl;
using ParametersType = TranslationTransformType::ParametersType;
ParametersType initialParameters(translationTransform->GetNumberOfParameters());
initialParameters[0] = 1.0; // Initial offset in mm along X
initialParameters[1] = 1.0; // Initial offset in mm along Y
translationTransform->SetParameters(initialParameters);
registration->SetInitialTransform(translationTransform);
registration->SetOptimizer(optimizer);
registration->SetMetric(multiMetric2);
optimizer->SetMetric(multiMetric2);
optimizer->SetNumberOfIterations(numberOfIterations);
optimizer->SetLearningRate(0.01);
optimizer->SetMaximumStepSizeInPhysicalUnits(1.0);
using CommandType = itkObjectToObjectMultiMetricv4RegistrationTestCommandIterationUpdate<OptimizerType>;
auto observer = CommandType::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
//optimizer->StartOptimization();
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;
system("pause");
return EXIT_FAILURE;
}
using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
OptimizerType::ParametersType finalParameters = translationTransform->GetParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
unsigned int numberOfIteration = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();
std::cout << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIteration << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
TranslationTransformType::Pointer finalTransform = TranslationTransformType::New();
finalTransform->SetParameters(finalParameters);
finalTransform->SetFixedParameters(translationTransform->GetFixedParameters());
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(finalTransform);
resample->SetInput(movingImageReader->GetOutput());
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType = itk::CastImageFilter<FixedImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName("out.png");
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
system("pause");
return EXIT_SUCCESS;
}
A post was split to a new topic: How to tune registration
This code will only give the values of the first metric during each iteration, right?
itk::ObjectToObjectMultiMetricv4, in this documentation its written that i can use GetWeightedValue() to get the weighted metric value. But multiMetric2->GetWeightedValue() will only give final metric value.
Is there any way to get metric values of both the metrics and weighted metric value during each iteration?