template <typename TRegistration>
{
public:
using Self = RegistrationInterfaceCommand;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() = default;
public:
using RegistrationType = TRegistration;
using RegistrationPointer = RegistrationType *;
using OptimizerPointer = OptimizerType *;
void
{
if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))
{
return;
}
auto registration = static_cast<RegistrationPointer>(object);
auto optimizer =
static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());
unsigned int currentLevel = registration->GetCurrentLevel();
typename RegistrationType::ShrinkFactorsPerDimensionContainerType
shrinkFactors =
registration->GetShrinkFactorsPerDimension(currentLevel);
typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
registration->GetSmoothingSigmasPerLevel();
std::cout << "-------------------------------------" << std::endl;
std::cout << " Current level = " << currentLevel << std::endl;
std::cout << " shrink factor = " << shrinkFactors << std::endl;
std::cout << " smoothing sigma = ";
std::cout << smoothingSigmas[currentLevel] << std::endl;
std::cout << std::endl;
if (registration->GetCurrentLevel() == 0)
{
optimizer->SetLearningRate(16.00);
optimizer->SetMinimumStepLength(2.5);
}
else
{
optimizer->SetLearningRate(optimizer->GetCurrentStepLength());
optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() *
0.2);
}
}
void
{
return;
}
};
{
public:
using Self = CommandIterationUpdate;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerPointer = const OptimizerType *;
void
{
}
void
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << " ";
std::cout << m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex{ 0 };
};
int
main(int argc, const char * 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;
const std::string fixedImageFile = argv[1];
const std::string movingImageFile = argv[2];
const std::string outImagefile = argv[3];
const PixelType backgroundGrayLevel = (argc > 4) ? std::stoi(argv[4]) : 100;
const std::string checkerBoardBefore = (argc > 5) ? argv[5] : "";
const std::string checkerBoardAfter = (argc > 6) ? argv[6] : "";
const int numberOfBins = (argc > 7) ? std::stoi(argv[7]) : 0;
using MetricType =
MovingImageType>;
using RegistrationType = itk::
ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
MetricType::Pointer metric = MetricType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetOptimizer(optimizer);
registration->SetMetric(metric);
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName(fixedImageFile);
movingImageReader->SetFileName(movingImageFile);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
using ParametersType = OptimizerType::ParametersType;
ParametersType initialParameters(transform->GetNumberOfParameters());
initialParameters[0] = 0.0;
initialParameters[1] = 0.0;
transform->SetParameters(initialParameters);
registration->SetInitialTransform(transform);
registration->InPlaceOn();
metric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
metric->SetNumberOfHistogramBins(numberOfBins);
}
optimizer->SetNumberOfIterations(200);
optimizer->SetRelaxationFactor(0.5);
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
constexpr unsigned int numberOfLevels = 3;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(3);
shrinkFactorsPerLevel[0] = 3;
shrinkFactorsPerLevel[1] = 2;
shrinkFactorsPerLevel[2] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(3);
smoothingSigmasPerLevel[0] = 0;
smoothingSigmasPerLevel[1] = 0;
smoothingSigmasPerLevel[2] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
using CommandType = RegistrationInterfaceCommand<RegistrationType>;
CommandType::Pointer command = CommandType::New();
registration->AddObserver(itk::MultiResolutionIterationEvent(), command);
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
ParametersType finalParameters = transform->GetParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
using ResampleFilterType =
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(transform);
resample->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
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 =
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName(outImagefile);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
resample->SetDefaultPixelValue(0);
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
for (int q = 0; q < argc; ++q)
{
std::cout << q << " " << argv[q] << std::endl;
}
if (checkerBoardBefore != std::string(""))
{
writer->SetFileName(checkerBoardBefore);
writer->Update();
}
resample->SetTransform(transform);
if (checkerBoardAfter != std::string(""))
{
writer->SetFileName(checkerBoardAfter);
writer->Update();
}
return EXIT_SUCCESS;
}