ITK  5.2.0
Insight Toolkit
Examples/RegistrationITKv4/BSplineWarping2.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.
*
*=========================================================================*/
// Software Guide : BeginLatex
//
// This example illustrates how to deform a 3D image using a
// BSplineTransform.
//
// \index{BSplineTransform}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
#include <fstream>
// The following section of code implements a Command observer
// used to monitor the evolution of the registration process.
//
#include "itkCommand.h"
class CommandProgressUpdate : public itk::Command
{
public:
using Self = CommandProgressUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandProgressUpdate() = 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
{
const auto * filter = static_cast<const itk::ProcessObject *>(object);
if (!itk::ProgressEvent().CheckEvent(&event))
{
return;
}
std::cout << filter->GetProgress() << std::endl;
}
};
int
main(int argc, char * argv[])
{
if (argc < 5)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " coefficientsFile fixedImage ";
std::cerr << "movingImage deformedMovingImage" << std::endl;
std::cerr << "[deformationField]" << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// Begin by creating the relevant types.
//
// Software Guide: EndLatex
// Software Guide : BeginCodeSnippet
constexpr unsigned int ImageDimension = 3;
using PixelType = unsigned char;
using FixedImageType = itk::Image<PixelType, ImageDimension>;
using MovingImageType = itk::Image<PixelType, ImageDimension>;
using FixedReaderType = itk::ImageFileReader<FixedImageType>;
using MovingReaderType = itk::ImageFileReader<MovingImageType>;
using MovingWriterType = itk::ImageFileWriter<MovingImageType>;
// Software Guide : EndCodeSnippet
FixedReaderType::Pointer fixedReader = FixedReaderType::New();
fixedReader->SetFileName(argv[2]);
try
{
fixedReader->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// Setup the moving reader and writer, and get the filenames from the
// command line arguments.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
MovingReaderType::Pointer movingReader = MovingReaderType::New();
MovingWriterType::Pointer movingWriter = MovingWriterType::New();
movingReader->SetFileName(argv[3]);
movingWriter->SetFileName(argv[4]);
// Software Guide : EndCodeSnippet
FixedImageType::ConstPointer fixedImage = fixedReader->GetOutput();
using FilterType =
FilterType::Pointer resampler = FilterType::New();
using InterpolatorType =
InterpolatorType::Pointer interpolator = InterpolatorType::New();
resampler->SetInterpolator(interpolator);
FixedImageType::SpacingType fixedSpacing = fixedImage->GetSpacing();
FixedImageType::PointType fixedOrigin = fixedImage->GetOrigin();
FixedImageType::DirectionType fixedDirection = fixedImage->GetDirection();
// Software Guide : BeginLatex
//
// Set the resampler spacing, origin, and direction to that of the fixed
// input image. Do the same with the size and output start index.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
resampler->SetOutputSpacing(fixedSpacing);
resampler->SetOutputOrigin(fixedOrigin);
resampler->SetOutputDirection(fixedDirection);
FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
FixedImageType::SizeType fixedSize = fixedRegion.GetSize();
resampler->SetSize(fixedSize);
resampler->SetOutputStartIndex(fixedRegion.GetIndex());
// Software Guide : EndCodeSnippet
resampler->SetInput(movingReader->GetOutput());
movingWriter->SetInput(resampler->GetOutput());
// Software Guide : BeginLatex
//
// We instantiate now the type of the \code{BSplineTransform} using
// as template parameters the type for coordinates representation, the
// dimension of the space, and the order of the B-spline.
//
// \index{BSplineTransform!New}
// \index{BSplineTransform!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int SpaceDimension = ImageDimension;
constexpr unsigned int SplineOrder = 3;
using CoordinateRepType = double;
using TransformType =
TransformType::Pointer bsplineTransform = TransformType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginCodeSnippet
constexpr unsigned int numberOfGridNodes = 8;
TransformType::PhysicalDimensionsType fixedPhysicalDimensions;
TransformType::MeshSizeType meshSize;
for (unsigned int i = 0; i < SpaceDimension; i++)
{
fixedPhysicalDimensions[i] =
fixedSpacing[i] * static_cast<double>(fixedSize[i] - 1);
}
meshSize.Fill(numberOfGridNodes - SplineOrder);
bsplineTransform->SetTransformDomainOrigin(fixedOrigin);
bsplineTransform->SetTransformDomainPhysicalDimensions(
fixedPhysicalDimensions);
bsplineTransform->SetTransformDomainMeshSize(meshSize);
bsplineTransform->SetTransformDomainDirection(fixedDirection);
using ParametersType = TransformType::ParametersType;
const unsigned int numberOfParameters =
bsplineTransform->GetNumberOfParameters();
const unsigned int numberOfNodes = numberOfParameters / SpaceDimension;
ParametersType parameters(numberOfParameters);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The B-spline grid should now be fed with coefficients at each node.
// Since this is a two-dimensional grid, each node should receive two
// coefficients. Each coefficient pair is representing a displacement
// vector at this node. The coefficients can be passed to the B-spline in
// the form of an array where the first set of elements are the first
// component of the displacements for all the nodes, and the second set of
// elemets is formed by the second component of the displacements for all
// the nodes.
//
// In this example we read such displacements from a file, but for
// convenience we have written this file using the pairs of $(x,y)$
// displacement for every node. The elements read from the file should
// therefore be reorganized when assigned to the elements of the array. We
// do this by storing all the odd elements from the file in the first block
// of the array, and all the even elements from the file in the second
// block of the array. Finally the array is passed to the B-spline
// transform using the \code{SetParameters()}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
std::ifstream infile;
infile.open(argv[1]);
for (unsigned int n = 0; n < numberOfNodes; n++)
{
infile >> parameters[n]; // X coordinate
infile >> parameters[n + numberOfNodes]; // Y coordinate
infile >> parameters[n + numberOfNodes * 2]; // Z coordinate
}
infile.close();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Finally the array is passed to the B-spline transform using the
// \code{SetParameters()}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
bsplineTransform->SetParameters(parameters);
// Software Guide : EndCodeSnippet
CommandProgressUpdate::Pointer observer = CommandProgressUpdate::New();
resampler->AddObserver(itk::ProgressEvent(), observer);
// Software Guide : BeginLatex
//
// At this point we are ready to use the transform as part of the resample
// filter. We trigger the execution of the pipeline by invoking
// \code{Update()} on the last filter of the pipeline, in this case the
// writer.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
resampler->SetTransform(bsplineTransform);
try
{
movingWriter->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Software Guide : EndCodeSnippet
using DisplacementFieldType = itk::Image<VectorType, ImageDimension>;
DisplacementFieldType::Pointer field = DisplacementFieldType::New();
field->SetRegions(fixedRegion);
field->SetOrigin(fixedOrigin);
field->SetSpacing(fixedSpacing);
field->SetDirection(fixedDirection);
field->Allocate();
FieldIterator fi(field, fixedRegion);
fi.GoToBegin();
TransformType::InputPointType fixedPoint;
TransformType::OutputPointType movingPoint;
VectorType displacement;
while (!fi.IsAtEnd())
{
index = fi.GetIndex();
field->TransformIndexToPhysicalPoint(index, fixedPoint);
movingPoint = bsplineTransform->TransformPoint(fixedPoint);
displacement = movingPoint - fixedPoint;
fi.Set(displacement);
++fi;
}
FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
fieldWriter->SetInput(field);
if (argc >= 6)
{
fieldWriter->SetFileName(argv[5]);
try
{
fieldWriter->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
}
if (argc >= 7)
{
fieldWriter->SetFileName(argv[6]);
try
{
using TransformWriterType = itk::TransformFileWriter;
TransformWriterType::Pointer transformWriter =
TransformWriterType::New();
transformWriter->AddTransform(bsplineTransform);
transformWriter->SetFileName(argv[6]);
transformWriter->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
}
return EXIT_SUCCESS;
}
itk::TransformFileWriter
itk::TransformFileWriterTemplate< double > TransformFileWriter
Definition: itkTransformFileWriter.h:129
itk::GTest::TypedefsAndConstructors::Dimension2::DirectionType
ImageBaseType::DirectionType DirectionType
Definition: itkGTestTypedefsAndConstructors.h:52
itk::GTest::TypedefsAndConstructors::Dimension2::VectorType
ImageBaseType::SpacingType VectorType
Definition: itkGTestTypedefsAndConstructors.h:53
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
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itk::SmartPointer< Self >
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::ImageRegionIterator
A multi-dimensional iterator templated over image type that walks a region of pixels.
Definition: itkImageRegionIterator.h:78
itk::GTest::TypedefsAndConstructors::Dimension2::IndexType
ImageBaseType::IndexType IndexType
Definition: itkGTestTypedefsAndConstructors.h:50
itk::LinearInterpolateImageFunction
Linearly interpolate an image at specified positions.
Definition: itkLinearInterpolateImageFunction.h:50
itk::Command
Superclass for callback/observer methods.
Definition: itkCommand.h:45
itk::BSplineTransform
Deformable transform using a BSpline representation.
Definition: itkBSplineTransform.h:103
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itkBSplineTransform.h
itk::Command
class ITK_FORWARD_EXPORT Command
Definition: itkObject.h:43
itk::GTest::TypedefsAndConstructors::Dimension2::RegionType
ImageBaseType::RegionType RegionType
Definition: itkGTestTypedefsAndConstructors.h:54
itk::Command::Execute
virtual void Execute(Object *caller, const EventObject &event)=0
itkImageFileWriter.h
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Object
Base class for most ITK classes.
Definition: itkObject.h:62
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::EventObject
Abstraction of the Events used to communicating among filters and with GUIs.
Definition: itkEventObject.h:57
itkResampleImageFilter.h
itkTransformFileWriter.h
itkCommand.h
itk::Size::GetSize
const SizeValueType * GetSize() const
Definition: itkSize.h:169