ITK  6.0.0
Insight Toolkit
SphinxExamples/src/Registration/Metricsv4/RegisterTwoPointSets/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.
*
*=========================================================================*/
// Adapted from ITK itkJensenHavrdaCharvatTsallisPointSetMetricRegistrationTest.cxx
#include "itkTransform.h"
#include "itkCommand.h"
#include <fstream>
template <typename TFilter>
class itkJensenHavrdaCharvatTsallisPointSetMetricRegistrationTestCommandIterationUpdate : public itk::Command
{
public:
using Self = itkJensenHavrdaCharvatTsallisPointSetMetricRegistrationTestCommandIterationUpdate;
itkNewMacro(Self);
protected:
itkJensenHavrdaCharvatTsallisPointSetMetricRegistrationTestCommandIterationUpdate() = 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
{
if (typeid(event) != typeid(itk::IterationEvent))
{
return;
}
const auto * optimizer = dynamic_cast<const TFilter *>(object);
if (!optimizer)
{
itkGenericExceptionMacro("Error dynamic_cast failed");
}
std::cout << "It: " << optimizer->GetCurrentIteration() << " metric value: " << optimizer->GetCurrentMetricValue();
std::cout << std::endl;
}
};
int
main(int argc, char * argv[])
{
constexpr unsigned int Dimension = 2;
unsigned int numberOfIterations = 10;
if (argc > 1)
{
numberOfIterations = std::stoi(argv[1]);
}
auto fixedPoints = PointSetType::New();
fixedPoints->Initialize();
auto movingPoints = PointSetType::New();
movingPoints->Initialize();
// two ellipses, one rotated slightly
/*
// Having trouble with these, as soon as there's a slight rotation added.
unsigned long count = 0;
for( float theta = 0; theta < 2.0 * itk::Math::pi; theta += 0.1 )
{
float radius = 100.0;
PointType fixedPoint;
fixedPoint[0] = 2 * radius * std::cos( theta );
fixedPoint[1] = radius * std::sin( theta );
fixedPoints->SetPoint( count, fixedPoint );
PointType movingPoint;
movingPoint[0] = 2 * radius * std::cos( theta + (0.02 * itk::Math::pi) ) + 2.0;
movingPoint[1] = radius * std::sin( theta + (0.02 * itk::Math::pi) ) + 2.0;
movingPoints->SetPoint( count, movingPoint );
count++;
}
*/
// two circles with a small offset
PointType offset;
for (unsigned int d = 0; d < Dimension; ++d)
{
offset[d] = 2.0;
}
unsigned long count = 0;
for (float theta = 0; theta < 2.0 * itk::Math::pi; theta += 0.1)
{
PointType fixedPoint;
float radius = 100.0;
fixedPoint[0] = radius * std::cos(theta);
fixedPoint[1] = radius * std::sin(theta);
if (Dimension > 2)
{
fixedPoint[2] = radius * std::sin(theta);
}
fixedPoints->SetPoint(count, fixedPoint);
PointType movingPoint;
movingPoint[0] = fixedPoint[0] + offset[0];
movingPoint[1] = fixedPoint[1] + offset[1];
if (Dimension > 2)
{
movingPoint[2] = fixedPoint[2] + offset[2];
}
movingPoints->SetPoint(count, movingPoint);
count++;
}
using AffineTransformType = itk::AffineTransform<double, Dimension>;
auto transform = AffineTransformType::New();
transform->SetIdentity();
// Instantiate the metric
auto metric = PointSetMetricType::New();
metric->SetFixedPointSet(fixedPoints);
metric->SetMovingPointSet(movingPoints);
metric->SetPointSetSigma(1.0);
metric->SetKernelSigma(10.0);
metric->SetUseAnisotropicCovariances(false);
metric->SetCovarianceKNeighborhood(5);
metric->SetEvaluationKNeighborhood(10);
metric->SetMovingTransform(transform);
metric->SetAlpha(1.1);
metric->Initialize();
// scales estimator
using RegistrationParameterScalesFromShiftType =
shiftScaleEstimator->SetMetric(metric);
// needed with pointset metrics
shiftScaleEstimator->SetVirtualDomainPointSet(metric->GetVirtualTransformedPointSet());
// optimizer
using OptimizerType = itk::GradientDescentOptimizerv4;
auto optimizer = OptimizerType::New();
optimizer->SetMetric(metric);
optimizer->SetNumberOfIterations(numberOfIterations);
optimizer->SetScalesEstimator(shiftScaleEstimator);
optimizer->SetMaximumStepSizeInPhysicalUnits(3.0);
using CommandType = itkJensenHavrdaCharvatTsallisPointSetMetricRegistrationTestCommandIterationUpdate<OptimizerType>;
auto observer = CommandType::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
optimizer->SetMinimumConvergenceValue(0.0);
optimizer->SetConvergenceWindowSize(10);
optimizer->StartOptimization();
std::cout << "numberOfIterations: " << numberOfIterations << std::endl;
std::cout << "Moving-source final value: " << optimizer->GetCurrentMetricValue() << std::endl;
std::cout << "Moving-source final position: " << optimizer->GetCurrentPosition() << std::endl;
std::cout << "Optimizer scales: " << optimizer->GetScales() << std::endl;
std::cout << "Optimizer learning rate: " << optimizer->GetLearningRate() << std::endl;
// applying the resultant transform to moving points and verify result
std::cout << "Fixed\tMoving\tMovingTransformed\tFixedTransformed\tDiff" << std::endl;
bool passed = true;
PointType::ValueType tolerance = 1e-2;
AffineTransformType::InverseTransformBasePointer movingInverse = metric->GetMovingTransform()->GetInverseTransform();
AffineTransformType::InverseTransformBasePointer fixedInverse = metric->GetFixedTransform()->GetInverseTransform();
for (unsigned int n = 0; n < metric->GetNumberOfComponents(); ++n)
{
// compare the points in virtual domain
PointType transformedMovingPoint = movingInverse->TransformPoint(movingPoints->GetPoint(n));
PointType transformedFixedPoint = fixedInverse->TransformPoint(fixedPoints->GetPoint(n));
PointType difference;
difference[0] = transformedMovingPoint[0] - transformedFixedPoint[0];
difference[1] = transformedMovingPoint[1] - transformedFixedPoint[1];
std::cout << fixedPoints->GetPoint(n) << "\t" << movingPoints->GetPoint(n) << "\t" << transformedMovingPoint << "\t"
<< transformedFixedPoint << "\t" << difference << std::endl;
if (fabs(difference[0]) > tolerance || fabs(difference[1]) > tolerance)
{
passed = false;
}
}
if (!passed)
{
std::cerr << "Results do not match truth within tolerance." << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
Pointer
SmartPointer< Self > Pointer
Definition: itkAddImageFilter.h:93
itk::PointSet
A superclass of the N-dimensional mesh structure; supports point (geometric coordinate and attribute)...
Definition: itkPointSet.h:81
itk::GTest::TypedefsAndConstructors::Dimension2::PointType
ImageBaseType::PointType PointType
Definition: itkGTestTypedefsAndConstructors.h:51
itkJensenHavrdaCharvatTsallisPointSetToPointSetMetricv4.h
itk::SmartPointer< Self >
itkAffineTransform.h
itk::AffineTransform
Definition: itkAffineTransform.h:101
itk::Command
Superclass for callback/observer methods.
Definition: itkCommand.h:45
itkRegistrationParameterScalesFromPhysicalShift.h
itk::Command
class ITK_FORWARD_EXPORT Command
Definition: itkObject.h:42
itkGradientDescentOptimizerv4.h
itk::JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4
Implementation of the Jensen Havrda Charvat Tsallis Point Set metric.
Definition: itkJensenHavrdaCharvatTsallisPointSetToPointSetMetricv4.h:72
itk::Command::Execute
virtual void Execute(Object *caller, const EventObject &event)=0
itk::GradientDescentOptimizerv4
GradientDescentOptimizerv4Template< double > GradientDescentOptimizerv4
Definition: itkGradientDescentOptimizerv4.h:255
itk::Object
Base class for most ITK classes.
Definition: itkObject.h:61
itk::Math::e
static constexpr double e
Definition: itkMath.h:56
itk::RegistrationParameterScalesFromPhysicalShift
Registration helper class for estimating scales of transform parameters a step sizes,...
Definition: itkRegistrationParameterScalesFromPhysicalShift.h:35
itk::EventObject
Abstraction of the Events used to communicating among filters and with GUIs.
Definition: itkEventObject.h:58
New
static Pointer New()
AddImageFilter
Definition: itkAddImageFilter.h:81
itkTransform.h
itk::Math::pi
static constexpr double pi
Definition: itkMath.h:66
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkCommand.h
Superclass
BinaryGeneratorImageFilter< TInputImage1, TInputImage2, TOutputImage > Superclass
Definition: itkAddImageFilter.h:90