template <typename TRegistration>
{
public:
using Self = RegistrationInterfaceCommand;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() = default;
public:
using RegistrationType = TRegistration;
void
{
}
void
{
if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << "\nObserving from class " << object->GetNameOfClass();
{
std::cout << " \"" << object->GetObjectName() << "\"" << std::endl;
}
const auto * registration = static_cast<const RegistrationType *>(object);
const unsigned int currentLevel = registration->GetCurrentLevel();
const typename RegistrationType::ShrinkFactorsPerDimensionContainerType
shrinkFactors =
registration->GetShrinkFactorsPerDimension(currentLevel);
typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
registration->GetSmoothingSigmasPerLevel();
std::cout << "-------------------------------------" << std::endl;
std::cout << " Current multi-resolution level = " << currentLevel
<< std::endl;
std::cout << " shrink factor = " << shrinkFactors << std::endl;
std::cout << " smoothing sigma = " << smoothingSigmas[currentLevel]
<< std::endl;
std::cout << std::endl;
}
};
{
public:
using Self = CommandIterationUpdate;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerPointer = const OptimizerType *;
void
{
}
void
Execute(
const itk::Object *
object,
const itk::EventObject & event)
override
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << " "
<< m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex{ 0 };
};
namespace
{
int
ExampleMain(int argc, const char * const argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile [backgroundGrayLevel]";
std::cerr << " [checkerboardbefore] [CheckerBoardAfter]";
std::cerr << " [numberOfBins] " << std::endl;
return EXIT_FAILURE;
}
using PixelType = float;
using MetricType =
MovingImageType>;
MovingImageType,
TTransformType>;
auto transOptimizer = TOptimizerType::New();
auto transMetric = MetricType::New();
auto transRegistration = TRegistrationType::New();
transRegistration->SetOptimizer(transOptimizer);
transRegistration->SetMetric(transMetric);
auto movingInitTx = TTransformType::New();
using ParametersType = TOptimizerType::ParametersType;
ParametersType initialParameters(movingInitTx->GetNumberOfParameters());
initialParameters[0] = 3.0;
initialParameters[1] = 5.0;
movingInitTx->SetParameters(initialParameters);
transRegistration->SetMovingInitialTransform(movingInitTx);
auto compositeTransform = CompositeTransformType::New();
compositeTransform->AddTransform(movingInitTx);
transRegistration->SetFixedImage(fixedImage);
transRegistration->SetMovingImage(movingImage);
transRegistration->SetObjectName("TranslationRegistration");
constexpr unsigned int numberOfLevels1 = 1;
TRegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel1;
shrinkFactorsPerLevel1.SetSize(numberOfLevels1);
shrinkFactorsPerLevel1[0] = 3;
TRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel1;
smoothingSigmasPerLevel1.SetSize(numberOfLevels1);
smoothingSigmasPerLevel1[0] = 2;
transRegistration->SetNumberOfLevels(numberOfLevels1);
transRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel1);
transRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel1);
transMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
transMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
transOptimizer->SetNumberOfIterations(200);
transOptimizer->SetRelaxationFactor(0.5);
transOptimizer->SetLearningRate(16);
transOptimizer->SetMinimumStepLength(1.5);
auto observer1 = CommandIterationUpdate::New();
transOptimizer->AddObserver(itk::IterationEvent(), observer1);
using TranslationCommandType =
RegistrationInterfaceCommand<TRegistrationType>;
auto command1 = TranslationCommandType::New();
transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command1);
transRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< transRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
compositeTransform->AddTransform(
transRegistration->GetModifiableTransform());
using AOptimizerType =
MovingImageType,
ATransformType>;
auto affineOptimizer = AOptimizerType::New();
auto affineMetric = MetricType::New();
auto affineRegistration = ARegistrationType::New();
affineRegistration->SetOptimizer(affineOptimizer);
affineRegistration->SetMetric(affineMetric);
affineRegistration->SetMovingInitialTransform(compositeTransform);
affineRegistration->SetFixedImage(fixedImage);
affineRegistration->SetMovingImage(movingImage);
affineRegistration->SetObjectName("AffineRegistration");
affineMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
affineMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
auto affineTx = ATransformType::New();
using SpacingType = FixedImageType::SpacingType;
using OriginType = FixedImageType::PointType;
using SizeType = FixedImageType::SizeType;
const SpacingType fixedSpacing = fixedImage->GetSpacing();
const OriginType fixedOrigin = fixedImage->GetOrigin();
const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
ATransformType::InputPointType centerFixed;
centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
const unsigned int numberOfFixedParameters =
affineTx->GetFixedParameters().Size();
ATransformType::ParametersType fixedParameters(numberOfFixedParameters);
for (unsigned int i = 0; i < numberOfFixedParameters; ++i)
{
fixedParameters[i] = centerFixed[i];
}
affineTx->SetFixedParameters(fixedParameters);
affineRegistration->SetInitialTransform(affineTx);
using ScalesEstimatorType =
auto scalesEstimator = ScalesEstimatorType::New();
scalesEstimator->SetMetric(affineMetric);
scalesEstimator->SetTransformForward(true);
affineOptimizer->SetScalesEstimator(scalesEstimator);
affineOptimizer->SetDoEstimateLearningRateOnce(true);
affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
affineOptimizer->SetLowerLimit(0);
affineOptimizer->SetUpperLimit(2);
affineOptimizer->SetEpsilon(0.2);
affineOptimizer->SetNumberOfIterations(200);
affineOptimizer->SetMinimumConvergenceValue(1e-6);
affineOptimizer->SetConvergenceWindowSize(5);
auto observer2 = CommandIterationUpdate::New();
affineOptimizer->AddObserver(itk::IterationEvent(), observer2);
constexpr unsigned int numberOfLevels2 = 2;
ARegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel2;
shrinkFactorsPerLevel2.SetSize(numberOfLevels2);
shrinkFactorsPerLevel2[0] = 2;
shrinkFactorsPerLevel2[1] = 1;
ARegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel2;
smoothingSigmasPerLevel2.SetSize(numberOfLevels2);
smoothingSigmasPerLevel2[0] = 1;
smoothingSigmasPerLevel2[1] = 0;
affineRegistration->SetNumberOfLevels(numberOfLevels2);
affineRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel2);
affineRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel2);
using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
auto command2 = AffineCommandType::New();
affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command2);
affineRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< affineRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
compositeTransform->AddTransform(
affineRegistration->GetModifiableTransform());
std::cout << "\nInitial parameters of the registration process: "
<< std::endl
<< movingInitTx->GetParameters() << std::endl;
std::cout << "\nTranslation parameters after registration: " << std::endl
<< transOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: "
<< transOptimizer->GetCurrentStepLength() << std::endl;
std::cout << "\nAffine parameters after registration: " << std::endl
<< affineOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: " << affineOptimizer->GetLearningRate()
<< std::endl;
using ResampleFilterType =
auto resample = ResampleFilterType::New();
resample->SetTransform(compositeTransform);
resample->SetInput(movingImage);
PixelType backgroundGrayLevel = 100;
if (argc > 4)
{
backgroundGrayLevel = std::stoi(argv[4]);
}
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(backgroundGrayLevel);
using OutputPixelType = unsigned char;
using CastFilterType =
auto caster = CastFilterType::New();
caster->SetInput(resample->GetOutput());
auto checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
resample->SetDefaultPixelValue(0);
TransformType::Pointer identityTransform;
identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 5)
{
}
resample->SetTransform(compositeTransform);
if (argc > 6)
{
}
return EXIT_SUCCESS;
}
}
int
main(int argc, char * argv[])
{
try
{
return ExampleMain(argc, argv);
}
{
std::cerr << "ITK exception caught:\n" << exceptionObject << '\n';
}
catch (const std::exception & stdException)
{
std::cerr << "std exception caught:\n" << stdException.what() << '\n';
}
catch (...)
{
std::cerr << "Unhandled exception!\n";
}
return EXIT_FAILURE;
}
Casts input pixels to output pixel type.
Combines two images in a checkerboard pattern.
Superclass for callback/observer methods.
virtual void Execute(Object *caller, const EventObject &event)=0
Conjugate gradient descent optimizer with a golden section line search for nonlinear optimization.
Abstraction of the Events used to communicating among filters and with GUIs.
Standard exception handling object.
Gradient descent optimizer.
Interface method for the current registration framework.
Templated n-dimensional image class.
Base class for most ITK classes.
virtual const std::string & GetObjectName() const
Registration helper class for estimating scales of transform parameters a step sizes,...
Regular Step Gradient descent optimizer.
Resample an image via a coordinate transform.
Implements transparent reference counting.
constexpr unsigned int Dimension
ImageBaseType::RegionType RegionType
ImageBaseType::SizeType SizeType
ITK_TEMPLATE_EXPORT void WriteImage(TImagePointer &&image, const std::string &filename, bool compress=false)
TOutputImage::Pointer ReadImage(const std::string &filename)
const SizeValueType * GetSize() const