diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index e97344c..dff2dd7 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -19,7 +19,7 @@ SET(SRCS itkImageProcessorHelpers.cpp vtkContourTopogramProjectionFilter.cxx DRTMetaInformation.cpp - + itkReg23.cpp ) SET(HDR @@ -30,6 +30,7 @@ SET(HDR itkgSiddonJacobsRayCastInterpolateImageFunction.hxx vtkContourTopogramProjectionFilter.h DRTMetaInformation.h + itkReg23.h ) ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR}) diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index ff6ae1c..1681aad 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -547,7 +547,9 @@ R23MetaInformation:: R23MetaInformation (){ m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; - + m_OptimizerType = tOptimizerTypeEnum::POWELL; + m_MetricType = tMetricTypeEnum::NCC; + m_MaxIterations = 6; } diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 0e6e9cd..2c59bfc 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -35,6 +35,19 @@ typedef enum eDegreeOfFreedomType { OTHER } tDegreeOfFreedomEnum; + +typedef enum eOptimizerType{ + POWELL, + AMOEBA, + EXHAUSTIVE +} tOptimizerTypeEnum; + +typedef enum eMetricType{ + NCC, + MI +} tMetricTypeEnum; + + inline int GetNumberOfDOF(eDegreeOfFreedomType dof) { switch (dof) { @@ -597,16 +610,35 @@ public: itkTypeMacro(R23MetaInformation, itk::Object); /** object information streaming **/ - void PrintSelf(std::ostream& os, itk::Indent indent) const; + void PrintSelf(std::ostream& os, itk::Indent indent) const override; itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + itkSetEnumMacro(OptimizerType, tOptimizerTypeEnum); + itkGetEnumMacro(OptimizerType, tOptimizerTypeEnum); + + + itkSetEnumMacro(MetricType, tMetricTypeEnum); + itkGetEnumMacro(MetricType, tMetricTypeEnum); + + itkSetEnumMacro(MaxIterations, int); + itkGetEnumMacro(MaxIterations, int); + protected: tDegreeOfFreedomEnum - m_DegreeOfFreedom; + m_DegreeOfFreedom; + + tOptimizerTypeEnum + m_OptimizerType; + + tMetricTypeEnum + m_MetricType; + + int + m_MaxIterations; /** Default Constructor **/ R23MetaInformation (); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 63c2628..030e36d 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -40,7 +40,7 @@ itkImageProcessor::itkImageProcessor() { iNWUnits = 1; - // std::cout<<"itkImageProcessor::itkImageProcessor contructor"<SetIdentity(); transform2 = TransformType::New(); @@ -53,11 +53,6 @@ itkImageProcessor::itkImageProcessor() interpolator1 = InterpolatorType::New(); interpolator2 = InterpolatorType::New(); - // TZero[0]=TZero[1]=TZero[2]=0.; - // RZero[0]=RZero[1]=RZero[2]=0.; - - - toVTK2D1 = ITKtoVTKFilterType::New(); toVTK2D2 = ITKtoVTKFilterType::New(); toVTKLocalizer1 = ITKtoVTKFilterType::New(); @@ -141,37 +136,14 @@ itkImageProcessor::itkImageProcessor() m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D); m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D); - // setup the Automatic Registration - NCCmetric = MetricType::New(); - MImetric = MIMetricType::New(); - PowellOptimizer = PowellOptimizerType::New(); - AmoebaOptimizer = AmoebaOptimizerType::New(); optimizerObserver = CommandIterationUpdate::New(); - - - ExhaustiveOptimizer = ExhaustiveOptimizerType::New(); ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New(); - registration = RegistrationType::New(); - optimizerObserver->SetProcess(registration); - - /* 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->SetTransform1(transform1); - registration->SetTransform2(transform2); - registration->SetInterpolator1(interpolator1); - registration->SetInterpolator2(interpolator2); - - m_UseFullROI = true; // if the full image ROI shall be used - m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file. + // TOOOOODOOOOO optimizerObserver->SetProcess(registration); + m_R23 = itkReg23::New(); + m_R23->SetOptimizerObserver(optimizerObserver); } @@ -262,11 +234,6 @@ double itkImageProcessor::GetPanelOffset2(){ -void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof) -{ - m_r23MetaInfo->SetDegreeOfFreedom(dof); -} - tDegreeOfFreedomEnum itkImageProcessor::GetDegreeOfFreedom() { @@ -1032,12 +999,12 @@ int itkImageProcessor::load2D(const char * pcFName){ m_Duplicator->SetInputImage(caster2DInput->GetOutput() ); m_Duplicator->Update(); - InternalImageType::IndexType pIndex; - pIndex[0]=200; - pIndex[1]=300; - InternalImageType::PixelType pPixel = - m_Duplicator->GetOutput()->GetPixel(pIndex); - std::cout<<"processro PIXEL: "<GetOutput()->GetPixel(pIndex); + //std::cout<<"processro PIXEL: "<GetOutput()->GetDirection()<ComputeGradientOff(); - MImetric->ComputeGradientOff(); - - NCCmetric->SetSubtractMean(true); - - NCCmetric->SetMaxTranslation(maxTranslation); - MImetric->SetMaxTranslation(maxTranslation); - - m_r23MetaInfo->SetDegreeOfFreedom(dof); - - - if (verbose) { - NCCmetric->DebugOn(); - MImetric->DebugOn(); - } - - - registration->SetFixedImage1(m_PASourceDupli->GetOutput()); - registration->SetFixedImage2(m_LATSourceDupli->GetOutput()); - registration->SetMovingImage(m_VolumeSourceDupli->GetOutput()); - - //TODO: Here we could set the ROI for image registartion - - - if (m_UseFullROI == false) { - - auto region1 = m_PASourceDupli->GetOutput()->GetBufferedRegion(); - auto region2 = m_LATSourceDupli->GetOutput()->GetBufferedRegion(); - std::cout << "region1 of opti - NO ROI " << region1 << std::endl; - std::cout << "region2 of opti - NO ROI " << region2 << std::endl; - - registration->SetFixedImageRegion1(roiAutoReg1); - registration->SetFixedImageRegion2(roiAutoReg2); - MImetric->SetFixedImageRegion1(roiAutoReg1); - MImetric->SetFixedImageRegion2(roiAutoReg2); - NCCmetric->SetFixedImageRegion1(roiAutoReg1); - NCCmetric->SetFixedImageRegion2(roiAutoReg2); - - std::cout << "using ROIs for Auto Reg " << std::endl; - - std::cout << "region1 of opti " << roiAutoReg1 << std::endl; - std::cout << "region2 of opti " << roiAutoReg2 << std::endl; - - }else{ - - auto region1 = m_PASourceDupli->GetOutput()->GetBufferedRegion(); - auto region2 = m_LATSourceDupli->GetOutput()->GetBufferedRegion(); - std::cout << "region1 of opti - NO ROI " << region1 << std::endl; - std::cout << "region2 of opti - NO ROI " << region2 << std::endl; - - registration->SetFixedImageRegion1(region1); - registration->SetFixedImageRegion2(region2); - MImetric->SetFixedImageRegion1(roiAutoReg1); - MImetric->SetFixedImageRegion2(roiAutoReg2); - NCCmetric->SetFixedImageRegion1(roiAutoReg1); - NCCmetric->SetFixedImageRegion2(roiAutoReg2); - } - - - // registration->SetTransformMetaInfo(m_TransformMetaInfo); - registration->SetTransformMetaInfo(m_r23MetaInfo); - registration->SetinternalMeta1(m_InternalTransf1); - registration->SetinternalMeta2(m_InternalTransf2); - - registration->SetFilter1(filter1); - registration->SetFilter2(filter2); - - - IsocTransf = TransformType::New(); - IsocTransf->SetComputeZYX(true); - IsocTransf->SetIdentity(); - - IsocTransf->SetRotation( - m_TransformMetaInfo->GetR()[0], - m_TransformMetaInfo->GetR()[1], - m_TransformMetaInfo->GetR()[2] - ); - - TransformType::OutputVectorType TranslV; - TranslV[0] = m_TransformMetaInfo->GetT()[0]; - TranslV[1] = m_TransformMetaInfo->GetT()[1]; - TranslV[2] = m_TransformMetaInfo->GetT()[2]; - - IsocTransf->SetTranslation(TranslV); - - registration->SetIsocIECTransform(IsocTransf); - - if (verbose) { - registration->DebugOn(); - registration->Print(std::cout); - } - - if (!m_UseExhaustiveOptmizer) { - if (!m_UseAmeobaOptimizer) { - - // The metric will return 0 or a "large" negative number in case of high correspondence of the images - // Thus, we need to minimze - 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) { - PowellOptimizer->DebugOn(); - PowellOptimizer->Print(std::cout); - } - - - PowellOptimizer->RemoveAllObservers(); - PowellOptimizer->AddObserver(itk::AnyEvent(), optimizerObserver); - - registration->SetOptimizer(PowellOptimizer); - registration->RemoveAllObservers(); - - } else { - - AmoebaOptimizer->SetMinimize(true); - AmoebaOptimizer->SetMaximumNumberOfIterations(100); - AmoebaOptimizer->SetParametersConvergenceTolerance(0.1); - AmoebaOptimizer->SetFunctionConvergenceTolerance(0.000002); - - AmoebaOptimizer->SetAutomaticInitialSimplex(false); - //Initial size of the simplex (otherwise it is a very small number and optimizer stops immeditaly) - // 2 mm / 2 degrees seems to be a good value which performs well. - AmoebaOptimizerType::ParametersType simplexDelta(3); - int numberOfDOF = GetNumberOfDOF(dof); - if (numberOfDOF == 0) { - itkExceptionMacro(<< "Unkown or unsupported degree of freedom"); - } - for (int i = 0; i < numberOfDOF; i++) { - simplexDelta[i] = 2.0; - } - AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta); - - AmoebaOptimizer->RemoveAllObservers(); - AmoebaOptimizer->AddObserver(itk::IterationEvent(), optimizerObserver); - - registration->SetOptimizer(AmoebaOptimizer); - } - - } else { - //only support for 3DOF - const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651. - const double stepLength = 0.1; - - //auto parameters = transform1->GetParameters(); - //transform1->SetParameters(parameters); - - ExhaustiveOptimizerType::StepsType steps(3); - steps[0] = numberOfSteps; - steps[1] = numberOfSteps; - steps[2] = numberOfSteps; - ExhaustiveOptimizer->SetNumberOfSteps(steps); - ExhaustiveOptimizer->SetStepLength(stepLength); - - ExhaustiveOptimizer->RemoveAllObservers(); - ExhaustiveOptimizer->AddObserver(itk::IterationEvent(), ExhaustiveOptimizerObserver); - - ExhaustiveOptimizer->SetStepLength(stepLength); - registration->SetOptimizer(ExhaustiveOptimizer); - } - - std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl; -} - -int itkImageProcessor::StartRegistration(std::string extraInfo) -{ - - std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; - // TODO: Check if the registartion pipeline has been initialized - using ParametersType = RegistrationType::ParametersType; - - //auto startParameters = transform1->GetParameters(); - - time_t t = time(0); // get time now - struct tm* now = localtime(&t); - - bool useDumpToFile = false; - - char buffer[255]; - strftime(buffer, 255, "test-%F-%H%M%S.txt", now); - - std::fstream fs; - - if (useDumpToFile || m_UseExhaustiveOptmizer) { - fs.open(buffer, std::fstream::out); - fs << extraInfo; - fs << "Value\tX\tY\tZ " << std::endl; - - if (m_UseExhaustiveOptmizer) { - ExhaustiveOptimizerObserver->set_stream(fs); - } - } - - // Start the registration - // ~~~~~~~~~~~~~~~~~~~~~~ - - // Create a timer to record calculation time. - itk::TimeProbesCollectorBase timer; - - if (verbose) { - std::cout << "Starting the registration now..." << std::endl; - } - - try { - // timer.Start("Registration"); - // Start the registration. - registration->StartRegistration(); - // timer.Stop("Registration"); - } catch (itk::ExceptionObject& err) { - - registration->ResetPipeline(); - std::cout << "ExceptionObject caught !" << std::endl; - std::cout << err << std::endl; - return -1; - } - - auto oldprecision = fs.precision(); - - if (m_UseExhaustiveOptmizer) { - fs.precision(std::numeric_limits::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: " << PowellOptimizer->GetValue() << std::endl; //same as GetCurrentCost - fs.precision(oldprecision); - fs << " FinalPosition: " << PowellOptimizer->GetCurrentPosition() << std::endl; - fs << " Iterations: " << PowellOptimizer->GetCurrentIteration() << std::endl; - fs << " StopConditionDescription: " << PowellOptimizer->GetStopConditionDescription() << std::endl; - } - - fs.close(); - - - - 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(); - // std::cout<<"-->startParameters "<finalParameters1 "<GetParameters()<finalParameters2 "<GetParameters()<GetCurrentPosition(); - // finalParameters = finalParameters - startParameters; - - // std::cout<<"startParameters "<GetCenter()<GetTranslation()<GetCenter()<GetTranslation()<SetParameters(finalParameters); - - // CalculateExternalUserTransform(offsetTransform, m_DRTImage1MetaInfo); - - // const double translationAlongX = finalParameters[3]; - // const double translationAlongY = finalParameters[4]; - // const double translationAlongZ = finalParameters[5]; - // const double rotationAlongX = finalParameters[0]; - // const double rotationAlongY = finalParameters[1]; - // const double rotationAlongZ = finalParameters[2]; - - // const int numberOfIterations = optimizer->GetCurrentIteration(); - - // std::cout << "Result = " << std::endl; - // std::cout << " Rotation Along X = " << rotationAlongX / dtr << " deg" << std::endl; - // std::cout << " Rotation Along Y = " << rotationAlongY / dtr << " deg" << std::endl; - // std::cout << " Rotation Along Z = " << rotationAlongZ / dtr << " deg" << std::endl; - // std::cout << " Translation X = " << translationAlongX << " mm" << std::endl; - // std::cout << " Translation Y = " << translationAlongY << " mm" << std::endl; - // std::cout << " Translation Z = " << translationAlongZ << " mm" << std::endl; - // std::cout << " Number Of Iterations = " << numberOfIterations << std::endl; - // std::cout << " Metric value = " << m_OptmizerValue << std::endl; - - std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl; - - return 0; -} Optimizer::ParametersType itkImageProcessor::GetFinalR23Parameters(){ @@ -1481,32 +1098,31 @@ itkImageProcessor::GetFinalR23Parameters(){ Optimizer::ParametersType pPars(6); pPars.fill(0.); - if(PowellOptimizer == nullptr){ - return pPars; - } + itk::Optimizer::ParametersType currentPosition = + m_R23->GetCurrentPosition(); switch (m_r23MetaInfo->GetDegreeOfFreedom()) { case THREE_DOF: - pPars[3] = PowellOptimizer->GetCurrentPosition()[0] + pPars[3] = currentPosition[0] - m_TransformMetaInfo->GetHiddenTranslations()[0]; - pPars[4] = PowellOptimizer->GetCurrentPosition()[1] + pPars[4] = currentPosition[1] - m_TransformMetaInfo->GetHiddenTranslations()[1]; - pPars[5] = PowellOptimizer->GetCurrentPosition()[2] + pPars[5] = currentPosition[2] - m_TransformMetaInfo->GetHiddenTranslations()[2]; break; case SIX_DOF: - pPars[0] = PowellOptimizer->GetCurrentPosition()[0] + pPars[0] = currentPosition[0] - m_TransformMetaInfo->GetHiddenRotations()[0]; - pPars[1] = PowellOptimizer->GetCurrentPosition()[1] + pPars[1] = currentPosition[1] - m_TransformMetaInfo->GetHiddenRotations()[1]; - pPars[2] = PowellOptimizer->GetCurrentPosition()[2] + pPars[2] = currentPosition[2] - m_TransformMetaInfo->GetHiddenRotations()[2]; - pPars[3] = PowellOptimizer->GetCurrentPosition()[3] + pPars[3] = currentPosition[3] - m_TransformMetaInfo->GetHiddenTranslations()[0]; - pPars[4] = PowellOptimizer->GetCurrentPosition()[4] + pPars[4] = currentPosition[4] - m_TransformMetaInfo->GetHiddenTranslations()[1]; - pPars[5] = PowellOptimizer->GetCurrentPosition()[5] + pPars[5] = currentPosition[5] - m_TransformMetaInfo->GetHiddenTranslations()[2]; break; @@ -1521,6 +1137,41 @@ itkImageProcessor::GetFinalR23Parameters(){ } + +void itkImageProcessor::SetOptimizer(std::string optimizer) +{ + + if (optimizer.compare("Powell") == 0) { + m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::POWELL); + } else if (optimizer.compare("Amoeba") == 0) { + m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::AMOEBA); + } else if (optimizer.compare("Exhaustive") == 0) { + m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::EXHAUSTIVE); + } else { + itkExceptionMacro(<< "Unkown optimzer string : " << optimizer); + } + + +} +void itkImageProcessor::SetMetric(std::string metric) +{ + + if (metric.compare("NCC") == 0) { + m_r23MetaInfo->SetMetricType(tMetricTypeEnum::NCC); + } else if (metric.compare("MI") == 0) { + m_r23MetaInfo->SetMetricType(tMetricTypeEnum::MI); + } else { + itkExceptionMacro(<< "Unkown metric string : " << metric); + } + + +} + +void itkImageProcessor::SetMaxNumberOfIterations(int iNum) +{ + m_r23MetaInfo->SetMaxIterations(iNum); +} + void itkImageProcessor::InitializeProjector() { @@ -2795,77 +2446,6 @@ itkImageProcessor::GetCompleteIsocentricTransform } -double itkImageProcessor::GetOptimizerValue() -{ - return m_OptmizerValue; -} - -double itkImageProcessor::GetCurrentMetricValue() -{ - //Note: this fails if the metric has not been propperly initilzed. - - if (!m_UseMutualInformation) { - if (NCCmetric == NULL) { - return nan(""); - } - return NCCmetric->GetValue(); - } else { - if (MImetric == NULL) { - return nan(""); - } - return MImetric->GetValue(); - } -} - - -void itkImageProcessor::SetOptimizer(std::string optimizer) -{ - if (optimizer.compare("Powell") == 0) { - m_UseAmeobaOptimizer = false; - m_UseExhaustiveOptmizer = false; - 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) -{ - if (metric.compare("NCC") == 0) { - m_UseMutualInformation = false; - } else if (metric.compare("MI") == 0) { - m_UseMutualInformation = true; - } else { - itkExceptionMacro(<< "Unkown metric string : " << metric); - } - - if (m_UseMutualInformation) { - registration->SetMetric(this->MImetric); - } else { - registration->SetMetric(this->NCCmetric); - } - -} -void itkImageProcessor::SetFullROI(bool fullROI) -{ - m_UseFullROI = fullROI; - // TODO: Remove this function when ROI functinalitz has been implemented. -} - -void itkImageProcessor::SetMaxNumberOfIterations(int iNum) -{ - m_MaxNumberOfIterations = iNum; - // TODO: Remove this function when ROI functinalitz has been implemented. -} void itkImageProcessor::SetRegionFixed1( int iIdx0, int iIdx1, int iSz0, int iSz1){ @@ -2888,6 +2468,8 @@ void itkImageProcessor::SetRegionFixed1( std::cout << "itkImageProcessor " << std::endl; std::cout << roiAutoReg1 << std::endl; + this->m_R23->SetroiAutoReg1(roiAutoReg1); + } void itkImageProcessor::SetRegionFixed2( @@ -2911,7 +2493,43 @@ void itkImageProcessor::SetRegionFixed2( std::cout << "itkImageProcessor " << std::endl; std::cout << roiAutoReg2 << std::endl; + this->m_R23->SetroiAutoReg2(roiAutoReg2); + } +void itkImageProcessor::ResetROIRegions(){ + + auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion(); + auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion(); + + this->m_R23->SetroiAutoReg1(region1temp); + this->m_R23->SetroiAutoReg2(region2temp); + +} + +void itkImageProcessor::InizializeRegistration(double stepLength, + double maxTranslation, + eDegreeOfFreedomType dof){ + + m_r23MetaInfo->SetDegreeOfFreedom(dof); + m_R23->SetPA(m_PASourceDupli->GetOutput()); + m_R23->SetLAT(m_LATSourceDupli->GetOutput()); + m_R23->SetVolume(m_VolumeSourceDupli->GetOutput()); + m_R23->Setr23Meta(m_r23MetaInfo); + m_R23->SetInternalTransf1(m_InternalTransf1); + m_R23->SetInternalTransf2(m_InternalTransf2); + m_R23->Setfilter1(filter1); + m_R23->Setfilter2(filter2); + m_R23->Setinterpolator1(interpolator1); + m_R23->Setinterpolator2(interpolator2); + m_R23->SetTransformMetaInfo(m_TransformMetaInfo); + + m_R23->InitializeRegistration( + stepLength,maxTranslation,dof); + +} + + + } // end namespace itk diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 42ac0ee..b4f2400 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -30,13 +30,6 @@ gfattori 08.11.2021 #include "itkObjectFactory.h" #include "itkSmartPointer.h" -#include "itkExhaustiveOptimizer.h" -#include "itkMutualInformationTwoImageToOneImageMetric.h" -#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h" -#include "itkPowellOptimizer.h" -#include "itkTwoProjectionImageRegistrationMethod.h" -#include "itkAmoebaOptimizer.h" - #include "itkImageToVTKImageFilter.h" #include "itkImageFileReader.h" @@ -53,6 +46,8 @@ gfattori 08.11.2021 #include "itkImageProcessorHelpers.h" +#include "itkReg23.h" + namespace itk { @@ -78,6 +73,33 @@ public: CommandIterationUpdate::Pointer optimizerObserver; + ExhaustiveCommandIterationUpdate::Pointer + ExhaustiveOptimizerObserver; + + /** Sets the degree of freedom for optimzation*/ + tDegreeOfFreedomEnum GetDegreeOfFreedom(); + /** Returns user components only of the autoReg + * */ + Optimizer::ParametersType + GetFinalR23Parameters(); + /** Define the region to be used as ROI on fixed images */ + void SetRegionFixed1(int,int,int,int); + void SetRegionFixed2(int,int,int,int); + void ResetROIRegions(); + void InizializeRegistration(double stepLength, + double maxTranslation, + eDegreeOfFreedomType dof); + + /** Maximum number of iterations for the optimizer */ + void SetMaxNumberOfIterations(int); + itkGetMacro(R23, itkReg23::Pointer); + + + /** Set Optimizer type */ + void SetOptimizer(std::string); + /** Set Metric type */ + void SetMetric(std::string); + /** Set number of logic CPU to be made available to interpolators*/ void SetNumberOfWorkingUnits(int iN); @@ -116,10 +138,6 @@ public: double GetPanelOffset1(); double GetPanelOffset2(); - /** Sets the degree of freedom for omptimzation*/ - void SetDegreeOfFreedom(tDegreeOfFreedomEnum); - tDegreeOfFreedomEnum GetDegreeOfFreedom(); - void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double); void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); @@ -147,16 +165,6 @@ public: void SetInitialTranslations(double, double, double); void SetInitialRotations(double, double, double); - /** Returns user components only of the autoReg - * */ - Optimizer::ParametersType - GetFinalR23Parameters(); - - /** Initialize the registration pipeline*/ - void InitializeRegistration(double, double, eDegreeOfFreedomType); - /** Start the registration process*/ - int StartRegistration(std::string extraInfo); - ///** Get transform parameters for 3D Volume */ /** Get the complete transform, likely for export to SRO. * This includes user and hidden contributions. @@ -215,25 +223,6 @@ public: void WriteProjectionImages(); void Write2DImages(); - /** Auto Reg23 methods */ - - /** Get the current cost function value from the optimizer*/ - double GetOptimizerValue(); - - /** Get the cost function value for the current transform*/ - double GetCurrentMetricValue(); - - /** Set Optimizer type */ - void SetOptimizer(std::string); - /** Set Metric type */ - void SetMetric(std::string); - /** Use full image as ROI */ - void SetFullROI(bool); - /** Maximum number of iterations for the optimizer */ - void SetMaxNumberOfIterations(int); - /** Define the region to be used as ROI on fixed images */ - void SetRegionFixed1(int,int,int,int); - void SetRegionFixed2(int,int,int,int); protected: /** Various pixel types */ @@ -268,17 +257,6 @@ private: using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; using ResampleFilterType = itk::ResampleImageFilter; - /** Optimizer which tries to find the minimun (Powell Optimizer)*/ - using PowellOptimizerType = itk::PowellOptimizer; - /** Optimizer which tries to find the minimn (Powell Optimizer)*/ - using AmoebaOptimizerType = itk::AmoebaOptimizer; - /** Optimizer which samples the whol space */ - using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer; - /** Metric to calculate similarites between images*/ - using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric; - using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric; - /** The thing which actuall does the image registration*/ - using RegistrationType = itk::TwoProjectionImageRegistrationMethod; using RoiForRegistration = itk::ImageRegion; @@ -318,20 +296,7 @@ private: imageDRT1In, imageDRT2In; - RegistrationType::Pointer - registration; - MetricType::Pointer - NCCmetric; - MIMetricType::Pointer - MImetric; - PowellOptimizerType::Pointer - PowellOptimizer; - AmoebaOptimizerType::Pointer - AmoebaOptimizer; - ExhaustiveOptimizerType::Pointer - ExhaustiveOptimizer; - ExhaustiveCommandIterationUpdate::Pointer - ExhaustiveOptimizerObserver; + DuplicatorType::Pointer m_LATSourceDupli, @@ -426,23 +391,19 @@ private: m_InternalTransf1, m_InternalTransf2; - double m_OptmizerValue; - int m_MaxNumberOfIterations; RoiForRegistration roiAutoReg1, roiAutoReg2; - bool m_UseExhaustiveOptmizer; - bool m_UseAmeobaOptimizer; - bool m_UseFullROI; - bool m_UseMutualInformation; - bool m_UseDumptoFile; - int iNWUnits; + + itk::itkReg23::Pointer + m_R23; + }; diff --git a/reg23Topograms/itkDTRrecon/itkReg23.cpp b/reg23Topograms/itkDTRrecon/itkReg23.cpp new file mode 100644 index 0000000..4f1c903 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkReg23.cpp @@ -0,0 +1,401 @@ +#include "itkReg23.h" + + + +#include "itkTimeProbesCollectorBase.h" + +namespace itk +{ + +itkReg23::itkReg23() +{ + // Metrics + NCCmetric = MetricType::New(); + MImetric = MIMetricType::New(); + + // Optimizers + PowellOptimizer = PowellOptimizerType::New(); + AmoebaOptimizer = AmoebaOptimizerType::New(); + ExhaustiveOptimizer = ExhaustiveOptimizerType::New(); + + // Registration + registration = RegistrationType::New(); + + //internalTransforms + Transform1 = TransformType::New(); + Transform2 = TransformType::New(); + IsocTransform = TransformType::New(); + + registration->SetTransform1(Transform1); + registration->SetTransform2(Transform2); + + + + // to be provided by the user + m_r23Meta = nullptr; + m_Volume = nullptr; + m_PA = nullptr; + m_LAT = nullptr; + m_InternalTransf1 = nullptr; + m_InternalTransf2 = nullptr; + m_filter1 = nullptr; + m_filter2 = nullptr; + m_TransformMetaInfo = nullptr; + m_OptimizerObserver = nullptr; + + // m_UseFullROI = true; // if the full image ROI shall be used + m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file. + +} + + +itkReg23::~itkReg23() +{ + +} + +void itkReg23::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + +void itkReg23::InitializeRegistration( + double stepLength, + double maxTranslation, + eDegreeOfFreedomType dof) +{ + std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; + + if (!m_PA) { + itkExceptionMacro(<< "PA topogram not present"); + } + + if (!m_LAT) { + itkExceptionMacro(<< "LAT topogram not present"); + } + + if (!m_Volume) { + itkExceptionMacro(<< "CT data not present"); + } + + if (!m_r23Meta) { + itkExceptionMacro(<< "r23Meta data not present"); + } + + if (!m_OptimizerObserver) { + itkExceptionMacro(<< "m_OptimizerObserver data not present"); + } + + if (!m_InternalTransf1) { + itkExceptionMacro(<< "m_InternalTransf1 data not present"); + } + + if (!m_InternalTransf2) { + itkExceptionMacro(<< "m_InternalTransf2 data not present"); + } + + if (!m_filter1) { + itkExceptionMacro(<< "m_filter1 data not present"); + } + + if (!m_filter2) { + itkExceptionMacro(<< "m_filter2 data not present"); + } + + if (!m_TransformMetaInfo) { + itkExceptionMacro(<< "m_TransformMetaInfo data not present"); + } + + if(!m_interpolator1){ + itkExceptionMacro(<< "m_interpolator1 data not present"); + } + + if(!m_interpolator2){ + itkExceptionMacro(<< "m_interpolator2 data not present"); + } + + AmoebaOptimizerType::ParametersType simplexDelta(3); + ExhaustiveOptimizerType::StepsType steps(3); + const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651. + //const double stepLength = 0.1; + + switch (m_r23Meta->GetOptimizerType()) { + + case tOptimizerTypeEnum::POWELL: + std::cout<< "Using POWELL Optimizer" <SetMaximize(false); // for NCC and MI + PowellOptimizer->SetMaximumIteration(m_r23Meta->GetMaxIterations()); + PowellOptimizer->SetMaximumLineIteration(4); // for Powell's method + PowellOptimizer->SetStepLength(stepLength); + PowellOptimizer->SetStepTolerance(0.01); + PowellOptimizer->SetValueTolerance(0.000002); + PowellOptimizer->RemoveAllObservers(); + PowellOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver); + + registration->SetOptimizer(PowellOptimizer); + + break; + + case tOptimizerTypeEnum::AMOEBA: + std::cout<< "Using AMOEBA Optimizer" <SetMinimize(true); + AmoebaOptimizer->SetMaximumNumberOfIterations(m_r23Meta->GetMaxIterations()); + AmoebaOptimizer->SetParametersConvergenceTolerance(0.1); + AmoebaOptimizer->SetFunctionConvergenceTolerance(0.000002); + AmoebaOptimizer->SetAutomaticInitialSimplex(false); + //Initial size of the simplex (otherwise it is a very small number and optimizer stops immeditaly) + // 2 mm / 2 degrees seems to be a good value which performs well. + if (GetNumberOfDOF(dof) == 0) { + itkExceptionMacro(<< "Unkown or unsupported degree of freedom"); + } + for (int i = 0; i < GetNumberOfDOF(dof); i++) { + simplexDelta[i] = 2.0; + } + AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta); + AmoebaOptimizer->RemoveAllObservers(); + AmoebaOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver); + + registration->SetOptimizer(AmoebaOptimizer); + + break; + + case tOptimizerTypeEnum::EXHAUSTIVE: + std::cout<< "Using Extensive Optimizer" <SetNumberOfSteps(steps); + ExhaustiveOptimizer->SetStepLength(stepLength); + + registration->SetOptimizer(ExhaustiveOptimizer); + + break; + + default: + break; + + } + + switch (m_r23Meta->GetMetricType()) { + + case tMetricTypeEnum::MI: + std::cout<< "Using MI Metric" <SetMetric(MImetric); + MImetric->ComputeGradientOff(); + MImetric->SetMaxTranslation(maxTranslation); + + break; + + case tMetricTypeEnum::NCC: + std::cout<< "Using NCC Metric" <SetMetric(NCCmetric); + NCCmetric->ComputeGradientOff(); + NCCmetric->SetSubtractMean(true); + NCCmetric->SetMaxTranslation(maxTranslation); + + break; + + default: + break; + + } + + registration->SetFixedImage1(m_PA); + registration->SetFixedImage2(m_LAT); + registration->SetMovingImage(m_Volume); + + registration->SetFixedImageRegion1(m_roiAutoReg1); + registration->SetFixedImageRegion2(m_roiAutoReg2); + std::cout<< "Using Region1:"<SetTransformMetaInfo(m_r23Meta); + + registration->SetinternalMeta1(m_InternalTransf1); + registration->SetinternalMeta2(m_InternalTransf2); + + registration->SetInterpolator1(m_interpolator1); + registration->SetInterpolator2(m_interpolator2); + + registration->SetFilter1(m_filter1); + registration->SetFilter2(m_filter2); + + + IsocTransform = TransformType::New(); + IsocTransform->SetComputeZYX(true); + IsocTransform->SetIdentity(); + + IsocTransform->SetRotation( + m_TransformMetaInfo->GetR()[0], + m_TransformMetaInfo->GetR()[1], + m_TransformMetaInfo->GetR()[2] + ); + + TransformType::OutputVectorType TranslV; + TranslV[0] = m_TransformMetaInfo->GetT()[0]; + TranslV[1] = m_TransformMetaInfo->GetT()[1]; + TranslV[2] = m_TransformMetaInfo->GetT()[2]; + + IsocTransform->SetTranslation(TranslV); + + registration->SetIsocIECTransform(IsocTransform); + + // if (verbose) { + // registration->DebugOn(); + // registration->Print(std::cout); + // } + + + std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl; +} + + + +int itkReg23::StartRegistration(std::string extraInfo) +{ + + std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; + // TODO: Check if the registartion pipeline has been initialized + using ParametersType = RegistrationType::ParametersType; + + //auto startParameters = transform1->GetParameters(); + + time_t t = time(0); // get time now + struct tm* now = localtime(&t); + + bool useDumpToFile = false; + + char buffer[255]; + strftime(buffer, 255, "test-%F-%H%M%S.txt", now); + + std::fstream fs; + + if (useDumpToFile) { + fs.open(buffer, std::fstream::out); + fs << extraInfo; + fs << "Value\tX\tY\tZ " << std::endl; + + // if (r23Meta->GetOptimizerType() == tOptimizerTypeEnum::EXHAUSTIVE) { + // ExhaustiveOptimizerObserver->set_stream(fs); + // } + } + + // Start the registration + // ~~~~~~~~~~~~~~~~~~~~~~ + + // Create a timer to record calculation time. + itk::TimeProbesCollectorBase timer; + + if (verbose) { + std::cout << "Starting the registration now..." << std::endl; + } + + try { + // timer.Start("Registration"); + // Start the registration. + registration->StartRegistration(); + // timer.Stop("Registration"); + } catch (itk::ExceptionObject& err) { + + registration->ResetPipeline(); + std::cout << "ExceptionObject caught !" << std::endl; + std::cout << err << std::endl; + return -1; + } + + fs.close(); + + + itk::Optimizer::ParametersType finalParameters = + this->GetCurrentPosition(); + std::cout<<" FinalPars: "<< finalParameters <GetOptimizerType()) { + + case tOptimizerTypeEnum::POWELL: + finalParameters = PowellOptimizer->GetCurrentPosition(); + m_OptmizerValue = PowellOptimizer->GetValue(); + break; + + case tOptimizerTypeEnum::AMOEBA: + finalParameters = AmoebaOptimizer->GetCurrentPosition(); + m_OptmizerValue = AmoebaOptimizer->GetValue(); + break; + + case tOptimizerTypeEnum::EXHAUSTIVE: + finalParameters = ExhaustiveOptimizer->GetCurrentPosition(); + break; + + default: + break; + + } + + return finalParameters; +} + +double itkReg23::GetCurrentMetricValue() +{ + + switch (m_r23Meta->GetMetricType()) { + + case tMetricTypeEnum::MI: + if(!MImetric){ + itkExceptionMacro(<< "MImetric not present"); + return -1.; + } else { + return MImetric->GetValue(); + } + + break; + + case tMetricTypeEnum::NCC: + if(!NCCmetric){ + itkExceptionMacro(<< "NCCmetric not present"); + return -1.; + } else { + return NCCmetric->GetValue(); + } + break; + + default: + break; + + } + + return -1.; + +} + + +//void itkReg23::SetFullROI(bool fullROI) +//{ +// m_UseFullROI = fullROI; +// // TODO: Remove this function when ROI functinalitz has been implemented. +//} + + +} diff --git a/reg23Topograms/itkDTRrecon/itkReg23.h b/reg23Topograms/itkDTRrecon/itkReg23.h new file mode 100644 index 0000000..0154672 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkReg23.h @@ -0,0 +1,197 @@ +#ifndef ITKREG23_H +#define ITKREG23_h + +#include "itkCommand.h" + +#include "itkExhaustiveOptimizer.h" +#include "itkMutualInformationTwoImageToOneImageMetric.h" +#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h" +#include "itkPowellOptimizer.h" +#include "itkTwoProjectionImageRegistrationMethod.h" +#include "itkAmoebaOptimizer.h" + +#include "itkCommand.h" +#include "itkObject.h" +#include "itkObjectFactory.h" +#include "itkSmartPointer.h" + +#include "itkImageProcessorHelpers.h" +#include "itkQtIterationUpdate.h" + + +namespace itk{ + +class ITK_EXPORT itkReg23 : public itk::Object +{ + +public: + + /** Standard "Self" typedef. */ + typedef itkReg23 Self; + /** Standard "Superclass" typedef. */ + typedef Object Superclass; + /** Smart pointer typedef support */ + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + /** Method of creation through the object factory. */ + itkNewMacro(Self); + /** Run-time type information (and related methods). */ + itkTypeMacro(itkReg23, Object); + + using ChangeInformationFilterType = itk::ChangeInformationImageFilter; + using RoiForRegistration = itk::ImageRegion; + using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; + + /* ---- User provided */ + + itkSetMacro(OptimizerObserver, CommandIterationUpdate::Pointer); + itkGetMacro(OptimizerObserver, CommandIterationUpdate::Pointer); + + itkSetMacro(r23Meta, R23MetaInformation::Pointer); + itkGetMacro(r23Meta, R23MetaInformation::Pointer); + + itkSetMacro(Volume, InternalImageType::Pointer); + itkGetMacro(Volume, InternalImageType::Pointer); + + itkSetMacro(PA, InternalImageType::Pointer); + itkGetMacro(PA, InternalImageType::Pointer); + + itkSetMacro(LAT, InternalImageType::Pointer); + itkGetMacro(LAT, InternalImageType::Pointer); + + itkSetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer); + itkGetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer); + + itkSetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer); + itkGetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer); + + itkSetMacro(filter1, ChangeInformationFilterType::Pointer); + itkGetMacro(filter1, ChangeInformationFilterType::Pointer); + + itkSetMacro(filter2, ChangeInformationFilterType::Pointer); + itkGetMacro(filter2, ChangeInformationFilterType::Pointer); + + itkSetMacro(interpolator1, InterpolatorType::Pointer); + itkGetMacro(interpolator1, InterpolatorType::Pointer); + + itkSetMacro(interpolator2, InterpolatorType::Pointer); + itkGetMacro(interpolator2, InterpolatorType::Pointer); + + itkSetMacro(TransformMetaInfo, TransformMetaInformation::Pointer); + itkGetMacro(TransformMetaInfo, TransformMetaInformation::Pointer); + + /* ---- User provided END */ + itkSetMacro(roiAutoReg1, RoiForRegistration); + itkGetMacro(roiAutoReg1, RoiForRegistration); + + itkSetMacro(roiAutoReg2, RoiForRegistration); + itkGetMacro(roiAutoReg2, RoiForRegistration); + + + + /** Auto Reg23 methods */ + + /** Initialize the registration pipeline*/ + void InitializeRegistration(double stepLength, + double maxTranslation, + eDegreeOfFreedomType dof); + /** Start the registration process*/ + int StartRegistration(std::string extraInfo); + /** Get the current cost function value from the optimizer*/ + double GetOptimizerValue(); + + /** Get the cost function value for the current transform*/ + double GetCurrentMetricValue(); + + /** Get the cost function value for the current transform*/ + Optimizer::ParametersType GetCurrentPosition(); + + /** Auto Reg23 methods END */ + +protected: + itkReg23(); + virtual ~itkReg23(); + void PrintSelf(std::ostream& os, itk::Indent indent) const; + + +private: + itkReg23(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + /** Optimizer which tries to find the minimun (Powell Optimizer)*/ + using PowellOptimizerType = itk::PowellOptimizer; + /** Optimizer which tries to find the minimn (Powell Optimizer)*/ + using AmoebaOptimizerType = itk::AmoebaOptimizer; + /** Optimizer which samples the whol space */ + using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer; + /** Metric to calculate similarites between images*/ + using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric; + using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric; + /** The thing which actuall does the image registration*/ + using RegistrationType = itk::TwoProjectionImageRegistrationMethod; + + + + RegistrationType::Pointer + registration; + MetricType::Pointer + NCCmetric; + MIMetricType::Pointer + MImetric; + PowellOptimizerType::Pointer + PowellOptimizer; + AmoebaOptimizerType::Pointer + AmoebaOptimizer; + ExhaustiveOptimizerType::Pointer + ExhaustiveOptimizer; + + /* ---- User provided */ + R23MetaInformation::Pointer + m_r23Meta; + + InternalImageType::Pointer + m_Volume, + m_PA, + m_LAT; + + InternalTransformMetaInformation::Pointer + m_InternalTransf1, + m_InternalTransf2; + + ChangeInformationFilterType::Pointer + m_filter1, + m_filter2; + + InterpolatorType::Pointer + m_interpolator1, + m_interpolator2; + + TransformMetaInformation::Pointer + m_TransformMetaInfo; + + CommandIterationUpdate::Pointer + m_OptimizerObserver; + + /* ---- User provided END */ + + + TransformType::Pointer + IsocTransform, + Transform1, + Transform2; + + RoiForRegistration + m_roiAutoReg1, + m_roiAutoReg2; + + bool m_UseDumptoFile; + + + double m_OptmizerValue; + + +}; + +} + +#endif