diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index c93893c..60032bc 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -2,6 +2,8 @@ cmake_minimum_required(VERSION 2.8.12) SET(LIB_NAME itkDTRrecon) +add_subdirectory(autoreg) + find_package(ITK REQUIRED) include(${ITK_USE_FILE}) INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS}) @@ -34,6 +36,7 @@ SET(LINK_LIBS ${VTK_LIBRARIES} ${ITK_LIBRARIES} gdcmMSFF + autoreg ) TARGET_LINK_LIBRARIES( diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index befd3b7..748cded 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -531,6 +531,8 @@ TransformMetaInformation (){ m_CurrentTransform.Fill(0.); + m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + } void diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 5ec3556..d79ffd1 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -11,15 +11,70 @@ namespace itk { -typedef enum eProjectionOrientationType -{NA = 0, LAT=2, PA=1} -tProjOrientationType; +typedef enum eProjectionOrientationType{ + NA = 0, + LAT=2, + PA=1 +} tProjOrientationType; /** Enum type for Orientation */ -typedef enum eImageOrientationType -{NotDefined = 0, HFS = 1, FFS = 2} -tPatOrientation; +typedef enum eImageOrientationType{ + NotDefined = 0, + HFS = 1, + FFS = 2 +} tPatOrientation; +typedef enum eDegreeOfFreedomType { + UNKNOWN = 0, + THREE_DOF, + SIX_DOF, + X_ONLY, + Y_ONLY, + Z_ONLY, + ROTATION_ONLY, + OTHER +} tDegreeOfFreedomEnum; + +inline int GetNumberOfDOF(eDegreeOfFreedomType dof) +{ + switch (dof) { + case SIX_DOF: + return 6; + break; + case THREE_DOF: + case ROTATION_ONLY: + return 3; + break; + case X_ONLY: + case Y_ONLY: + case Z_ONLY: + return 1; + break; + default: + return 0; + } +} + + +class DegreeOfFreedom : public itk::Object { + +public: + typedef DegreeOfFreedom Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + + itkNewMacro(Self); + +protected: + bool + m_TranslationX, + m_TranslationY, + m_TranslationZ, + m_RotationX, + m_RotationY, + m_RotationZ; + tDegreeOfFreedomEnum m_dof; +}; class TopogramImageMetaInformation : public itk::Object{ @@ -490,6 +545,8 @@ public: itkSetMacro(CurrentTransform,TransformMatrixType); itkGetMacro(CurrentTransform,TransformMatrixType); + itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); protected: @@ -503,6 +560,8 @@ protected: m_ReferenceTransform, m_CurrentTransform; + tDegreeOfFreedomEnum + m_DegreeOfFreedom; /** Default Constructor **/ TransformMetaInformation (); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 54ab2f0..e613dd6 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -14,6 +14,7 @@ gfattori 08.11.2021 #include "itkImageProcessor.h" #include "itkCreateObjectFunction.h" +#include "itkDRTHelpers.h" #include "itkVersion.h" #include #include @@ -24,6 +25,8 @@ gfattori 08.11.2021 #include "itkOrientImageFilter.h" #include +#include "itkTimeProbesCollectorBase.h" + #include @@ -31,34 +34,6 @@ gfattori 08.11.2021 namespace itk { -/* this is in the end a IEC to HFS... - * but we keep this for ourselves... - */ -static double Standard_DRT2LPS [9] = { - 1,0,0, - 0,0,-1, - 0,1,0 }; - -static double PAElementsIEC[9] = { - 1, 0, 0, - 0, -1, 0, - 0, 0, -1 }; - -static double LATElementsIEC[9] = { - 0, 0, -1, - 0, -1, 0, - -1, 0, 0}; - -static double HFS2IEC[9] = { - 1, 0, 0, - 0, 0, 1, - 0, -1, 0}; - -static double FFS2IEC[9] = { - -1, 0, 0, - 0, 0, -1, - 0, -1, 0}; - //FUNCTION DECLARATION NECESSARY FOR COPY/PASTE FROM vtkGDCMImageReader // const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds); @@ -97,6 +72,8 @@ itkImageProcessor::itkImageProcessor() // TZero[0]=TZero[1]=TZero[2]=0.; // RZero[0]=RZero[1]=RZero[2]=0.; + + toVTK2D1 = ITKtoVTKFilterType::New(); toVTK2D2 = ITKtoVTKFilterType::New(); toVTKLocalizer1 = ITKtoVTKFilterType::New(); @@ -217,7 +194,37 @@ itkImageProcessor::itkImageProcessor() // }; + // setup the Automatic Registration + metric = MetricType::New(); + mimetric = MIMetricType::New(); + optimizer = OptimizerType::New(); + amoebaoptimizer = AmoebaOptimizerType::New(); + + optimizerObserver = CommandIterationUpdate::New(); + + exhaustiveOptimizer = ExhaustiveOptimizerType::New(); + exhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New(); + + registration = RegistrationType::New(); + + 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() @@ -257,6 +264,17 @@ double itkImageProcessor::GetSCD2(){ m_DRTImage1MetaInfo ->GetSCD(); } +void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof) +{ + m_TransformMetaInfo->SetDegreeOfFreedom(dof); +} + +tDegreeOfFreedomEnum +itkImageProcessor::GetDegreeOfFreedom() +{ + return m_TransformMetaInfo->GetDegreeOfFreedom(); +} + void itkImageProcessor::SetCustom_ImportTransform(double dTx, double dTy, double dTz, @@ -384,8 +402,6 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2) void itkImageProcessor::SetCustom_UpdateMetaInfo(){ - - this->UpdateProjectionGeometryMeta(); } @@ -1254,6 +1270,35 @@ itkImageProcessor::MapTransformToNewOrigin( return m_OutputTransform; } + +// This External User Transform thing need to be checked out +void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo) +{ + //crude hack to get from internal transform in LPS to external transform in IEC. + //result is stored in m_TransformMetaInfo + + //TODO: support for projectionTransformCenter and userTransformCenter which are not the same + //currentlz we support only pFakeIsoLPS situations. + + InternalImageType::DirectionType LPStoIEC_Directions; + LPStoIEC_Directions = m_CTMetaInfo->GetLPS2IECDirections(); + + auto parameters = transform->GetParameters(); + + ImageType3D::PointType translationUser; + translationUser[0] = parameters[3]; + translationUser[1] = parameters[4]; + translationUser[2] = parameters[5]; + ImageType3D::PointType rotationUser; + rotationUser[0] = parameters[0]; + rotationUser[1] = parameters[1]; + rotationUser[2] = parameters[2]; + + m_TransformMetaInfo->SetUserTranslations(LPStoIEC_Directions * translationUser); + m_TransformMetaInfo->SetUserRotations(LPStoIEC_Directions * rotationUser); +} + + itkImageProcessor::TransformType::Pointer itkImageProcessor::CalculateInternalTransform( ImageType3D::PointType m_TranslationOffset, //IEC @@ -1323,6 +1368,362 @@ itkImageProcessor::CalculateInternalTransform( m_OutputTransform; } +void itkImageProcessor::InitializeRegistration(int maximumIteration, 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); + + 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 + + auto region1 = m_PASourceDupli->GetOutput()->GetBufferedRegion(); + auto region2 = m_LATSourceDupli->GetOutput()->GetBufferedRegion(); + + // bool setROI = false; + // if (setROI) { + // auto index1 = region1.GetIndex(); + // auto size1 = region1.GetSize(); + // auto point1 = m_PASourceDupli->GetOutput()->GetOrigin(); + + // index1[0] = size1[0] / 2 - 50; + // index1[1] = size1[1] / 2 - 50; + + // size1[0] = 100; + // size1[1] = 100; + // size1[2] = 1; + + // region1.SetIndex(index1); + // region1.SetSize(size1); + // } + + registration->SetFixedImageRegion1(region1); + registration->SetFixedImageRegion2(region2); + + registration->SetTransformMetaInfo(m_TransformMetaInfo); + + registration->SetFilter1(filter1); + registration->SetFilter2(filter2); + + //CalculateInternalTransform(transform1, m_DRTImage1MetaInfo); + /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ1 *********/ + ImageType3D::PointType ZeroPoint; + ZeroPoint.Fill(0.); + InternalImageType::DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + TransformType::Pointer CurrTransform; + if(m_RTMetaInfo == NULL) + { + + //AAAAAA TODOOOOO CERCA + ImageType3D::PointType pFakeIsoLPS = + m_DRTImage1MetaInfo->GetProjectionOriginLPS(); + //ImageType3D::PointType pFakeIsoLPS = + // m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter1()); + CurrTransform = + CalculateInternalTransform( + ZeroPoint, + ZeroPoint, + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + pFakeIsoLPS, + ZeroPoint, + IECtoLPS_Directions + ); + } else { + + CurrTransform = + CalculateInternalTransform( + ZeroPoint , + m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetImportOffset(), + IECtoLPS_Directions + ); + + } + + + transform1->SetComputeZYX(true); + transform1->SetIdentity(); + transform1->SetTranslation( + CurrTransform->GetTranslation()); + transform1->SetRotation( + CurrTransform->GetAngleX(), + CurrTransform->GetAngleY(), + CurrTransform->GetAngleZ() + ); + transform1->SetCenter( + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); + /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ + + + //CalculateInternalTransform(transform2, m_DRTImage2MetaInfo); + /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/ + if(m_RTMetaInfo == NULL) + { + ImageType3D::PointType pFakeIsoLPS = + m_DRTImage2MetaInfo->GetProjectionOriginLPS(); + + //ImageType3D::PointType pFakeIsoLPS = + // m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter2()); + + + + CurrTransform = + CalculateInternalTransform( + ZeroPoint, + ZeroPoint, + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + pFakeIsoLPS, + ZeroPoint, + IECtoLPS_Directions + ); + } else { + + CurrTransform = + CalculateInternalTransform( + ZeroPoint , + m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetImportOffset(), + IECtoLPS_Directions + ); + + } + transform2->SetComputeZYX(true); + transform2->SetIdentity(); + transform2->SetTranslation( + CurrTransform->GetTranslation()); + transform2->SetRotation( + CurrTransform->GetAngleX(), + CurrTransform->GetAngleY(), + CurrTransform->GetAngleZ() + ); + transform2->SetCenter( + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ); + /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ + + + 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(maximumIteration); + 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::IterationEvent(), optimizerObserver); + + registration->SetOptimizer(optimizer); + } 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) { + 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::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::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 = transform1->GetParameters(); + optimizer->GetCurrentPosition(); + finalParameters = finalParameters - startParameters; + + 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; +} void itkImageProcessor::InitializeProjector() { @@ -2643,6 +3044,85 @@ double* itkImageProcessor::GetTransformParameters(){ return dTransfParam; } +// Doubts about this user thing +void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) +{ + ImageType3D::PointType translations = m_TransformMetaInfo->GetTranslations(); + ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations(); + + dX = translations[0] + userTranslations[0]; + dY = translations[1] + userTranslations[1]; + dZ = translations[2] + userTranslations[2]; +} + +void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ) +{ + ImageType3D::PointType rotations = m_TransformMetaInfo->GetRotations(); + ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations(); + + dRX = rotations[0] + userRotations[0]; + dRY = rotations[1] + userRotations[1]; + dRZ = rotations[2] + userRotations[2]; +} + + +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::SetROI(double, double, double, double) +{ + //TODO: Not implemanted yet. ROI are hard coded and can be enabled using SetFullROI. +} + +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. +} } // end namespace itk diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 5e2de92..19ad4d8 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -30,6 +30,13 @@ 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" @@ -46,13 +53,104 @@ gfattori 08.11.2021 namespace itk { -/* reference string required for comparison with tag values */ -static const char *ImageOrientationStrings[] = { - "NotDefined", - "HFS", - "FFS" + +class CommandIterationUpdate : public itk::Command { + // TODO: Move to own files. + +public: + using Self = CommandIterationUpdate; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer; + 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 { @@ -97,6 +195,10 @@ public: double GetSCD1(); double GetSCD2(); + /** 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); @@ -117,6 +219,19 @@ public: void SetInitialTranslations(double, double, double); void SetInitialRotations(double, double, double); + /** Get transform parameters for 3D Volume */ + void GetFinalTranslations(double&, double&, double&); + void GetFinalRotations(double&, double&, double&); + + /** Set user transform parameters for 3D Volume in LPS*/ + void SetUserTranslationsLPS(double, double, double); + void SetUserRotationsLPS(double, double, double); + + /** Initialize the registration pipeline*/ + void InitializeRegistration(int, double, double, eDegreeOfFreedomType); + /** Start the registration process*/ + int StartRegistration(std::string extraInfo); + /** Get transform parameters for 3D Volume */ double* GetTransformParameters(); @@ -129,8 +244,6 @@ public: /** Get Projection origin in LPS coordinates*/ const std::vector GetRTImportOffset(); - - /** Get the Localizer intensity window and * center for rendering */ double GetLocalizerDisplayWindowLevel(int ); @@ -153,6 +266,17 @@ public: void WriteProjectionImages(); void Write2DImages(); + /** Get the current cost function value from the optimizer*/ + double GetOptimizerValue(); + + /** Get the cost function value for the current transform*/ + double GetCurrentMetricValue(); + + void SetROI(double, double, double, double); + + void SetOptimizer(std::string); + void SetMetric(std::string); + void SetFullROI(bool); protected: @@ -189,6 +313,16 @@ private: using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; using ResampleFilterType = itk::ResampleImageFilter; + /** Optimizer which tries to find the minimn (Powell Optimizer)*/ + using OptimizerType = itk::PowellOptimizer; + using AmoebaOptimizerType = itk::AmoebaOptimizer; + /** Optimizer which samples the whol space */ + using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer; + /** Metric to calculate similarites between images*/ + using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric; + using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric; + /** The thing which actuall does the image registration*/ + using RegistrationType = itk::TwoProjectionImageRegistrationMethod; /** Image reader types */ using ImageReaderType2D = itk::ImageFileReader; @@ -206,6 +340,8 @@ private: /** ITK to VTK filter */ using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; + void + CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo); TransformType::Pointer CalculateInternalTransform( @@ -239,6 +375,23 @@ private: imageDRT1In, imageDRT2In; + RegistrationType::Pointer + registration; + MetricType::Pointer + metric; + MIMetricType::Pointer + mimetric; + OptimizerType::Pointer + optimizer; + AmoebaOptimizerType::Pointer + amoebaoptimizer; + CommandIterationUpdate::Pointer + optimizerObserver; + ExhaustiveOptimizerType::Pointer + exhaustiveOptimizer; + ExhaustiveCommandIterationUpdate::Pointer + exhaustiveOptimizerObserver; + DuplicatorType::Pointer m_LATSourceDupli, m_PASourceDupli, @@ -335,6 +488,15 @@ private: TransformMetaInformation::Pointer m_TransformMetaInfo; + + double m_OptmizerValue; + + bool m_UseExhaustiveOptmizer; + bool m_UseAmeobaOptimizer; + bool m_UseFullROI; + bool m_UseMutualInformation; + bool m_UseDumptoFile; + };