autoReg23 code restructure
- all automatic registration stuff moved to another itk object - most of input are passed as pointers from the parent
This commit is contained in:
@ -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})
|
||||
|
@ -547,7 +547,9 @@ R23MetaInformation::
|
||||
R23MetaInformation (){
|
||||
|
||||
m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
|
||||
|
||||
m_OptimizerType = tOptimizerTypeEnum::POWELL;
|
||||
m_MetricType = tMetricTypeEnum::NCC;
|
||||
m_MaxIterations = 6;
|
||||
}
|
||||
|
||||
|
||||
|
@ -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 ();
|
||||
|
@ -40,7 +40,7 @@ itkImageProcessor::itkImageProcessor()
|
||||
{
|
||||
iNWUnits = 1;
|
||||
|
||||
// std::cout<<"itkImageProcessor::itkImageProcessor contructor"<<std::endl;
|
||||
|
||||
transform1 = TransformType::New();
|
||||
transform1->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: "<<pPixel <<std::endl;
|
||||
// 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;
|
||||
@ -1054,27 +1021,6 @@ int itkImageProcessor::load2D(const char * pcFName){
|
||||
|
||||
|
||||
|
||||
//int itkImageProcessor::query2DimageReconstructionDiameter(const char *pcFName){
|
||||
|
||||
// /* Check if we can open the file */
|
||||
// gdcm::Reader R;
|
||||
// R.SetFileName(pcFName);
|
||||
// if(!R.Read())
|
||||
// { cerr<<"ERROR: cannot read file: "<<pcFName<<endl;
|
||||
// return -1;
|
||||
// }
|
||||
|
||||
// /* Check if it's a localizer image */
|
||||
// const gdcm::File &file = R.GetFile();
|
||||
// const gdcm::DataSet &ds = file.GetDataSet();
|
||||
|
||||
// char* sTmpString = new char [255];
|
||||
// strcpy( sTmpString,gGetStringValueFromTag(gdcm::Tag(0x0018,0x1100), ds));
|
||||
|
||||
// return std::atoi(sTmpString);
|
||||
//}
|
||||
|
||||
|
||||
double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
|
||||
{
|
||||
|
||||
@ -1140,335 +1086,6 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
|
||||
|
||||
|
||||
|
||||
void itkImageProcessor::InitializeRegistration(
|
||||
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");
|
||||
}
|
||||
|
||||
NCCmetric->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<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: " << 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: "<<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;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
// auto finalParameters = transform1->GetParameters();
|
||||
// std::cout<<"-->startParameters "<<startParameters<<std::endl;
|
||||
// std::cout<<"-->finalParameters1 "<<transform1->GetParameters()<<std::endl;
|
||||
// std::cout<<"-->finalParameters2 "<<transform2->GetParameters()<<std::endl;
|
||||
|
||||
// //optimizer->GetCurrentPosition();
|
||||
// finalParameters = finalParameters - startParameters;
|
||||
|
||||
// std::cout<<"startParameters "<<startParameters<<std::endl;
|
||||
// std::cout<<"finalParameters "<<finalParameters<<std::endl;
|
||||
// std::cout<<"finalParametersDelta "<<finalParameters<<std::endl;
|
||||
// std::cout<<"transform1 "<<transform1->GetCenter()<<std::endl;
|
||||
// std::cout<<"transform1 "<<transform1->GetTranslation()<<std::endl;
|
||||
// std::cout<<"transform2 "<<transform2->GetCenter()<<std::endl;
|
||||
// std::cout<<"transform2 "<<transform2->GetTranslation()<<std::endl;
|
||||
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
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
|
||||
|
||||
|
@ -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<InternalImageType, double>;
|
||||
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
|
||||
|
||||
/** 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<InternalImageType, InternalImageType>;
|
||||
using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
|
||||
/** The thing which actuall does the image registration*/
|
||||
using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
|
||||
|
||||
using RoiForRegistration = itk::ImageRegion<Dimension>;
|
||||
|
||||
@ -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;
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
401
reg23Topograms/itkDTRrecon/itkReg23.cpp
Normal file
401
reg23Topograms/itkDTRrecon/itkReg23.cpp
Normal file
@ -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" <<std::endl;
|
||||
|
||||
PowellOptimizer->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" <<std::endl;
|
||||
|
||||
AmoebaOptimizer->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" <<std::endl;
|
||||
|
||||
steps[0] = numberOfSteps;
|
||||
steps[1] = numberOfSteps;
|
||||
steps[2] = numberOfSteps;
|
||||
ExhaustiveOptimizer->SetNumberOfSteps(steps);
|
||||
ExhaustiveOptimizer->SetStepLength(stepLength);
|
||||
|
||||
registration->SetOptimizer(ExhaustiveOptimizer);
|
||||
|
||||
break;
|
||||
|
||||
default:
|
||||
break;
|
||||
|
||||
}
|
||||
|
||||
switch (m_r23Meta->GetMetricType()) {
|
||||
|
||||
case tMetricTypeEnum::MI:
|
||||
std::cout<< "Using MI Metric" <<std::endl;
|
||||
|
||||
registration->SetMetric(MImetric);
|
||||
MImetric->ComputeGradientOff();
|
||||
MImetric->SetMaxTranslation(maxTranslation);
|
||||
|
||||
break;
|
||||
|
||||
case tMetricTypeEnum::NCC:
|
||||
std::cout<< "Using NCC Metric" <<std::endl;
|
||||
registration->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:"<<m_roiAutoReg1 <<std::endl;
|
||||
std::cout<< "Using Region2:"<<m_roiAutoReg2 <<std::endl;
|
||||
|
||||
|
||||
registration->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 <<std::endl;
|
||||
|
||||
|
||||
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
double itkReg23::GetOptimizerValue()
|
||||
{
|
||||
return m_OptmizerValue;
|
||||
}
|
||||
|
||||
|
||||
/** Get the cost function value for the current transform*/
|
||||
Optimizer::ParametersType itkReg23::GetCurrentPosition(){
|
||||
|
||||
|
||||
itk::Optimizer::ParametersType finalParameters;
|
||||
|
||||
switch (m_r23Meta->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.
|
||||
//}
|
||||
|
||||
|
||||
}
|
197
reg23Topograms/itkDTRrecon/itkReg23.h
Normal file
197
reg23Topograms/itkDTRrecon/itkReg23.h
Normal file
@ -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<Self> Pointer;
|
||||
typedef itk::SmartPointer<const Self> ConstPointer;
|
||||
/** Method of creation through the object factory. */
|
||||
itkNewMacro(Self);
|
||||
/** Run-time type information (and related methods). */
|
||||
itkTypeMacro(itkReg23, Object);
|
||||
|
||||
using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
|
||||
using RoiForRegistration = itk::ImageRegion<Dimension>;
|
||||
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
|
||||
|
||||
/* ---- 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<InternalImageType, InternalImageType>;
|
||||
using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
|
||||
/** The thing which actuall does the image registration*/
|
||||
using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
|
||||
|
||||
|
||||
|
||||
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
|
Reference in New Issue
Block a user