394 lines
11 KiB
C++
394 lines
11 KiB
C++
#include "itkReg23.h"
|
|
|
|
|
|
|
|
#include "itkTimeProbesCollectorBase.h"
|
|
|
|
namespace itk
|
|
{
|
|
|
|
itkReg23::itkReg23()
|
|
{
|
|
// Metrics
|
|
NCCmetric = MetricType::New();
|
|
MImetric = MIMetricType::New();
|
|
|
|
|
|
// Optimizers
|
|
PowellOptimizer = PowellOptimizerType::New();
|
|
AmoebaOptimizer = AmoebaOptimizerType::New();
|
|
ExhaustiveOptimizer = ExhaustiveOptimizerType::New();
|
|
|
|
// Registration
|
|
registration = RegistrationType::New();
|
|
|
|
//internalTransforms
|
|
Transform1 = TransformType::New();
|
|
Transform2 = TransformType::New();
|
|
IsocTransform = TransformType::New();
|
|
|
|
registration->SetTransform1(Transform1);
|
|
registration->SetTransform2(Transform2);
|
|
|
|
// to be provided by the user
|
|
m_r23Meta = nullptr;
|
|
m_Volume = nullptr;
|
|
m_PA = nullptr;
|
|
m_LAT = nullptr;
|
|
m_InternalTransf1 = nullptr;
|
|
m_InternalTransf2 = nullptr;
|
|
m_filter1 = nullptr;
|
|
m_filter2 = nullptr;
|
|
m_TransformMetaInfo = nullptr;
|
|
m_OptimizerObserver = nullptr;
|
|
|
|
m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file.
|
|
|
|
this->setROISizeToZero();
|
|
}
|
|
|
|
itkReg23::~itkReg23()
|
|
{
|
|
|
|
}
|
|
|
|
void itkReg23::setROISizeToZero(){
|
|
|
|
ImageType3D::SizeType zeroSize;
|
|
zeroSize.Fill(0);
|
|
|
|
m_roiAutoReg1.SetSize(zeroSize);
|
|
m_roiAutoReg2.SetSize(zeroSize);
|
|
|
|
}
|
|
|
|
void itkReg23::PrintSelf(std::ostream& os, Indent indent) const
|
|
{
|
|
Superclass::PrintSelf(os, indent);
|
|
}
|
|
|
|
void itkReg23::InitializeRegistration()
|
|
{
|
|
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
|
|
|
|
if (!m_PA) {
|
|
itkExceptionMacro(<< "PA topogram not present");
|
|
}
|
|
|
|
if (!m_LAT) {
|
|
itkExceptionMacro(<< "LAT topogram not present");
|
|
}
|
|
|
|
if (!m_Volume) {
|
|
itkExceptionMacro(<< "CT data not present");
|
|
}
|
|
|
|
if (!m_r23Meta) {
|
|
itkExceptionMacro(<< "r23Meta data not present");
|
|
}
|
|
|
|
if (!m_OptimizerObserver) {
|
|
itkExceptionMacro(<< "m_OptimizerObserver data not present");
|
|
}
|
|
|
|
if (!m_InternalTransf1) {
|
|
itkExceptionMacro(<< "m_InternalTransf1 data not present");
|
|
}
|
|
|
|
if (!m_InternalTransf2) {
|
|
itkExceptionMacro(<< "m_InternalTransf2 data not present");
|
|
}
|
|
|
|
if (!m_filter1) {
|
|
itkExceptionMacro(<< "m_filter1 data not present");
|
|
}
|
|
|
|
if (!m_filter2) {
|
|
itkExceptionMacro(<< "m_filter2 data not present");
|
|
}
|
|
|
|
if (!m_TransformMetaInfo) {
|
|
itkExceptionMacro(<< "m_TransformMetaInfo data not present");
|
|
}
|
|
|
|
if(!m_interpolator1){
|
|
itkExceptionMacro(<< "m_interpolator1 data not present");
|
|
}
|
|
|
|
if(!m_interpolator2){
|
|
itkExceptionMacro(<< "m_interpolator2 data not present");
|
|
}
|
|
|
|
switch (m_r23Meta->GetOptimizerType()) {
|
|
|
|
case tOptimizerTypeEnum::POWELL:
|
|
std::cout<< "Using POWELL Optimizer" <<std::endl;
|
|
|
|
PowellOptimizer->SetMaximize(false); // for NCC and MI
|
|
PowellOptimizer->SetMaximumIteration(m_PowellMeta->GetMaxIterations());
|
|
PowellOptimizer->SetMaximumLineIteration(m_PowellMeta->GetMaximumLineInteration()); // for Powell's method
|
|
PowellOptimizer->SetStepLength(m_PowellMeta->GetStepLength());
|
|
PowellOptimizer->SetStepTolerance(m_PowellMeta->GetStepTolerance());
|
|
PowellOptimizer->SetValueTolerance(m_PowellMeta->GetValueTolerance());
|
|
PowellOptimizer->RemoveAllObservers();
|
|
PowellOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
|
|
m_OptimizerObserver->setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum::POWELL);
|
|
|
|
registration->SetOptimizer(PowellOptimizer);
|
|
|
|
break;
|
|
|
|
case tOptimizerTypeEnum::AMOEBA:
|
|
std::cout<< "Using AMOEBA Optimizer" <<std::endl;
|
|
|
|
AmoebaOptimizer->SetMinimize(true);
|
|
AmoebaOptimizer->SetMaximumNumberOfIterations(m_AmoebaMeta->GetMaxIterations());
|
|
AmoebaOptimizer->SetParametersConvergenceTolerance(m_AmoebaMeta->GetParametersConvergenceTolerance());
|
|
AmoebaOptimizer->SetFunctionConvergenceTolerance(m_AmoebaMeta->GetFunctionConvergenceTolerance());
|
|
AmoebaOptimizer->SetAutomaticInitialSimplex(false);
|
|
//Initial size of the simplex (otherwise it is a very small number and optimizer stops immeditaly)
|
|
// 2 mm / 2 degrees seems to be a good value which performs well.
|
|
if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 0) {
|
|
itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
|
|
return;
|
|
}else if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 3){
|
|
AmoebaOptimizerType::ParametersType simplexDelta3dof(3);
|
|
for (int i = 0; i < GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()); i++) {
|
|
simplexDelta3dof[i] = m_AmoebaMeta->GetSimplexDelta();
|
|
}
|
|
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta3dof);
|
|
|
|
}else if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 6){
|
|
AmoebaOptimizerType::ParametersType simplexDelta6dof(6);
|
|
for (int i = 0; i < GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()); i++) {
|
|
simplexDelta6dof[i] = m_AmoebaMeta->GetSimplexDelta();
|
|
}
|
|
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta6dof);
|
|
}
|
|
|
|
AmoebaOptimizer->RemoveAllObservers();
|
|
AmoebaOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
|
|
m_OptimizerObserver->setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum::AMOEBA);
|
|
|
|
registration->SetOptimizer(AmoebaOptimizer);
|
|
|
|
break;
|
|
|
|
case tOptimizerTypeEnum::EXHAUSTIVE:
|
|
std::cout<< "Using Extensive Optimizer" <<std::endl;
|
|
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
|
|
}
|
|
|
|
switch (m_r23Meta->GetMetricType()) {
|
|
|
|
case tMetricTypeEnum::MI:
|
|
std::cout<< "Using MI Metric" <<std::endl;
|
|
|
|
registration->SetMetric(MImetric);
|
|
MImetric->SetNumberOfHistogramBins(m_MIMeta->GetNumberOfHistogramBins());
|
|
MImetric->ComputeGradientOff();
|
|
MImetric->SetMaxTranslation(m_MIMeta->GetMaxTranslation());
|
|
|
|
break;
|
|
|
|
case tMetricTypeEnum::NCC:
|
|
std::cout<< "Using NCC Metric" <<std::endl;
|
|
registration->SetMetric(NCCmetric);
|
|
NCCmetric->ComputeGradientOff();
|
|
NCCmetric->SetSubtractMean(true);
|
|
NCCmetric->SetMaxTranslation(m_NCCMeta->GetMaxTranslation());
|
|
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
|
|
}
|
|
|
|
registration->SetFixedImage1(m_PA);
|
|
registration->SetFixedImage2(m_LAT);
|
|
registration->SetMovingImage(m_Volume);
|
|
|
|
|
|
if( this->GetroiAutoReg1().GetNumberOfPixels() == 0 ){
|
|
auto defaultRegion = m_PA->GetBufferedRegion();
|
|
this->SetroiAutoReg1(defaultRegion);
|
|
std::cout<< "Image1 region not defined. Using the whole image buffered region" <<std::endl;
|
|
}
|
|
if( this->GetroiAutoReg2().GetNumberOfPixels() == 0 ){
|
|
auto defaultRegion = m_LAT->GetBufferedRegion();
|
|
this->SetroiAutoReg2(defaultRegion);
|
|
std::cout<< "Image2 region not defined. Using the whole image buffered region" <<std::endl;
|
|
}
|
|
|
|
registration->SetFixedImageRegion1(m_roiAutoReg1);
|
|
registration->SetFixedImageRegion2(m_roiAutoReg2);
|
|
|
|
registration->SetTransformMetaInfo(m_r23Meta);
|
|
|
|
registration->SetinternalMeta1(m_InternalTransf1);
|
|
registration->SetinternalMeta2(m_InternalTransf2);
|
|
|
|
registration->SetInterpolator1(m_interpolator1);
|
|
registration->SetInterpolator2(m_interpolator2);
|
|
|
|
registration->SetFilter1(m_filter1);
|
|
registration->SetFilter2(m_filter2);
|
|
|
|
|
|
IsocTransform = TransformType::New();
|
|
IsocTransform->SetComputeZYX(true);
|
|
IsocTransform->SetIdentity();
|
|
|
|
IsocTransform->SetRotation(
|
|
m_TransformMetaInfo->GetR()[0],
|
|
m_TransformMetaInfo->GetR()[1],
|
|
m_TransformMetaInfo->GetR()[2]
|
|
);
|
|
|
|
TransformType::OutputVectorType TranslV;
|
|
TranslV[0] = m_TransformMetaInfo->GetT()[0];
|
|
TranslV[1] = m_TransformMetaInfo->GetT()[1];
|
|
TranslV[2] = m_TransformMetaInfo->GetT()[2];
|
|
|
|
IsocTransform->SetTranslation(TranslV);
|
|
|
|
registration->SetIsocIECTransform(IsocTransform);
|
|
|
|
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
|
|
}
|
|
|
|
int itkReg23::StartRegistration(std::string extraInfo)
|
|
{
|
|
|
|
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
|
|
// TODO: Check if the registartion pipeline has been initialized
|
|
using ParametersType = RegistrationType::ParametersType;
|
|
|
|
time_t t = time(0); // get time now
|
|
struct tm* now = localtime(&t);
|
|
|
|
bool useDumpToFile = false;
|
|
|
|
char buffer[255];
|
|
strftime(buffer, 255, "test-%F-%H%M%S.txt", now);
|
|
|
|
std::fstream fs;
|
|
|
|
if (useDumpToFile) {
|
|
fs.open(buffer, std::fstream::out);
|
|
fs << extraInfo;
|
|
fs << "Value\tX\tY\tZ " << std::endl;
|
|
|
|
}
|
|
|
|
// Start the registration
|
|
// ~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Create a timer to record calculation time.
|
|
itk::TimeProbesCollectorBase timer;
|
|
|
|
if (verbose) {
|
|
std::cout << "Starting the registration now..." << std::endl;
|
|
}
|
|
|
|
try {
|
|
// timer.Start("Registration");
|
|
// Start the registration.
|
|
registration->StartRegistration();
|
|
// timer.Stop("Registration");
|
|
} catch (itk::ExceptionObject& err) {
|
|
|
|
registration->ResetPipeline();
|
|
std::cout << "ExceptionObject caught !" << std::endl;
|
|
std::cout << err << std::endl;
|
|
return -1;
|
|
}
|
|
|
|
fs.close();
|
|
|
|
|
|
itk::Optimizer::ParametersType finalParameters =
|
|
this->GetCurrentPosition();
|
|
std::cout<<" FinalPars: "<< finalParameters <<std::endl;
|
|
|
|
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
double itkReg23::GetOptimizerValue()
|
|
{
|
|
return m_OptmizerValue;
|
|
}
|
|
|
|
/** Get the cost function value for the current transform*/
|
|
Optimizer::ParametersType itkReg23::GetCurrentPosition(){
|
|
|
|
itk::Optimizer::ParametersType finalParameters;
|
|
|
|
switch (m_r23Meta->GetOptimizerType()) {
|
|
|
|
case tOptimizerTypeEnum::POWELL:
|
|
finalParameters = PowellOptimizer->GetCurrentPosition();
|
|
m_OptmizerValue = PowellOptimizer->GetValue();
|
|
break;
|
|
|
|
case tOptimizerTypeEnum::AMOEBA:
|
|
finalParameters = AmoebaOptimizer->GetCurrentPosition();
|
|
m_OptmizerValue = AmoebaOptimizer->GetValue();
|
|
break;
|
|
|
|
case tOptimizerTypeEnum::EXHAUSTIVE:
|
|
finalParameters = ExhaustiveOptimizer->GetCurrentPosition();
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
|
|
}
|
|
|
|
return finalParameters;
|
|
}
|
|
|
|
double itkReg23::GetCurrentMetricValue()
|
|
{
|
|
|
|
switch (m_r23Meta->GetMetricType()) {
|
|
|
|
case tMetricTypeEnum::MI:
|
|
if(!MImetric){
|
|
itkExceptionMacro(<< "MImetric not present");
|
|
return -1.;
|
|
} else {
|
|
return MImetric->GetValue();
|
|
}
|
|
|
|
break;
|
|
|
|
case tMetricTypeEnum::NCC:
|
|
if(!NCCmetric){
|
|
itkExceptionMacro(<< "NCCmetric not present");
|
|
return -1.;
|
|
} else {
|
|
return NCCmetric->GetValue();
|
|
}
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
|
|
}
|
|
|
|
return -1.;
|
|
|
|
}
|
|
|
|
|
|
}
|