Merge branch 'ScoutRTUIDevelOptim' into 'ScoutRTUIDevel'

First integration of Automatic Registration is running (on Approve button)....

See merge request cpt_bioeng/drt!20

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(

reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp

@@ -531,6 +531,8 @@ TransformMetaInformation (){
 
     m_CurrentTransform.Fill(0.);
 
+    m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
+
 }
 
 void

reg23Topograms/itkDTRrecon/DRTMetaInformation.h

@@ -11,15 +11,70 @@
 namespace itk
 {
 
-typedef enum eProjectionOrientationType
+typedef enum eProjectionOrientationType{
-{NA = 0, LAT=2, PA=1}
+    NA = 0,
-tProjOrientationType;
+    LAT=2,
+    PA=1
+} tProjOrientationType;
 
 /** Enum type for Orientation */
-typedef enum eImageOrientationType
+typedef enum eImageOrientationType{
-{NotDefined = 0, HFS = 1, FFS = 2}
+    NotDefined = 0,
-tPatOrientation;
+    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<Self> 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 ();

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 <gdcmReader.h>
 #include <itkMatrix.h>
@@ -24,6 +25,8 @@ gfattori 08.11.2021
 #include "itkOrientImageFilter.h"
 #include <gdcmIPPSorter.h>
 
+#include "itkTimeProbesCollectorBase.h"
+
 #include <string>
 
 
@@ -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<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 = 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
 

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[] = {
+class CommandIterationUpdate : public itk::Command {
-    "NotDefined",
+    // TODO: Move to own files.
-    "HFS",
+
-    "FFS"
+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
 {
 
@@ -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 <double> 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<InternalImageType, double>;
     using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
 
+    /** 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<InternalImageType, InternalImageType>;
+    using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
+    /** The thing which actuall does the image registration*/
+    using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
 
     /** Image reader types */
     using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
@@ -206,6 +340,8 @@ private:
     /** ITK to VTK filter */
     using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
 
+    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;
+
 };