ITK  5.3.0
Insight Toolkit
Examples/RegistrationITKv4/IterativeClosestPoint2.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 perform Iterative Closest Point (ICP)
// registration in ITK using sets of 3D points.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// The first step is to include the relevant headers.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include <iostream>
#include <fstream>
// Software Guide : EndCodeSnippet
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerType = itk::LevenbergMarquardtOptimizer;
using OptimizerPointer = const OptimizerType *;
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
{
auto optimizer = dynamic_cast<OptimizerPointer>(object);
if (optimizer == nullptr)
{
itkExceptionMacro("Could not cast optimizer.");
}
if (!itk::IterationEvent().CheckEvent(&event))
{
return;
}
std::cout << "Value = " << optimizer->GetCachedValue() << std::endl;
std::cout << "Position = " << optimizer->GetCachedCurrentPosition();
std::cout << std::endl << std::endl;
}
};
int
main(int argc, char * argv[])
{
if (argc < 3)
{
std::cerr << "Arguments Missing. " << std::endl;
std::cerr
<< "Usage: IterativeClosestPoint2 fixedPointsFile movingPointsFile "
<< std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 3;
// Software Guide : BeginLatex
//
// First, define the necessary types for the moving and fixed point sets.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using PointSetType = itk::PointSet<float, Dimension>;
PointSetType::Pointer fixedPointSet = PointSetType::New();
PointSetType::Pointer movingPointSet = PointSetType::New();
using PointsContainer = PointSetType::PointsContainer;
PointsContainer::Pointer fixedPointContainer = PointsContainer::New();
PointsContainer::Pointer movingPointContainer = PointsContainer::New();
PointType fixedPoint;
PointType movingPoint;
// Software Guide : EndCodeSnippet
// Read the file containing coordinates of fixed points.
std::ifstream fixedFile;
fixedFile.open(argv[1]);
if (fixedFile.fail())
{
std::cerr << "Error opening points file with name : " << std::endl;
std::cerr << argv[1] << std::endl;
return EXIT_FAILURE;
}
unsigned int pointId = 0;
fixedFile >> fixedPoint;
while (!fixedFile.eof())
{
fixedPointContainer->InsertElement(pointId, fixedPoint);
fixedFile >> fixedPoint;
pointId++;
}
fixedPointSet->SetPoints(fixedPointContainer);
std::cout << "Number of fixed Points = "
<< fixedPointSet->GetNumberOfPoints() << std::endl;
// Read the file containing coordinates of moving points.
std::ifstream movingFile;
movingFile.open(argv[2]);
if (movingFile.fail())
{
std::cerr << "Error opening points file with name : " << std::endl;
std::cerr << argv[2] << std::endl;
return EXIT_FAILURE;
}
pointId = 0;
movingFile >> movingPoint;
while (!movingFile.eof())
{
movingPointContainer->InsertElement(pointId, movingPoint);
movingFile >> movingPoint;
pointId++;
}
movingPointSet->SetPoints(movingPointContainer);
std::cout << "Number of moving Points = "
<< movingPointSet->GetNumberOfPoints() << std::endl;
// Software Guide : BeginLatex
//
// After the points are read in from files, setup the metric to be used
// later by the registration.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using MetricType =
MetricType::Pointer metric = MetricType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Next, setup the tranform, optimizers, and registration.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using TransformType = itk::Euler3DTransform<double>;
TransformType::Pointer transform = TransformType::New();
// Optimizer Type
using OptimizerType = itk::LevenbergMarquardtOptimizer;
OptimizerType::Pointer optimizer = OptimizerType::New();
optimizer->SetUseCostFunctionGradient(false);
// Registration Method
using RegistrationType =
RegistrationType::Pointer registration = RegistrationType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Scale the translation components of the Transform in the Optimizer
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
OptimizerType::ScalesType scales(transform->GetNumberOfParameters());
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Next, set the scales and ranges for translations and rotations in the
// transform. Also, set the convergence criteria and number of iterations
// to be used by the optimizer.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
constexpr double translationScale = 1000.0; // dynamic range of translations
constexpr double rotationScale = 1.0; // dynamic range of rotations
scales[0] = 1.0 / rotationScale;
scales[1] = 1.0 / rotationScale;
scales[2] = 1.0 / rotationScale;
scales[3] = 1.0 / translationScale;
scales[4] = 1.0 / translationScale;
scales[5] = 1.0 / translationScale;
unsigned long numberOfIterations = 2000;
double gradientTolerance = 1e-4; // convergence criterion
double valueTolerance = 1e-4; // convergence criterion
double epsilonFunction = 1e-5; // convergence criterion
optimizer->SetScales(scales);
optimizer->SetNumberOfIterations(numberOfIterations);
optimizer->SetValueTolerance(valueTolerance);
optimizer->SetGradientTolerance(gradientTolerance);
optimizer->SetEpsilonFunction(epsilonFunction);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Here we start with an identity transform, although the user will usually
// be able to provide a better guess than this.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
transform->SetIdentity();
// Software Guide : EndCodeSnippet
registration->SetInitialTransformParameters(transform->GetParameters());
// Software Guide : BeginLatex
//
// Connect all the components required for the registration.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
registration->SetTransform(transform);
registration->SetFixedPointSet(fixedPointSet);
registration->SetMovingPointSet(movingPointSet);
// Software Guide : EndCodeSnippet
//
// Connect an observer
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
try
{
registration->Update();
}
catch (const itk::ExceptionObject & e)
{
std::cerr << e << std::endl;
return EXIT_FAILURE;
}
std::cout << "Solution = " << transform->GetParameters() << std::endl;
std::cout << "Stopping condition: "
<< optimizer->GetStopConditionDescription() << std::endl;
return EXIT_SUCCESS;
}
itk::Euler3DTransform
Euler3DTransform of a vector space (e.g. space coordinates)
Definition: itkEuler3DTransform.h:47
itk::PointSet
A superclass of the N-dimensional mesh structure; supports point (geometric coordinate and attribute)...
Definition: itkPointSet.h:82
itkLevenbergMarquardtOptimizer.h
itk::EuclideanDistancePointMetric
Computes the minimum distance between a moving point-set and a fixed point-set. A vector of minimum c...
Definition: itkEuclideanDistancePointMetric.h:47
itk::GTest::TypedefsAndConstructors::Dimension2::PointType
ImageBaseType::PointType PointType
Definition: itkGTestTypedefsAndConstructors.h:51
itkPointSetToPointSetRegistrationMethod.h
itk::SmartPointer< Self >
itk::Command
Superclass for callback/observer methods.
Definition: itkCommand.h:45
itkEuclideanDistancePointMetric.h
itk::Command
class ITK_FORWARD_EXPORT Command
Definition: itkObject.h:43
itkEuler3DTransform.h
itk::Command::Execute
virtual void Execute(Object *caller, const EventObject &event)=0
itk::Object
Base class for most ITK classes.
Definition: itkObject.h:62
itk::Math::e
static constexpr double e
Definition: itkMath.h:54
itk::EventObject
Abstraction of the Events used to communicating among filters and with GUIs.
Definition: itkEventObject.h:57
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itk::PointSetToPointSetRegistrationMethod
Base class for PointSet to PointSet Registration Methods.
Definition: itkPointSetToPointSetRegistrationMethod.h:66
itk::LevenbergMarquardtOptimizer
Wrap of the vnl_levenberg_marquardt algorithm.
Definition: itkLevenbergMarquardtOptimizer.h:37