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

Merge branch 'scoutRTUIDevelREG' into 'ScoutRTUIDevel'

Browse Source 
Restructuring of reg23 code and ROI fix

See merge request cpt_bioeng/drt!29
This commit is contained in:
via_r
2023-05-23 15:06:08 +02:00
parent 8d2976ebbc 47fe531aae
commit 205f2f9db1
11 changed files with 843 additions and 612 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
reg23Topograms/itkDTRrecon/CMakeLists.txt
View File

@ -19,7 +19,7 @@ SET(SRCS
itkImageProcessorHelpers.cpp
vtkContourTopogramProjectionFilter.cxx
DRTMetaInformation.cpp
itkReg23.cpp
)
SET(HDR
@ -30,6 +30,7 @@ SET(HDR
itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
vtkContourTopogramProjectionFilter.h
DRTMetaInformation.h
itkReg23.h
)
ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR})

4
reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
View File

@ -547,7 +547,9 @@ R23MetaInformation::
R23MetaInformation (){
m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
m_OptimizerType = tOptimizerTypeEnum::POWELL;
m_MetricType = tMetricTypeEnum::NCC;
m_MaxIterations = 6;
}

36
reg23Topograms/itkDTRrecon/DRTMetaInformation.h
View File

@ -35,6 +35,19 @@ typedef enum eDegreeOfFreedomType {
OTHER
} tDegreeOfFreedomEnum;
typedef enum eOptimizerType{
POWELL,
AMOEBA,
EXHAUSTIVE
} tOptimizerTypeEnum;
typedef enum eMetricType{
NCC,
MI
} tMetricTypeEnum;
inline int GetNumberOfDOF(eDegreeOfFreedomType dof)
{
switch (dof) {
@ -597,16 +610,35 @@ public:
itkTypeMacro(R23MetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
void PrintSelf(std::ostream& os, itk::Indent indent) const override;
itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
itkSetEnumMacro(OptimizerType, tOptimizerTypeEnum);
itkGetEnumMacro(OptimizerType, tOptimizerTypeEnum);
itkSetEnumMacro(MetricType, tMetricTypeEnum);
itkGetEnumMacro(MetricType, tMetricTypeEnum);
itkSetEnumMacro(MaxIterations, int);
itkGetEnumMacro(MaxIterations, int);
protected:
tDegreeOfFreedomEnum
m_DegreeOfFreedom;
m_DegreeOfFreedom;
tOptimizerTypeEnum
m_OptimizerType;
tMetricTypeEnum
m_MetricType;
int
m_MaxIterations;
/** Default Constructor **/
R23MetaInformation ();

3
reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt
View File

@ -1,6 +1,5 @@
SET(LIB_NAME autoreg)
SET(CMAKE_CXX_FLAGS "-std=c++11 -fPIC")
#SET(CMAKE_CXX_FLAGS "-std=c++11 -fPIC")
find_package(ITK REQUIRED)
include(${ITK_USE_FILE})

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

@ -149,6 +149,12 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
//metric1->UseAllPixelsOn();
metric1->SetNumberOfHistogramBins(50);
// InternalImageType::IndexType pIndex;
// pIndex[0]=200;
// pIndex[1]=300;
// InternalImageType::PixelType pPixel = fixedImage1->GetPixel(pIndex);
// std::cout<<"MATTES PIXEL: "<<pPixel <<std::endl;
metric1->SetFixedImage(fixedImage1);
metric1->SetMovingImage(this->m_Filter1->GetOutput());

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

@ -76,7 +76,6 @@ NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Calcul
this->m_NumberOfPixelsCounted = 0;
int numberOfPixelsVisited = 0;
//filter->Update();
filter->Update();
MovingImageConstPointer imageDRTIn = filter->GetOutput();
//this->WriteImage(imageDRTIn, "", "");
@ -189,15 +188,6 @@ NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Calcul
++ti;
}
// sfm = sfm - m_meanCurr_f * sm - m_meanCurr_m * sf + this->m_NumberOfPixelsCounted * m_meanCurr_m * m_meanCurr_f;
// const RealType denom = -1 * sqrt(m_varianceCurr_f * m_varianceCurr_m) * this->m_NumberOfPixelsCounted;
// if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
// measure = sfm / denom;
// } else {
// measure = NumericTraits<MeasureType>::ZeroValue();
// }
if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0) {
sff -= (sf * sf / this->m_NumberOfPixelsCounted);
smm -= (sm * sm / this->m_NumberOfPixelsCounted);
@ -274,6 +264,11 @@ NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetVal
// Calculate the measure value between fixed image 1 and the moving image
auto oldprecision = std::cout.precision();
// std::cout<<"Region " <<this->GetFixedImageRegion1() <<std::endl;
// std::cout<<"Buffered Region " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
// std::cout<<"Buffered Region " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
MeasureType measure1 = CalculateMeasure(1);
// std::cout.precision(std::numeric_limits<double>::digits10 + 2);
// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;

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

@ -250,6 +250,12 @@ void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
m_FixedImage2->GetSource()->Update();
}
//std::cout<<"FixedImageRegion1 " <<this->GetFixedImageRegion1() <<std::endl;
//std::cout<<"m_FixedImage1 buffered " <<this->m_FixedImage1->GetBufferedRegion() <<std::endl;
//std::cout<<"Moving buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
//std::cout<<"GetFixedImage1 buffered " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
//std::cout<<"Filter1 buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
// Make sure the FixedImageRegion is within the FixedImage buffered region
if (!m_FixedImageRegion1.Crop(m_FixedImage1->GetBufferedRegion())) {
itkExceptionMacro(<< "FixedImageRegion1 does not overlap the fixed image buffered region");

685
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -40,7 +40,7 @@ itkImageProcessor::itkImageProcessor()
{
iNWUnits = 1;
// std::cout<<"itkImageProcessor::itkImageProcessor contructor"<<std::endl;
transform1 = TransformType::New();
transform1->SetIdentity();
transform2 = TransformType::New();
@ -53,11 +53,6 @@ itkImageProcessor::itkImageProcessor()
interpolator1 = InterpolatorType::New();
interpolator2 = InterpolatorType::New();
// TZero[0]=TZero[1]=TZero[2]=0.;
// RZero[0]=RZero[1]=RZero[2]=0.;
toVTK2D1 = ITKtoVTKFilterType::New();
toVTK2D2 = ITKtoVTKFilterType::New();
toVTKLocalizer1 = ITKtoVTKFilterType::New();
@ -141,38 +136,16 @@ itkImageProcessor::itkImageProcessor()
m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D);
m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D);
// setup the Automatic Registration
metric = MetricType::New();
mimetric = MIMetricType::New();
optimizer = OptimizerType::New();
amoebaoptimizer = AmoebaOptimizerType::New();
optimizerObserver = CommandIterationUpdate::New();
ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();
// TOOOOODOOOOO optimizerObserver->SetProcess(registration);
m_R23 = itkReg23::New();
m_R23->SetOptimizerObserver(optimizerObserver);
exhaustiveOptimizer = ExhaustiveOptimizerType::New();
exhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();
registration = RegistrationType::New();
optimizerObserver->SetProcess(registration);
if (m_UseMutualInformation) {
registration->SetMetric(mimetric);
} else {
registration->SetMetric(metric);
}
registration->SetOptimizer(optimizer);
registration->SetTransform1(transform1);
registration->SetTransform2(transform2);
registration->SetInterpolator1(interpolator1);
registration->SetInterpolator2(interpolator2);
m_UseExhaustiveOptmizer = false; //if the exhaustive optimizer shal be used, ovverrides ameoba
m_UseAmeobaOptimizer = false; //if the ameoba optimzer shall be used, otherwise Powell. Overriden by Exhaustive option.
m_UseFullROI = true; // if the full image ROI shall be used
m_UseMutualInformation = false; //oterwise NCC is used
m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file.
}
itkImageProcessor::~itkImageProcessor()
@ -261,11 +234,6 @@ double itkImageProcessor::GetPanelOffset2(){
void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof)
{
m_r23MetaInfo->SetDegreeOfFreedom(dof);
}
tDegreeOfFreedomEnum
itkImageProcessor::GetDegreeOfFreedom()
{
@ -1031,7 +999,12 @@ int itkImageProcessor::load2D(const char * pcFName){
m_Duplicator->SetInputImage(caster2DInput->GetOutput() );
m_Duplicator->Update();
// InternalImageType::IndexType pIndex;
// pIndex[0]=200;
// pIndex[1]=300;
// InternalImageType::PixelType pPixel =
// m_Duplicator->GetOutput()->GetPixel(pIndex);
//std::cout<<"processro PIXEL: "<<pPixel <<std::endl;
// std::cout<< " ////////////// 2D Topo BEG " <<std::endl;
// std::cout<<"GetDirection()"<<m_Duplicator->GetOutput()->GetDirection()<<std::endl;
@ -1048,27 +1021,6 @@ int itkImageProcessor::load2D(const char * pcFName){
//int itkImageProcessor::query2DimageReconstructionDiameter(const char *pcFName){
// /* Check if we can open the file */
// gdcm::Reader R;
// R.SetFileName(pcFName);
// if(!R.Read())
// { cerr<<"ERROR: cannot read file: "<<pcFName<<endl;
// return -1;
// }
// /* Check if it's a localizer image */
// const gdcm::File &file = R.GetFile();
// const gdcm::DataSet &ds = file.GetDataSet();
// char* sTmpString = new char [255];
// strcpy( sTmpString,gGetStringValueFromTag(gdcm::Tag(0x0018,0x1100), ds));
// return std::atoi(sTmpString);
//}
double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
{
@ -1134,312 +1086,6 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
void itkImageProcessor::InitializeRegistration(
double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof)
{
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
if (!m_PASourceDupli) {
itkExceptionMacro(<< "PA topogram not present");
}
if (!m_LATSourceDupli) {
itkExceptionMacro(<< "LAT topogram not present");
}
if (!m_VolumeSourceDupli) {
itkExceptionMacro(<< "CT data not present");
}
metric->ComputeGradientOff();
metric->SetSubtractMean(true);
metric->SetMaxTranslation(maxTranslation);
// m_TransformMetaInfo->SetDegreeOfFreedom(dof);
m_r23MetaInfo->SetDegreeOfFreedom(dof);
mimetric->ComputeGradientOff();
if (verbose) {
metric->DebugOn();
}
registration->SetFixedImage1(m_PASourceDupli->GetOutput());
registration->SetFixedImage2(m_LATSourceDupli->GetOutput());
registration->SetMovingImage(m_VolumeSourceDupli->GetOutput());
//TODO: Here we could set the ROI for image registartion
if (m_UseFullROI == false) {
auto region1 = m_PASourceDupli->GetOutput()->GetBufferedRegion();
auto region2 = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
std::cout << "region1 of opti - NO ROI " << region1 << std::endl;
std::cout << "region2 of opti - NO ROI " << region2 << std::endl;
registration->SetFixedImageRegion1(roiAutoReg1);
registration->SetFixedImageRegion2(roiAutoReg2);
mimetric->SetFixedImageRegion1(roiAutoReg1);
mimetric->SetFixedImageRegion2(roiAutoReg2);
std::cout << "using ROIs for Auto Reg " << std::endl;
std::cout << "region1 of opti " << roiAutoReg1 << std::endl;
std::cout << "region2 of opti " << roiAutoReg2 << std::endl;
}else{
auto region1 = m_PASourceDupli->GetOutput()->GetBufferedRegion();
auto region2 = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
std::cout << "region1 of opti - NO ROI " << region1 << std::endl;
std::cout << "region2 of opti - NO ROI " << region2 << std::endl;
registration->SetFixedImageRegion1(region1);
registration->SetFixedImageRegion2(region2);
mimetric->SetFixedImageRegion1(roiAutoReg1);
mimetric->SetFixedImageRegion2(roiAutoReg2);
}
// registration->SetTransformMetaInfo(m_TransformMetaInfo);
registration->SetTransformMetaInfo(m_r23MetaInfo);
registration->SetinternalMeta1(m_InternalTransf1);
registration->SetinternalMeta2(m_InternalTransf2);
registration->SetFilter1(filter1);
registration->SetFilter2(filter2);
IsocTransf = TransformType::New();
IsocTransf->SetComputeZYX(true);
IsocTransf->SetIdentity();
IsocTransf->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];
IsocTransf->SetTranslation(TranslV);
registration->SetIsocIECTransform(IsocTransf);
if (verbose) {
registration->DebugOn();
registration->Print(std::cout);
}
if (!m_UseExhaustiveOptmizer) {
if (!m_UseAmeobaOptimizer) {
// 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(m_MaxNumberOfIterations);
optimizer->SetMaximumLineIteration(4); // for Powell's method
optimizer->SetStepLength(stepLength);
optimizer->SetStepTolerance(0.01);
optimizer->SetValueTolerance(0.000002);
if (verbose) {
optimizer->DebugOn();
optimizer->Print(std::cout);
}
optimizer->RemoveAllObservers();
optimizer->AddObserver(itk::AnyEvent(), optimizerObserver);
registration->SetOptimizer(optimizer);
registration->RemoveAllObservers();
} else {
amoebaoptimizer->SetMinimize(true);
amoebaoptimizer->SetMaximumNumberOfIterations(100);
amoebaoptimizer->SetParametersConvergenceTolerance(0.1);
amoebaoptimizer->SetFunctionConvergenceTolerance(0.000002);
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.
OptimizerType::ParametersType simplexDelta(3);
int numberOfDOF = GetNumberOfDOF(dof);
if (numberOfDOF == 0) {
itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
}
for (int i = 0; i < numberOfDOF; i++) {
simplexDelta[i] = 2.0;
}
amoebaoptimizer->SetInitialSimplexDelta(simplexDelta);
amoebaoptimizer->RemoveAllObservers();
amoebaoptimizer->AddObserver(itk::IterationEvent(), optimizerObserver);
registration->SetOptimizer(amoebaoptimizer);
}
} else {
//only support for 3DOF
const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
const double stepLength = 0.1;
//auto parameters = transform1->GetParameters();
//transform1->SetParameters(parameters);
ExhaustiveOptimizerType::StepsType steps(3);
steps[0] = numberOfSteps;
steps[1] = numberOfSteps;
steps[2] = numberOfSteps;
exhaustiveOptimizer->SetNumberOfSteps(steps);
exhaustiveOptimizer->SetStepLength(stepLength);
exhaustiveOptimizer->RemoveAllObservers();
exhaustiveOptimizer->AddObserver(itk::IterationEvent(), exhaustiveOptimizerObserver);
exhaustiveOptimizer->SetStepLength(stepLength);
registration->SetOptimizer(exhaustiveOptimizer);
}
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
}
int itkImageProcessor::StartRegistration(std::string extraInfo)
{
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
// TODO: Check if the registartion pipeline has been initialized
using ParametersType = RegistrationType::ParametersType;
//auto startParameters = transform1->GetParameters();
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 || m_UseExhaustiveOptmizer) {
fs.open(buffer, std::fstream::out);
fs << extraInfo;
fs << "Value\tX\tY\tZ " << std::endl;
if (m_UseExhaustiveOptmizer) {
exhaustiveOptimizerObserver->set_stream(fs);
}
}
// 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;
}
auto oldprecision = fs.precision();
if (m_UseExhaustiveOptmizer) {
fs.precision(std::numeric_limits<double>::digits10 + 2);
fs << " MinimumMetricValue: " << exhaustiveOptimizer->GetMinimumMetricValue() << std::endl;
fs << " MaximumMetricValue: " << exhaustiveOptimizer->GetMaximumMetricValue() << std::endl;
fs.precision(oldprecision);
fs << " MinimumMetricValuePosition: " << exhaustiveOptimizer->GetMinimumMetricValuePosition() << std::endl;
fs << " MaximumMetricValuePosition: " << exhaustiveOptimizer->GetMaximumMetricValuePosition() << std::endl;
fs << " StopConditionDescription: " << exhaustiveOptimizer->GetStopConditionDescription() << std::endl;
fs << " MaximumMetricValuePosition: " << exhaustiveOptimizer->GetMaximumMetricValuePosition() << std::endl;
} else if (useDumpToFile) {
fs.precision(std::numeric_limits<double>::digits10 + 2);
fs << " FinalMetricValue: " << optimizer->GetValue() << std::endl; //same as GetCurrentCost
fs.precision(oldprecision);
fs << " FinalPosition: " << optimizer->GetCurrentPosition() << std::endl;
fs << " Iterations: " << optimizer->GetCurrentIteration() << std::endl;
fs << " StopConditionDescription: " << optimizer->GetStopConditionDescription() << std::endl;
}
fs.close();
m_OptmizerValue = optimizer->GetValue();
auto finalParameters = optimizer->GetCurrentPosition();
std::cout<<"REG FINAL PARS: "<<finalParameters<<std::endl;
// auto finalParameters = transform1->GetParameters();
// std::cout<<"-->startParameters "<<startParameters<<std::endl;
// std::cout<<"-->finalParameters1 "<<transform1->GetParameters()<<std::endl;
// std::cout<<"-->finalParameters2 "<<transform2->GetParameters()<<std::endl;
// //optimizer->GetCurrentPosition();
// finalParameters = finalParameters - startParameters;
// std::cout<<"startParameters "<<startParameters<<std::endl;
// std::cout<<"finalParameters "<<finalParameters<<std::endl;
// std::cout<<"finalParametersDelta "<<finalParameters<<std::endl;
// std::cout<<"transform1 "<<transform1->GetCenter()<<std::endl;
// std::cout<<"transform1 "<<transform1->GetTranslation()<<std::endl;
// std::cout<<"transform2 "<<transform2->GetCenter()<<std::endl;
// std::cout<<"transform2 "<<transform2->GetTranslation()<<std::endl;
// TransformType::Pointer offsetTransform = TransformType::New();
// offsetTransform->SetParameters(finalParameters);
// CalculateExternalUserTransform(offsetTransform, m_DRTImage1MetaInfo);
// const double translationAlongX = finalParameters[3];
// const double translationAlongY = finalParameters[4];
// const double translationAlongZ = finalParameters[5];
// const double rotationAlongX = finalParameters[0];
// const double rotationAlongY = finalParameters[1];
// const double rotationAlongZ = finalParameters[2];
// const int numberOfIterations = optimizer->GetCurrentIteration();
// std::cout << "Result = " << 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 << " Translation X = " << translationAlongX << " mm" << std::endl;
// std::cout << " Translation Y = " << translationAlongY << " mm" << std::endl;
// std::cout << " Translation Z = " << translationAlongZ << " mm" << std::endl;
// std::cout << " Number Of Iterations = " << numberOfIterations << std::endl;
// std::cout << " Metric value = " << m_OptmizerValue << std::endl;
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
return 0;
}
Optimizer::ParametersType
itkImageProcessor::GetFinalR23Parameters(){
@ -1452,32 +1098,31 @@ itkImageProcessor::GetFinalR23Parameters(){
Optimizer::ParametersType pPars(6);
pPars.fill(0.);
if(optimizer == nullptr){
return pPars;
}
itk::Optimizer::ParametersType currentPosition =
m_R23->GetCurrentPosition();
switch (m_r23MetaInfo->GetDegreeOfFreedom()) {
case THREE_DOF:
pPars[3] = optimizer->GetCurrentPosition()[0]
pPars[3] = currentPosition[0]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = optimizer->GetCurrentPosition()[1]
pPars[4] = currentPosition[1]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = optimizer->GetCurrentPosition()[2]
pPars[5] = currentPosition[2]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
break;
case SIX_DOF:
pPars[0] = optimizer->GetCurrentPosition()[0]
pPars[0] = currentPosition[0]
- m_TransformMetaInfo->GetHiddenRotations()[0];
pPars[1] = optimizer->GetCurrentPosition()[1]
pPars[1] = currentPosition[1]
- m_TransformMetaInfo->GetHiddenRotations()[1];
pPars[2] = optimizer->GetCurrentPosition()[2]
pPars[2] = currentPosition[2]
- m_TransformMetaInfo->GetHiddenRotations()[2];
pPars[3] = optimizer->GetCurrentPosition()[3]
pPars[3] = currentPosition[3]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = optimizer->GetCurrentPosition()[4]
pPars[4] = currentPosition[4]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = optimizer->GetCurrentPosition()[5]
pPars[5] = currentPosition[5]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
break;
@ -1492,6 +1137,41 @@ itkImageProcessor::GetFinalR23Parameters(){
}
void itkImageProcessor::SetOptimizer(std::string optimizer)
{
if (optimizer.compare("Powell") == 0) {
m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::POWELL);
} else if (optimizer.compare("Amoeba") == 0) {
m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::AMOEBA);
} else if (optimizer.compare("Exhaustive") == 0) {
m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::EXHAUSTIVE);
} else {
itkExceptionMacro(<< "Unkown optimzer string : " << optimizer);
}
}
void itkImageProcessor::SetMetric(std::string metric)
{
if (metric.compare("NCC") == 0) {
m_r23MetaInfo->SetMetricType(tMetricTypeEnum::NCC);
} else if (metric.compare("MI") == 0) {
m_r23MetaInfo->SetMetricType(tMetricTypeEnum::MI);
} else {
itkExceptionMacro(<< "Unkown metric string : " << metric);
}
}
void itkImageProcessor::SetMaxNumberOfIterations(int iNum)
{
m_r23MetaInfo->SetMaxIterations(iNum);
}
void itkImageProcessor::InitializeProjector()
{
@ -2288,44 +1968,45 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK()
toVTKLocalizer1->SetInput(rescaler1->GetOutput());
toVTKLocalizer1->Update();
if(false) {
using ImageRegionType3D = ImageType3D::RegionType;
using SizeType3D = ImageRegionType3D::SizeType;
using ImageDirectionType3D = ImageType3D::DirectionType;
// if(false) {
//// using ImageRegionType3D = ImageType3D::RegionType;
//// using SizeType3D = ImageRegionType3D::SizeType;
//// using ImageDirectionType3D = ImageType3D::DirectionType;
ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion();
SizeType3D size3D = region3D.GetSize();
ImageType3D::PointType origin3D = rescaler1->GetOutput()->GetOrigin();
const itk::Vector<double, 3> resolution3D = rescaler1->GetOutput()->GetSpacing();
ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection();
//// ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion();
//// SizeType3D size3D = region3D.GetSize();
//// //ImageType3D::PointType origin3D = rescaler1->GetOutput()->GetOrigin();
//// const itk::Vector<double, 3> resolution3D = rescaler1->GetOutput()->GetSpacing();
//// ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection();
/* calculate image size in 3D Space */
using VectorType = itk::Vector<double, 3>;
VectorType Dest;
Dest[0]=(size3D[0]-1) * resolution3D[0];
Dest[1]=(size3D[1]-1) * resolution3D[1];
Dest[2]=(size3D[2]-1) * resolution3D[2];
//// /* calculate image size in 3D Space */
//// using VectorType = itk::Vector<double, 3>;
//// VectorType Dest;
//// Dest[0]=(size3D[0]-1) * resolution3D[0];
//// Dest[1]=(size3D[1]-1) * resolution3D[1];
//// Dest[2]=(size3D[2]-1) * resolution3D[2];
ImageType3D::PointType LastVoxelPosition =
origin3D + (imagDirection * Dest);
// std::cout<<" -------- Localizer 1 ITK --------"<<std::endl;
// std::cout<<"size3D " <<size3D <<std::endl;
// std::cout<<"resolution3D " <<resolution3D <<std::endl;
// std::cout<<"imagDirection " <<imagDirection <<std::endl;
// std::cout<<"First Voxel position " <<origin3D <<std::endl;
// std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
// // ImageType3D::PointType LastVoxelPosition =
// // origin3D + (imagDirection * Dest);
// // std::cout<<" -------- Localizer 1 ITK --------"<<std::endl;
// // std::cout<<"size3D " <<size3D <<std::endl;
// // std::cout<<"resolution3D " <<resolution3D <<std::endl;
// // std::cout<<"imagDirection " <<imagDirection <<std::endl;
// // std::cout<<"First Voxel position " <<origin3D <<std::endl;
// // std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
double* dBounds = toVTKLocalizer1->GetOutput()->GetBounds();
// // double* dBounds = toVTKLocalizer1->GetOutput()->GetBounds();
// std::cout<< "-------- Localizer 1 VTK --------" <<std::endl;
// std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
// <<dBounds[2]<<" "<<dBounds[3]<<" "
// <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
// std::cout<< "-------- -------- --------" <<std::endl;
}
// // std::cout<< "-------- Localizer 1 VTK --------" <<std::endl;
// // std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
// // <<dBounds[2]<<" "<<dBounds[3]<<" "
// // <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
// // std::cout<< "-------- -------- --------" <<std::endl;
// }
return
toVTKLocalizer1->GetOutput();
@ -2350,42 +2031,42 @@ vtkImageData* itkImageProcessor::GetLocalizer2VTK()
if(true) {
using ImageRegionType3D = ImageType3D::RegionType;
using SizeType3D = ImageRegionType3D::SizeType;
using ImageDirectionType3D = ImageType3D::DirectionType;
// if(true) {
// using ImageRegionType3D = ImageType3D::RegionType;
// using SizeType3D = ImageRegionType3D::SizeType;
// using ImageDirectionType3D = ImageType3D::DirectionType;
ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion();
SizeType3D size3D = region3D.GetSize();
ImageType3D::PointType origin3D = rescaler2->GetOutput()->GetOrigin();
const itk::Vector<double, 3> resolution3D = rescaler2->GetOutput()->GetSpacing();
ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection();
// ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion();
// SizeType3D size3D = region3D.GetSize();
// ImageType3D::PointType origin3D = rescaler2->GetOutput()->GetOrigin();
// const itk::Vector<double, 3> resolution3D = rescaler2->GetOutput()->GetSpacing();
// ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection();
/* calculate image size in 3D Space */
using VectorType = itk::Vector<double, 3>;
VectorType Dest;
Dest[0]=(size3D[0]-1) * resolution3D[0];
Dest[1]=(size3D[1]-1) * resolution3D[1];
Dest[2]=(size3D[2]-1) * resolution3D[2];
// /* calculate image size in 3D Space */
// using VectorType = itk::Vector<double, 3>;
// VectorType Dest;
// Dest[0]=(size3D[0]-1) * resolution3D[0];
// Dest[1]=(size3D[1]-1) * resolution3D[1];
// Dest[2]=(size3D[2]-1) * resolution3D[2];
ImageType3D::PointType LastVoxelPosition =
origin3D + (imagDirection * Dest);
// ImageType3D::PointType LastVoxelPosition =
// origin3D + (imagDirection * Dest);
// std::cout<<" -------- Localizer 2 ITK --------"<<std::endl;
// std::cout<<"size3D " <<size3D <<std::endl;
// std::cout<<"resolution3D " <<resolution3D <<std::endl;
// std::cout<<"imagDirection " <<imagDirection <<std::endl;
// std::cout<<"First Voxel position " <<origin3D <<std::endl;
// std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
// // std::cout<<" -------- Localizer 2 ITK --------"<<std::endl;
// // std::cout<<"size3D " <<size3D <<std::endl;
// // std::cout<<"resolution3D " <<resolution3D <<std::endl;
// // std::cout<<"imagDirection " <<imagDirection <<std::endl;
// // std::cout<<"First Voxel position " <<origin3D <<std::endl;
// // std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
// double* dBounds = toVTKLocalizer2->GetOutput()->GetBounds();
// std::cout<< "-------- Localizer 2 VTK --------" <<std::endl;
// std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
// <<dBounds[2]<<" "<<dBounds[3]<<" "
// <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
// std::cout<< "-------- -------- --------" <<std::endl;
// // double* dBounds = toVTKLocalizer2->GetOutput()->GetBounds();
// // std::cout<< "-------- Localizer 2 VTK --------" <<std::endl;
// // std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
// // <<dBounds[2]<<" "<<dBounds[3]<<" "
// // <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
// // std::cout<< "-------- -------- --------" <<std::endl;
}
// }
@ -2474,11 +2155,6 @@ vtkImageData* itkImageProcessor::GetProjection1VTK()
vtkMatrix4x4 * itkImageProcessor::GetProjection1VTKTransform()
{
// std::cout<< "Composed " <<std::endl;
// interpolator1->GetComposedTransform()->Print(std::cout);
// std::cout<< "Inverse " <<std::endl;
// interpolator1->GetInverseTransform()->Print(std::cout);
m_Projection1VTKTransform->Identity();
for(int iIdx = 0; iIdx<3 ; iIdx++){
@ -2718,22 +2394,6 @@ itkImageProcessor::GetCompleteIsocentricTransform
return nullptr;
}
// std::cout<<" --- GetCompleteIsocentricTransform --- "<<std::endl;
// std::cout<< "R: "<<
// m_TransformMetaInfo->GetR()[0] <<" "<<
// m_TransformMetaInfo->GetR()[1] <<" "<<
// m_TransformMetaInfo->GetR()[2] <<" "<<std::endl;
// std::cout<< "T: "<<
// m_TransformMetaInfo->GetT()[0] <<" "<<
// m_TransformMetaInfo->GetT()[1] <<" "<<
// m_TransformMetaInfo->GetT()[2] <<" "<<std::endl;
// std::cout<< "Import Offset: "<<
// m_CTMetaInfo->GetImportOffset()[0]<<" "<<
// m_CTMetaInfo->GetImportOffset()[1]<<" "<<
// m_CTMetaInfo->GetImportOffset()[2]<<" "<<std::endl;
// std::cout<<" --- --- --- "<<std::endl;
ImageType3D::PointType mISORTIEC =
m_CTMetaInfo->GetLPS2IECDirections() * m_RTMetaInfo->GetIsocenterLPS();
@ -2766,65 +2426,6 @@ itkImageProcessor::GetCompleteIsocentricTransform
}
double itkImageProcessor::GetOptimizerValue()
{
return m_OptmizerValue;
}
double itkImageProcessor::GetCurrentMetricValue()
{
//Note: this fails if the metric has not been propperly initilzed.
if (!m_UseMutualInformation) {
if (metric == NULL) {
return nan("");
}
return metric->GetValue();
} else {
if (mimetric == NULL) {
return nan("");
}
return mimetric->GetValue();
}
}
void itkImageProcessor::SetOptimizer(std::string optimizer)
{
if (optimizer.compare("Powell") == 0) {
m_UseAmeobaOptimizer = false;
m_UseExhaustiveOptmizer = false;
} else if (optimizer.compare("Amoeba") == 0) {
m_UseExhaustiveOptmizer = false;
m_UseAmeobaOptimizer = true;
} else if (optimizer.compare("Exhaustive") == 0) {
m_UseExhaustiveOptmizer = true;
m_UseAmeobaOptimizer = false;
} else {
itkExceptionMacro(<< "Unkown optimzer string : " << optimizer);
}
}
void itkImageProcessor::SetMetric(std::string metric)
{
if (metric.compare("NCC") == 0) {
m_UseMutualInformation = false;
} else if (metric.compare("MI") == 0) {
m_UseMutualInformation = true;
} else {
itkExceptionMacro(<< "Unkown metric string : " << optimizer);
}
}
void itkImageProcessor::SetFullROI(bool fullROI)
{
m_UseFullROI = 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.
}
void itkImageProcessor::SetRegionFixed1(
int iIdx0, int iIdx1, int iSz0, int iSz1){
@ -2844,8 +2445,10 @@ void itkImageProcessor::SetRegionFixed1(
roiAutoReg1.SetIndex(index1);
roiAutoReg1.SetSize(size1);
std::cout << "itkImageProcessor " << std::endl;
std::cout << roiAutoReg1 << std::endl;
//std::cout << "itkImageProcessor " << std::endl;
//std::cout << roiAutoReg1 << std::endl;
this->m_R23->SetroiAutoReg1(roiAutoReg1);
}
@ -2867,10 +2470,46 @@ void itkImageProcessor::SetRegionFixed2(
roiAutoReg2.SetIndex(index2);
roiAutoReg2.SetSize(size2);
std::cout << "itkImageProcessor " << std::endl;
std::cout << roiAutoReg2 << std::endl;
//std::cout << "itkImageProcessor " << std::endl;
//std::cout << roiAutoReg2 << std::endl;
this->m_R23->SetroiAutoReg2(roiAutoReg2);
}
void itkImageProcessor::ResetROIRegions(){
auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion();
auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
this->m_R23->SetroiAutoReg1(region1temp);
this->m_R23->SetroiAutoReg2(region2temp);
}
void itkImageProcessor::InizializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof){
m_r23MetaInfo->SetDegreeOfFreedom(dof);
m_R23->SetPA(m_PASourceDupli->GetOutput());
m_R23->SetLAT(m_LATSourceDupli->GetOutput());
m_R23->SetVolume(m_VolumeSourceDupli->GetOutput());
m_R23->Setr23Meta(m_r23MetaInfo);
m_R23->SetInternalTransf1(m_InternalTransf1);
m_R23->SetInternalTransf2(m_InternalTransf2);
m_R23->Setfilter1(filter1);
m_R23->Setfilter2(filter2);
m_R23->Setinterpolator1(interpolator1);
m_R23->Setinterpolator2(interpolator2);
m_R23->SetTransformMetaInfo(m_TransformMetaInfo);
m_R23->InitializeRegistration(
stepLength,maxTranslation,dof);
}
} // end namespace itk

107
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -30,13 +30,6 @@ gfattori 08.11.2021
#include "itkObjectFactory.h"
#include "itkSmartPointer.h"
#include "itkExhaustiveOptimizer.h"
#include "itkMutualInformationTwoImageToOneImageMetric.h"
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
#include "itkPowellOptimizer.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
#include "itkAmoebaOptimizer.h"
#include "itkImageToVTKImageFilter.h"
#include "itkImageFileReader.h"
@ -53,6 +46,8 @@ gfattori 08.11.2021
#include "itkImageProcessorHelpers.h"
#include "itkReg23.h"
namespace itk
{
@ -78,6 +73,33 @@ public:
CommandIterationUpdate::Pointer
optimizerObserver;
ExhaustiveCommandIterationUpdate::Pointer
ExhaustiveOptimizerObserver;
/** Sets the degree of freedom for optimzation*/
tDegreeOfFreedomEnum GetDegreeOfFreedom();
/** Returns user components only of the autoReg
* */
Optimizer::ParametersType
GetFinalR23Parameters();
/** Define the region to be used as ROI on fixed images */
void SetRegionFixed1(int,int,int,int);
void SetRegionFixed2(int,int,int,int);
void ResetROIRegions();
void InizializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof);
/** Maximum number of iterations for the optimizer */
void SetMaxNumberOfIterations(int);
itkGetMacro(R23, itkReg23::Pointer);
/** Set Optimizer type */
void SetOptimizer(std::string);
/** Set Metric type */
void SetMetric(std::string);
/** Set number of logic CPU to be made available to interpolators*/
void SetNumberOfWorkingUnits(int iN);
@ -116,10 +138,6 @@ public:
double GetPanelOffset1();
double GetPanelOffset2();
/** Sets the degree of freedom for omptimzation*/
void SetDegreeOfFreedom(tDegreeOfFreedomEnum);
tDegreeOfFreedomEnum GetDegreeOfFreedom();
void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double);
void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double);
@ -147,16 +165,6 @@ public:
void SetInitialTranslations(double, double, double);
void SetInitialRotations(double, double, double);
/** Returns user components only of the autoReg
* */
Optimizer::ParametersType
GetFinalR23Parameters();
/** Initialize the registration pipeline*/
void InitializeRegistration(double, double, eDegreeOfFreedomType);
/** Start the registration process*/
int StartRegistration(std::string extraInfo);
///** Get transform parameters for 3D Volume */
/** Get the complete transform, likely for export to SRO.
* This includes user and hidden contributions.
@ -215,25 +223,6 @@ public:
void WriteProjectionImages();
void Write2DImages();
/** Auto Reg23 methods */
/** Get the current cost function value from the optimizer*/
double GetOptimizerValue();
/** Get the cost function value for the current transform*/
double GetCurrentMetricValue();
/** Set Optimizer type */
void SetOptimizer(std::string);
/** Set Metric type */
void SetMetric(std::string);
/** Use full image as ROI */
void SetFullROI(bool);
/** Maximum number of iterations for the optimizer */
void SetMaxNumberOfIterations(int);
/** Define the region to be used as ROI on fixed images */
void SetRegionFixed1(int,int,int,int);
void SetRegionFixed2(int,int,int,int);
protected:
/** Various pixel types */
@ -268,17 +257,6 @@ private:
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
/** Optimizer which tries to find the minimun (Powell Optimizer)*/
using OptimizerType = itk::PowellOptimizer;
/** Optimizer which tries to find the minimn (Powell Optimizer)*/
using AmoebaOptimizerType = itk::AmoebaOptimizer;
/** Optimizer which samples the whol space */
using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer;
/** Metric to calculate similarites between images*/
using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
/** The thing which actuall does the image registration*/
using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
using RoiForRegistration = itk::ImageRegion<Dimension>;
@ -318,20 +296,7 @@ private:
imageDRT1In,
imageDRT2In;
RegistrationType::Pointer
registration;
MetricType::Pointer
metric;
MIMetricType::Pointer
mimetric;
OptimizerType::Pointer
optimizer;
AmoebaOptimizerType::Pointer
amoebaoptimizer;
ExhaustiveOptimizerType::Pointer
exhaustiveOptimizer;
ExhaustiveCommandIterationUpdate::Pointer
exhaustiveOptimizerObserver;
DuplicatorType::Pointer
m_LATSourceDupli,
@ -426,23 +391,19 @@ private:
m_InternalTransf1,
m_InternalTransf2;
double m_OptmizerValue;
int m_MaxNumberOfIterations;
RoiForRegistration
roiAutoReg1,
roiAutoReg2;
bool m_UseExhaustiveOptmizer;
bool m_UseAmeobaOptimizer;
bool m_UseFullROI;
bool m_UseMutualInformation;
bool m_UseDumptoFile;
int
iNWUnits;
itk::itkReg23::Pointer
m_R23;
};

393
reg23Topograms/itkDTRrecon/itkReg23.cpp Normal file
View File

@ -0,0 +1,393 @@
#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_UseFullROI = true; // if the full image ROI shall be used
m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file.
}
itkReg23::~itkReg23()
{
}
void itkReg23::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
void itkReg23::InitializeRegistration(
double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof)
{
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");
}
AmoebaOptimizerType::ParametersType simplexDelta(3);
ExhaustiveOptimizerType::StepsType steps(3);
const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
//const double stepLength = 0.1;
switch (m_r23Meta->GetOptimizerType()) {
case tOptimizerTypeEnum::POWELL:
std::cout<< "Using POWELL Optimizer" <<std::endl;
PowellOptimizer->SetMaximize(false); // for NCC and MI
PowellOptimizer->SetMaximumIteration(m_r23Meta->GetMaxIterations());
PowellOptimizer->SetMaximumLineIteration(4); // for Powell's method
PowellOptimizer->SetStepLength(stepLength);
PowellOptimizer->SetStepTolerance(0.01);
PowellOptimizer->SetValueTolerance(0.000002);
PowellOptimizer->RemoveAllObservers();
PowellOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
registration->SetOptimizer(PowellOptimizer);
break;
case tOptimizerTypeEnum::AMOEBA:
std::cout<< "Using AMOEBA Optimizer" <<std::endl;
AmoebaOptimizer->SetMinimize(true);
AmoebaOptimizer->SetMaximumNumberOfIterations(m_r23Meta->GetMaxIterations());
AmoebaOptimizer->SetParametersConvergenceTolerance(0.1);
AmoebaOptimizer->SetFunctionConvergenceTolerance(0.000002);
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(dof) == 0) {
itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
}
for (int i = 0; i < GetNumberOfDOF(dof); i++) {
simplexDelta[i] = 2.0;
}
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta);
AmoebaOptimizer->RemoveAllObservers();
AmoebaOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
registration->SetOptimizer(AmoebaOptimizer);
break;
case tOptimizerTypeEnum::EXHAUSTIVE:
std::cout<< "Using Extensive Optimizer" <<std::endl;
steps[0] = numberOfSteps;
steps[1] = numberOfSteps;
steps[2] = numberOfSteps;
ExhaustiveOptimizer->SetNumberOfSteps(steps);
ExhaustiveOptimizer->SetStepLength(stepLength);
registration->SetOptimizer(ExhaustiveOptimizer);
break;
default:
break;
}
switch (m_r23Meta->GetMetricType()) {
case tMetricTypeEnum::MI:
std::cout<< "Using MI Metric" <<std::endl;
registration->SetMetric(MImetric);
MImetric->ComputeGradientOff();
MImetric->SetMaxTranslation(maxTranslation);
break;
case tMetricTypeEnum::NCC:
std::cout<< "Using NCC Metric" <<std::endl;
registration->SetMetric(NCCmetric);
NCCmetric->ComputeGradientOff();
NCCmetric->SetSubtractMean(true);
NCCmetric->SetMaxTranslation(maxTranslation);
break;
default:
break;
}
registration->SetFixedImage1(m_PA);
registration->SetFixedImage2(m_LAT);
registration->SetMovingImage(m_Volume);
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);
// if (verbose) {
// registration->DebugOn();
// registration->Print(std::cout);
// }
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;
//auto startParameters = transform1->GetParameters();
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;
// if (r23Meta->GetOptimizerType() == tOptimizerTypeEnum::EXHAUSTIVE) {
// ExhaustiveOptimizerObserver->set_stream(fs);
// }
}
// 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.;
}
}

197
reg23Topograms/itkDTRrecon/itkReg23.h Normal file
View File

@ -0,0 +1,197 @@
#ifndef ITKREG23_H
#define ITKREG23_h
#include "itkCommand.h"
#include "itkExhaustiveOptimizer.h"
#include "itkMutualInformationTwoImageToOneImageMetric.h"
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
#include "itkPowellOptimizer.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
#include "itkAmoebaOptimizer.h"
#include "itkCommand.h"
#include "itkObject.h"
#include "itkObjectFactory.h"
#include "itkSmartPointer.h"
#include "itkImageProcessorHelpers.h"
#include "itkQtIterationUpdate.h"
namespace itk{
class ITK_EXPORT itkReg23 : public itk::Object
{
public:
/** Standard "Self" typedef. */
typedef itkReg23 Self;
/** Standard "Superclass" typedef. */
typedef Object Superclass;
/** Smart pointer typedef support */
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Method of creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(itkReg23, Object);
using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
using RoiForRegistration = itk::ImageRegion<Dimension>;
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
/* ---- User provided */
itkSetMacro(OptimizerObserver, CommandIterationUpdate::Pointer);
itkGetMacro(OptimizerObserver, CommandIterationUpdate::Pointer);
itkSetMacro(r23Meta, R23MetaInformation::Pointer);
itkGetMacro(r23Meta, R23MetaInformation::Pointer);
itkSetMacro(Volume, InternalImageType::Pointer);
itkGetMacro(Volume, InternalImageType::Pointer);
itkSetMacro(PA, InternalImageType::Pointer);
itkGetMacro(PA, InternalImageType::Pointer);
itkSetMacro(LAT, InternalImageType::Pointer);
itkGetMacro(LAT, InternalImageType::Pointer);
itkSetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer);
itkGetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer);
itkSetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer);
itkGetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer);
itkSetMacro(filter1, ChangeInformationFilterType::Pointer);
itkGetMacro(filter1, ChangeInformationFilterType::Pointer);
itkSetMacro(filter2, ChangeInformationFilterType::Pointer);
itkGetMacro(filter2, ChangeInformationFilterType::Pointer);
itkSetMacro(interpolator1, InterpolatorType::Pointer);
itkGetMacro(interpolator1, InterpolatorType::Pointer);
itkSetMacro(interpolator2, InterpolatorType::Pointer);
itkGetMacro(interpolator2, InterpolatorType::Pointer);
itkSetMacro(TransformMetaInfo, TransformMetaInformation::Pointer);
itkGetMacro(TransformMetaInfo, TransformMetaInformation::Pointer);
/* ---- User provided END */
itkSetMacro(roiAutoReg1, RoiForRegistration);
itkGetMacro(roiAutoReg1, RoiForRegistration);
itkSetMacro(roiAutoReg2, RoiForRegistration);
itkGetMacro(roiAutoReg2, RoiForRegistration);
/** Auto Reg23 methods */
/** Initialize the registration pipeline*/
void InitializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof);
/** Start the registration process*/
int StartRegistration(std::string extraInfo);
/** Get the current cost function value from the optimizer*/
double GetOptimizerValue();
/** Get the cost function value for the current transform*/
double GetCurrentMetricValue();
/** Get the cost function value for the current transform*/
Optimizer::ParametersType GetCurrentPosition();
/** Auto Reg23 methods END */
protected:
itkReg23();
virtual ~itkReg23();
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
itkReg23(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
/** Optimizer which tries to find the minimun (Powell Optimizer)*/
using PowellOptimizerType = itk::PowellOptimizer;
/** Optimizer which tries to find the minimn (Powell Optimizer)*/
using AmoebaOptimizerType = itk::AmoebaOptimizer;
/** Optimizer which samples the whol space */
using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer;
/** Metric to calculate similarites between images*/
using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
/** The thing which actuall does the image registration*/
using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
RegistrationType::Pointer
registration;
MetricType::Pointer
NCCmetric;
MIMetricType::Pointer
MImetric;
PowellOptimizerType::Pointer
PowellOptimizer;
AmoebaOptimizerType::Pointer
AmoebaOptimizer;
ExhaustiveOptimizerType::Pointer
ExhaustiveOptimizer;
/* ---- User provided */
R23MetaInformation::Pointer
m_r23Meta;
InternalImageType::Pointer
m_Volume,
m_PA,
m_LAT;
InternalTransformMetaInformation::Pointer
m_InternalTransf1,
m_InternalTransf2;
ChangeInformationFilterType::Pointer
m_filter1,
m_filter2;
InterpolatorType::Pointer
m_interpolator1,
m_interpolator2;
TransformMetaInformation::Pointer
m_TransformMetaInfo;
CommandIterationUpdate::Pointer
m_OptimizerObserver;
/* ---- User provided END */
TransformType::Pointer
IsocTransform,
Transform1,
Transform2;
RoiForRegistration
m_roiAutoReg1,
m_roiAutoReg2;
bool m_UseDumptoFile;
double m_OptmizerValue;
};
}
#endif