diff --git a/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt b/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt index adee1d8..4921214 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt +++ b/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}) diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx index 2e18525..b7c03a7 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx @@ -149,6 +149,12 @@ MutualInformationTwoImageToOneImageMetric::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: "<SetFixedImage(fixedImage1); metric1->SetMovingImage(this->m_Filter1->GetOutput()); diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx index 06e9c34..cd14a02 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx @@ -76,7 +76,6 @@ NormalizedCorrelationTwoImageToOneImageMetric::Calcul this->m_NumberOfPixelsCounted = 0; int numberOfPixelsVisited = 0; - //filter->Update(); filter->Update(); MovingImageConstPointer imageDRTIn = filter->GetOutput(); //this->WriteImage(imageDRTIn, "", ""); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index ecaf8ed..63c2628 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -142,37 +142,38 @@ itkImageProcessor::itkImageProcessor() m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D); // setup the Automatic Registration - metric = MetricType::New(); - mimetric = MIMetricType::New(); - optimizer = OptimizerType::New(); + NCCmetric = MetricType::New(); + MImetric = MIMetricType::New(); + PowellOptimizer = PowellOptimizerType::New(); - amoebaoptimizer = AmoebaOptimizerType::New(); + AmoebaOptimizer = AmoebaOptimizerType::New(); optimizerObserver = CommandIterationUpdate::New(); - exhaustiveOptimizer = ExhaustiveOptimizerType::New(); - exhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New(); + ExhaustiveOptimizer = ExhaustiveOptimizerType::New(); + ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New(); registration = RegistrationType::New(); optimizerObserver->SetProcess(registration); - if (m_UseMutualInformation) { - registration->SetMetric(mimetric); - } else { - registration->SetMetric(metric); - } + /* default is MI */ + m_UseMutualInformation = true; //oterwise NCC is used + registration->SetMetric(MImetric); + /* 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: "<GetOutput()->GetDirection()<ComputeGradientOff(); - metric->SetSubtractMean(true); - metric->SetMaxTranslation(maxTranslation); + NCCmetric->ComputeGradientOff(); + MImetric->ComputeGradientOff(); + 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->SetFixedImageRegion2(roiAutoReg2); + MImetric->SetFixedImageRegion1(roiAutoReg1); + 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->SetFixedImageRegion2(roiAutoReg2); + MImetric->SetFixedImageRegion1(roiAutoReg1); + 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 - optimizer->SetMaximumIteration(m_MaxNumberOfIterations); - optimizer->SetMaximumLineIteration(4); // for Powell's method - optimizer->SetStepLength(stepLength); - optimizer->SetStepTolerance(0.01); - optimizer->SetValueTolerance(0.000002); + PowellOptimizer->SetMaximize(false); // for NCC and MI + PowellOptimizer->SetMaximumIteration(m_MaxNumberOfIterations); + PowellOptimizer->SetMaximumLineIteration(4); // for Powell's method + PowellOptimizer->SetStepLength(stepLength); + PowellOptimizer->SetStepTolerance(0.01); + PowellOptimizer->SetValueTolerance(0.000002); if (verbose) { - optimizer->DebugOn(); - optimizer->Print(std::cout); + PowellOptimizer->DebugOn(); + PowellOptimizer->Print(std::cout); } - optimizer->RemoveAllObservers(); - optimizer->AddObserver(itk::AnyEvent(), optimizerObserver); + PowellOptimizer->RemoveAllObservers(); + PowellOptimizer->AddObserver(itk::AnyEvent(), optimizerObserver); - registration->SetOptimizer(optimizer); + registration->SetOptimizer(PowellOptimizer); registration->RemoveAllObservers(); } else { - amoebaoptimizer->SetMinimize(true); - amoebaoptimizer->SetMaximumNumberOfIterations(100); - amoebaoptimizer->SetParametersConvergenceTolerance(0.1); - amoebaoptimizer->SetFunctionConvergenceTolerance(0.000002); + AmoebaOptimizer->SetMinimize(true); + AmoebaOptimizer->SetMaximumNumberOfIterations(100); + AmoebaOptimizer->SetParametersConvergenceTolerance(0.1); + 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->AddObserver(itk::IterationEvent(), optimizerObserver); + AmoebaOptimizer->RemoveAllObservers(); + 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->SetStepLength(stepLength); + ExhaustiveOptimizer->SetNumberOfSteps(steps); + ExhaustiveOptimizer->SetStepLength(stepLength); - exhaustiveOptimizer->RemoveAllObservers(); - exhaustiveOptimizer->AddObserver(itk::IterationEvent(), exhaustiveOptimizerObserver); + ExhaustiveOptimizer->RemoveAllObservers(); + ExhaustiveOptimizer->AddObserver(itk::IterationEvent(), ExhaustiveOptimizerObserver); - exhaustiveOptimizer->SetStepLength(stepLength); - registration->SetOptimizer(exhaustiveOptimizer); + ExhaustiveOptimizer->SetStepLength(stepLength); + 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::digits10 + 2); - fs << " MinimumMetricValue: " << exhaustiveOptimizer->GetMinimumMetricValue() << std::endl; - fs << " MaximumMetricValue: " << exhaustiveOptimizer->GetMaximumMetricValue() << std::endl; + 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; + 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 << " FinalMetricValue: " << PowellOptimizer->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 << " FinalPosition: " << PowellOptimizer->GetCurrentPosition() << std::endl; + fs << " Iterations: " << PowellOptimizer->GetCurrentIteration() << std::endl; + fs << " StopConditionDescription: " << PowellOptimizer->GetStopConditionDescription() << std::endl; } fs.close(); - m_OptmizerValue = optimizer->GetValue(); - auto finalParameters = optimizer->GetCurrentPosition(); + m_OptmizerValue = PowellOptimizer->GetValue(); + + if(!m_UseExhaustiveOptmizer){ + if(!m_UseAmeobaOptimizer){ + //Powell + auto finalParameters = PowellOptimizer->GetCurrentPosition(); + std::cout<<"REG FINAL PARS: "<GetCurrentPosition(); + std::cout<<"REG FINAL PARS: "<GetCurrentPosition(); + std::cout<<"REG FINAL PARS: "<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) { diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index df9f45b..42ac0ee 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -269,7 +269,7 @@ private: using ResampleFilterType = itk::ResampleImageFilter; /** 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; - OptimizerType::Pointer - optimizer; + MImetric; + PowellOptimizerType::Pointer + PowellOptimizer; AmoebaOptimizerType::Pointer - amoebaoptimizer; + AmoebaOptimizer; ExhaustiveOptimizerType::Pointer - exhaustiveOptimizer; + ExhaustiveOptimizer; ExhaustiveCommandIterationUpdate::Pointer - exhaustiveOptimizerObserver; + ExhaustiveOptimizerObserver; DuplicatorType::Pointer m_LATSourceDupli,