ITK  6.0.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
*
* https://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;
InputImageType::Pointer fixedImage = itk::ReadImage<InputImageType>(fixedImageFile);
InputImageType::Pointer movingImage = itk::ReadImage<InputImageType>(movingImageFile);
auto affineTransform = AffineTransformType::New();
using DisplacementTransformType =
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
auto identityTransform = IdentityTransformType::New();
// The metric
using VirtualImageType = itk::Image<double, Dimension>;
using MetricTraitsType =
using MetricType =
using PointSetType = MetricType::FixedSampledPointSetType;
auto metric = MetricType::New();
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)
{
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>;
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;
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
auto resampler = ResampleFilterType::New();
resampler->SetTransform(compositeTransform);
resampler->SetInput(movingImage);
resampler->SetReferenceImage(fixedImage);
resampler->SetUseReferenceImage(true);
// Write the warped image into a file
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
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;
}
Pointer
SmartPointer< Self > Pointer
Definition: itkAddImageFilter.h:93
itk::CastImageFilter
Casts input pixels to output pixel type.
Definition: itkCastImageFilter.h:100
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::GaussianSmoothingOnUpdateDisplacementFieldTransform
Modifies the UpdateTransformParameters method to perform a Gaussian smoothing of the displacement fie...
Definition: itkGaussianSmoothingOnUpdateDisplacementFieldTransform.h:46
itk::point
*par Constraints *The filter requires an image with at least two dimensions and a vector *length of at least The theory supports extension to scalar but *the implementation of the itk vector classes do not **The template parameter TRealType must be floating point(float or double) or *a user-defined "real" numerical type with arithmetic operations defined *sufficient to compute derivatives. **\par Performance *This filter will automatically multithread if run with *SetUsePrincipleComponents
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:88
New
static Pointer New()
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkCommand.h
itk::WriteImage
ITK_TEMPLATE_EXPORT void WriteImage(TImagePointer &&image, const std::string &filename, bool compress=false)
Definition: itkImageFileWriter.h:256
itk::SizeValueType
unsigned long SizeValueType
Definition: itkIntTypes.h:86