Perform Registration on Vector Images#
Synopsis#
Register images that have multi-component pixels like an itk::VectorImage or an itk::Image< itk::Vector, Dimension >.
Results#
Fixed image. |
Moving image. |
Resampled, registered moving image. |
Code#
C++#
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkGradientDescentOptimizerv4.h"
#include "itkRegistrationParameterScalesFromPhysicalShift.h"
#include "itkVectorImageToImageMetricTraitsv4.h"
#include "itkGaussianSmoothingOnUpdateDisplacementFieldTransform.h"
#include "itkCastImageFilter.h"
#include "itkCommand.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <iomanip>
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " resampledMovingImageFile ";
std::cerr << " [numberOfAffineIterations numberOfDisplacementIterations] ";
std::cerr << std::endl;
return EXIT_FAILURE;
}
const char * fixedImageFile = argv[1];
const char * movingImageFile = argv[2];
const char * resampledMovingImageFile = argv[3];
unsigned int numberOfAffineIterations = 2;
unsigned int numberOfDisplacementIterations = 2;
if (argc >= 5)
{
numberOfAffineIterations = std::stoi(argv[4]);
}
if (argc >= 6)
{
numberOfDisplacementIterations = std::stoi(argv[5]);
}
constexpr unsigned int Dimension = 2;
// The input images have red, blue, and green pixel components.
constexpr unsigned int NumberOfPixelComponents = 3;
using PixelComponentType = float;
using InputPixelType = itk::Vector<PixelComponentType, NumberOfPixelComponents>;
using InputImageType = itk::Image<InputPixelType, Dimension>;
using ParametersValueType = double;
InputImageType::Pointer fixedImage = itk::ReadImage<InputImageType>(fixedImageFile);
InputImageType::Pointer movingImage = itk::ReadImage<InputImageType>(movingImageFile);
using AffineTransformType = itk::AffineTransform<ParametersValueType, Dimension>;
auto affineTransform = AffineTransformType::New();
using DisplacementTransformType =
itk::GaussianSmoothingOnUpdateDisplacementFieldTransform<ParametersValueType, Dimension>;
auto displacementTransform = DisplacementTransformType::New();
using DisplacementFieldType = DisplacementTransformType::DisplacementFieldType;
auto displacementField = DisplacementFieldType::New();
// Set the displacement field to be the same as the fixed image region, which will
// act by default as the virtual domain in this example.
displacementField->SetRegions(fixedImage->GetLargestPossibleRegion());
// make sure the displacement field has the same spatial information as the image
displacementField->CopyInformation(fixedImage);
std::cout << "fixedImage->GetLargestPossibleRegion(): " << fixedImage->GetLargestPossibleRegion() << std::endl;
displacementField->Allocate();
// Fill it with 0's
DisplacementTransformType::OutputVectorType zeroVector;
zeroVector.Fill(0.0);
displacementField->FillBuffer(zeroVector);
// Assign to transform
displacementTransform->SetDisplacementField(displacementField);
displacementTransform->SetGaussianSmoothingVarianceForTheUpdateField(5.0);
displacementTransform->SetGaussianSmoothingVarianceForTheTotalField(6.0);
// Identity transform for fixed image
using IdentityTransformType = itk::IdentityTransform<ParametersValueType, Dimension>;
auto identityTransform = IdentityTransformType::New();
// The metric
using VirtualImageType = itk::Image<double, Dimension>;
using MetricTraitsType =
itk::VectorImageToImageMetricTraitsv4<InputImageType, InputImageType, VirtualImageType, InputPixelType::Length>;
using MetricType =
itk::MeanSquaresImageToImageMetricv4<InputImageType, InputImageType, VirtualImageType, double, MetricTraitsType>;
using PointSetType = MetricType::FixedSampledPointSetType;
auto metric = MetricType::New();
using PointType = PointSetType::PointType;
auto pointSet = PointSetType::New();
itk::ImageRegionIteratorWithIndex<InputImageType> fixedImageIt(fixedImage, fixedImage->GetLargestPossibleRegion());
itk::SizeValueType index = 0;
itk::SizeValueType count = 0;
for (fixedImageIt.GoToBegin(); !fixedImageIt.IsAtEnd(); ++fixedImageIt)
{
// take every N^th point
if (count % 2 == 0)
{
PointType point;
fixedImage->TransformIndexToPhysicalPoint(fixedImageIt.GetIndex(), point);
pointSet->SetPoint(index, point);
++index;
}
++count;
}
metric->SetFixedSampledPointSet(pointSet);
metric->SetUseSampledPointSet(true);
// Assign images and transforms.
// By not setting a virtual domain image or virtual domain settings,
// the metric will use the fixed image for the virtual domain.
metric->SetFixedImage(fixedImage);
metric->SetMovingImage(movingImage);
metric->SetFixedTransform(identityTransform);
metric->SetMovingTransform(affineTransform);
const bool gaussian = false;
metric->SetUseMovingImageGradientFilter(gaussian);
metric->SetUseFixedImageGradientFilter(gaussian);
metric->Initialize();
using RegistrationParameterScalesFromShiftType = itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
RegistrationParameterScalesFromShiftType::Pointer shiftScaleEstimator =
RegistrationParameterScalesFromShiftType::New();
shiftScaleEstimator->SetMetric(metric);
//
// Affine registration
//
std::cout << "First do an affine registration..." << std::endl;
using OptimizerType = itk::GradientDescentOptimizerv4;
auto optimizer = OptimizerType::New();
optimizer->SetMetric(metric);
optimizer->SetNumberOfIterations(numberOfAffineIterations);
optimizer->SetScalesEstimator(shiftScaleEstimator);
optimizer->StartOptimization();
std::cout << "After optimization affine parameters are: " << affineTransform->GetParameters() << std::endl;
//
// Deformable registration
//
std::cout << "Now, add the displacement field to the composite transform..." << std::endl;
using CompositeType = itk::CompositeTransform<ParametersValueType, Dimension>;
auto compositeTransform = CompositeType::New();
compositeTransform->AddTransform(affineTransform);
compositeTransform->AddTransform(displacementTransform);
// Optimize only the displacement field, but still using the
// previously computed affine transformation.
compositeTransform->SetOnlyMostRecentTransformToOptimizeOn();
metric->SetMovingTransform(compositeTransform);
metric->SetUseSampledPointSet(false);
metric->Initialize();
// Optimizer
optimizer->SetScalesEstimator(shiftScaleEstimator);
optimizer->SetMetric(metric);
optimizer->SetNumberOfIterations(numberOfDisplacementIterations);
try
{
if (numberOfDisplacementIterations > 0)
{
optimizer->StartOptimization();
}
else
{
std::cout << "*** Skipping displacement field optimization.\n";
}
}
catch (const itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
// Warp the image with the displacement field
using ResampleFilterType = itk::ResampleImageFilter<InputImageType, InputImageType>;
auto resampler = ResampleFilterType::New();
resampler->SetTransform(compositeTransform);
resampler->SetInput(movingImage);
resampler->SetReferenceImage(fixedImage);
resampler->SetUseReferenceImage(true);
// Write the warped image into a file
using OutputPixelType = itk::Vector<unsigned char, NumberOfPixelComponents>;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType = itk::CastImageFilter<InputImageType, OutputImageType>;
auto caster = CastFilterType::New();
caster->SetInput(resampler->GetOutput());
try
{
itk::WriteImage(caster->GetOutput(), resampledMovingImageFile);
}
catch (const itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
Classes demonstrated#
-
template<typename TFixedImageType, typename TMovingImageType, typename TVirtualImageType, unsigned int NumberOfComponents, typename TCoordRep = typename ObjectToObjectMetricBase::CoordinateRepresentationType>
class VectorImageToImageMetricTraitsv4 A simple structure holding type information for ImageToImageMetricv4 classes.
This class provides type information for class members and methods used in gradient calculation. This class is used for images with vector pixel types, including VectorImage. For images with scalar pixel types, see itkDefaultImageToImageMetricTraitsv4.
- See
itkDefaultImageToImageMetricTraitsv4
