Review of multiple optimizers - metrics

- not usuable util the opti and metric types are better defined using enums
- ROI are passed to both metrics
- Full functionalty only for Powell(NCC or MI), mainly due to GetFinalR23Parameters which get those from Powell....

reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt

@@ -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})

reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx

@@ -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());
 

reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx

@@ -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, "", "");

reg23Topograms/itkDTRrecon/itkImageProcessor.cpp

@@ -142,37 +142,38 @@ itkImageProcessor::itkImageProcessor()
     m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D);
 
     // setup the Automatic Registration
-    metric = MetricType::New();
+    NCCmetric = MetricType::New();
-    mimetric = MIMetricType::New();
+    MImetric = MIMetricType::New();
-    optimizer = OptimizerType::New();
+    PowellOptimizer = PowellOptimizerType::New();
 
-    amoebaoptimizer = AmoebaOptimizerType::New();
+    AmoebaOptimizer = AmoebaOptimizerType::New();
     optimizerObserver = CommandIterationUpdate::New();
 
 
-    exhaustiveOptimizer = ExhaustiveOptimizerType::New();
+    ExhaustiveOptimizer = ExhaustiveOptimizerType::New();
-    exhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();
+    ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();
 
     registration = RegistrationType::New();
     optimizerObserver->SetProcess(registration);
 
-    if (m_UseMutualInformation) {
+    /* default is MI */
-        registration->SetMetric(mimetric);
+    m_UseMutualInformation = true; //oterwise NCC is used
-    } else {
+    registration->SetMetric(MImetric);
-        registration->SetMetric(metric);
+    /* default is Powell */
-    }
+    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.
+    registration->SetOptimizer(PowellOptimizer);
 
-    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()
@@ -1031,7 +1032,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;
@@ -1153,18 +1159,20 @@ void itkImageProcessor::InitializeRegistration(
         itkExceptionMacro(<< "CT data not present");
     }
 
-    metric->ComputeGradientOff();
+    NCCmetric->ComputeGradientOff();
-    metric->SetSubtractMean(true);
+    MImetric->ComputeGradientOff();
-    metric->SetMaxTranslation(maxTranslation);
 
+    NCCmetric->SetSubtractMean(true);
+
+    NCCmetric->SetMaxTranslation(maxTranslation);
+    MImetric->SetMaxTranslation(maxTranslation);
 
-    //    m_TransformMetaInfo->SetDegreeOfFreedom(dof);
     m_r23MetaInfo->SetDegreeOfFreedom(dof);
 
-    mimetric->ComputeGradientOff();
 
     if (verbose) {
-        metric->DebugOn();
+        NCCmetric->DebugOn();
+        MImetric->DebugOn();
     }
 
 
@@ -1184,8 +1192,10 @@ void itkImageProcessor::InitializeRegistration(
 
         registration->SetFixedImageRegion1(roiAutoReg1);
         registration->SetFixedImageRegion2(roiAutoReg2);
-        mimetric->SetFixedImageRegion1(roiAutoReg1);
+        MImetric->SetFixedImageRegion1(roiAutoReg1);
-        mimetric->SetFixedImageRegion2(roiAutoReg2);
+        MImetric->SetFixedImageRegion2(roiAutoReg2);
+        NCCmetric->SetFixedImageRegion1(roiAutoReg1);
+        NCCmetric->SetFixedImageRegion2(roiAutoReg2);
 
         std::cout << "using ROIs for Auto Reg " << std::endl;
 
@@ -1201,8 +1211,10 @@ void itkImageProcessor::InitializeRegistration(
 
         registration->SetFixedImageRegion1(region1);
         registration->SetFixedImageRegion2(region2);
-        mimetric->SetFixedImageRegion1(roiAutoReg1);
+        MImetric->SetFixedImageRegion1(roiAutoReg1);
-        mimetric->SetFixedImageRegion2(roiAutoReg2);
+        MImetric->SetFixedImageRegion2(roiAutoReg2);
+        NCCmetric->SetFixedImageRegion1(roiAutoReg1);
+        NCCmetric->SetFixedImageRegion2(roiAutoReg2);
     }
 
 
@@ -1244,35 +1256,35 @@ void itkImageProcessor::InitializeRegistration(
 
             // 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
+            PowellOptimizer->SetMaximize(false); // for NCC and MI
-            optimizer->SetMaximumIteration(m_MaxNumberOfIterations);
+            PowellOptimizer->SetMaximumIteration(m_MaxNumberOfIterations);
-            optimizer->SetMaximumLineIteration(4); // for Powell's method
+            PowellOptimizer->SetMaximumLineIteration(4); // for Powell's method
-            optimizer->SetStepLength(stepLength);
+            PowellOptimizer->SetStepLength(stepLength);
-            optimizer->SetStepTolerance(0.01);
+            PowellOptimizer->SetStepTolerance(0.01);
-            optimizer->SetValueTolerance(0.000002);
+            PowellOptimizer->SetValueTolerance(0.000002);
             if (verbose) {
-                optimizer->DebugOn();
+                PowellOptimizer->DebugOn();
-                optimizer->Print(std::cout);
+                PowellOptimizer->Print(std::cout);
             }
 
 
-            optimizer->RemoveAllObservers();
+            PowellOptimizer->RemoveAllObservers();
-            optimizer->AddObserver(itk::AnyEvent(), optimizerObserver);
+            PowellOptimizer->AddObserver(itk::AnyEvent(), optimizerObserver);
 
-            registration->SetOptimizer(optimizer);
+            registration->SetOptimizer(PowellOptimizer);
             registration->RemoveAllObservers();
 
         } else {
 
-            amoebaoptimizer->SetMinimize(true);
+            AmoebaOptimizer->SetMinimize(true);
-            amoebaoptimizer->SetMaximumNumberOfIterations(100);
+            AmoebaOptimizer->SetMaximumNumberOfIterations(100);
-            amoebaoptimizer->SetParametersConvergenceTolerance(0.1);
+            AmoebaOptimizer->SetParametersConvergenceTolerance(0.1);
-            amoebaoptimizer->SetFunctionConvergenceTolerance(0.000002);
+            AmoebaOptimizer->SetFunctionConvergenceTolerance(0.000002);
 
-            amoebaoptimizer->SetAutomaticInitialSimplex(false);
+            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);
+            AmoebaOptimizerType::ParametersType simplexDelta(3);
             int numberOfDOF = GetNumberOfDOF(dof);
             if (numberOfDOF == 0) {
                 itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
@@ -1280,12 +1292,12 @@ void itkImageProcessor::InitializeRegistration(
             for (int i = 0; i < numberOfDOF; i++) {
                 simplexDelta[i] = 2.0;
             }
-            amoebaoptimizer->SetInitialSimplexDelta(simplexDelta);
+            AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta);
 
-            amoebaoptimizer->RemoveAllObservers();
+            AmoebaOptimizer->RemoveAllObservers();
-            amoebaoptimizer->AddObserver(itk::IterationEvent(), optimizerObserver);
+            AmoebaOptimizer->AddObserver(itk::IterationEvent(), optimizerObserver);
 
-            registration->SetOptimizer(amoebaoptimizer);
+            registration->SetOptimizer(AmoebaOptimizer);
         }
 
     } else {
@@ -1300,14 +1312,14 @@ void itkImageProcessor::InitializeRegistration(
         steps[0] = numberOfSteps;
         steps[1] = numberOfSteps;
         steps[2] = numberOfSteps;
-        exhaustiveOptimizer->SetNumberOfSteps(steps);
+        ExhaustiveOptimizer->SetNumberOfSteps(steps);
-        exhaustiveOptimizer->SetStepLength(stepLength);
+        ExhaustiveOptimizer->SetStepLength(stepLength);
 
-        exhaustiveOptimizer->RemoveAllObservers();
+        ExhaustiveOptimizer->RemoveAllObservers();
-        exhaustiveOptimizer->AddObserver(itk::IterationEvent(), exhaustiveOptimizerObserver);
+        ExhaustiveOptimizer->AddObserver(itk::IterationEvent(), ExhaustiveOptimizerObserver);
 
-        exhaustiveOptimizer->SetStepLength(stepLength);
+        ExhaustiveOptimizer->SetStepLength(stepLength);
-        registration->SetOptimizer(exhaustiveOptimizer);
+        registration->SetOptimizer(ExhaustiveOptimizer);
     }
 
     std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
@@ -1338,7 +1350,7 @@ int itkImageProcessor::StartRegistration(std::string extraInfo)
         fs << "Value\tX\tY\tZ " << std::endl;
 
         if (m_UseExhaustiveOptmizer) {
-            exhaustiveOptimizerObserver->set_stream(fs);
+            ExhaustiveOptimizerObserver->set_stream(fs);
         }
     }
 
@@ -1369,30 +1381,47 @@ int itkImageProcessor::StartRegistration(std::string extraInfo)
 
     if (m_UseExhaustiveOptmizer) {
         fs.precision(std::numeric_limits<double>::digits10 + 2);
-        fs << "  MinimumMetricValue: " << exhaustiveOptimizer->GetMinimumMetricValue() << std::endl;
+        fs << "  MinimumMetricValue: " << ExhaustiveOptimizer->GetMinimumMetricValue() << std::endl;
-        fs << "  MaximumMetricValue: " << exhaustiveOptimizer->GetMaximumMetricValue() << std::endl;
+        fs << "  MaximumMetricValue: " << ExhaustiveOptimizer->GetMaximumMetricValue() << std::endl;
         fs.precision(oldprecision);
-        fs << "  MinimumMetricValuePosition: " << exhaustiveOptimizer->GetMinimumMetricValuePosition() << std::endl;
+        fs << "  MinimumMetricValuePosition: " << ExhaustiveOptimizer->GetMinimumMetricValuePosition() << std::endl;
-        fs << "  MaximumMetricValuePosition: " << exhaustiveOptimizer->GetMaximumMetricValuePosition() << std::endl;
+        fs << "  MaximumMetricValuePosition: " << ExhaustiveOptimizer->GetMaximumMetricValuePosition() << std::endl;
-        fs << "  StopConditionDescription: " << exhaustiveOptimizer->GetStopConditionDescription() << std::endl;
+        fs << "  StopConditionDescription: " << ExhaustiveOptimizer->GetStopConditionDescription() << std::endl;
-        fs << "  MaximumMetricValuePosition: " << exhaustiveOptimizer->GetMaximumMetricValuePosition() << 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 << "  FinalMetricValue: " << PowellOptimizer->GetValue() << std::endl; //same as GetCurrentCost
         fs.precision(oldprecision);
-        fs << "  FinalPosition: " << optimizer->GetCurrentPosition() << std::endl;
+        fs << "  FinalPosition: " << PowellOptimizer->GetCurrentPosition() << std::endl;
-        fs << "  Iterations: " << optimizer->GetCurrentIteration() << std::endl;
+        fs << "  Iterations: " << PowellOptimizer->GetCurrentIteration() << std::endl;
-        fs << "  StopConditionDescription: " << optimizer->GetStopConditionDescription() << std::endl;
+        fs << "  StopConditionDescription: " << PowellOptimizer->GetStopConditionDescription() << std::endl;
     }
 
     fs.close();
 
 
 
-    m_OptmizerValue = optimizer->GetValue();
+    m_OptmizerValue = PowellOptimizer->GetValue();
-    auto finalParameters = optimizer->GetCurrentPosition();
+
+    if(!m_UseExhaustiveOptmizer){
+        if(!m_UseAmeobaOptimizer){
+            //Powell
+            auto finalParameters = PowellOptimizer->GetCurrentPosition();
+            std::cout<<"REG FINAL PARS: "<<finalParameters<<std::endl;
+
+        }else{
+            //Amoeba
+            auto finalParameters = AmoebaOptimizer->GetCurrentPosition();
+            std::cout<<"REG FINAL PARS: "<<finalParameters<<std::endl;
+
+        }
+    }else{
+        //use exhaustive
+        auto finalParameters = ExhaustiveOptimizer->GetCurrentPosition();
+        std::cout<<"REG FINAL PARS: "<<finalParameters<<std::endl;
+
+    }
 
-    std::cout<<"REG FINAL PARS: "<<finalParameters<<std::endl;
 
 
     //    auto finalParameters = transform1->GetParameters();
@@ -1452,32 +1481,32 @@ itkImageProcessor::GetFinalR23Parameters(){
     Optimizer::ParametersType pPars(6);
     pPars.fill(0.);
 
-    if(optimizer == nullptr){
+    if(PowellOptimizer == nullptr){
         return pPars;
     }
 
     switch (m_r23MetaInfo->GetDegreeOfFreedom()) {
     case THREE_DOF:
-        pPars[3] = optimizer->GetCurrentPosition()[0]
+        pPars[3] = PowellOptimizer->GetCurrentPosition()[0]
                 - m_TransformMetaInfo->GetHiddenTranslations()[0];
-        pPars[4] = optimizer->GetCurrentPosition()[1]
+        pPars[4] = PowellOptimizer->GetCurrentPosition()[1]
                 - m_TransformMetaInfo->GetHiddenTranslations()[1];
-        pPars[5] = optimizer->GetCurrentPosition()[2]
+        pPars[5] = PowellOptimizer->GetCurrentPosition()[2]
                 - m_TransformMetaInfo->GetHiddenTranslations()[2];
         break;
 
     case SIX_DOF:
-        pPars[0] = optimizer->GetCurrentPosition()[0]
+        pPars[0] = PowellOptimizer->GetCurrentPosition()[0]
                 - m_TransformMetaInfo->GetHiddenRotations()[0];
-        pPars[1] = optimizer->GetCurrentPosition()[1]
+        pPars[1] = PowellOptimizer->GetCurrentPosition()[1]
                 - m_TransformMetaInfo->GetHiddenRotations()[1];
-        pPars[2] = optimizer->GetCurrentPosition()[2]
+        pPars[2] = PowellOptimizer->GetCurrentPosition()[2]
                 - m_TransformMetaInfo->GetHiddenRotations()[2];
-        pPars[3] = optimizer->GetCurrentPosition()[3]
+        pPars[3] = PowellOptimizer->GetCurrentPosition()[3]
                 - m_TransformMetaInfo->GetHiddenTranslations()[0];
-        pPars[4] = optimizer->GetCurrentPosition()[4]
+        pPars[4] = PowellOptimizer->GetCurrentPosition()[4]
                 - m_TransformMetaInfo->GetHiddenTranslations()[1];
-        pPars[5] = optimizer->GetCurrentPosition()[5]
+        pPars[5] = PowellOptimizer->GetCurrentPosition()[5]
                 - m_TransformMetaInfo->GetHiddenTranslations()[2];
         break;
 
@@ -2776,15 +2805,15 @@ double itkImageProcessor::GetCurrentMetricValue()
     //Note: this fails if the metric has not been propperly initilzed.
 
     if (!m_UseMutualInformation) {
-        if (metric == NULL) {
+        if (NCCmetric == NULL) {
             return nan("");
         }
-        return metric->GetValue();
+        return NCCmetric->GetValue();
     } else {
-        if (mimetric == NULL) {
+        if (MImetric == NULL) {
             return nan("");
         }
-        return mimetric->GetValue();
+        return MImetric->GetValue();
     }
 }
 
@@ -2794,15 +2823,20 @@ void itkImageProcessor::SetOptimizer(std::string optimizer)
     if (optimizer.compare("Powell") == 0) {
         m_UseAmeobaOptimizer = false;
         m_UseExhaustiveOptmizer = false;
+        registration->SetOptimizer(PowellOptimizer);
     } else if (optimizer.compare("Amoeba") == 0) {
         m_UseExhaustiveOptmizer = false;
         m_UseAmeobaOptimizer = true;
+        registration->SetOptimizer(AmoebaOptimizer);
     } else if (optimizer.compare("Exhaustive") == 0) {
         m_UseExhaustiveOptmizer = true;
         m_UseAmeobaOptimizer = false;
+        registration->SetOptimizer(ExhaustiveOptimizer);
     } else {
         itkExceptionMacro(<< "Unkown optimzer string : " << optimizer);
     }
+
+
 }
 void itkImageProcessor::SetMetric(std::string metric)
 {
@@ -2811,8 +2845,15 @@ void itkImageProcessor::SetMetric(std::string metric)
     } else if (metric.compare("MI") == 0) {
         m_UseMutualInformation = true;
     } else {
-        itkExceptionMacro(<< "Unkown metric string : " << optimizer);
+        itkExceptionMacro(<< "Unkown metric string : " << metric);
     }
+
+    if (m_UseMutualInformation) {
+         registration->SetMetric(this->MImetric);
+     } else {
+         registration->SetMetric(this->NCCmetric);
+     }
+
 }
 void itkImageProcessor::SetFullROI(bool fullROI)
 {

reg23Topograms/itkDTRrecon/itkImageProcessor.h

@@ -269,7 +269,7 @@ private:
     using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
 
     /** Optimizer which tries to find the minimun (Powell Optimizer)*/
-    using OptimizerType = itk::PowellOptimizer;
+    using PowellOptimizerType = itk::PowellOptimizer;
     /** Optimizer which tries to find the minimn (Powell Optimizer)*/
     using AmoebaOptimizerType = itk::AmoebaOptimizer;
     /** Optimizer which samples the whol space */
@@ -321,17 +321,17 @@ private:
     RegistrationType::Pointer
     registration;
     MetricType::Pointer
-    metric;
+    NCCmetric;
     MIMetricType::Pointer
-    mimetric;
+    MImetric;
-    OptimizerType::Pointer
+    PowellOptimizerType::Pointer
-    optimizer;
+    PowellOptimizer;
     AmoebaOptimizerType::Pointer
-    amoebaoptimizer;
+    AmoebaOptimizer;
     ExhaustiveOptimizerType::Pointer
-    exhaustiveOptimizer;
+    ExhaustiveOptimizer;
     ExhaustiveCommandIterationUpdate::Pointer
-    exhaustiveOptimizerObserver;
+    ExhaustiveOptimizerObserver;
 
     DuplicatorType::Pointer
     m_LATSourceDupli,