ITK  6.0.0
Insight Toolkit
SphinxExamples/src/Registration/Common/GlobalRegistrationOfTwoImages/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 "itkImage.h"
constexpr unsigned int Dimension = 2;
using PixelType = unsigned char;
static void
CreateEllipseImage(ImageType::Pointer image);
static void
CreateSphereImage(ImageType::Pointer image);
int
main()
{
// The transform that will map the fixed image into the moving image.
// An optimizer is required to explore the parameter space of the transform
// in search of optimal values of the metric.
// The metric will compare how well the two images match each other. Metric
// types are usually parameterized by the image types as it can be seen in
// the following type declaration.
// Finally, the type of the interpolator is declared. The interpolator will
// evaluate the intensities of the moving image at non-grid positions.
// The registration method type is instantiated using the types of the
// fixed and moving images. This class is responsible for interconnecting
// all the components that we have described so far.
// Create components
auto metric = MetricType::New();
auto transform = TransformType::New();
auto optimizer = OptimizerType::New();
auto interpolator = InterpolatorType::New();
auto registration = RegistrationType::New();
// Each component is now connected to the instance of the registration method.
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
registration->SetTransform(transform);
registration->SetInterpolator(interpolator);
// Get the two images
auto fixedImage = ImageType::New();
auto movingImage = ImageType::New();
CreateSphereImage(fixedImage);
CreateEllipseImage(movingImage);
// Write the two synthetic inputs
itk::WriteImage(fixedImage, "fixed.png");
itk::WriteImage(movingImage, "moving.png");
// Set the registration inputs
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
// Initialize the transform
using ParametersType = RegistrationType::ParametersType;
ParametersType initialParameters(transform->GetNumberOfParameters());
initialParameters[0] = 0.0; // Initial offset along X
initialParameters[1] = 0.0; // Initial offset along Y
registration->SetInitialTransformParameters(initialParameters);
optimizer->SetMaximumStepLength(4.00);
optimizer->SetMinimumStepLength(0.01);
// Set a stopping criterion
optimizer->SetNumberOfIterations(200);
// Connect an observer
// auto observer = CommandIterationUpdate::New();
// optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
registration->Update();
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// The result of the registration process is an array of parameters that
// defines the spatial transformation in an unique way. This final result is
// obtained using the \code{GetLastTransformParameters()} method.
ParametersType finalParameters = registration->GetLastTransformParameters();
// In the case of the \doxygen{TranslationTransform}, there is a
// straightforward interpretation of the parameters. Each element of the
// array corresponds to a translation along one spatial dimension.
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
// The optimizer can be queried for the actual number of iterations
// performed to reach convergence. The \code{GetCurrentIteration()}
// method returns this value. A large number of iterations may be an
// indication that the maximum step length has been set too small, which
// is undesirable since it results in long computational times.
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
// The value of the image metric corresponding to the last set of parameters
// can be obtained with the \code{GetValue()} method of the optimizer.
const double bestValue = optimizer->GetValue();
// Print out results
//
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;
// It is common, as the last step of a registration task, to use the
// resulting transform to map the moving image into the fixed image space.
// This is easily done with the \doxygen{ResampleImageFilter}. Please
// refer to Section~\ref{sec:ResampleImageFilter} for details on the use
// of this filter. First, a ResampleImageFilter type is instantiated
// using the image types. It is convenient to use the fixed image type as
// the output type since it is likely that the transformed moving image
// will be compared with the fixed image.
// A resampling filter is created and the moving image is connected as its input.
auto resampler = ResampleFilterType::New();
resampler->SetInput(movingImage);
// The Transform that is produced as output of the Registration method is
// also passed as input to the resampling filter. Note the use of the
// methods \code{GetOutput()} and \code{Get()}. This combination is needed
// here because the registration method acts as a filter whose output is a
// transform decorated in the form of a \doxygen{DataObject}. For details in
// this construction you may want to read the documentation of the
// \doxygen{DataObjectDecorator}.
resampler->SetTransform(registration->GetOutput()->Get());
// As described in Section \ref{sec:ResampleImageFilter}, the
// ResampleImageFilter requires additional parameters to be specified, in
// particular, the spacing, origin and size of the output image. The default
// pixel value is also set to a distinct gray level in order to highlight
// the regions that are mapped outside of the moving image.
resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resampler->SetOutputOrigin(fixedImage->GetOrigin());
resampler->SetOutputSpacing(fixedImage->GetSpacing());
resampler->SetOutputDirection(fixedImage->GetDirection());
resampler->SetDefaultPixelValue(100);
// The output of the filter is passed to a writer that will store the
// image in a file. An \doxygen{CastImageFilter} is used to convert the
// pixel type of the resampled image to the final type used by the
// writer. The cast and writer filters are instantiated below.
auto caster = CastFilterType::New();
caster->SetInput(resampler->GetOutput());
itk::WriteImage(caster->GetOutput(), "output.png");
/*
// The fixed image and the transformed moving image can easily be compared
// using the \doxygen{SubtractImageFilter}. This pixel-wise filter computes
// the difference between homologous pixels of its two input images.
using DifferenceFilterType = itk::SubtractImageFilter<
FixedImageType,
FixedImageType,
FixedImageType >;
auto difference = DifferenceFilterType::New();
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resampler->GetOutput() );
*/
return EXIT_SUCCESS;
}
void
CreateEllipseImage(ImageType::Pointer image)
{
using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;
size[0] = 100;
size[1] = 100;
imageFilter->SetSize(size);
ImageType::SpacingType spacing;
spacing.Fill(1);
imageFilter->SetSpacing(spacing);
auto ellipse = EllipseType::New();
EllipseType::ArrayType radiusArray;
radiusArray[0] = 10;
radiusArray[1] = 20;
ellipse->SetRadiusInObjectSpace(radiusArray);
using TransformType = EllipseType::TransformType;
auto transform = TransformType::New();
transform->SetIdentity();
TransformType::OutputVectorType translation;
translation[0] = 65;
translation[1] = 45;
transform->Translate(translation, false);
ellipse->SetObjectToParentTransform(transform);
imageFilter->SetInput(ellipse);
ellipse->SetDefaultInsideValue(255);
ellipse->SetDefaultOutsideValue(0);
imageFilter->SetUseObjectValue(true);
imageFilter->SetOutsideValue(0);
imageFilter->Update();
image->Graft(imageFilter->GetOutput());
}
void
CreateSphereImage(ImageType::Pointer image)
{
using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;
size[0] = 100;
size[1] = 100;
imageFilter->SetSize(size);
ImageType::SpacingType spacing;
spacing.Fill(1);
imageFilter->SetSpacing(spacing);
auto ellipse = EllipseType::New();
EllipseType::ArrayType radiusArray;
radiusArray[0] = 10;
radiusArray[1] = 10;
ellipse->SetRadiusInObjectSpace(radiusArray);
using TransformType = EllipseType::TransformType;
auto transform = TransformType::New();
transform->SetIdentity();
TransformType::OutputVectorType translation;
translation[0] = 50;
translation[1] = 50;
transform->Translate(translation, false);
ellipse->SetObjectToParentTransform(transform);
imageFilter->SetInput(ellipse);
ellipse->SetDefaultInsideValue(255);
ellipse->SetDefaultOutsideValue(0);
imageFilter->SetUseObjectValue(true);
imageFilter->SetOutsideValue(0);
imageFilter->Update();
image->Graft(imageFilter->GetOutput());
}
Pointer
SmartPointer< Self > Pointer
Definition: itkAddImageFilter.h:93
itk::CastImageFilter
Casts input pixels to output pixel type.
Definition: itkCastImageFilter.h:100
itkLinearInterpolateImageFunction.h
itkEllipseSpatialObject.h
itkRegularStepGradientDescentOptimizer.h
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itkImage.h
itk::ImageRegistrationMethod
Base class for Image Registration Methods.
Definition: itkImageRegistrationMethod.h:70
itk::SpatialObjectToImageFilter
Base class for filters that take a SpatialObject as input and produce an image as output....
Definition: itkSpatialObjectToImageFilter.h:41
itkCastImageFilter.h
itk::EllipseSpatialObject
Definition: itkEllipseSpatialObject.h:38
itk::RegularStepGradientDescentOptimizer
Implement a gradient descent optimizer.
Definition: itkRegularStepGradientDescentOptimizer.h:33
itkTranslationTransform.h
itkMeanSquaresImageToImageMetric.h
itk::LinearInterpolateImageFunction
Linearly interpolate an image at specified positions.
Definition: itkLinearInterpolateImageFunction.h:51
itkSpatialObjectToImageFilter.h
itk::TranslationTransform
Translation transformation of a vector space (e.g. space coordinates)
Definition: itkTranslationTransform.h:43
itkImageRegistrationMethod.h
itkRescaleIntensityImageFilter.h
itk::MeanSquaresImageToImageMetric
TODO.
Definition: itkMeanSquaresImageToImageMetric.h:41
itkImageFileWriter.h
itk::Size::SetSize
void SetSize(const SizeValueType val[VDimension])
Definition: itkSize.h:179
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
New
static Pointer New()
itkResampleImageFilter.h
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itk::WriteImage
ITK_TEMPLATE_EXPORT void WriteImage(TImagePointer &&image, const std::string &filename, bool compress=false)
Definition: itkImageFileWriter.h:256