From 085f8e090067b4f68f791af65fdda9e3b23a0d00 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Wed, 22 Feb 2023 15:09:26 +0100 Subject: [PATCH] Added dialogs for approval and Automatic Registration. Added abort of auto registration --- reg23Topograms/itkDTRrecon/CMakeLists.txt | 1 + .../itkDTRrecon/autoreg/itkDRTHelpers.h | 2 +- ...ualInformationTwoImageToOneImageMetric.hxx | 2 + ...zedCorrelationTwoImageToOneImageMetric.hxx | 25 +-- .../itkTwoProjectionImageRegistrationMethod.h | 9 + ...tkTwoProjectionImageRegistrationMethod.hxx | 45 +++++ .../autoreg/itkgTwoImageToOneImageMetric.hxx | 16 +- .../itkDTRrecon/itkImageProcessor.cpp | 32 ++- .../itkDTRrecon/itkImageProcessor.h | 130 ++---------- .../itkDTRrecon/itkQtIterationUpdate.h | 185 ++++++++++++++++++ 10 files changed, 308 insertions(+), 139 deletions(-) create mode 100644 reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index 60032bc..fd530ce 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -23,6 +23,7 @@ SET(SRCS SET(HDR itkImageProcessor.h + itkQtIterationUpdate.h itkgSiddonJacobsRayCastInterpolateImageFunction.h itkgSiddonJacobsRayCastInterpolateImageFunction.hxx vtkContourTopogramProjectionFilter.h diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h b/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h index e8e563e..36580ed 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h @@ -35,7 +35,7 @@ static const char* ImageOrientationStrings[] = { // constant for converting degrees into radians const float dtr = itk::Math::pi_over_180; -static const bool verbose = true; +static const bool verbose = false; /* this is in the end a IEC to HFS... * but we keep this for ourselves... diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx index 748df79..089c7f9 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx @@ -25,6 +25,8 @@ #include "itkNormalizeImageFilter.h" #include "itkTranslationTransform.h" +#include "itkProgressReporter.h" + #include namespace itk { diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx index f76c15f..06e9c34 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx @@ -211,11 +211,11 @@ NormalizedCorrelationTwoImageToOneImageMetric::Calcul } else { measure = NumericTraits::Zero; } - std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl; - std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl; - std::cout << "max coordinates: " << maxX << " " << maxX << " " << maxZ << std::endl; - std::cout << "min coordinates: " << minX << " " << minY << " " << minZ << std::endl; +// std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl; +// std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl; +// std::cout << "max coordinates: " << maxX << " " << maxX << " " << maxZ << std::endl; +// std::cout << "min coordinates: " << minX << " " << minY << " " << minZ << std::endl; return measure; } @@ -275,18 +275,19 @@ NormalizedCorrelationTwoImageToOneImageMetric::GetVal auto oldprecision = std::cout.precision(); MeasureType measure1 = CalculateMeasure(1); - std::cout.precision(std::numeric_limits::digits10 + 2); - std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; - std::cout.precision(oldprecision); +// std::cout.precision(std::numeric_limits::digits10 + 2); +// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout.precision(oldprecision); MeasureType measure2 = CalculateMeasure(2); - std::cout.precision(std::numeric_limits::digits10 + 2); - std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout.precision(std::numeric_limits::digits10 + 2); +// std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; auto measure = (measure1 + measure2) / 2.0; - std::cout << "Measure: " << measure << std::endl; - std::cout << "=======================================================" << std::endl; - std::cout.precision(oldprecision); + +// std::cout << "Measure: " << measure << std::endl; +// std::cout << "=======================================================" << std::endl; +// std::cout.precision(oldprecision); return measure; } diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h index de7a291..3dd6241 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -25,6 +25,10 @@ #include "itkProcessObject.h" #include "itkSingleValuedNonLinearOptimizer.h" #include "itkgTwoImageToOneImageMetric.h" + +#include "itkObject.h" +#include "itkEventObject.h" + #include namespace itk { @@ -231,6 +235,11 @@ protected: /** Provides derived classes with the ability to set this private var */ itkSetMacro(LastOptimizerParameters, ParametersType); + /** Entry-point for internal ITK filter observer. **/ + void OnInternalFilterProgressReceptor(itk::Object *caller, + const itk::EventObject &event); + + private: MetricPointer m_Metric; OptimizerType::Pointer m_Optimizer; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx index fa77cd7..b486ea0 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -18,6 +18,8 @@ #ifndef itkTwoProjectionImageRegistrationMethod_hxx #define itkTwoProjectionImageRegistrationMethod_hxx +#include "itkCommand.h" + #include "itkTwoProjectionImageRegistrationMethod.h" namespace itk { @@ -175,6 +177,32 @@ void TwoProjectionImageRegistrationMethod::Initialize m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters); } +template +void TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor( + itk::Object *caller, const itk::EventObject &event) +{ + itk::ProcessObject *filter = dynamic_cast(caller); + + if(!itk::ProgressEvent().CheckEvent(&event) || !filter) + return; + + if (filter) + { +// double p = m_CurrentProgressStart; +// // filter->GetProgress() in [0;1] +// p += m_CurrentProgressSpan * filter->GetProgress(); +// // NOTE: filter is buggy - not in [0;1] if multi-threading is activated! +// if (p > (m_CurrentProgressStart + m_CurrentProgressSpan)) +// p = m_CurrentProgressStart + m_CurrentProgressSpan; +// TaskProgressInfo(m_CurrentProgressDirection, p); + +// if (m_CancelRequest) +// filter->SetAbortGenerateData(true); // may be handled by filter + + std::cout << "Porca Madonna " << std::endl; + } +} + /* * Starts the Registration Process */ @@ -182,6 +210,22 @@ template void TwoProjectionImageRegistrationMethod::StartRegistration() { +// typedef itk::MemberCommand ITKCommandType; +// typedef ITKCommandType::Pointer MemberPointer ; + + +// ITKCommandType::Pointer cmd = ITKCommandType::New(); +// itk::Command cmd= itk::Command::Pointer(this); + +// cmd->SetCallbackFunction(this, +// (ITKCommandType::Pointer)&TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor); + + itk::ProgressReporter progress( + this, + 1, 1000,100); + +// this->SetAbortGenerateData(true); + ParametersType empty(1); empty.Fill(0.0); try { @@ -253,6 +297,7 @@ void TwoProjectionImageRegistrationMethod::PrintSelf( template void TwoProjectionImageRegistrationMethod::GenerateData() { + this->StartRegistration(); } diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index 0e9c600..90a9548 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -97,14 +97,14 @@ bool gTwoImageToOneImageMetric::SetTransformParameter break; } - std::cout << "New Transform Parameters = " << std::endl; - std::cout << " Translation X = " << TranslationAlongX << " mm" << std::endl; - std::cout << " Translation Y = " << TranslationAlongY << " mm" << std::endl; - std::cout << " Translation Z = " << TranslationAlongZ << " mm" << std::endl; - std::cout << " Rotation Along X = " << RotationAlongX / dtr << " deg" << std::endl; - std::cout << " Rotation Along Y = " << RotationAlongY / dtr << " deg" << std::endl; - std::cout << " Rotation Along Z = " << RotationAlongZ / dtr << " deg" << std::endl; - std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; +// std::cout << "New Transform Parameters = " << std::endl; +// std::cout << " Translation X = " << TranslationAlongX << " mm" << std::endl; +// std::cout << " Translation Y = " << TranslationAlongY << " mm" << std::endl; +// std::cout << " Translation Z = " << TranslationAlongZ << " mm" << std::endl; +// std::cout << " Rotation Along X = " << RotationAlongX / dtr << " deg" << std::endl; +// std::cout << " Rotation Along Y = " << RotationAlongY / dtr << " deg" << std::endl; +// std::cout << " Rotation Along Z = " << RotationAlongZ / dtr << " deg" << std::endl; +// std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; transformParameters[0] = RotationAlongX; transformParameters[1] = RotationAlongY; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index e613dd6..aa2c7b8 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -1,4 +1,4 @@ - + /* @@ -25,6 +25,8 @@ gfattori 08.11.2021 #include "itkOrientImageFilter.h" #include +#include "itkCommand.h" + #include "itkTimeProbesCollectorBase.h" #include @@ -200,13 +202,14 @@ itkImageProcessor::itkImageProcessor() optimizer = OptimizerType::New(); amoebaoptimizer = AmoebaOptimizerType::New(); - optimizerObserver = CommandIterationUpdate::New(); + exhaustiveOptimizer = ExhaustiveOptimizerType::New(); exhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New(); registration = RegistrationType::New(); + optimizerObserver->SetProcess(registration); if (m_UseMutualInformation) { registration->SetMetric(mimetric); @@ -1368,7 +1371,7 @@ itkImageProcessor::CalculateInternalTransform( m_OutputTransform; } -void itkImageProcessor::InitializeRegistration(int maximumIteration, double stepLength, double maxTranslation, eDegreeOfFreedomType dof) +void itkImageProcessor::InitializeRegistration(double stepLength, double maxTranslation, eDegreeOfFreedomType dof) { std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; @@ -1396,6 +1399,7 @@ void itkImageProcessor::InitializeRegistration(int maximumIteration, double step metric->DebugOn(); } + registration->SetFixedImage1(m_PASourceDupli->GetOutput()); registration->SetFixedImage2(m_LATSourceDupli->GetOutput()); registration->SetMovingImage(m_VolumeSourceDupli->GetOutput()); @@ -1553,7 +1557,7 @@ void itkImageProcessor::InitializeRegistration(int maximumIteration, double step // The metric will return 0 or a "large" negative number in case of high correspondence of the images // Thus, we need to minimze optimizer->SetMaximize(false); // for NCC and MI - optimizer->SetMaximumIteration(maximumIteration); + optimizer->SetMaximumIteration(m_MaxNumberOfIterations); optimizer->SetMaximumLineIteration(4); // for Powell's method optimizer->SetStepLength(stepLength); optimizer->SetStepTolerance(0.01); @@ -1563,10 +1567,13 @@ void itkImageProcessor::InitializeRegistration(int maximumIteration, double step optimizer->Print(std::cout); } + optimizer->RemoveAllObservers(); - optimizer->AddObserver(itk::IterationEvent(), optimizerObserver); + optimizer->AddObserver(itk::AnyEvent(), optimizerObserver); registration->SetOptimizer(optimizer); + registration->RemoveAllObservers(); + } else { amoebaoptimizer->SetMinimize(true); @@ -1658,11 +1665,13 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) } try { - timer.Start("Registration"); +// timer.Start("Registration"); // Start the registration. registration->StartRegistration(); - timer.Stop("Registration"); +// timer.Stop("Registration"); } catch (itk::ExceptionObject& err) { + + registration->ResetPipeline(); std::cout << "ExceptionObject caught !" << std::endl; std::cout << err << std::endl; return -1; @@ -1728,7 +1737,8 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) void itkImageProcessor::InitializeProjector() { -// std::cout<< "itkImageProcessor::InitializeProjector()" <; - itkNewMacro(Self); - -protected: - CommandIterationUpdate() = default; - -public: - using OptimizerType = itk::PowellOptimizer; - using OptimizerPointer = const OptimizerType*; - - using AmoebaOptimizerType = itk::AmoebaOptimizer; - ; - using AmoebaOptimizerPointer = const OptimizerType*; - - 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 - { - auto optimizer = dynamic_cast(object); - if (typeid(event) == typeid(itk::IterationEvent)) { - - //Feedback from the optimizer executed at the end of every itteration - // currently just print the result into the cout. We might add - // functionality to register and emit signals to update the UI. - - std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl; - auto oldprecision = std::cout.precision(); - std::cout.precision(std::numeric_limits::digits10 + 2); - std::cout << "Similarity: " << optimizer->GetValue() << std::endl; - std::cout.precision(oldprecision); - std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl; - } - return; - } -}; - -class ExhaustiveCommandIterationUpdate : public itk::Command { - // TODO: Move to own files. - -public: - using Self = ExhaustiveCommandIterationUpdate; - using Superclass = itk::Command; - using Pointer = itk::SmartPointer; - itkNewMacro(Self); - - std::ostream* os; - -protected: - ExhaustiveCommandIterationUpdate() = default; - -public: - using OptimizerType = itk::ExhaustiveOptimizer; - ; - using OptimizerPointer = const OptimizerType*; - - void set_stream(std::ostream& stream) - { - os = &stream; - } - - 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 - { - auto optimizer = dynamic_cast(object); - if (typeid(event) == typeid(itk::IterationEvent)) { - - //crude LPS to IEC transform for HFS. - auto position = optimizer->GetCurrentPosition(); - - auto oldprecision = os->precision(); - os->precision(std::numeric_limits::digits10 + 2); - *os << optimizer->GetCurrentValue(); - os->precision(oldprecision); - *os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl; - } - return; - } -}; - - class ITK_EXPORT itkImageProcessor : public itk::Object { @@ -169,6 +72,11 @@ public: /** Run-time type information (and related methods). */ itkTypeMacro(itkImageProcessor, Object); + + CommandIterationUpdate::Pointer + optimizerObserver; + + /** Input data load methods*/ int load3DSerieFromFolder(const char* ); int load3DSerieFromFiles( std::vector ); @@ -228,7 +136,7 @@ public: void SetUserRotationsLPS(double, double, double); /** Initialize the registration pipeline*/ - void InitializeRegistration(int, double, double, eDegreeOfFreedomType); + void InitializeRegistration(double, double, eDegreeOfFreedomType); /** Start the registration process*/ int StartRegistration(std::string extraInfo); @@ -277,7 +185,10 @@ public: void SetOptimizer(std::string); void SetMetric(std::string); void SetFullROI(bool); + void SetMaxNumberOfIterations(int); + /** Optimizer which tries to find the minimn (Powell Optimizer)*/ + using OptimizerType = itk::PowellOptimizer; protected: /** Various pixel types */ @@ -299,7 +210,7 @@ private: /** Fill Meta after 3D volume load */ int fill3DVolumeMeta(InternalImageType::Pointer, - tPatOrientation); + tPatOrientation); /** Image types */ using ImageType2D = itk::Image; @@ -314,7 +225,7 @@ private: using ResampleFilterType = itk::ResampleImageFilter; /** Optimizer which tries to find the minimn (Powell Optimizer)*/ - using OptimizerType = itk::PowellOptimizer; + // using OptimizerType = itk::PowellOptimizer; using AmoebaOptimizerType = itk::AmoebaOptimizer; /** Optimizer which samples the whol space */ using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer; @@ -376,21 +287,19 @@ private: imageDRT2In; RegistrationType::Pointer - registration; + registration; MetricType::Pointer - metric; + metric; MIMetricType::Pointer - mimetric; + mimetric; OptimizerType::Pointer - optimizer; + optimizer; AmoebaOptimizerType::Pointer - amoebaoptimizer; - CommandIterationUpdate::Pointer - optimizerObserver; + amoebaoptimizer; ExhaustiveOptimizerType::Pointer - exhaustiveOptimizer; + exhaustiveOptimizer; ExhaustiveCommandIterationUpdate::Pointer - exhaustiveOptimizerObserver; + exhaustiveOptimizerObserver; DuplicatorType::Pointer m_LATSourceDupli, @@ -490,6 +399,7 @@ private: double m_OptmizerValue; + int m_MaxNumberOfIterations; bool m_UseExhaustiveOptmizer; bool m_UseAmeobaOptimizer; diff --git a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h new file mode 100644 index 0000000..5c055d4 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h @@ -0,0 +1,185 @@ +#ifndef ITKQTITERATIONUPDATE_H +#define ITKQTITERATIONUPDATE_H + +#include "itkCommand.h" +#include +#include "itkPowellOptimizer.h" +#include "itkAmoebaOptimizer.h" +#include "itkExhaustiveOptimizer.h" + +#include "itkTwoProjectionImageRegistrationMethod.h" + + +class QObjectIterationUpdate : public QObject{ + Q_OBJECT +public: + QObjectIterationUpdate(){ + bAbortProcessCommand = false; + }; + + bool getAbortFlag(){ + return bAbortProcessCommand; + }; + void setAbortFlag(bool bVal){ + bAbortProcessCommand = bVal; + }; + +public slots: + void onAbortProcessCommand(){ + std::cout << " Abort STOCAZZO " << std::endl; + bAbortProcessCommand = true; + }; + void onIteration(double dI,double dX,double dY,double dZ){ + emit + sendRegistrationProgress(dI,dX,dY,dZ); + + } +private: + bool + bAbortProcessCommand; + + + + + +signals: +void sendRegistrationProgress(double,double,double,double); + +}; + +namespace itk +{ + + +class CommandIterationUpdate : public itk::Command { + // TODO: Move to own files. + + constexpr static unsigned int Dimension = 3; + +public: + QObjectIterationUpdate *objIterUpdate; + + using Self = CommandIterationUpdate; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer; + itkNewMacro(CommandIterationUpdate); + + using InternalPixelType = float; + using InternalImageType = itk::Image; + using ProcesssType = typename itk::TwoProjectionImageRegistrationMethod; + using ProcesssPointer = typename ProcesssType::Pointer; + + /** Set/Get the Process. */ + itkSetObjectMacro(Process, ProcesssType); + itkGetConstObjectMacro(Process, ProcesssType); + +private: + +protected: + CommandIterationUpdate() { + objIterUpdate = new QObjectIterationUpdate; + } + ProcesssPointer m_Process; + +public: + using OptimizerType = itk::PowellOptimizer; + using OptimizerPointer = const OptimizerType*; + + using AmoebaOptimizerType = itk::AmoebaOptimizer; + using AmoebaOptimizerPointer = const OptimizerType*; + + 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(objIterUpdate->getAbortFlag()){ + std::cout << "Abortisci per piacere" << std::endl; + objIterUpdate->setAbortFlag(false); + throw itk::ProcessAborted(); + } +// std::cout << "Progress: " << this->m_Process->GetAbortGenerateData() << std::endl; + + auto optimizer = dynamic_cast(object); + + if (typeid(event) == typeid(itk::IterationEvent)) { + + //Feedback from the optimizer executed at the end of every itteration + // currently just print the result into the cout. We might add + // functionality to register and emit signals to update the UI. + + std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl; + auto oldprecision = std::cout.precision(); + std::cout.precision(std::numeric_limits::digits10 + 2); + std::cout << "Similarity: " << optimizer->GetValue() << std::endl; + std::cout.precision(oldprecision); + std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl; + + objIterUpdate->onIteration( + optimizer->GetCurrentIteration()+1, + optimizer->GetCurrentPosition()[0], + optimizer->GetCurrentPosition()[2], + -optimizer->GetCurrentPosition()[1] + ); + } + return; + } +}; + +class ExhaustiveCommandIterationUpdate : public itk::Command { + // TODO: Move to own files. + +public: + using Self = ExhaustiveCommandIterationUpdate; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer; + itkNewMacro(Self); + + std::ostream* os; + +protected: + ExhaustiveCommandIterationUpdate() = default; + +public: + using OptimizerType = itk::ExhaustiveOptimizer; + ; + using OptimizerPointer = const OptimizerType*; + + void set_stream(std::ostream& stream) + { + os = &stream; + } + + 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 + { + auto optimizer = dynamic_cast(object); + if (typeid(event) == typeid(itk::IterationEvent)) { + + //crude LPS to IEC transform for HFS. + auto position = optimizer->GetCurrentPosition(); + + auto oldprecision = os->precision(); + os->precision(std::numeric_limits::digits10 + 2); + *os << optimizer->GetCurrentValue(); + os->precision(oldprecision); + *os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl; + } + return; + } +}; + +} + +#endif +