#include <iomanip>
#include <cstdio>
template <class TInput>
class RescaleDynamicRangeFunctor
{
public:
using OutputPixelType = unsigned char;
RescaleDynamicRangeFunctor() = default;
~RescaleDynamicRangeFunctor() = default;
inline OutputPixelType
operator()(const TInput & A)
{
if ((A > 0.0))
{
if (-(30.0 * std::log(A)) > 255)
{
return static_cast<OutputPixelType>(255);
}
else
{
return itk::Math::Round<OutputPixelType>(-(30.0 * std::log(A)));
}
}
else
{
return static_cast<OutputPixelType>(255);
}
}
};
namespace
{
class HistogramWriter
{
public:
using InternalPixelType = float;
using MetricType =
InternalImageType>;
using HistogramType = MetricType::HistogramType;
using HistogramToEntropyImageFilterType =
using HistogramToImageFilterPointer =
using OutputImageType = HistogramToEntropyImageFilterType::OutputImageType;
using OutputPixelType = HistogramToEntropyImageFilterType::OutputPixelType;
HistogramWriter()
: m_Metric(nullptr)
{
this->m_HistogramFileWriter->SetInput(this->m_Filter->GetOutput());
}
~HistogramWriter() = default;
void
SetMetric(MetricPointer metric)
{
this->m_Metric = metric;
}
MetricPointer
GetMetric() const
{
return this->m_Metric;
}
void
WriteHistogramFile(unsigned int iterationNumber)
{
const std::string outputFileBase = "JointHistogram";
std::ostringstream outputFilename;
outputFilename << outputFileBase << "." << std::setfill('0')
<< std::setw(3) << iterationNumber << "."
<< "mhd";
m_HistogramFileWriter->SetFileName(outputFilename.str());
this->m_Filter->SetInput(this->GetMetric()->GetHistogram());
this->m_Filter->Modified();
try
{
m_Filter->Update();
}
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
try
{
m_HistogramFileWriter->Update();
}
{
std::cerr << "Exception thrown " << excp << std::endl;
}
std::cout << "Joint Histogram file: ";
std::cout << outputFilename.str() << " written" << std::endl;
}
void
WriteHistogramFile(std::string & outputFilename)
{
this->m_Filter->SetInput(this->GetMetric()->GetHistogram());
this->m_Filter->Modified();
using RescaleDynamicRangeFunctorType =
RescaleDynamicRangeFunctor<OutputPixelType>;
using RescaleDynamicRangeFilterType =
RescaledOutputImageType,
RescaleDynamicRangeFunctorType>;
rescaler->SetInput(m_Filter->GetOutput());
rescaledWriter->SetInput(rescaler->GetOutput());
rescaledWriter->SetFileName(outputFilename);
try
{
m_Filter->Update();
}
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
try
{
rescaledWriter->Update();
}
{
std::cerr << "Exception thrown " << excp << std::endl;
}
std::cout << "Joint Histogram file: " << outputFilename << " written"
<< std::endl;
}
private:
MetricPointer m_Metric;
HistogramToImageFilterPointer m_Filter;
HistogramFileWriterPointer m_HistogramFileWriter;
std::string m_OutputFile;
};
{
public:
using Self = CommandIterationUpdate;
protected:
CommandIterationUpdate() { m_WriteHistogramsAfterEveryIteration = false; }
public:
using OptimizerPointer = const OptimizerType *;
void
{
}
void
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!itk::IterationEvent().CheckEvent(&event) || optimizer == nullptr)
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
if (optimizer->GetCurrentIteration() == 0)
{
m_JointHistogramWriter.WriteHistogramFile(m_InitialHistogramFile);
}
if (m_WriteHistogramsAfterEveryIteration)
{
m_JointHistogramWriter.WriteHistogramFile(
optimizer->GetCurrentIteration());
}
}
void
SetWriteHistogramsAfterEveryIteration(bool value)
{
m_WriteHistogramsAfterEveryIteration = value;
}
void
SetInitialHistogramFile(const char * filename)
{
m_InitialHistogramFile = filename;
}
HistogramWriter m_JointHistogramWriter;
private:
bool m_WriteHistogramsAfterEveryIteration;
std::string m_InitialHistogramFile;
};
}
int
main(int argc, char * argv[])
{
if (argc < 8)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << "outputImagefile WriteJointHistogramsAfterEveryIteration ";
std::cerr << "JointHistogramPriorToRegistrationFile ";
std::cerr << "JointHistogramAfterRegistrationFile ";
std::cerr
<< "NumberOfHistogramBinsForWritingTheMutualInformationHistogramMetric";
std::cerr << std::endl;
return EXIT_FAILURE;
}
using PixelType = unsigned char;
using InternalPixelType = float;
using InterpolatorType =
using RegistrationType =
using MetricType =
InternalImageType>;
registration->SetOptimizer(optimizer);
registration->SetTransform(transform);
registration->SetInterpolator(interpolator);
const unsigned int numberOfHistogramBins = std::stoi(argv[7]);
histogramSize[0] = numberOfHistogramBins;
histogramSize[1] = numberOfHistogramBins;
metric->SetHistogramSize(histogramSize);
const unsigned int numberOfParameters = transform->GetNumberOfParameters();
using ScalesType = MetricType::ScalesType;
ScalesType scales(numberOfParameters);
scales.Fill(1.0);
metric->SetDerivativeStepLengthScales(scales);
observer->m_JointHistogramWriter.SetMetric(metric);
registration->SetMetric(metric);
fixedImageReader->SetFileName(argv[1]);
movingImageReader->SetFileName(argv[2]);
using FixedNormalizeFilterType =
using MovingNormalizeFilterType =
using GaussianFilterType =
fixedSmoother->SetVariance(2.0);
movingSmoother->SetVariance(2.0);
fixedNormalizer->SetInput(fixedImageReader->GetOutput());
movingNormalizer->SetInput(movingImageReader->GetOutput());
fixedSmoother->SetInput(fixedNormalizer->GetOutput());
movingSmoother->SetInput(movingNormalizer->GetOutput());
registration->SetFixedImage(fixedSmoother->GetOutput());
registration->SetMovingImage(movingSmoother->GetOutput());
try
{
fixedNormalizer->Update();
}
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
registration->SetFixedImageRegion(
fixedNormalizer->GetOutput()->GetBufferedRegion());
using ParametersType = RegistrationType::ParametersType;
ParametersType initialParameters(transform->GetNumberOfParameters());
initialParameters[0] = 0.0;
initialParameters[1] = 0.0;
registration->SetInitialTransformParameters(initialParameters);
optimizer->SetMaximumStepLength(4.00);
optimizer->SetMinimumStepLength(0.01);
optimizer->SetRelaxationFactor(0.90);
optimizer->SetNumberOfIterations(200);
optimizer->MaximizeOn();
optimizer->AddObserver(itk::IterationEvent(), observer);
observer->SetInitialHistogramFile(argv[5]);
if (std::stoi(argv[4]))
{
observer->SetWriteHistogramsAfterEveryIteration(true);
}
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
ParametersType finalParameters = registration->GetLastTransformParameters();
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const 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;
std::string histogramAfter(argv[6]);
try
{
observer->m_JointHistogramWriter.WriteHistogramFile(histogramAfter);
}
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
using ResampleFilterType =
finalTransform->SetParameters(finalParameters);
finalTransform->SetFixedParameters(transform->GetFixedParameters());
resample->SetTransform(finalTransform);
resample->SetInput(movingImageReader->GetOutput());
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
using OutputPixelType = unsigned char;
using CastFilterType =
writer->SetFileName(argv[3]);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
try
{
writer->Update();
}
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
return EXIT_SUCCESS;
}