ITK  5.2.0
Insight Toolkit
Examples/RegistrationITKv4/ThinPlateSplineWarp.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.
*
*=========================================================================*/
// Command Line Arguments:
// Insight/Testing/Data/LandmarkWarping3Landmarks1.txt
// inputImage deformedImage deformationField
//
// Software Guide : BeginLatex
//
// This example deforms a 3D volume with the Thin plate spline.
// \index{ThinPlateSplineKernelTransform}
//
// Software Guide : EndLatex
#include "itkImage.h"
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
#include "itkPointSet.h"
#include <fstream>
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " landmarksFile inputImage ";
std::cerr << "DeformedImage " << std::endl;
std::cerr << "deformationField" << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int ImageDimension = 3;
using PixelType = unsigned char;
using InputImageType = itk::Image<PixelType, ImageDimension>;
using DeformedImageWriterType = itk::ImageFileWriter<InputImageType>;
using FieldVectorType = itk::Vector<float, ImageDimension>;
using DisplacementFieldType = itk::Image<FieldVectorType, ImageDimension>;
using CoordinateRepType = double;
using TransformType =
using PointSetType = TransformType::PointSetType;
using PointIdType = PointSetType::PointIdentifier;
using ResamplerType =
using InterpolatorType =
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(argv[2]);
try
{
reader->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// Landmarks correspondances may be associated with the
// SplineKernelTransforms via Point Set containers. Let us define containers
// for the landmarks.
//
// Software Guide : EndLatex
// Define container for landmarks
// Software Guide : BeginCodeSnippet
PointSetType::Pointer sourceLandMarks = PointSetType::New();
PointSetType::Pointer targetLandMarks = PointSetType::New();
PointSetType::PointsContainer::Pointer sourceLandMarkContainer =
sourceLandMarks->GetPoints();
PointSetType::PointsContainer::Pointer targetLandMarkContainer =
targetLandMarks->GetPoints();
// Software Guide : EndCodeSnippet
// Read in the list of landmarks
std::ifstream infile;
infile.open(argv[1]);
while (!infile.eof())
{
infile >> p1[0] >> p1[1] >> p1[2] >> p2[0] >> p2[1] >> p2[2];
// Software Guide : BeginCodeSnippet
sourceLandMarkContainer->InsertElement(id, p1);
targetLandMarkContainer->InsertElement(id++, p2);
// Software Guide : EndCodeSnippet
}
infile.close();
// Software Guide : BeginCodeSnippet
TransformType::Pointer tps = TransformType::New();
tps->SetSourceLandmarks(sourceLandMarks);
tps->SetTargetLandmarks(targetLandMarks);
tps->ComputeWMatrix();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The image is then resampled to produce an output image as defined by the
// transform. Here we use a LinearInterpolator.
//
// Software Guide : EndLatex
// Set the resampler params
InputImageType::ConstPointer inputImage = reader->GetOutput();
ResamplerType::Pointer resampler = ResamplerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
resampler->SetInterpolator(interpolator);
InputImageType::SpacingType spacing = inputImage->GetSpacing();
InputImageType::PointType origin = inputImage->GetOrigin();
InputImageType::DirectionType direction = inputImage->GetDirection();
InputImageType::RegionType region = inputImage->GetBufferedRegion();
// Software Guide : BeginCodeSnippet
resampler->SetOutputSpacing(spacing);
resampler->SetOutputDirection(direction);
resampler->SetOutputOrigin(origin);
resampler->SetSize(size);
resampler->SetTransform(tps);
// Software Guide : EndCodeSnippet
resampler->SetOutputStartIndex(region.GetIndex());
resampler->SetInput(reader->GetOutput());
// Set and write deformed image
DeformedImageWriterType::Pointer deformedImageWriter =
DeformedImageWriterType::New();
deformedImageWriter->SetInput(resampler->GetOutput());
deformedImageWriter->SetFileName(argv[3]);
try
{
deformedImageWriter->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// The deformation field is computed as the difference between the input and
// the deformed image by using an iterator.
//
// Software Guide : EndLatex
// Compute the deformation field
DisplacementFieldType::Pointer field = DisplacementFieldType::New();
field->SetRegions(region);
field->SetOrigin(origin);
field->SetSpacing(spacing);
field->Allocate();
FieldIterator fi(field, region);
fi.GoToBegin();
TransformType::InputPointType point1;
TransformType::OutputPointType point2;
FieldVectorType displacement;
while (!fi.IsAtEnd())
{
index = fi.GetIndex();
field->TransformIndexToPhysicalPoint(index, point1);
point2 = tps->TransformPoint(point1);
for (unsigned int i = 0; i < ImageDimension; i++)
{
displacement[i] = point2[i] - point1[i];
}
fi.Set(displacement);
++fi;
}
// Write computed deformation field
FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
fieldWriter->SetFileName(argv[4]);
fieldWriter->SetInput(field);
try
{
fieldWriter->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
itk::GTest::TypedefsAndConstructors::Dimension2::DirectionType
ImageBaseType::DirectionType DirectionType
Definition: itkGTestTypedefsAndConstructors.h:52
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
itkImage.h
itkThinPlateSplineKernelTransform.h
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::ThinPlateSplineKernelTransform
Definition: itkThinPlateSplineKernelTransform.h:35
itk::ImageRegionIterator
A multi-dimensional iterator templated over image type that walks a region of pixels.
Definition: itkImageRegionIterator.h:80
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::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:88
itk::Index::GetIndex
const IndexValueType * GetIndex() const
Definition: itkIndex.h:228
itk::GTest::TypedefsAndConstructors::Dimension2::RegionType
ImageBaseType::RegionType RegionType
Definition: itkGTestTypedefsAndConstructors.h:54
itkImageFileWriter.h
itk::NumericTraits::ZeroValue
static T ZeroValue()
Definition: itkNumericTraits.h:148
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Point
A templated class holding a geometric point in n-Dimensional space.
Definition: itkPoint.h:53
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itkPointSet.h
itkResampleImageFilter.h
itk::Size::GetSize
const SizeValueType * GetSize() const
Definition: itkSize.h:169