fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

Added dialogs for approval and Automatic Registration. Added abort of auto registration

Browse Source
This commit is contained in:
Proton local user
2023-02-22 15:09:26 +01:00
parent d82bd12ef9
commit 085f8e0900
10 changed files with 308 additions and 139 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
reg23Topograms/itkDTRrecon/CMakeLists.txt
View File

@ -23,6 +23,7 @@ SET(SRCS
SET(HDR
itkImageProcessor.h
itkQtIterationUpdate.h
itkgSiddonJacobsRayCastInterpolateImageFunction.h
itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
vtkContourTopogramProjectionFilter.h

2
reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h
View File

@ -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...

2
reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx
View File

@ -25,6 +25,8 @@
#include "itkNormalizeImageFilter.h"
#include "itkTranslationTransform.h"
#include "itkProgressReporter.h"
#include <limits>
namespace itk {

25
reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx
View File

@ -211,11 +211,11 @@ NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Calcul
} else {
measure = NumericTraits<MeasureType>::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<TFixedImage, TMovingImage>::GetVal
auto oldprecision = std::cout.precision();
MeasureType measure1 = CalculateMeasure(1);
std::cout.precision(std::numeric_limits<double>::digits10 + 2);
std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
std::cout.precision(oldprecision);
// std::cout.precision(std::numeric_limits<double>::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<double>::digits10 + 2);
std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
// std::cout.precision(std::numeric_limits<double>::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;
}

9
reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h
View File

@ -25,6 +25,10 @@
#include "itkProcessObject.h"
#include "itkSingleValuedNonLinearOptimizer.h"
#include "itkgTwoImageToOneImageMetric.h"
#include "itkObject.h"
#include "itkEventObject.h"
#include <iostream>
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;

45
reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx
View File

@ -18,6 +18,8 @@
#ifndef itkTwoProjectionImageRegistrationMethod_hxx
#define itkTwoProjectionImageRegistrationMethod_hxx
#include "itkCommand.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
namespace itk {
@ -175,6 +177,32 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize
m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters);
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
itk::Object *caller, const itk::EventObject &event)
{
itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(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 <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
{
// typedef itk::MemberCommand<TwoProjectionImageRegistrationMethod> 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<TFixedImage, TMovingImage>::PrintSelf(
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GenerateData()
{
this->StartRegistration();
}

16
reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx
View File

@ -97,14 +97,14 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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;

32
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -1,4 +1,4 @@

/*
@ -25,6 +25,8 @@ gfattori 08.11.2021
#include "itkOrientImageFilter.h"
#include <gdcmIPPSorter.h>
#include "itkCommand.h"
#include "itkTimeProbesCollectorBase.h"
#include <string>
@ -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()" <<std::endl;
if(m_DRTImage1MetaInfo == NULL ||
m_DRTImage2MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
@ -3124,5 +3134,11 @@ void itkImageProcessor::SetFullROI(bool fullROI)
// TODO: Remove this function when ROI functinalitz has been implemented.
}
void itkImageProcessor::SetMaxNumberOfIterations(int iNum)
{
m_MaxNumberOfIterations = iNum;
// TODO: Remove this function when ROI functinalitz has been implemented.
}
} // end namespace itk

130
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -45,112 +45,15 @@ gfattori 08.11.2021
#include "gdcmGlobal.h"
#include "vtkImageData.h"
#include "vtkTransform.h"
#include "vtkSmartPointer.h"
#include "itkQtIterationUpdate.h"
namespace itk
{
class CommandIterationUpdate : public itk::Command {
// TODO: Move to own files.
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
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<OptimizerPointer>(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<double>::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<Self>;
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<OptimizerPointer>(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<double>::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<std::string> );
@ -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<PixelType3D, Dimension>;
@ -314,7 +225,7 @@ private:
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
/** 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;

185
reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h Normal file
View File

@ -0,0 +1,185 @@
#ifndef ITKQTITERATIONUPDATE_H
#define ITKQTITERATIONUPDATE_H
#include "itkCommand.h"
#include <QObject>
#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<Self>;
itkNewMacro(CommandIterationUpdate);
using InternalPixelType = float;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
using ProcesssType = typename itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
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<OptimizerPointer>(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<double>::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<Self>;
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<OptimizerPointer>(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<double>::digits10 + 2);
*os << optimizer->GetCurrentValue();
os->precision(oldprecision);
*os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl;
}
return;
}
};
}
#endif