ITK  5.2.0
Insight Toolkit
SphinxExamples/src/Registration/Metricsv4/PerformRegistrationOnVectorImages/Code.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "itkCommand.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 InputImageType = itk::Image<InputPixelType, Dimension>;
using ParametersValueType = double;
ReaderType::Pointer fixedImageReader = ReaderType::New();
fixedImageReader->SetFileName(fixedImageFile);
ReaderType::Pointer movingImageReader = ReaderType::New();
movingImageReader->SetFileName(movingImageFile);
// Get the input images
try
{
fixedImageReader->Update();
movingImageReader->Update();
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error while reading images: " << error << std::endl;
return EXIT_FAILURE;
}
InputImageType::Pointer fixedImage = fixedImageReader->GetOutput();
InputImageType::Pointer movingImage = movingImageReader->GetOutput();
AffineTransformType::Pointer affineTransform = AffineTransformType::New();
using DisplacementTransformType =
DisplacementTransformType::Pointer displacementTransform = DisplacementTransformType::New();
using DisplacementFieldType = DisplacementTransformType::DisplacementFieldType;
DisplacementFieldType::Pointer 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
IdentityTransformType::Pointer identityTransform = IdentityTransformType::New();
// The metric
using VirtualImageType = itk::Image<double, Dimension>;
using MetricTraitsType =
using MetricType =
using PointSetType = MetricType::FixedSampledPointSetType;
MetricType::Pointer metric = MetricType::New();
PointSetType::Pointer 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;
OptimizerType::Pointer 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;
CompositeType::Pointer 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 (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
// Warp the image with the displacement field
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform(compositeTransform);
resampler->SetInput(movingImageReader->GetOutput());
resampler->SetReferenceImage(fixedImage);
resampler->SetUseReferenceImage(true);
// Write the warped image into a file
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput(resampler->GetOutput());
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(resampledMovingImageFile);
writer->SetInput(caster->GetOutput());
try
{
writer->UpdateLargestPossibleRegion();
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
itk::CastImageFilter
Casts input pixels to output pixel type.
Definition: itkCastImageFilter.h:104
itk::CompositeTransform
This class contains a list of transforms and concatenates them by composition.
Definition: itkCompositeTransform.h:87
itk::IdentityTransform
Implementation of an Identity Transform.
Definition: itkIdentityTransform.h:50
itk::GTest::TypedefsAndConstructors::Dimension2::PointType
ImageBaseType::PointType PointType
Definition: itkGTestTypedefsAndConstructors.h:51
itk::Vector
A templated class holding a n-Dimensional vector.
Definition: itkVector.h:62
itkImageFileReader.h
itkCastImageFilter.h
itk::AffineTransform
Definition: itkAffineTransform.h:101
itkMeanSquaresImageToImageMetricv4.h
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::GaussianSmoothingOnUpdateDisplacementFieldTransform
Modifies the UpdateTransformParameters method to peform a Gaussian smoothing of the displacement fiel...
Definition: itkGaussianSmoothingOnUpdateDisplacementFieldTransform.h:46
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itkRegistrationParameterScalesFromPhysicalShift.h
itkGradientDescentOptimizerv4.h
itk::GradientDescentOptimizerv4
GradientDescentOptimizerv4Template< double > GradientDescentOptimizerv4
Definition: itkGradientDescentOptimizerv4.h:255
itkVectorImageToImageMetricTraitsv4.h
itk::ImageRegionIteratorWithIndex
A multi-dimensional iterator templated over image type that walks pixels within a region and is speci...
Definition: itkImageRegionIteratorWithIndex.h:73
itkImageFileWriter.h
itkGaussianSmoothingOnUpdateDisplacementFieldTransform.h
itk::MeanSquaresImageToImageMetricv4
Class implementing a mean squares metric.
Definition: itkMeanSquaresImageToImageMetricv4.h:46
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::RegistrationParameterScalesFromPhysicalShift
Registration helper class for estimating scales of transform parameters a step sizes,...
Definition: itkRegistrationParameterScalesFromPhysicalShift.h:35
itk::VectorImageToImageMetricTraitsv4
A simple structure holding type information for ImageToImageMetricv4 classes.
Definition: itkVectorImageToImageMetricTraitsv4.h:47
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkCommand.h
itk::SizeValueType
unsigned long SizeValueType
Definition: itkIntTypes.h:83