diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index c93893c..dff2dd7 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -2,6 +2,8 @@ cmake_minimum_required(VERSION 2.8.12) SET(LIB_NAME itkDTRrecon) +add_subdirectory(autoreg) + find_package(ITK REQUIRED) include(${ITK_USE_FILE}) INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS}) @@ -14,17 +16,21 @@ INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS}) SET(SRCS itkImageProcessor.cpp + itkImageProcessorHelpers.cpp vtkContourTopogramProjectionFilter.cxx DRTMetaInformation.cpp - + itkReg23.cpp ) SET(HDR itkImageProcessor.h + itkImageProcessorHelpers.h + itkQtIterationUpdate.h itkgSiddonJacobsRayCastInterpolateImageFunction.h itkgSiddonJacobsRayCastInterpolateImageFunction.hxx vtkContourTopogramProjectionFilter.h DRTMetaInformation.h + itkReg23.h ) ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR}) @@ -34,6 +40,7 @@ SET(LINK_LIBS ${VTK_LIBRARIES} ${ITK_LIBRARIES} gdcmMSFF + autoreg ) TARGET_LINK_LIBRARIES( diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 8368cc7..1681aad 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -36,7 +36,7 @@ void TopogramImageMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const { - Superclass::PrintSelf(os, indent); + Superclass::PrintSelf(os, indent); } @@ -53,21 +53,21 @@ SetPatientOrientation(tPatOrientation m_orient) this->m_PatientOrientation = m_orient; for(int iIdx = 0 ; iIdx < 3; iIdx++){ - for(int jIdx = 0 ; jIdx < 3; jIdx++){ + for(int jIdx = 0 ; jIdx < 3; jIdx++){ - switch (this->m_PatientOrientation) { - case tPatOrientation::HFS: - m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3]; - break; - case tPatOrientation ::FFS: - m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3]; - break; - default: - m_LPS2IECDirections.SetIdentity(); - break; - } + switch (this->m_PatientOrientation) { + case tPatOrientation::HFS: + m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3]; + break; + case tPatOrientation ::FFS: + m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3]; + break; + default: + m_LPS2IECDirections.SetIdentity(); + break; } } + } } @@ -100,7 +100,7 @@ void DRTImageMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const { - Superclass::PrintSelf(os, indent); + Superclass::PrintSelf(os, indent); } @@ -137,69 +137,67 @@ void DRTImageMetaInformation::SetSizeWithBounds( SpacingType Spacing) { - double dXmin = dBounds[0]; - double dYmin = dBounds[1]; - double dZmin = dBounds[2]; - double dXmax = dBounds[3]; - double dYmax = dBounds[4]; - double dZmax = dBounds[5]; +// double dXmin = dBounds[0]; +// double dYmin = dBounds[1]; +// double dZmin = dBounds[2]; +// double dXmax = dBounds[3]; +// double dYmax = dBounds[4]; +// double dZmax = dBounds[5]; - double - dNPixMin1 = 0., - dNPixMin2 = 0., - dNPixMax1 = 0., - dNPixMax2 = 0.; +// double +// dNPixMin1 = 0., +// dNPixMin2 = 0., +// dNPixMax1 = 0., +// dNPixMax2 = 0.; - dNPixMin1 = abs(dXmin / Spacing[0]); - dNPixMin2 = abs(dYmin / Spacing[0]); - dNPixMax1 = abs(dXmax / Spacing[0]); - dNPixMax2 = abs(dYmax / Spacing[0]); +// dNPixMin1 = abs(dXmin / Spacing[0]); +// dNPixMin2 = abs(dYmin / Spacing[0]); +// dNPixMax1 = abs(dXmax / Spacing[0]); +// dNPixMax2 = abs(dYmax / Spacing[0]); - double dPixMaxOverlap =dNPixMin1; +// double dPixMaxOverlap =dNPixMin1; - if(dPixMaxOverlap > dNPixMin2){ - dPixMaxOverlap = dNPixMin2; - } +// if(dPixMaxOverlap > dNPixMin2){ +// dPixMaxOverlap = dNPixMin2; +// } - if(dPixMaxOverlap > dNPixMax1){ - dPixMaxOverlap = dNPixMax1; - } +// if(dPixMaxOverlap > dNPixMax1){ +// dPixMaxOverlap = dNPixMax1; +// } - if(dPixMaxOverlap > dNPixMax2){ - dPixMaxOverlap = dNPixMax2; - } +// if(dPixMaxOverlap > dNPixMax2){ +// dPixMaxOverlap = dNPixMax2; +// } - if(floor(dPixMaxOverlap*2) > MaxSize[0]) { - m_Size[0] = MaxSize[0]; - } else { - m_Size[0] = floor(dPixMaxOverlap*2); - } +// if(floor(dPixMaxOverlap*2) > MaxSize[0]) { +// m_Size[0] = MaxSize[0]; +// } else { +// m_Size[0] = floor(dPixMaxOverlap*2); +// } - dNPixMin1 = abs(dZmin / Spacing[1]); - dNPixMax1 = abs(dZmax / Spacing[1]); +// dNPixMin1 = abs(dZmin / Spacing[1]); +// dNPixMax1 = abs(dZmax / Spacing[1]); - dPixMaxOverlap =dNPixMin1; - if(dPixMaxOverlap > dNPixMax1){ - dPixMaxOverlap = dNPixMax1; - } - if(floor(dPixMaxOverlap*2) > MaxSize[1]) { - m_Size[1] = MaxSize[1]; - } else { - m_Size[1] = floor(dPixMaxOverlap*2); - } +// dPixMaxOverlap =dNPixMin1; +// if(dPixMaxOverlap > dNPixMax1){ +// dPixMaxOverlap = dNPixMax1; +// } +// if(floor(dPixMaxOverlap*2) > MaxSize[1]) { +// m_Size[1] = MaxSize[1]; +// } else { +// m_Size[1] = floor(dPixMaxOverlap*2); +// } + + m_Size[0] = MaxSize[0]; + m_Size[1] = MaxSize[1]; m_Size[2] = 1; m_Spacing = Spacing; - - - return; - - - + return; } @@ -263,9 +261,9 @@ CTVolumeImageMetaInformation::GetProjectionOriginLPS( PointType ProjectionCenter) { -// PointType Origin; -// Origin = this->m_OriginLPS; -// Origin = Origin - this->m_ImportOffset; + // PointType Origin; + // Origin = this->m_OriginLPS; + // Origin = Origin - this->m_ImportOffset; DirectionType IECtoLPS_Directions; IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose(); @@ -279,8 +277,8 @@ CTVolumeImageMetaInformation::GetProjectionOriginLPS( * middle of the volume */ m_ProjectionOriginLPS[2] = this->GetCOV()[2] - this->m_ImportOffset[2]; - return - m_ProjectionOriginLPS; + return + m_ProjectionOriginLPS; } @@ -305,8 +303,8 @@ CTVolumeImageMetaInformation::GetProjectionOriginLPSZero( m_ProjectionOriginLPS[2] = this->GetCOV()[2] - this->m_ImportOffset[2]; - return - m_ProjectionOriginLPS - Origin; + return + m_ProjectionOriginLPS - Origin; } @@ -315,29 +313,29 @@ CTVolumeImageMetaInformation::CalculateRegions(PointType pProjectionOriginLPS) { /* calculate image size in 3D Space by finding the last voxel position */ - PointType Dest; - Dest[0]=(m_Size[0]-1) * m_Spacing [0]; - Dest[1]=(m_Size[1]-1) * m_Spacing [1]; - Dest[2]=(m_Size[2]-1) * m_Spacing [2]; + PointType Dest; + Dest[0]=(m_Size[0]-1) * m_Spacing [0]; + Dest[1]=(m_Size[1]-1) * m_Spacing [1]; + Dest[2]=(m_Size[2]-1) * m_Spacing [2]; - PointType LastVoxelPosition = - (m_ImageDirections * Dest); + PointType LastVoxelPosition = + (m_ImageDirections * Dest); - LastVoxelPosition [0] += m_OriginLPS[0]; - LastVoxelPosition [1] += m_OriginLPS[1]; - LastVoxelPosition [2] += m_OriginLPS[2]; + LastVoxelPosition [0] += m_OriginLPS[0]; + LastVoxelPosition [1] += m_OriginLPS[1]; + LastVoxelPosition [2] += m_OriginLPS[2]; - std::vector vBounds; - vBounds.clear(); - vBounds.push_back( pProjectionOriginLPS[0] - LastVoxelPosition [0] ); - vBounds.push_back( pProjectionOriginLPS[1] - LastVoxelPosition [1] ); - vBounds.push_back( pProjectionOriginLPS[2] - LastVoxelPosition [2] ); - vBounds.push_back( pProjectionOriginLPS[0] - m_OriginLPS [0] ); - vBounds.push_back( pProjectionOriginLPS[1] - m_OriginLPS [1] ); - vBounds.push_back( pProjectionOriginLPS[2] - m_OriginLPS [2] ); + std::vector vBounds; + vBounds.clear(); + vBounds.push_back( pProjectionOriginLPS[0] - LastVoxelPosition [0] ); + vBounds.push_back( pProjectionOriginLPS[1] - LastVoxelPosition [1] ); + vBounds.push_back( pProjectionOriginLPS[2] - LastVoxelPosition [2] ); + vBounds.push_back( pProjectionOriginLPS[0] - m_OriginLPS [0] ); + vBounds.push_back( pProjectionOriginLPS[1] - m_OriginLPS [1] ); + vBounds.push_back( pProjectionOriginLPS[2] - m_OriginLPS [2] ); - return - vBounds; + return + vBounds; } @@ -365,7 +363,7 @@ void CTVolumeImageMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const { - Superclass::PrintSelf(os, indent); + Superclass::PrintSelf(os, indent); } @@ -381,21 +379,21 @@ SetPatientOrientation(tPatOrientation m_orient) this->m_PatientOrientation = m_orient; for(int iIdx = 0 ; iIdx < 3; iIdx++){ - for(int jIdx = 0 ; jIdx < 3; jIdx++){ + for(int jIdx = 0 ; jIdx < 3; jIdx++){ - switch (this->m_PatientOrientation) { - case tPatOrientation::HFS: - m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3]; - break; - case tPatOrientation ::FFS: - m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3]; - break; - default: - m_LPS2IECDirections.SetIdentity(); - break; - } + switch (this->m_PatientOrientation) { + case tPatOrientation::HFS: + m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3]; + break; + case tPatOrientation ::FFS: + m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3]; + break; + default: + m_LPS2IECDirections.SetIdentity(); + break; } } + } } @@ -452,6 +450,10 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_SCD2Offset = 0.; + this->m_Panel1Offset = 0.; + + this->m_Panel2Offset = 0.; + this->m_IntensityThreshold=0.; this->m_DRT1Size.Fill(0.); @@ -472,13 +474,15 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_ProjectionCenter.Fill(0.); + this->m_UseRotationsForHiddenTransform = true; + } void DRTProjectionGeometryImageMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const { - Superclass::PrintSelf(os, indent); + Superclass::PrintSelf(os, indent); } @@ -505,7 +509,7 @@ void RTGeometryMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const { - Superclass::PrintSelf(os, indent); + Superclass::PrintSelf(os, indent); } @@ -515,23 +519,100 @@ RTGeometryMetaInformation } +InternalTransformMetaInformation::InternalTransformMetaInformation(){ + m_pCalProjCenter.Fill(0.); + m_pRTIsocenter.Fill(0.); + + m_IECtoLPSDirs.SetIdentity(); + +} + +void +InternalTransformMetaInformation +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + +InternalTransformMetaInformation +::~InternalTransformMetaInformation () +{ + +} + + +R23MetaInformation:: +R23MetaInformation (){ + + m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + m_OptimizerType = tOptimizerTypeEnum::POWELL; + m_MetricType = tMetricTypeEnum::NCC; + m_MaxIterations = 6; +} + + +void +R23MetaInformation +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + + +R23MetaInformation +::~R23MetaInformation () +{ + +} TransformMetaInformation :: TransformMetaInformation (){ - m_Translations.Fill(0.); + m_HiddenTranslations.Fill(0.); - m_Rotations.Fill(0.); + m_HiddenRotations.Fill(0.); - m_UserRotations.Fill(0.); + m_UserRotations.Fill(0.); - m_UserTranslations.Fill(0.); + m_UserTranslations.Fill(0.); - m_ReferenceTransform.Fill(0.); + //m_ReferenceTransform.Fill(0.); - m_CurrentTransform.Fill(0.); + //m_CurrentTransform.Fill(0.); + + // m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + +} + +TransformMetaInformation::PointType +TransformMetaInformation::GetT(){ + + PointType + pT; + + pT[0] = m_UserTranslations[0] + m_HiddenTranslations[0]; + pT[1] = m_UserTranslations[1] + m_HiddenTranslations[1]; + pT[2] = m_UserTranslations[2] + m_HiddenTranslations[2]; + + return + pT; + +} + +TransformMetaInformation::PointType +TransformMetaInformation::TransformMetaInformation::GetR(){ + + PointType + pR; + + pR[0] = m_UserRotations[0] + m_HiddenRotations[0]; + pR[1] = m_UserRotations[1] + m_HiddenRotations[1]; + pR[2] = m_UserRotations[2] + m_HiddenRotations[2]; + + return + pR; } @@ -539,7 +620,7 @@ void TransformMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const { - Superclass::PrintSelf(os, indent); + Superclass::PrintSelf(os, indent); } diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 5ec3556..2c59bfc 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -11,15 +11,83 @@ namespace itk { -typedef enum eProjectionOrientationType -{NA = 0, LAT=2, PA=1} -tProjOrientationType; +typedef enum eProjectionOrientationType{ + NA = 0, + LAT=2, + PA=1 +} tProjOrientationType; /** Enum type for Orientation */ -typedef enum eImageOrientationType -{NotDefined = 0, HFS = 1, FFS = 2} -tPatOrientation; +typedef enum eImageOrientationType{ + NotDefined = 0, + HFS = 1, + FFS = 2 +} tPatOrientation; +typedef enum eDegreeOfFreedomType { + UNKNOWN = 0, + THREE_DOF, + SIX_DOF, + X_ONLY, + Y_ONLY, + Z_ONLY, + ROTATION_ONLY, + OTHER +} tDegreeOfFreedomEnum; + + +typedef enum eOptimizerType{ + POWELL, + AMOEBA, + EXHAUSTIVE +} tOptimizerTypeEnum; + +typedef enum eMetricType{ + NCC, + MI +} tMetricTypeEnum; + + +inline int GetNumberOfDOF(eDegreeOfFreedomType dof) +{ + switch (dof) { + case SIX_DOF: + return 6; + break; + case THREE_DOF: + case ROTATION_ONLY: + return 3; + break; + case X_ONLY: + case Y_ONLY: + case Z_ONLY: + return 1; + break; + default: + return 0; + } +} + + +class DegreeOfFreedom : public itk::Object { + +public: + typedef DegreeOfFreedom Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + + itkNewMacro(Self); + +protected: + bool + m_TranslationX, + m_TranslationY, + m_TranslationZ, + m_RotationX, + m_RotationY, + m_RotationZ; + tDegreeOfFreedomEnum m_dof; +}; class TopogramImageMetaInformation : public itk::Object{ @@ -317,6 +385,12 @@ public: itkSetMacro(SCD2Offset, double); itkGetMacro(SCD2Offset, double); + itkSetMacro(Panel1Offset, double); + itkGetMacro(Panel1Offset, double); + + itkSetMacro(Panel2Offset, double); + itkGetMacro(Panel2Offset, double); + itkSetMacro(IntensityThreshold, double); itkGetMacro(IntensityThreshold, double); @@ -347,6 +421,10 @@ public: itkSetMacro(ProjectionCenter, PointType); itkGetMacro(ProjectionCenter, PointType); + itkSetMacro(UseRotationsForHiddenTransform, bool); + itkGetMacro(UseRotationsForHiddenTransform, bool); + + protected: double @@ -358,6 +436,8 @@ protected: m_SCD2, m_SCD1Offset, m_SCD2Offset, + m_Panel1Offset, + m_Panel2Offset, m_IntensityThreshold; @@ -384,6 +464,9 @@ protected: m_ProjectionCenterOffset2; + bool + m_UseRotationsForHiddenTransform; + /** Default Constructor **/ DRTProjectionGeometryImageMetaInformation (); /** Default Destructor **/ @@ -451,6 +534,126 @@ private: }; + + + +class InternalTransformMetaInformation : + public itk::Object{ + +public: + /** standard typedefs **/ + typedef InternalTransformMetaInformation Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + + typedef itk::Point PointType; + typedef itk::Matrix TransformMatrixType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(InternalTransformMetaInformation, itk::Object); + + /** object information streaming **/ + void PrintSelf(std::ostream& os, itk::Indent indent) const; + + itkSetMacro(pCalProjCenter,PointType); + itkSetMacro(pRTIsocenter,PointType); + itkSetMacro(IECtoLPSDirs,TransformMatrixType); + itkGetMacro(pCalProjCenter,PointType); + itkGetMacro(pRTIsocenter,PointType); + itkGetMacro(IECtoLPSDirs,TransformMatrixType); + + +protected: + + PointType + m_pCalProjCenter, + m_pRTIsocenter; + + + TransformMatrixType + m_IECtoLPSDirs; + + + /** Default Constructor **/ + InternalTransformMetaInformation (); + /** Default Destructor **/ + virtual ~InternalTransformMetaInformation (); + +private: + /** purposely not implemented **/ + InternalTransformMetaInformation (const Self&); + + /** purposely not implemented **/ + void operator=(const Self&); + +}; + + + +class R23MetaInformation : + public itk::Object{ + +public: + /** standard typedefs **/ + typedef R23MetaInformation Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(R23MetaInformation, itk::Object); + + /** object information streaming **/ + 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; + + tOptimizerTypeEnum + m_OptimizerType; + + tMetricTypeEnum + m_MetricType; + + int + m_MaxIterations; + + /** Default Constructor **/ + R23MetaInformation (); + /** Default Destructor **/ + virtual ~R23MetaInformation (); + +private: + /** purposely not implemented **/ + R23MetaInformation (const Self&); + + /** purposely not implemented **/ + void operator=(const Self&); + +}; + class TransformMetaInformation : public itk::Object{ @@ -472,11 +675,11 @@ public: void PrintSelf(std::ostream& os, itk::Indent indent) const; - itkSetMacro(Translations,PointType); - itkGetMacro(Translations,PointType); + itkSetMacro(HiddenTranslations,PointType); + itkGetMacro(HiddenTranslations,PointType); - itkSetMacro(Rotations,PointType); - itkGetMacro(Rotations,PointType); + itkSetMacro(HiddenRotations,PointType); + itkGetMacro(HiddenRotations,PointType); itkSetMacro(UserTranslations,PointType); itkGetMacro(UserTranslations,PointType); @@ -484,25 +687,17 @@ public: itkSetMacro(UserRotations,PointType); itkGetMacro(UserRotations,PointType); - itkSetMacro(ReferenceTransform,TransformMatrixType); - itkGetMacro(ReferenceTransform,TransformMatrixType); - - itkSetMacro(CurrentTransform,TransformMatrixType); - itkGetMacro(CurrentTransform,TransformMatrixType); - + PointType GetT(); + PointType GetR(); protected: PointType - m_Translations, - m_Rotations, + m_HiddenTranslations, + m_HiddenRotations, m_UserTranslations, m_UserRotations; - TransformMatrixType - m_ReferenceTransform, - m_CurrentTransform; - /** Default Constructor **/ TransformMetaInformation (); diff --git a/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt b/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt new file mode 100644 index 0000000..4921214 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt @@ -0,0 +1,31 @@ +SET(LIB_NAME autoreg) +#SET(CMAKE_CXX_FLAGS "-std=c++11 -fPIC") + +find_package(ITK REQUIRED) +include(${ITK_USE_FILE}) +INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS}) + + +SET(HDR +itkDRTHelpers.h +itkgTwoImageToOneImageMetric.h +itkgTwoImageToOneImageMetric.hxx +itkMutualInformationTwoImageToOneImageMetric.h +itkMutualInformationTwoImageToOneImageMetric.hxx +itkNormalizedCorrelationTwoImageToOneImageMetric.h +itkNormalizedCorrelationTwoImageToOneImageMetric.hxx +itkTwoProjectionImageRegistrationMethod.h +itkTwoProjectionImageRegistrationMethod.hxx +) + +ADD_LIBRARY(${LIB_NAME} ${HDR}) + +target_link_libraries( + ${LIB_NAME} + ${ITK_LIBRARIES} +) + +target_include_directories (${LIB_NAME} + PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) + + diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h b/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h new file mode 100644 index 0000000..36580ed --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h @@ -0,0 +1,93 @@ +#ifndef ITKDRTHELPERS_H +#define ITKDRTHELPERS_H + +#include "itkMath.h" +#include +#include +#include + +//Function to get method and class name for logging purposes. +template +const char* computeMethodName(const char (&function)[FL], const char (&prettyFunction)[PFL]) +{ + using reverse_ptr = std::reverse_iterator; + thread_local static char result[PFL]; + const char* locFuncName = std::search(prettyFunction, prettyFunction + PFL - 1, function, function + FL - 1); + const char* locClassName = std::find(reverse_ptr(locFuncName), reverse_ptr(prettyFunction), ' ').base(); + const char* endFuncName = std::find(locFuncName, prettyFunction + PFL - 1, '('); + result[0] = '\0'; + std::strncat(result, locClassName, endFuncName - locClassName + 1); + std::strcat(result, ")"); + return result; +} +#define __COMPACT_PRETTY_FUNCTION__ computeMethodName(__FUNCTION__, __PRETTY_FUNCTION__) + +namespace itk { +/* reference string required for comparison with tag values */ +static const char* ImageOrientationStrings[] = { + "NotDefined", + "HFS", + "FFS", + "HFP", + "FFP", +}; + +// constant for converting degrees into radians +const float dtr = itk::Math::pi_over_180; + +static const bool verbose = false; + +/* this is in the end a IEC to HFS... + * but we keep this for ourselves... + */ +static double Standard_DRT2LPS[9] = { + 1, 0, 0, + 0, 0, -1, + 0, 1, 0 +}; + +static double PAElementsIEC[9] = { + 1, 0, 0, + 0, -1, 0, + 0, 0, -1 +}; + +static double LATElementsIEC[9] = { + 0, 0, -1, + 0, -1, 0, + -1, 0, 0 +}; + +static double HFS2IEC[9] = { + 1, 0, 0, + 0, 0, 1, + 0, -1, 0 +}; + +static double FFS2IEC[9] = { + -1, 0, 0, + 0, 0, -1, + 0, -1, 0 +}; + +static double HFP2IEC[9] = { + -1, 0, 0, + 0, 0, 1, + 0, 1, 0 +}; + +static double FFP2IEC[9] = { + 1, 0, 0, + 0, 0, -1, + 0, 1, 0 +}; + +static double PAT2IEC[9] = { + 1, 0, 0, + 0, 0, 1, + 0, -1, 0 +}; + +} + +#endif // ITKDRTHELPERS_H diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h new file mode 100644 index 0000000..678d7c0 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h @@ -0,0 +1,117 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMutualInformationTwoImageToOneImageMetric_h +#define itkMutualInformationTwoImageToOneImageMetric_h + +#include "itkCovariantVector.h" +#include "itkMacro.h" +#include "itkPoint.h" +#include "itkgTwoImageToOneImageMetric.h" + +namespace itk { +/** \class NormalizedCorrelationTwoImageToOneImageMetric + * \brief Computes similarity between two fixed images and one moving image + * + * This metric computes the correlation between pixels in the two fixed images + * and pixels in the moving image. The spatial correspondance between + * two fixed images and the moving image is established through a Transform. Pixel values are + * taken from the fixed images, their positions are mapped to the moving + * image and result in general in non-grid position on it. Values at these + * non-grid position of the moving image are interpolated using user-selected + * Interpolators. The correlation is normalized by the autocorrelations of both + * the fixed and moving images. + * + * \ingroup RegistrationMetrics + * \ingroup TwoProjectionRegistration + */ +template +class MutualInformationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric { +public: + ITK_DISALLOW_COPY_AND_ASSIGN(MutualInformationTwoImageToOneImageMetric); + + /** Standard class type alias. */ + using Self = MutualInformationTwoImageToOneImageMetric; + using Superclass = gTwoImageToOneImageMetric; + + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(MutualInformationTwoImageToOneImageMetric, Object); + + /** Types transferred from the base class */ + using RealType = typename Superclass::RealType; + using TransformType = typename Superclass::TransformType; + using TransformPointer = typename Superclass::TransformPointer; + using TransformParametersType = typename Superclass::TransformParametersType; + using TransformJacobianType = typename Superclass::TransformJacobianType; + using GradientPixelType = typename Superclass::GradientPixelType; + + using MeasureType = typename Superclass::MeasureType; + using DerivativeType = typename Superclass::DerivativeType; + using FixedImageType = typename Superclass::FixedImageType; + using MovingImageType = typename Superclass::MovingImageType; + using FixedImageConstPointer = typename Superclass::FixedImageConstPointer; + using MovingImageConstPointer = typename Superclass::MovingImageConstPointer; + + /** Get the derivatives of the match measure. */ + void + GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override; + + /** Get the value for single valued optimizers. */ + MeasureType + GetValue(const TransformParametersType& parameters) const override; + + /** Get the value using the current transforms. */ + MeasureType + GetValue() const; + + /** Get value and derivatives for multiple valued optimizers. */ + void + GetValueAndDerivative(const TransformParametersType& parameters, + MeasureType& Value, + DerivativeType& Derivative) const override; + + /** Set/Get SubtractMean boolean. If true, the sample mean is subtracted + * from the sample values in the cross-correlation formula and + * typically results in narrower valleys in the cost fucntion. + * Default value is false. */ + itkSetMacro(SubtractMean, bool); + itkGetConstReferenceMacro(SubtractMean, bool); + itkBooleanMacro(SubtractMean); + +protected: + MutualInformationTwoImageToOneImageMetric(); + ~MutualInformationTwoImageToOneImageMetric() override = default; + void + PrintSelf(std::ostream& os, Indent indent) const override; + +private: + bool m_SubtractMean; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkMutualInformationTwoImageToOneImageMetric.hxx" +#endif + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx new file mode 100644 index 0000000..b7c03a7 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx @@ -0,0 +1,237 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMutualInformationTwoImageToOneImageMetric_hxx +#define itkMutualInformationTwoImageToOneImageMetric_hxx + +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkMattesMutualInformationImageToImageMetric.h" +#include "itkMutualInformationTwoImageToOneImageMetric.h" +#include "itkNearestNeighborInterpolateImageFunction.h" +#include "itkNormalizeImageFilter.h" +#include "itkTranslationTransform.h" + +#include "itkProgressReporter.h" + +#include + +namespace itk { + +template +MutualInformationTwoImageToOneImageMetric::MutualInformationTwoImageToOneImageMetric() +{ + m_SubtractMean = false; +} + +template +typename MutualInformationTwoImageToOneImageMetric::MeasureType +MutualInformationTwoImageToOneImageMetric::GetValue( + const TransformParametersType& parameters) const +{ + FixedImageConstPointer fixedImage1 = this->m_FixedImage1; + + if (!fixedImage1) { + itkExceptionMacro(<< "Fixed image1 has not been assigned"); + } + + FixedImageConstPointer fixedImage2 = this->m_FixedImage2; + + if (!fixedImage2) { + itkExceptionMacro(<< "Fixed image2 has not been assigned"); + } + + if (!this->m_Filter1) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + if (!this->m_Filter2) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + this->SetTransformParameters(parameters); + + return GetValue(); +} + +template +typename MutualInformationTwoImageToOneImageMetric::MeasureType +MutualInformationTwoImageToOneImageMetric::GetValue() const +{ + + FixedImageConstPointer fixedImage1 = this->m_FixedImage1; + + if (!fixedImage1) { + itkExceptionMacro(<< "Fixed image1 has not been assigned"); + } + + FixedImageConstPointer fixedImage2 = this->m_FixedImage2; + + if (!fixedImage2) { + itkExceptionMacro(<< "Fixed image2 has not been assigned"); + } + + if (!this->m_Filter1) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + if (!this->m_Filter2) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + if (!this->m_Transform1) { + itkExceptionMacro(<< "m_Transform1 has not been assigned"); + } + + if (!this->m_Transform2) { + itkExceptionMacro(<< "m_Transform2 has not been assigned"); + } + + +// std::cout<< "----- itkMutualInformationTwoImage GetValue -----" <m_Transform1.GetPointer()->GetCenter() <m_Transform1.GetPointer()->GetParameters() <m_Interpolator1->SetTransform(this->m_Transform1); + this->m_Interpolator1->Initialize(); + this->m_Filter1->Update(); + + using InternalMetricType = itk::MattesMutualInformationImageToImageMetric; + /* + * NearestNeighborType in our images results in not resolving for the cranio-caudal direction. + * in LPS-Z we are very much bound to the CT resolution, therefore if we take the Nearest Neighbor + * we simply do not explore that direction... and will always have an uncertainty of 1 slice thickness. + * */ + //using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction; + /* + * LinearInterpolateImageFunction performs better for us, performance degradation not noticeable. + * */ + using LinearInterpType = itk::LinearInterpolateImageFunction; + + //We need to set transform and parameters. + auto dummyTransform = TransformType::New(); + auto dummyParameters = dummyTransform->GetParameters(); + auto dummyInterpolator1 = LinearInterpType::New(); + + auto metric1 = InternalMetricType::New(); + metric1->SetTransform(dummyTransform); + metric1->SetTransformParameters(dummyParameters); + metric1->SetInterpolator(dummyInterpolator1); + + //auto fixedImageRegion1 = fixedImage1->GetBufferedRegion(); + //metric1->SetFixedImageRegion(fixedImageRegion1); + /* We get the fixed image region from the parent class... */ + metric1->SetFixedImageRegion(Superclass::GetFixedImageRegion1()); + //std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl; + + auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion(); + unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels(); + + //auto numberOfSamples = static_cast(numberOfPixels * 0.50); //100% + // since we have a ROI, then we should not set allPixels to TRUE. + //metric1->UseAllPixelsOn(); + metric1->SetNumberOfHistogramBins(50); + +// InternalImageType::IndexType pIndex; +// pIndex[0]=200; +// pIndex[1]=300; +// InternalImageType::PixelType pPixel = fixedImage1->GetPixel(pIndex); +// std::cout<<"MATTES PIXEL: "<SetFixedImage(fixedImage1); + metric1->SetMovingImage(this->m_Filter1->GetOutput()); + + metric1->Initialize(); + + auto measure1 = metric1->GetValue(dummyParameters); + + this->m_Interpolator2->SetTransform(this->m_Transform2); + this->m_Interpolator2->Initialize(); + this->m_Filter2->Update(); + + auto dummyInterpolator2 = LinearInterpType::New(); + auto metric2 = InternalMetricType::New(); + metric2->SetTransform(dummyTransform); + metric2->SetTransformParameters(dummyParameters); + metric2->SetInterpolator(dummyInterpolator2); + +// auto fixedImageRegion2 = fixedImage2->GetBufferedRegion(); +// metric2->SetFixedImageRegion(fixedImageRegion2); + metric2->SetFixedImageRegion(Superclass::GetFixedImageRegion2()); + //std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl; + + movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion(); + numberOfPixels = movingImageRegion.GetNumberOfPixels(); + + //numberOfSamples = static_cast(numberOfPixels * 0.50); //100% + //metric2->SetNumberOfSpatialSamples(numberOfSamples); + //metric2->UseAllPixelsOn(); + metric2->SetNumberOfHistogramBins(50); + metric2->SetFixedImage(fixedImage2); + metric2->SetMovingImage(this->m_Filter2->GetOutput()); + metric2->Initialize(); + + auto measure2 = metric2->GetValue(dummyParameters); + + //std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl; + //std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl; + + // Calculate the measure value between fixed image 2 and the moving image + + auto measure = (measure1 + measure2) / 2.0; + +// auto oldprecision = std::cout.precision(); +// std::cout.precision(std::numeric_limits::digits10 + 2); +// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout << "Measure: " << measure << std::endl; +// std::cout << "=======================================================" << std::endl; +// std::cout.precision(oldprecision); + return measure; +} + +template +void MutualInformationTwoImageToOneImageMetric::GetDerivative( + const TransformParametersType& itkNotUsed(parameters), + DerivativeType& itkNotUsed(derivative)) const +{ + // under construction +} + +template +void MutualInformationTwoImageToOneImageMetric::GetValueAndDerivative( + const TransformParametersType& itkNotUsed(parameters), + MeasureType& itkNotUsed(value), + DerivativeType& itkNotUsed(derivative)) const +{ + // under construction +} + +template +void MutualInformationTwoImageToOneImageMetric::PrintSelf(std::ostream& os, + Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "SubtractMean: " << m_SubtractMean << std::endl; +} + +} // end namespace itk + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h new file mode 100644 index 0000000..0257965 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h @@ -0,0 +1,127 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkNormalizedCorrelationTwoImageToOneImageMetric_h +#define itkNormalizedCorrelationTwoImageToOneImageMetric_h + +#include "itkCovariantVector.h" +#include "itkMacro.h" +#include "itkNormalizedCorrelationImageToImageMetric.hxx" +#include "itkPoint.h" +#include "itkgTwoImageToOneImageMetric.h" + +namespace itk { +/** \class NormalizedCorrelationTwoImageToOneImageMetric + * \brief Computes similarity between two fixed images and one moving image + * + * This metric computes the correlation between pixels in the two fixed images + * and pixels in the moving image. The spatial correspondance between + * two fixed images and the moving image is established through a Transform. Pixel values are + * taken from the fixed images, their positions are mapped to the moving + * image and result in general in non-grid position on it. Values at these + * non-grid position of the moving image are interpolated using user-selected + * Interpolators. The correlation is normalized by the autocorrelations of both + * the fixed and moving images. + * + * \ingroup RegistrationMetrics + * \ingroup TwoProjectionRegistration + */ +template +class NormalizedCorrelationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric { +public: + ITK_DISALLOW_COPY_AND_ASSIGN(NormalizedCorrelationTwoImageToOneImageMetric); + + /** Standard class type alias. */ + using Self = NormalizedCorrelationTwoImageToOneImageMetric; + using Superclass = gTwoImageToOneImageMetric; + + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(NormalizedCorrelationTwoImageToOneImageMetric, Object); + + /** Types transferred from the base class */ + using RealType = typename Superclass::RealType; + using TransformType = typename Superclass::TransformType; + using TransformPointer = typename Superclass::TransformPointer; + using TransformParametersType = typename Superclass::TransformParametersType; + using TransformJacobianType = typename Superclass::TransformJacobianType; + using GradientPixelType = typename Superclass::GradientPixelType; + + using MeasureType = typename Superclass::MeasureType; + using DerivativeType = typename Superclass::DerivativeType; + + using FixedImageType = typename Superclass::FixedImageType; + using FixedImageConstPointer = typename Superclass::FixedImageConstPointer; + using FixedImageMaskPointer = typename Superclass::FixedImageMaskPointer; + using FixedImageRegionType = typename Superclass::FixedImageRegionType; + + using MovingImageType = typename Superclass::MovingImageType; + using MovingImageConstPointer = typename Superclass::MovingImageConstPointer; + using MovingImageMaskPointer = typename Superclass::MovingImageMaskPointer; + + using ChangeInformationFilterPointer = typename Superclass::ChangeInformationFilterPointer; + + /** Get the derivatives of the match measure. */ + void + GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override; + + /** Get the value for single valued optimizers, after applying paramters to the transforms. */ + MeasureType + GetValue(const TransformParametersType& parameters) const override; + + /** Get the value using the current transforms. */ + MeasureType + GetValue() const; + + /** Get value and derivatives for multiple valued optimizers. */ + void + GetValueAndDerivative(const TransformParametersType& parameters, + MeasureType& Value, + DerivativeType& Derivative) const override; + + /** Set/Get SubtractMean boolean. If true, the sample mean is subtracted + * from the sample values in the cross-correlation formula and + * typically results in narrower valleys in the cost fucntion. + * Default value is false. */ + itkSetMacro(SubtractMean, bool); + itkGetConstReferenceMacro(SubtractMean, bool); + itkBooleanMacro(SubtractMean); + +protected: + NormalizedCorrelationTwoImageToOneImageMetric(); + ~NormalizedCorrelationTwoImageToOneImageMetric() override = default; + void + PrintSelf(std::ostream& os, Indent indent) const override; + +private: + bool m_SubtractMean; + + MeasureType CalculateMeasure(int imageId) const; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkNormalizedCorrelationTwoImageToOneImageMetric.hxx" +#endif + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx new file mode 100644 index 0000000..f88d034 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx @@ -0,0 +1,317 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkNormalizedCorrelationTwoImageToOneImageMetric_hxx +#define itkNormalizedCorrelationTwoImageToOneImageMetric_hxx + +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h" + +# + +namespace itk { + +template +NormalizedCorrelationTwoImageToOneImageMetric::NormalizedCorrelationTwoImageToOneImageMetric() +{ + m_SubtractMean = false; +} + +template +typename NormalizedCorrelationTwoImageToOneImageMetric::MeasureType +NormalizedCorrelationTwoImageToOneImageMetric::CalculateMeasure(int imageId) const +{ + + using FixedIteratorType = itk::ImageRegionConstIteratorWithIndex; + using MovingIteratorType = itk::ImageRegionConstIteratorWithIndex; + using AccumulateType = typename NumericTraits::AccumulateType; + + FixedImageConstPointer fixedImage; + ChangeInformationFilterPointer filter; + FixedImageMaskPointer fixedImageMask; + FixedImageRegionType fixedImageRegion; + MeasureType measure; + + typename FixedImageType::IndexType indexFixed; + typename MovingImageType::IndexType indexMoving; + typename Superclass::InputPointType inputPoint; + + if (imageId == 1) { + this->m_Interpolator1->SetTransform(this->m_Transform1); + this->m_Interpolator1->Initialize(); + + fixedImage = this->m_FixedImage1; + fixedImageMask = this->m_FixedImageMask1; + fixedImageRegion = this->GetFixedImageRegion1(); + filter = this->m_Filter1; + } else if (imageId == 2) { + this->m_Interpolator2->SetTransform(this->m_Transform2); + this->m_Interpolator2->Initialize(); + + fixedImage = this->m_FixedImage2; + fixedImageMask = this->m_FixedImageMask2; + fixedImageRegion = this->GetFixedImageRegion2(); + filter = this->m_Filter2; + } else { + itkExceptionMacro(<< "Non supported image id."); + } + + FixedIteratorType ti(fixedImage, fixedImageRegion); + + this->m_NumberOfPixelsCounted = 0; + int numberOfPixelsVisited = 0; + + filter->Update(); + MovingImageConstPointer imageDRTIn = filter->GetOutput(); + //this->WriteImage(imageDRTIn, "", ""); + + AccumulateType sff = NumericTraits::ZeroValue(); + AccumulateType smm = NumericTraits::ZeroValue(); + AccumulateType sfm = NumericTraits::ZeroValue(); + AccumulateType sf = NumericTraits::ZeroValue(); + AccumulateType sm = NumericTraits::ZeroValue(); + + RealType movingValueMin = NumericTraits::max(); + RealType fixedValueMin = NumericTraits::max(); + RealType movingValueMax = NumericTraits::NonpositiveMin(); + RealType fixedValueMax = NumericTraits::NonpositiveMin(); + + RealType maxX = NumericTraits::NonpositiveMin(); + RealType maxY = NumericTraits::NonpositiveMin(); + RealType maxZ = NumericTraits::NonpositiveMin(); + RealType minX = NumericTraits::max(); + RealType minY = NumericTraits::max(); + RealType minZ = NumericTraits::max(); + + AccumulateType m_meanCurr_f = NumericTraits::ZeroValue(); + AccumulateType m_meanCurr_m = NumericTraits::ZeroValue(); + AccumulateType m_meanPrev_f = NumericTraits::ZeroValue(); + AccumulateType m_meanPrev_m = NumericTraits::ZeroValue(); + AccumulateType m_sPrev_f = NumericTraits::ZeroValue(); + AccumulateType m_sPrev_m = NumericTraits::ZeroValue(); + AccumulateType m_sCurr_f = NumericTraits::ZeroValue(); + AccumulateType m_sCurr_m = NumericTraits::ZeroValue(); + AccumulateType m_varianceCurr_f = NumericTraits::ZeroValue(); + AccumulateType m_varianceCurr_m = NumericTraits::ZeroValue(); + + while (!ti.IsAtEnd()) { + + indexFixed = ti.GetIndex(); + ++numberOfPixelsVisited; + + fixedImage->TransformIndexToPhysicalPoint(indexFixed, inputPoint); + + if (fixedImageMask && !fixedImageMask->IsInsideInWorldSpace(inputPoint)) { + ++ti; + continue; + } + + if (this->m_MovingImageMask && !this->m_MovingImageMask->IsInsideInWorldSpace(inputPoint)) { + ++ti; + continue; + } + + const RealType fixedValue = ti.Get(); + + if (fixedValue == -1024) { + ++ti; + continue; + } + + if (imageDRTIn->TransformPhysicalPointToIndex(inputPoint, indexMoving)) { + + double x = inputPoint[0]; + double y = inputPoint[1]; + double z = inputPoint[2]; + maxX = maxX > x ? maxX : x; + minX = minX < x ? minX : x; + maxY = maxY > y ? maxY : y; + minY = minY < y ? minY : y; + maxZ = maxZ > z ? maxZ : z; + minZ = minZ < z ? minZ : z; + + const RealType movingValue = imageDRTIn->GetPixel(indexMoving); + + if (NumericTraits::NonpositiveMin() == movingValue) { + ++ti; + continue; + } + + movingValueMin = movingValue < movingValueMin ? movingValue : movingValueMin; + movingValueMax = movingValue > movingValueMax ? movingValue : movingValueMax; + fixedValueMin = fixedValue < fixedValueMin ? fixedValue : fixedValueMin; + fixedValueMax = fixedValue > fixedValueMax ? fixedValue : fixedValueMax; + + sff += fixedValue * fixedValue; + smm += movingValue * movingValue; + sfm += fixedValue * movingValue; + if (this->m_SubtractMean) { + sf += fixedValue; + sm += movingValue; + } + this->m_NumberOfPixelsCounted++; + + if (this->m_NumberOfPixelsCounted == 1) { + m_meanCurr_f = fixedValue; + m_meanCurr_m = movingValue; + } else { + m_meanPrev_f = m_meanCurr_f; + m_meanPrev_m = m_meanCurr_m; + m_sPrev_f = m_sCurr_f; + m_sPrev_m = m_sCurr_m; + + m_meanCurr_f = m_meanPrev_f + (fixedValue - m_meanPrev_f) / this->m_NumberOfPixelsCounted; + m_sCurr_f = m_sPrev_f + (fixedValue - m_meanPrev_f) * (fixedValue - m_meanCurr_f); + m_varianceCurr_f = m_sCurr_f / (this->m_NumberOfPixelsCounted - 1); + + m_meanCurr_m = m_meanPrev_m + (movingValue - m_meanPrev_m) / this->m_NumberOfPixelsCounted; + m_sCurr_m = m_sPrev_m + (movingValue - m_meanPrev_m) * (movingValue - m_meanCurr_m); + m_varianceCurr_m = m_sCurr_m / (this->m_NumberOfPixelsCounted - 1); + } + } + + ++ti; + } + + if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0) { + sff -= (sf * sf / this->m_NumberOfPixelsCounted); + smm -= (sm * sm / this->m_NumberOfPixelsCounted); + sfm -= (sf * sm / this->m_NumberOfPixelsCounted); + } + + RealType denom = -1.0 * sqrt(sff * smm); + + if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { + measure = sfm / denom; + } else { + measure = NumericTraits::Zero; + } + +// std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl; +// std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl; +// std::cout << "max coordinates: " << maxX << " " << maxX << " " << maxZ << std::endl; +// std::cout << "min coordinates: " << minX << " " << minY << " " << minZ << std::endl; + + return measure; +} + +template +typename NormalizedCorrelationTwoImageToOneImageMetric::MeasureType +NormalizedCorrelationTwoImageToOneImageMetric::GetValue( + const TransformParametersType& parameters) const +{ + if (!this->m_FixedImage1) { + itkExceptionMacro(<< "Fixed image1 has not been assigned"); + } + + if (!this->m_FixedImage2) { + itkExceptionMacro(<< "Fixed image2 has not been assigned"); + } + + if (!this->m_Filter1) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + if (!this->m_Filter2) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + if (!this->SetTransformParameters(parameters)) { + std::cout << "Measure: 0.0, transform Invalid!" << std::endl; + std::cout << "=======================================================" << std::endl; + + return 0.0; + } + return GetValue(); +} + +template +typename NormalizedCorrelationTwoImageToOneImageMetric::MeasureType +NormalizedCorrelationTwoImageToOneImageMetric::GetValue() const +{ + + if (!this->m_FixedImage1) { + itkExceptionMacro(<< "Fixed image1 has not been assigned"); + } + + if (!this->m_FixedImage2) { + itkExceptionMacro(<< "Fixed image2 has not been assigned"); + } + + if (!this->m_Filter1) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + if (!this->m_Filter2) { + itkExceptionMacro(<< "Fixed filter1 has not been assigned"); + } + + // Calculate the measure value between fixed image 1 and the moving image + auto oldprecision = std::cout.precision(); + + +// std::cout<<"Region " <GetFixedImageRegion1() <GetFixedImage1()->GetBufferedRegion() <m_Filter1->GetOutput()->GetBufferedRegion() <::digits10 + 2); +// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; +// std::cout.precision(oldprecision); + + MeasureType measure2 = CalculateMeasure(2); +// std::cout.precision(std::numeric_limits::digits10 + 2); +// std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl; + + auto measure = (measure1 + measure2) / 2.0; + +// std::cout << "Measure: " << measure << std::endl; +// std::cout << "=======================================================" << std::endl; +// std::cout.precision(oldprecision); + + return measure; +} + +template +void NormalizedCorrelationTwoImageToOneImageMetric::GetDerivative( + const TransformParametersType& itkNotUsed(parameters), + DerivativeType& itkNotUsed(derivative)) const +{ + // under construction +} + +template +void NormalizedCorrelationTwoImageToOneImageMetric::GetValueAndDerivative( + const TransformParametersType& itkNotUsed(parameters), + MeasureType& itkNotUsed(value), + DerivativeType& itkNotUsed(derivative)) const +{ + // under construction +} + +template +void NormalizedCorrelationTwoImageToOneImageMetric::PrintSelf(std::ostream& os, + Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "SubtractMean: " << m_SubtractMean << std::endl; +} + +} // end namespace itk + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h new file mode 100644 index 0000000..9ddf934 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -0,0 +1,294 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkTwoProjectionImageRegistrationMethod_h +#define itkTwoProjectionImageRegistrationMethod_h + +#include "itkChangeInformationImageFilter.h" +#include "itkDataObjectDecorator.h" +#include "itkImage.h" +#include "itkMacro.h" +#include "itkProcessObject.h" +#include "itkSingleValuedNonLinearOptimizer.h" +#include "itkgTwoImageToOneImageMetric.h" + +#include "itkObject.h" +#include "itkEventObject.h" + +#include + +namespace itk { + +/** \class TwoProjectionImageRegistrationMethod + * \brief Base class for Image Registration Methods + * + * This Class define the generic interface for a registration method. + * + * This class is templated over the type of the two image to be + * registered. A generic Transform is used by this class. That allows + * to select at run time the particular type of transformation that + * is to be applied for registering the images. + * + * This method use a generic Metric in order to compare the two images. + * the final goal of the registration method is to find the set of + * parameters of the Transformation that optimizes the metric. + * + * The registration method also support a generic optimizer that can + * be selected at run-time. The only restriction for the optimizer is + * that it should be able to operate in single-valued cost functions + * given that the metrics used to compare images provide a single + * value as output. + * + * The terms : Fixed image and Moving image are used in this class + * to indicate what image is being mapped by the transform. + * + * This class uses the coordinate system of the Fixed image as a reference + * and searchs for a Transform that will map points from the space of the + * Fixed image to the space of the Moving image. + * + * For doing so, a Metric will be continously applied to compare the Fixed + * image with the Transformed Moving image. This process also requires to + * interpolate values from the Moving image. + * + * \ingroup RegistrationFilters + * \ingroup TwoProjectionRegistration + */ +template +class TwoProjectionImageRegistrationMethod : public ProcessObject { +public: + ITK_DISALLOW_COPY_AND_ASSIGN(TwoProjectionImageRegistrationMethod); + + /** Standard class type alias. */ + using Self = TwoProjectionImageRegistrationMethod; + using Superclass = ProcessObject; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(TwoProjectionImageRegistrationMethod, ProcessObject); + + /** Type of the Fixed image. */ + using FixedImageType = TFixedImage; + using FixedImageConstPointer = typename FixedImageType::ConstPointer; + + /** Type of the Moving image. */ + using MovingImageType = TMovingImage; + using MovingImageConstPointer = typename MovingImageType::ConstPointer; + using ChangeInformationFilterType = ChangeInformationImageFilter; + using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer; + + /** Type of the metric. */ + using MetricType = gTwoImageToOneImageMetric; + using MetricPointer = typename MetricType::Pointer; + using FixedImageRegionType = typename MetricType::FixedImageRegionType; + + /** Type of the Transform . */ + using TransformType = typename MetricType::TransformType; + using TransformPointer = typename TransformType::Pointer; + + /** Type for the output: Using Decorator pattern for enabling + * the Transform to be passed in the data pipeline */ + using TransformOutputType = DataObjectDecorator; + using TransformOutputPointer = typename TransformOutputType::Pointer; + using TransformOutputConstPointer = typename TransformOutputType::ConstPointer; + + /** Type of the Interpolator. */ + using InterpolatorType = typename MetricType::InterpolatorType; + using InterpolatorPointer = typename InterpolatorType::Pointer; + + /** Type of the optimizer. */ + using OptimizerType = SingleValuedNonLinearOptimizer; + + /** Type of the Transformation parameters This is the same type used to + * represent the search space of the optimization algorithm */ + using ParametersType = typename MetricType::TransformParametersType; + + /** Smart Pointer type to a DataObject. */ + using DataObjectPointer = typename DataObject::Pointer; + + /** Method that initiates the registration. This will Initialize and ensure + * that all inputs the registration needs are in place, via a call to + * Initialize() will then start the optimization process via a call to + * StartOptimization() */ + void + StartRegistration(); + + /** Method that initiates the optimization process. */ + void + StartOptimization(); + + /** Set/Get the Fixed images. */ + void + SetFixedImage1(const FixedImageType* fixedImage1); + void + SetFixedImage2(const FixedImageType* fixedImage2); + itkGetConstObjectMacro(FixedImage1, FixedImageType); + itkGetConstObjectMacro(FixedImage2, FixedImageType); + + /** Set/Get the Moving image. */ + void + SetMovingImage(const MovingImageType* movingImage); + itkGetConstObjectMacro(MovingImage, MovingImageType); + + /** Set/Get the Optimizer. */ + itkSetObjectMacro(Optimizer, OptimizerType); + itkGetConstObjectMacro(Optimizer, OptimizerType); + + /** Set/Get the Metric. */ + itkSetObjectMacro(Metric, MetricType); + itkGetConstObjectMacro(Metric, MetricType); + + /** Set/Get the Transfrom1. */ + itkSetObjectMacro(Transform1, TransformType); + itkGetConstObjectMacro(Transform1, TransformType); + + /** Set/Get the Transfrom2. */ + itkSetObjectMacro(Transform2, TransformType); + itkGetConstObjectMacro(Transform2, TransformType); + + /** Set/Get the IsocIECTransform. */ + itkSetObjectMacro(IsocIECTransform, TransformType); + itkGetConstObjectMacro(IsocIECTransform, TransformType); + + /** Set/Get the Interpolators. */ + itkSetObjectMacro(Interpolator1, InterpolatorType); + itkSetObjectMacro(Interpolator2, InterpolatorType); + itkGetConstObjectMacro(Interpolator1, InterpolatorType); + itkGetConstObjectMacro(Interpolator2, InterpolatorType); + + /** Set/Get the Meta informations. */ + itkSetObjectMacro(TransformMetaInfo, R23MetaInformation); + itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation); + + itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation); + itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation); + + itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation); + itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation); + + /** Set/Get the output filters. */ + itkSetObjectMacro(Filter1, ChangeInformationFilterType); + itkGetConstObjectMacro(Filter1, ChangeInformationFilterType); + itkSetObjectMacro(Filter2, ChangeInformationFilterType); + itkGetConstObjectMacro(Filter2, ChangeInformationFilterType); + + /** Set the region of the fixed image to be considered as region of + interest during the registration. This region will be passed to + the ImageMetric in order to restrict the metric computation to + consider only this region. + \warning The same region can also be set directly into the metric. + please avoid to set the region in both places since this can lead + to inconsistent configurations. */ + void + SetFixedImageRegion1(const FixedImageRegionType& region1); + void + SetFixedImageRegion2(const FixedImageRegionType& region2); + /** Get the region of the fixed image to be considered as region of + interest during the registration. This region will be passed to + the ImageMetric in order to restrict the metric computation to + consider only this region. */ + itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType); + itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType); + /** True if a region has been defined for the fixed image to which + the ImageMetric will limit its computation */ + itkGetMacro(FixedImageRegionDefined1, bool); + itkGetMacro(FixedImageRegionDefined2, bool); + /** Turn on/off the use of a fixed image region to which + the ImageMetric will limit its computation. + \warning The region must have been previously defined using the + SetFixedImageRegion member function */ + itkSetMacro(FixedImageRegionDefined1, bool); + itkSetMacro(FixedImageRegionDefined2, bool); + + /** Initialize by setting the interconnects between the components. */ + virtual void + Initialize(); + + /** Returns the transform resulting from the registration process */ + const TransformOutputType* + GetOutput() const; + + /** Make a DataObject of the correct type to be used as the specified + * output. */ + using Superclass::MakeOutput; + DataObjectPointer + MakeOutput(DataObjectPointerArraySizeType idx) override; + +protected: + TwoProjectionImageRegistrationMethod(); + ~TwoProjectionImageRegistrationMethod() override = default; + void + PrintSelf(std::ostream& os, Indent indent) const override; + + /** Method invoked by the pipeline in order to trigger the computation of + * the registration. */ + void + GenerateData() override; + + /** Provides derived classes with the ability to set this private var */ + itkSetMacro(LastOptimizerParameters, ParametersType); + + /** Entry-point for internal ITK filter observer. **/ +// void OnInternalFilterProgressReceptor(itk::Object *caller, +// const itk::EventObject &event); + + +private: + MetricPointer m_Metric; + OptimizerType::Pointer m_Optimizer; + + MovingImageConstPointer m_MovingImage; + FixedImageConstPointer m_FixedImage1; + FixedImageConstPointer m_FixedImage2; + + TransformPointer m_IsocIECTransform; + + TransformPointer m_Transform1; + TransformPointer m_Transform2; + InterpolatorPointer m_Interpolator1; + InterpolatorPointer m_Interpolator2; + + InternalTransformMetaInformation::Pointer + m_internalMeta1, + m_internalMeta2; + + R23MetaInformation::Pointer + m_TransformMetaInfo; + + ChangeInformationFilterPointer + m_Filter1, + m_Filter2; + + ParametersType m_InitialOptimizerParameters; + ParametersType m_LastOptimizerParameters; + + bool m_FixedImageRegionDefined1; + bool m_FixedImageRegionDefined2; + FixedImageRegionType m_FixedImageRegion1; + FixedImageRegionType m_FixedImageRegion2; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkTwoProjectionImageRegistrationMethod.hxx" +#endif + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx new file mode 100644 index 0000000..2c9b6e0 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -0,0 +1,381 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkTwoProjectionImageRegistrationMethod_hxx +#define itkTwoProjectionImageRegistrationMethod_hxx + +#include "itkCommand.h" + +#include "itkTwoProjectionImageRegistrationMethod.h" + +namespace itk { + +// TwoProjectionImageRegistrationMethod from Jian Wu. +// This part is not really needed and functionality should be moved to +// to the itkImageProcessor.cpp. This is just abstraction with an additional step. + +template +TwoProjectionImageRegistrationMethod::TwoProjectionImageRegistrationMethod() +{ + this->SetNumberOfRequiredOutputs(1); // for the Transform + + m_FixedImage1 = nullptr; // has to be provided by the user. + m_FixedImage2 = nullptr; // has to be provided by the user. + m_MovingImage = nullptr; // has to be provided by the user. + m_Transform1 = nullptr; // has to be provided by the user. + m_Transform2 = nullptr; // has to be provided by the user. + m_IsocIECTransform = nullptr; + m_TransformMetaInfo = nullptr; // has to be provided by the user. + m_internalMeta1 = nullptr; + m_internalMeta2 = nullptr; + m_Interpolator1 = nullptr; // has to be provided by the user. + m_Interpolator2 = nullptr; // has to be provided by the user. + m_Metric = nullptr; // has to be provided by the user. + m_Optimizer = nullptr; // has to be provided by the user. + + m_InitialOptimizerParameters = ParametersType(1); + m_LastOptimizerParameters = ParametersType(1); + + m_InitialOptimizerParameters.Fill(0.0f); + m_LastOptimizerParameters.Fill(0.0f); + + m_FixedImageRegionDefined1 = false; + m_FixedImageRegionDefined2 = false; + + TransformOutputPointer transformDecorator = static_cast(this->MakeOutput(0).GetPointer()); + + this->ProcessObject::SetNthOutput(0, transformDecorator.GetPointer()); +} + +/* + * Set the region of the fixed image 1 to be considered for registration + */ +template +void TwoProjectionImageRegistrationMethod::SetFixedImageRegion1( + const FixedImageRegionType& region1) +{ + m_FixedImageRegion1 = region1; + m_FixedImageRegionDefined1 = true; +} + +/* + * Set the region of the fixed image 2 to be considered for registration + */ +template +void TwoProjectionImageRegistrationMethod::SetFixedImageRegion2( + const FixedImageRegionType& region2) +{ + m_FixedImageRegion2 = region2; + m_FixedImageRegionDefined2 = true; +} + + +/* + * Initialize by setting the interconnects between components. + */ +template +void TwoProjectionImageRegistrationMethod::Initialize() +{ + + if (!m_FixedImage1) { + itkExceptionMacro(<< "FixedImage1 is not present"); + } + + if (!m_FixedImage2) { + itkExceptionMacro(<< "FixedImage2 is not present"); + } + + if (!m_MovingImage) { + itkExceptionMacro(<< "MovingImage is not present"); + } + + if (!m_Metric) { + itkExceptionMacro(<< "Metric is not present"); + } + + if (!m_Optimizer) { + itkExceptionMacro(<< "Optimizer is not present"); + } + + if (!m_Transform1) { + itkExceptionMacro(<< "Transform1 is not present"); + } + + if (!m_Transform2) { + itkExceptionMacro(<< "Transform2 is not present"); + } + + // + // Connect the transform to the Decorator. + auto* transformOutput = static_cast(this->ProcessObject::GetOutput(0)); + + //transformOutput->Set(m_Transform1.GetPointer()); + transformOutput->Set(m_IsocIECTransform.GetPointer()); + + if (!m_Interpolator1) { + itkExceptionMacro(<< "Interpolator1 is not present"); + } + + if (!m_Interpolator2) { + itkExceptionMacro(<< "Interpolator2 is not present"); + } + + if (!m_Filter1) { + itkExceptionMacro(<< "Filter1 is not present"); + } + + if (!m_Filter1) { + itkExceptionMacro(<< "Filter2 is not present"); + } + + if (!m_TransformMetaInfo) { + itkExceptionMacro(<< "TransformMetaInfo is not present"); + } + + if (!m_internalMeta1) { + itkExceptionMacro(<< "internalMeta1 is not present"); + } + + if (!m_internalMeta2) { + itkExceptionMacro(<< "internalMeta1 is not present"); + } + + // Setup the metric + m_Metric->SetMovingImage(m_MovingImage); + m_Metric->SetFixedImage1(m_FixedImage1); + m_Metric->SetFixedImage2(m_FixedImage2); + m_Metric->SetTransform1(m_Transform1); + m_Metric->SetTransform2(m_Transform2); + + m_Metric->SetIsocTransf(m_IsocIECTransform); + m_Metric->SetTransformMetaInfo(m_TransformMetaInfo); + m_Metric->SetinternalMeta1(m_internalMeta1); + m_Metric->SetinternalMeta2(m_internalMeta2); + + m_Metric->SetFilter1(m_Filter1); + m_Metric->SetFilter2(m_Filter2); + + m_Metric->SetInterpolator1(m_Interpolator1); + m_Metric->SetInterpolator2(m_Interpolator2); + + if (m_FixedImageRegionDefined1) { + m_Metric->SetFixedImageRegion1(m_FixedImageRegion1); + std::cout << "m_FixedImageRegionDefined1 is defined " << std::endl; + } else { + m_Metric->SetFixedImageRegion1(m_FixedImage1->GetBufferedRegion()); + } + + if (m_FixedImageRegionDefined2) { + m_Metric->SetFixedImageRegion2(m_FixedImageRegion2); + std::cout << "m_FixedImageRegionDefined2 is defined " << std::endl; + } else { + m_Metric->SetFixedImageRegion2(m_FixedImage2->GetBufferedRegion()); + } + + m_Metric->Initialize(); + + // Setup the optimizer + m_Optimizer->SetScales(m_Metric->GetWeightings()); + m_Optimizer->SetCostFunction(m_Metric); + + m_InitialOptimizerParameters = + m_Metric->GetParameters(); + + m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters); +} + +//template +//void TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor( +// itk::Object *caller, const itk::EventObject &event) +//{ +// itk::ProcessObject *filter = dynamic_cast(caller); + +// if(!itk::ProgressEvent().CheckEvent(&event) || !filter) +// return; + +// if (filter) +// { +//// double p = m_CurrentProgressStart; +//// // filter->GetProgress() in [0;1] +//// p += m_CurrentProgressSpan * filter->GetProgress(); +//// // NOTE: filter is buggy - not in [0;1] if multi-threading is activated! +//// if (p > (m_CurrentProgressStart + m_CurrentProgressSpan)) +//// p = m_CurrentProgressStart + m_CurrentProgressSpan; +//// TaskProgressInfo(m_CurrentProgressDirection, p); + +//// if (m_CancelRequest) +//// filter->SetAbortGenerateData(true); // may be handled by filter + +// std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl; +// } +//} + +/* + * Starts the Registration Process + */ +template +void TwoProjectionImageRegistrationMethod::StartRegistration() +{ + + itk::ProgressReporter progress( + this, + 1, 1000,100); + + ParametersType empty(1); + empty.Fill(0.0); + try { + // initialize the interconnects between components + this->Initialize(); + if (GetDebug()) { + m_Metric->Print(std::cout); + m_Optimizer->Print(std::cout); + this->Print(std::cout); + } + } catch (ExceptionObject& err) { + m_LastOptimizerParameters = empty; + + // pass exception to caller + throw err; + } + + this->StartOptimization(); +} + +/* + * Starts the Optimization process + */ +template +void TwoProjectionImageRegistrationMethod::StartOptimization() +{ + try { + // do the optimization + m_Optimizer->StartOptimization(); + + } catch (ExceptionObject& err) { + // An error has occurred in the optimization. + // Update the parameters + m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition(); + + // Pass exception to caller + throw err; + } + + // get the results + m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition(); + + //Update transform1 and transform2 with the last transform parameters from the optimizer. + m_Metric->SetTransformParameters(m_LastOptimizerParameters); +} + +template +void TwoProjectionImageRegistrationMethod::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Metric: " << m_Metric.GetPointer() << std::endl; + os << indent << "Optimizer: " << m_Optimizer.GetPointer() << std::endl; + os << indent << "Transform1: " << m_Transform1.GetPointer() << std::endl; + os << indent << "Transform2: " << m_Transform2.GetPointer() << std::endl; + os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl; + os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl; + os << indent << "Fixed Image 1: " << m_FixedImage1.GetPointer() << std::endl; + os << indent << "Fixed Image 2: " << m_FixedImage2.GetPointer() << std::endl; + os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl; + os << indent << "Transform MetaInfo: " << m_TransformMetaInfo.GetPointer() << std::endl; + os << indent << "Fixed Image 1 Region Defined: " << m_FixedImageRegionDefined1 << std::endl; + os << indent << "Fixed Image 2 Region Defined: " << m_FixedImageRegionDefined2 << std::endl; + os << indent << "Fixed Image 1 Region: " << m_FixedImageRegion1 << std::endl; + os << indent << "Fixed Image 2 Region: " << m_FixedImageRegion2 << std::endl; + os << indent << "Initial Optimizer Parameters: " << m_InitialOptimizerParameters << std::endl; + os << indent << "Last Optimizer Parameters: " << m_LastOptimizerParameters << std::endl; +} + +template +void TwoProjectionImageRegistrationMethod::GenerateData() +{ + + this->StartRegistration(); +} + +template +const typename TwoProjectionImageRegistrationMethod::TransformOutputType* +TwoProjectionImageRegistrationMethod::GetOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(0)); +} + +template +DataObject::Pointer +TwoProjectionImageRegistrationMethod::MakeOutput(DataObjectPointerArraySizeType output) +{ + switch (output) { + case 0: + return static_cast(TransformOutputType::New().GetPointer()); + break; + default: + itkExceptionMacro("MakeOutput request for an output number larger than the expected number of outputs"); + return nullptr; + } +} + +template +void TwoProjectionImageRegistrationMethod::SetFixedImage1(const FixedImageType* fixedImage1) +{ + itkDebugMacro("setting Fixed Image 1 to " << fixedImage1); + + if (this->m_FixedImage1.GetPointer() != fixedImage1) { + this->m_FixedImage1 = fixedImage1; + + // Process object is not const-correct so the const_cast is required here + this->ProcessObject::SetNthInput(0, const_cast(fixedImage1)); + + this->Modified(); + } +} + +template +void TwoProjectionImageRegistrationMethod::SetFixedImage2(const FixedImageType* fixedImage2) +{ + itkDebugMacro("setting Fixed Image 2 to " << fixedImage2); + + if (this->m_FixedImage2.GetPointer() != fixedImage2) { + this->m_FixedImage2 = fixedImage2; + + // Process object is not const-correct so the const_cast is required here + this->ProcessObject::SetNthInput(0, const_cast(fixedImage2)); + + this->Modified(); + } +} + +template +void TwoProjectionImageRegistrationMethod::SetMovingImage(const MovingImageType* movingImage) +{ + itkDebugMacro("setting Moving Image to " << movingImage); + + if (this->m_MovingImage.GetPointer() != movingImage) { + this->m_MovingImage = movingImage; + + // Process object is not const-correct so the const_cast is required here + this->ProcessObject::SetNthInput(1, const_cast(movingImage)); + + this->Modified(); + } +} + +} // end namespace itk + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h new file mode 100644 index 0000000..2c08287 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h @@ -0,0 +1,319 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkgTwoImageToOneImageMetric_h +#define itkgTwoImageToOneImageMetric_h + +#include "DRTMetaInformation.h" +#include "itkChangeInformationImageFilter.h" +#include "itkDRTHelpers.h" +#include "itkGradientRecursiveGaussianImageFilter.h" +#include "itkImageBase.h" +#include "itkInterpolateImageFunction.h" +#include "itkOptimizer.h" +#include "itkResampleImageFilter.h" +#include "itkSingleValuedCostFunction.h" +#include "itkSpatialObject.h" +#include "itkTransform.h" +#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h" + +#include "itkImageFileWriter.h" +#include "itkRescaleIntensityImageFilter.h" + +#include +#include +#include +#include +#include + +#include + +namespace itk { + +/** \class gTwoImageToOneImageMetric + * \brief Computes similarity between two fixed images and one fixed image. + * + * This Class is templated over the type of the two input images. + * It expects a Transform and two Interpolators to be plugged in. + * This particular class is the base class for a hierarchy of + * similarity metrics. + * + * This class computes a value that measures the similarity + * between two Fixed image and the transformed Moving images. + * The Interpolators are used to compute intensity values on + * non-grid positions resulting from mapping points through + * the Transform. + * + * + * \ingroup RegistrationMetrics + * \ingroup TwoProjectionRegistration + * + */ + +template +class gTwoImageToOneImageMetric : public SingleValuedCostFunction { +public: + ITK_DISALLOW_COPY_AND_ASSIGN(gTwoImageToOneImageMetric); + + /** Standard class type alias. */ + using Self = gTwoImageToOneImageMetric; + using Superclass = SingleValuedCostFunction; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Type used for representing point components */ + using CoordinateRepresentationType = Superclass::ParametersValueType; + + /** Run-time type information (and related methods). */ + itkTypeMacro(gTwoImageToOneImageMetric, SingleValuedCostFunction); + + /** Type of the moving Image. */ + using MovingImageType = TMovingImage; + using MovingImagePixelType = typename TMovingImage::PixelType; + using MovingImageConstPointer = typename MovingImageType::ConstPointer; + + /** Type of the fixed Image. */ + using FixedImageType = TFixedImage; + using FixedImageConstPointer = typename FixedImageType::ConstPointer; + using FixedImageRegionType = typename FixedImageType::RegionType; + + /** Constants for the image dimensions */ + static constexpr unsigned int MovingImageDimension = TMovingImage::ImageDimension; + static constexpr unsigned int FixedImageDimension = TFixedImage::ImageDimension; + + /** Type of the Transform Base class */ + using TransformType = itk::Euler3DTransform; + + using TransformPointer = typename TransformType::Pointer; + using InputPointType = typename TransformType::InputPointType; + using OutputPointType = typename TransformType::OutputPointType; + using TransformParametersType = typename TransformType::ParametersType; + using TransformJacobianType = typename TransformType::JacobianType; + + /** Type of the Interpolator Base class */ + using InterpolatorType = gSiddonJacobsRayCastInterpolateImageFunction; + + /** Gaussian filter to compute the gradient of the Moving Image */ + using RealType = typename NumericTraits::RealType; + using GradientPixelType = CovariantVector; + using GradientImageType = Image; + using GradientImagePointer = SmartPointer; + using GradientImageFilterType = GradientRecursiveGaussianImageFilter; + using GradientImageFilterPointer = typename GradientImageFilterType::Pointer; + using InterpolatorPointer = typename InterpolatorType::Pointer; + + using ResampleImageFilterType = ResampleImageFilter; + using ResampleImageFilterPointer = typename ResampleImageFilterType::Pointer; + using ChangeInformationFilterType = ChangeInformationImageFilter; + using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer; + + /** Type for the mask of the fixed image. Only pixels that are "inside" + this mask will be considered for the computation of the metric */ + typedef SpatialObject FixedImageMaskType; + using FixedImageMaskPointer = typename FixedImageMaskType::Pointer; + + /** Type for the mask of the moving image. Only pixels that are "inside" + this mask will be considered for the computation of the metric */ + typedef SpatialObject MovingImageMaskType; + using MovingImageMaskPointer = typename MovingImageMaskType::Pointer; + + /** Type of the measure. */ + using MeasureType = Superclass::MeasureType; + + /** Type of the derivative. */ + using DerivativeType = Superclass::DerivativeType; + + /** Type of the Transformation parameters This is the same type used to + * represent the search space of the optimization algorithm */ + using ParametersType = Superclass::ParametersType; + + /** Connect the Fixed Image. */ + itkSetConstObjectMacro(FixedImage1, FixedImageType); + + /** Connect the Fixed Image. */ + itkSetConstObjectMacro(FixedImage2, FixedImageType); + + /** Get the Fixed Image. */ + itkGetConstObjectMacro(FixedImage1, FixedImageType); + + /** Get the Fixed Image. */ + itkGetConstObjectMacro(FixedImage2, FixedImageType); + + /** Connect the Moving Image. */ + itkSetConstObjectMacro(MovingImage, MovingImageType); + + /** Get the Moving Image. */ + itkGetConstObjectMacro(MovingImage, MovingImageType); + + /** Connect the Transform1. */ + itkSetObjectMacro(Transform1, TransformType); + + /** Get a pointer to the Transform1. */ + itkGetConstObjectMacro(Transform1, TransformType); + + /** Connect the Transform2. */ + itkSetObjectMacro(Transform2, TransformType); + + /** Get a pointer to the Transform2. */ + itkGetConstObjectMacro(Transform2, TransformType); + + itkSetObjectMacro(IsocTransf, TransformType); + itkGetObjectMacro(IsocTransf, TransformType); + + + /** Connect the Interpolator. */ + itkSetObjectMacro(Interpolator1, InterpolatorType); + + /** Get a pointer to the Interpolator. */ + itkGetConstObjectMacro(Interpolator1, InterpolatorType); + + /** Connect the Interpolator. */ + itkSetObjectMacro(Interpolator2, InterpolatorType); + + /** Get a pointer to the Interpolator. */ + itkGetConstObjectMacro(Interpolator2, InterpolatorType); + + /** Get the number of pixels considered in the computation. */ + itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long); + + /** Connect the DRTGeometryMetaInfo. */ + itkSetObjectMacro(TransformMetaInfo, R23MetaInformation); + + /** Get a pointer to the DRTGeometryMetaInfo. */ + itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation); + + + itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation); + itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation); + + itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation); + itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation); + + /** Set the region over which the metric will be computed */ + itkSetMacro(FixedImageRegion1, FixedImageRegionType); + + /** Set the region over which the metric will be computed */ + itkSetMacro(FixedImageRegion2, FixedImageRegionType); + + /** Get the region over which the metric will be computed */ + itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType); + + /** Get the region over which the metric will be computed */ + itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType); + + /** Set/Get the output filters. */ + itkSetObjectMacro(Filter1, ChangeInformationFilterType); + itkGetConstObjectMacro(Filter1, ChangeInformationFilterType); + itkSetObjectMacro(Filter2, ChangeInformationFilterType); + itkGetConstObjectMacro(Filter2, ChangeInformationFilterType); + + /** Set/Get the moving image mask. */ + itkSetObjectMacro(MovingImageMask, MovingImageMaskType); + itkGetConstObjectMacro(MovingImageMask, MovingImageMaskType); + + /** Set/Get the fixed image mask. */ + itkSetObjectMacro(FixedImageMask1, FixedImageMaskType); + itkSetObjectMacro(FixedImageMask2, FixedImageMaskType); + itkGetConstObjectMacro(FixedImageMask1, FixedImageMaskType); + itkGetConstObjectMacro(FixedImageMask2, FixedImageMaskType); + + /** Set/Get gradient computation. */ + itkSetMacro(ComputeGradient, bool); + itkGetConstReferenceMacro(ComputeGradient, bool); + itkBooleanMacro(ComputeGradient); + + itkSetMacro(MaxTranslation, double); + itkGetConstReferenceMacro(MaxTranslation, double); + + /** Get Gradient Image. */ + itkGetConstObjectMacro(GradientImage, GradientImageType); + + /** Set the parameters defining the Transform. */ + bool + SetTransformParameters(const ParametersType& parameters) const; + + /** Return the number of parameters required by the Transform */ + + unsigned int GetNumberOfParameters() const; + itk::Optimizer::ScalesType GetWeightings() const; + ParametersType GetParameters() const; + + /** Initialize the Metric by making sure that all the components + * are present and plugged together correctly */ + virtual void + Initialize(); + +protected: + gTwoImageToOneImageMetric(); + ~gTwoImageToOneImageMetric() override = default; + void + PrintSelf(std::ostream& os, Indent indent) const override; + + constexpr static unsigned int Dimension = 3; + + mutable unsigned long m_NumberOfPixelsCounted; + mutable double m_bestMeassure; + + FixedImageConstPointer m_FixedImage1; + FixedImageConstPointer m_FixedImage2; + MovingImageConstPointer m_MovingImage; + + mutable TransformPointer m_Transform1; + mutable TransformPointer m_Transform2; + mutable TransformPointer m_IsocTransf; + + InterpolatorPointer m_Interpolator1; + InterpolatorPointer m_Interpolator2; + + bool m_ComputeGradient; + GradientImagePointer m_GradientImage; + + using OutputPixelType = unsigned char; + using OutputImageType = itk::Image; + + mutable FixedImageMaskPointer m_FixedImageMask1; + mutable FixedImageMaskPointer m_FixedImageMask2; + mutable MovingImageMaskPointer m_MovingImageMask; + + R23MetaInformation::Pointer + m_TransformMetaInfo; + + InternalTransformMetaInformation::Pointer + m_internalMeta1, + m_internalMeta2; + + ChangeInformationFilterPointer + m_Filter1, + m_Filter2; + + double m_MaxTranslation; + + /* Writes the image to the disk. If path is empty it will be saved to /img/out/*/ + void WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const; + +private: + FixedImageRegionType m_FixedImageRegion1; + FixedImageRegionType m_FixedImageRegion2; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkgTwoImageToOneImageMetric.hxx" +#endif + +#endif diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx new file mode 100644 index 0000000..46a0b0f --- /dev/null +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -0,0 +1,476 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkgTwoImageToOneImageMetric_hxx +#define itkgTwoImageToOneImageMetric_hxx + +#include "itkgTwoImageToOneImageMetric.h" +#include + +namespace itk { + +template +gTwoImageToOneImageMetric::gTwoImageToOneImageMetric() +{ + m_FixedImage1 = nullptr; // has to be provided by the user. + m_FixedImage2 = nullptr; // has to be provided by the user. + m_MovingImage = nullptr; // has to be provided by the user. + m_Transform1 = nullptr; // has to be provided by the user. + m_Transform2 = nullptr; // has to be provided by the user. + m_TransformMetaInfo = nullptr; // has to be provided by the user. + m_internalMeta1 = nullptr; + m_internalMeta2 = nullptr; + m_Interpolator1 = nullptr; // has to be provided by the user. + m_Interpolator2 = nullptr; // has to be provided by the user. + m_GradientImage = nullptr; // will receive the output of the filter; + m_ComputeGradient = true; // metric computes gradient by default + m_NumberOfPixelsCounted = 0; // initialize to zero + m_bestMeassure = std::numeric_limits::max(); + m_MaxTranslation = std::numeric_limits::max(); +} + +/* + * Set the parameters that define a unique transform + */ +template +bool gTwoImageToOneImageMetric::SetTransformParameters(const ParametersType& parameters) const +{ + + if (!m_Transform1) { + itkExceptionMacro(<< "Transform 1 has not been assigned"); + } + if (!m_Transform2) { + itkExceptionMacro(<< "Transform 2 has not been assigned"); + } + if (!m_internalMeta1){ + itkExceptionMacro(<< "m_internalMeta1 has not been assigned"); + } + if (!m_internalMeta2){ + itkExceptionMacro(<< "m_internalMeta2 has not been assigned"); + } + if (!m_IsocTransf){ + itkExceptionMacro(<< "m_IsocTransf has not been assigned"); + } + + //auto transformParameters = m_Transform1->GetParameters(); + auto transformParameters = m_IsocTransf->GetParameters(); + double TranslationAlongX = transformParameters[3]; + double TranslationAlongY = transformParameters[4]; + double TranslationAlongZ = transformParameters[5]; + double RotationAlongX = transformParameters[0]; + double RotationAlongY = transformParameters[1]; + double RotationAlongZ = transformParameters[2]; + + //std::cout << "m_IsocTransf PARS: "<< transformParameters <GetDegreeOfFreedom()) { + case X_ONLY: + TranslationAlongX = parameters[0]; + break; + case Y_ONLY: + TranslationAlongY = parameters[0]; + break; + case Z_ONLY: + TranslationAlongZ = parameters[0]; + break; + case THREE_DOF: + TranslationAlongX = parameters[0]; + TranslationAlongY = parameters[1]; + TranslationAlongZ = parameters[2]; + break; + case ROTATION_ONLY: + RotationAlongX = parameters[0]; + RotationAlongY = parameters[1]; + RotationAlongZ = parameters[2]; + break; + case SIX_DOF: + TranslationAlongX = parameters[0]; + TranslationAlongY = parameters[1]; + TranslationAlongZ = parameters[2]; + RotationAlongX = parameters[3]; + RotationAlongY = parameters[4]; + RotationAlongZ = parameters[5]; + break; + + default: + itkExceptionMacro(<< "Unsupported degree of freedeom"); + break; + } + +// std::cout << "New Transform Parameters ISOC IEC = " << 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 << " 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 << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; + + + // individual R and T components of the IsoC transform + ImageType3D::PointType + pT; + pT[0] = TranslationAlongX; + pT[1] = TranslationAlongY; + pT[2] = TranslationAlongZ; + + ImageType3D::PointType + pR; + pR[0] = RotationAlongX; + pR[1] = RotationAlongY; + pR[2] = RotationAlongZ; + + // check if we are within tolerance + bool transformValid = + (std::abs(pT[0]) < m_MaxTranslation) && + (std::abs(pT[1]) < m_MaxTranslation) && + (std::abs(pT[2]) < m_MaxTranslation); + + if (transformValid) { + /* + Here we calculate the transform for eack projection corresponding to the isocentric + one we are optimizing. + This calculation takes into account calibrated center of projection for each view. + */ + + // Calculate internal transform1 + TransformType::Pointer CurrTransform; + CurrTransform = CalculateInternalTransformV3( + pT, + pR, + m_internalMeta1->GetpCalProjCenter(), + m_internalMeta1->GetpRTIsocenter(), + m_internalMeta1->GetIECtoLPSDirs()); + + m_Transform1->SetIdentity(); + m_Transform1->SetParameters( + CurrTransform->GetParameters()); + m_Transform1->SetCenter( + m_internalMeta1->GetpCalProjCenter()); + +// std::cout<<"----- itkgTwoImageToOne SetTransformParameters -----"<GetCenter()<GetParameters()<GetpCalProjCenter(), + m_internalMeta2->GetpRTIsocenter(), + m_internalMeta2->GetIECtoLPSDirs()); + + m_Transform2->SetIdentity(); + m_Transform2->SetParameters( + CurrTransform->GetParameters()); + m_Transform2->SetCenter( + m_internalMeta2->GetpCalProjCenter()); + + } else { + std::cout << " Transform invalid, out of range!" << std::endl; + std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; + } + + return transformValid; +} + +template +void gTwoImageToOneImageMetric::Initialize() +{ + + if (!m_IsocTransf) { + itkExceptionMacro(<< "Transform is not present"); + } + + if (!m_Transform1) { + itkExceptionMacro(<< "Transform1 is not present"); + } + + if (!m_Transform2) { + itkExceptionMacro(<< "Transform2 is not present"); + } + + if(!m_internalMeta1) { + itkExceptionMacro(<< "m_internalMeta1 is not present"); + } + if(!m_internalMeta2) { + itkExceptionMacro(<< "m_internalMeta2 is not present"); + } + + if (!m_Interpolator1) { + itkExceptionMacro(<< "Interpolator1 is not present"); + } + if (!m_Interpolator2) { + itkExceptionMacro(<< "Interpolator2 is not present"); + } + + if (!m_FixedImage1) { + itkExceptionMacro(<< "FixedImage1 is not present"); + } + + if (!m_FixedImage2) { + itkExceptionMacro(<< "FixedImage2 is not present"); + } + + if (m_FixedImageRegion1.GetNumberOfPixels() == 0) { + itkExceptionMacro(<< "FixedImageRegion1 is empty"); + } + + if (m_FixedImageRegion2.GetNumberOfPixels() == 0) { + itkExceptionMacro(<< "FixedImageRegion2 is empty"); + } + + // If the image is provided by a source, update the source. + if (m_MovingImage->GetSource()) { + m_MovingImage->GetSource()->Update(); + } + + // If the image is provided by a source, update the source. + if (m_FixedImage1->GetSource()) { + m_FixedImage1->GetSource()->Update(); + } + + if (m_FixedImage2->GetSource()) { + m_FixedImage2->GetSource()->Update(); + } + + //std::cout<<"FixedImageRegion1 " <GetFixedImageRegion1() <m_FixedImage1->GetBufferedRegion() <m_Filter1->GetOutput()->GetBufferedRegion() <GetFixedImage1()->GetBufferedRegion() <m_Filter1->GetOutput()->GetBufferedRegion() <GetBufferedRegion())) { + itkExceptionMacro(<< "FixedImageRegion1 does not overlap the fixed image buffered region"); + } + + if (!m_FixedImageRegion2.Crop(m_FixedImage2->GetBufferedRegion())) { + itkExceptionMacro(<< "FixedImageRegion2 does not overlap the fixed image buffered region"); + } + + this->SetTransformParameters(this->GetParameters()); + + if (m_ComputeGradient) { + + GradientImageFilterPointer gradientFilter = GradientImageFilterType::New(); + + gradientFilter->SetInput(m_MovingImage); + + const typename MovingImageType::SpacingType& spacing = m_MovingImage->GetSpacing(); + double maximumSpacing = 0.0; + for (unsigned int i = 0; i < MovingImageDimension; i++) { + if (spacing[i] > maximumSpacing) { + maximumSpacing = spacing[i]; + } + } + gradientFilter->SetSigma(maximumSpacing); + gradientFilter->SetNormalizeAcrossScale(true); + + gradientFilter->Update(); + + m_GradientImage = gradientFilter->GetOutput(); + } + + // If there are any observers on the metric, call them to give the + // user code a chance to set parameters on the metric + this->InvokeEvent(InitializeEvent()); +} + +template +unsigned int gTwoImageToOneImageMetric::GetNumberOfParameters() const +{ + if (!m_TransformMetaInfo) { + itkExceptionMacro(<< "TransformMetaInfo is missing"); + } + + switch (m_TransformMetaInfo->GetDegreeOfFreedom()) { + case X_ONLY: + case Y_ONLY: + case Z_ONLY: + return 1; + case THREE_DOF: + case ROTATION_ONLY: + return 3; + case SIX_DOF: + return 6; + default: + itkExceptionMacro(<< "Unsupported degree of freedeom"); + break; + } + + return m_IsocTransf->GetNumberOfParameters(); +} + +template +itk::Optimizer::ScalesType gTwoImageToOneImageMetric::GetWeightings() const +{ + itk::Optimizer::ScalesType weightings(this->GetNumberOfParameters()); + + switch (m_TransformMetaInfo->GetDegreeOfFreedom()) { + case X_ONLY: + case Y_ONLY: + case Z_ONLY: + weightings[0] = 1.0; + break; + ; + case THREE_DOF: + weightings[0] = 1.; + weightings[1] = 1.; + weightings[2] = 1.; + break; + case ROTATION_ONLY: + weightings[0] = 1. / dtr; + weightings[1] = 1. / dtr; + weightings[2] = 1. / dtr; + break; + case SIX_DOF: + //G: ITHINKTHOSEARE FLIPPED!! + weightings[0] = 1.; + weightings[1] = 1.; + weightings[2] = 1.; + weightings[3] = 1. / dtr; + weightings[4] = 1. / dtr; + weightings[5] = 1. / dtr; + break; + + default: + itkExceptionMacro(<< "Unsupported degree of freedeom"); + break; + } + return weightings; +} + +template +typename gTwoImageToOneImageMetric:: +ParametersType gTwoImageToOneImageMetric::GetParameters() const +{ + + // ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ + ParametersType parametersTransform = m_IsocTransf->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ + ParametersType parameters(this->GetNumberOfParameters()); + switch (m_TransformMetaInfo->GetDegreeOfFreedom()) { + case X_ONLY: + parameters[0] = parametersTransform[3]; + break; + case Y_ONLY: + parameters[0] = parametersTransform[4]; + break; + case Z_ONLY: + parameters[0] = parametersTransform[5]; + break; + case THREE_DOF: + parameters[0] = parametersTransform[3]; + parameters[1] = parametersTransform[4]; + parameters[2] = parametersTransform[5]; + break; + case ROTATION_ONLY: + parameters[0] = parametersTransform[0]; + parameters[1] = parametersTransform[1]; + parameters[2] = parametersTransform[2]; + break; + case SIX_DOF: + parameters[0] = parametersTransform[3]; + parameters[1] = parametersTransform[4]; + parameters[2] = parametersTransform[5]; + parameters[3] = parametersTransform[0]; + parameters[4] = parametersTransform[1]; + parameters[5] = parametersTransform[2]; + break; + + default: + itkExceptionMacro(<< "Unsupported degree of freedeom"); + break; + } + return parameters; +} + +template +void gTwoImageToOneImageMetric::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "ComputeGradient: " << static_cast::PrintType>(m_ComputeGradient) + << std::endl; + os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl; + os << indent << "Fixed Image 1: " << m_FixedImage1.GetPointer() << std::endl; + os << indent << "Fixed Image 2: " << m_FixedImage2.GetPointer() << std::endl; + os << indent << "Gradient Image: " << m_GradientImage.GetPointer() << std::endl; + os << indent << "Transform1: " << m_Transform1.GetPointer() << std::endl; + os << indent << "Transform2: " << m_Transform2.GetPointer() << std::endl; + os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl; + os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl; + os << indent << "Filter 1: " << m_Filter1.GetPointer() << std::endl; + os << indent << "Filter 2: " << m_Filter2.GetPointer() << std::endl; + os << indent << "Transform MetaInfoe: " << m_TransformMetaInfo.GetPointer() << std::endl; + os << indent << "FixedImageRegion 1: " << m_FixedImageRegion1 << std::endl; + os << indent << "FixedImageRegion 2: " << m_FixedImageRegion2 << std::endl; + os << indent << "Moving Image Mask: " << m_MovingImageMask.GetPointer() << std::endl; + os << indent << "Fixed Image Mask 1: " << m_FixedImageMask1.GetPointer() << std::endl; + os << indent << "Fixed Image Mask 2: " << m_FixedImageMask2.GetPointer() << std::endl; + os << indent << "Number of Pixels Counted: " << m_NumberOfPixelsCounted << std::endl; +} + +template +void gTwoImageToOneImageMetric::WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const +{ + + using RescaleFilterType = itk::RescaleIntensityImageFilter; + using WriterType = itk::ImageFileWriter; + + if (img) { + + // Rescale the intensity of the projection images to 0-255 for output. + typename RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New(); + rescaler1->SetOutputMinimum(0); + rescaler1->SetOutputMaximum(255); + rescaler1->SetInput(img); + + WriterType::Pointer writer1 = WriterType::New(); + + if (fileName.empty()) { + + char result[PATH_MAX]; + ssize_t count = readlink("/proc/self/exe", result, PATH_MAX); + if (count != -1) { + path = dirname(result); + path += "/img/out/"; + } + } + + if (fileName.empty()) { + + fileName = path; + fileName += "2D1.tif"; + } + writer1->SetFileName(fileName); + writer1->SetInput(rescaler1->GetOutput()); + + try { + std::cout << "Writing image: " << fileName << std::endl; + writer1->Update(); + } catch (itk::ExceptionObject& err) { + std::cerr << "ERROR: ExceptionObject caught !" << std::endl; + std::cerr << err << std::endl; + } + } +} + +} // end namespace itk + +#endif diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 5fb3ad1..35f6aad 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -1,4 +1,4 @@ - + /* @@ -14,6 +14,7 @@ gfattori 08.11.2021 #include "itkImageProcessor.h" #include "itkCreateObjectFunction.h" +#include "itkDRTHelpers.h" #include "itkVersion.h" #include #include @@ -24,6 +25,10 @@ gfattori 08.11.2021 #include "itkOrientImageFilter.h" #include +#include "itkCommand.h" + +#include "itkTimeProbesCollectorBase.h" + #include @@ -31,72 +36,23 @@ gfattori 08.11.2021 namespace itk { -/* this is in the end a IEC to HFS... - * but we keep this for ourselves... - */ -static double Standard_DRT2LPS [9] = { - 1,0,0, - 0,0,-1, - 0,1,0 }; - -static double PAElementsIEC[9] = { - 1, 0, 0, - 0, -1, 0, - 0, 0, -1 }; - -static double LATElementsIEC[9] = { - 0, 0, -1, - 0, -1, 0, - -1, 0, 0}; - -static double HFS2IEC[9] = { - 1, 0, 0, - 0, 0, 1, - 0, -1, 0}; - -static double FFS2IEC[9] = { - -1, 0, 0, - 0, 0, -1, - 0, -1, 0}; - - -//FUNCTION DECLARATION NECESSARY FOR COPY/PASTE FROM vtkGDCMImageReader -// const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds); - -const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds) -{ - static std::string buffer; - buffer = ""; // cleanup previous call - if( ds.FindDataElement( t ) ) - { - const gdcm::DataElement& de = ds.GetDataElement( t ); - const gdcm::ByteValue *bv = de.GetByteValue(); - if( bv ) // Can be Type 2 - { - buffer = std::string( bv->GetPointer(), bv->GetLength() ); - // Will be padded with at least one \0 - } - } - // Since return is a const char* the very first \0 will be considered - return buffer.c_str(); -}; - - itkImageProcessor::itkImageProcessor() { + iNWUnits = 1; + - // std::cout<<"itkImageProcessor::itkImageProcessor contructor"<SetIdentity(); transform2 = TransformType::New(); transform2->SetIdentity(); + m_InternalTransf1 = InternalTransformMetaInformation::New(); + m_InternalTransf2 = InternalTransformMetaInformation::New(); + + 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(); @@ -136,14 +92,14 @@ itkImageProcessor::itkImageProcessor() m_RTMetaInfo = NULL; - m_DRTGeometryMetaInfo = NULL; - m_TransformMetaInfo = NULL; - - m_TransformMetaInfo = TransformMetaInformation::New(); + m_r23MetaInfo = NULL; + m_r23MetaInfo = R23MetaInformation::New(); + /* Initialise the projection geoemtry with defaults */ + m_DRTGeometryMetaInfo = NULL; m_DRTGeometryMetaInfo = DRTProjectionGeometryImageMetaInformation::New(); m_DRTGeometryMetaInfo->SetSCD1(570.); @@ -157,7 +113,7 @@ itkImageProcessor::itkImageProcessor() Point3D[1]=0.; Point3D[2]=175.; m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); - m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); + // m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); ImageType3D::SizeType ImageSize; ImageSize[0]=512; @@ -181,41 +137,13 @@ itkImageProcessor::itkImageProcessor() m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D); - // // TEST MAPPING - // { + optimizerObserver = CommandIterationUpdate::New(); + ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New(); - // ImageType3D::PointType m_COR; - // ImageType3D::PointType m_TransformOrigin; // Origin of the user defined trasnform. likely isocenter. - // ImageType3D::PointType m_Translations; - // ImageType3D::PointType m_Rotations; + // TOOOOODOOOOO optimizerObserver->SetProcess(registration); - // TransformType::Pointer m_OutputTransform ; - - // m_COR[0] = 4.70067; - // m_COR[1] = -775.183; - // m_COR[2] = 202.321 ; - - // m_TransformOrigin.Fill(0.); - - // m_Translations[0] = 85.33068531; - // m_Translations[1] = 61.87811729; - // m_Translations[2] = -1.13774466 ; - - // m_Rotations[0] = -0.09535216; - // m_Rotations[1] = 0.15113885; - // m_Rotations[2] = -5.19765723; - - // m_OutputTransform= - // this->CalculateIsocentricTransform( - // m_COR, - // m_TransformOrigin, - // m_Translations, - // m_Rotations - // ); - - // m_OutputTransform->Print(std::cout); - - // }; + m_R23 = itkReg23::New(); + m_R23->SetOptimizerObserver(optimizerObserver); } @@ -225,6 +153,37 @@ itkImageProcessor::~itkImageProcessor() } +void itkImageProcessor::SetNumberOfWorkingUnits(int iN){ + if(iN>0 && iN<50){ + iNWUnits = iN; + } +} + + +const CTVolumeImageMetaInformation::Pointer +itkImageProcessor::GetCTMetaInfo( + ){ + + return m_CTMetaInfo; + +} + +const RTGeometryMetaInformation::Pointer +itkImageProcessor::GetRTMetaInfo( + ){ + + return m_RTMetaInfo; + +} + +const DRTProjectionGeometryImageMetaInformation::Pointer +itkImageProcessor::GetDRTGeoMetaInfo(){ + + return m_DRTGeometryMetaInfo; + +} + + void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os, indent); @@ -239,6 +198,7 @@ void itkImageProcessor::SetSCDOffset(double dOff1, double dOff2) { m_DRTGeometryMetaInfo->SetSCD1Offset(dOff1); m_DRTGeometryMetaInfo->SetSCD2Offset(dOff2); + } void itkImageProcessor::SetSCD(double dDist1, double dDist2) @@ -253,9 +213,45 @@ double itkImageProcessor::GetSCD1(){ } double itkImageProcessor::GetSCD2(){ return - m_DRTImage1MetaInfo ->GetSCD(); + m_DRTImage2MetaInfo ->GetSCD(); } +void itkImageProcessor::SetPanelOffsets(double dOff1, double dOff2) +{ + m_DRTGeometryMetaInfo->SetPanel1Offset(dOff1); + m_DRTGeometryMetaInfo->SetPanel2Offset(dOff2); +} + +double itkImageProcessor::GetPanelOffset1(){ + return + m_DRTGeometryMetaInfo->GetPanel1Offset(); +} + +double itkImageProcessor::GetPanelOffset2(){ + return + m_DRTGeometryMetaInfo ->GetPanel2Offset(); +} + + + +tDegreeOfFreedomEnum +itkImageProcessor::GetDegreeOfFreedom() +{ + return + m_r23MetaInfo->GetDegreeOfFreedom(); +} + +void itkImageProcessor::SetUseRotationsForImportOffset(bool bVal){ + + if(m_DRTGeometryMetaInfo == NULL){ + return; + } + m_DRTGeometryMetaInfo->SetUseRotationsForHiddenTransform(bVal); + + return; +} + + void itkImageProcessor::SetCustom_ImportTransform(double dTx, double dTy, double dTz, @@ -331,7 +327,6 @@ void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC( void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double nY2) { - if(m_DRTGeometryMetaInfo == NULL) { // todo } @@ -370,6 +365,33 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2) //TODO UPDATE TO FOLLOW } +void itkImageProcessor::SetCustom_UpdateMetaInfo(){ + + this->UpdateProjectionGeometryMeta(); +} + +int itkImageProcessor::unload3DVolumeAndMeta(){ + + std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; + + if (m_VolumeSourceDupli) { + m_VolumeSourceDupli = NULL; + m_VolumeSourceDupli = DuplicatorType::New(); + } + + m_CTMetaInfo = NULL; + + m_DRTImage1MetaInfo = NULL; + m_DRTImage2MetaInfo = NULL; + + m_TransformMetaInfo = NULL; + m_TransformMetaInfo = TransformMetaInformation::New(); + + + return 1; +} + + int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ tPatOrientation m_PatOrientation @@ -409,7 +431,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ imageSeriesReader3D->SetImageIO(dicomIO); imageSeriesReader3D->SetFileNames(v_sortedFnames); - //imageSeriesReader3D->SetNumberOfWorkUnits(8); + imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits); imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt //imageSeriesReader3D->SetSpacingWarningRelThreshold(); @@ -451,6 +473,15 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ caster3D->SetInput(imageSeriesReader3D->GetOutput()); caster3D->Update(); + if (m_VolumeSourceDupli == NULL) + std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl; + + // if (m_VolumeSourceDupli) { + // std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl; + // m_VolumeSourceDupli = NULL; + // m_VolumeSourceDupli = DuplicatorType::New(); + // } + m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput()); m_VolumeSourceDupli->Update(); @@ -460,137 +491,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ } catch (itk::ExceptionObject & ex) { - std::cout << ex << std::endl; - return EXIT_FAILURE; - } - - InternalImageType::Pointer m_InputImage = - m_VolumeSourceDupli->GetOutput(); - return - this->fill3DVolumeMeta(m_InputImage, m_PatOrientation); - -} - -int itkImageProcessor::load3DSerieFromFolder(const char * pcDirName) -{ - - using NamesGeneratorType = itk::GDCMSeriesFileNames; - NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New(); - - nameGenerator->SetUseSeriesDetails(true); - nameGenerator->AddSeriesRestriction("0008|0021"); - nameGenerator->SetGlobalWarningDisplay(false); - nameGenerator->SetDirectory(pcDirName); - - - tPatOrientation m_PatOrientation - = tPatOrientation::NotDefined; - - try - { - using SeriesIdContainer = std::vector; - const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs(); - auto seriesItr = seriesUID.begin(); - auto seriesEnd = seriesUID.end(); - - if (seriesItr != seriesEnd) - { - std::cout << "The directory: "; - std::cout << pcDirName << std::endl; - std::cout << "Contains the following DICOM Series: "; - std::cout << std::endl; - } - else - { - std::cout << "No DICOMs in: " << pcDirName << std::endl; - return EXIT_FAILURE; - } - - while (seriesItr != seriesEnd) - { - std::cout << seriesItr->c_str() << std::endl; - ++seriesItr; - } - - /* we expect only one serie at this point... */ - if(seriesUID.size() != 1){ - std::cout<<"More than one serie in the selected folder. returning..."<c_str(); - seriesItr++; - - std::cout << "\nReading: "; - std::cout << seriesIdentifier << std::endl; - using FileNamesContainer = std::vector; - FileNamesContainer fileNames = nameGenerator->GetFileNames(seriesIdentifier); - - using ImageIOType = itk::GDCMImageIO; - ImageIOType::Pointer dicomIO = ImageIOType::New(); - ImageSeriesReaderType3D::Pointer - imageSeriesReader3D = ImageSeriesReaderType3D::New(); - - imageSeriesReader3D->SetImageIO(dicomIO); - imageSeriesReader3D->SetFileNames(fileNames); - imageSeriesReader3D->SetNumberOfWorkUnits(8); - imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt - - imageSeriesReader3D->Update(); - - if(fileNames.size()>0){ - - gdcm::Reader R; - R.SetFileName(fileNames.at(0).c_str()); - if(!R.Read()) - { - cerr<<"ERROR: cannot read file: "<SetInput(imageSeriesReader3D->GetOutput()); - caster3D->Update(); - - m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput()); - m_VolumeSourceDupli->Update(); - - caster3D = NULL; - - } - - } - catch (itk::ExceptionObject & ex) - { - std::cout << ex << std::endl; + // std::cout << ex << std::endl; return EXIT_FAILURE; } @@ -608,19 +509,23 @@ int itkImageProcessor::fill3DVolumeMeta( if(m_CTMetaInfo != NULL){ - // TODO UNLOAD + std::cout << "NEVER HERE m_CTMetaInfo" << std::endl; + m_CTMetaInfo = NULL; } if(m_DRTImage1MetaInfo != NULL){ - // TODO UNLOAD + std::cout << "NEVER HERE m_DRTImage1MetaInfo" << std::endl; + m_DRTImage1MetaInfo = NULL; } if(m_DRTImage2MetaInfo != NULL){ - // TODO UNLOAD + std::cout << "NEVER HERE m_DRTImage2MetaInfo" << std::endl; + m_DRTImage2MetaInfo = NULL; } if(m_RTMetaInfo != NULL){ - // TODO UNLOAD + std::cout << "NEVER HERE m_RTMetaInfo" << std::endl; + m_RTMetaInfo = NULL; } @@ -668,22 +573,22 @@ int itkImageProcessor::fill3DVolumeMeta( m_DRTGeometryMetaInfo->GetSCD2Offset() ); - std::cout<< "///////////////// 3D VOLUME BEG ///////////////" <GetOriginLPS() <GetCOV() <GetOriginLPS() <GetCOV() <GetImportOffset() <GetLPS2IECDirections() <GetSize() <GetSpacing() <GetLPS2IECDirections() <GetSize() <GetSpacing() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <UpdateProjectionGeometryMeta(); @@ -696,6 +601,11 @@ int itkImageProcessor::fill3DVolumeMeta( pZeroOrigin[1] = 0.; pZeroOrigin[2] = 0.; + if (m_3DInputChangeInformationToZero) { + m_3DInputChangeInformationToZero = NULL; + m_3DInputChangeInformationToZero = ChangeInformationFilterType::New(); + } + m_3DInputChangeInformationToZero->SetInput( m_VolumeSourceDupli->GetOutput()); m_3DInputChangeInformationToZero->SetOutputOrigin( @@ -818,6 +728,8 @@ void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2) { m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1); m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2); + + std::cout << "SetProjectionAngleOffsetIEC " << dOff1 << " " << dOff2 << std::endl; } @@ -842,6 +754,38 @@ double itkImageProcessor::GetProjectionAngle2LPS(){ m_DRTImage2MetaInfo->GetProjectionAngleLPS(); } +int itkImageProcessor::unload2DAndMeta(int iImgType){ + + + switch (iImgType) { + case eProjectionOrientationType::PA: + + std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << " PA "<< std::endl; + m_TImage1MetaInfo = NULL; + + if (m_PASourceDupli) { + m_PASourceDupli = NULL; + m_PASourceDupli = DuplicatorType::New(); + } + break; + + case eProjectionOrientationType::LAT: + + std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << " LAT "<< std::endl; + m_TImage2MetaInfo = NULL; + + if (m_LATSourceDupli) { + m_LATSourceDupli = NULL; + m_LATSourceDupli = DuplicatorType::New(); + } + + break; + } + + return 1; + +} + int itkImageProcessor::load2D(const char * pcFName){ /* Check if we can open the file */ @@ -860,7 +804,7 @@ int itkImageProcessor::load2D(const char * pcFName){ strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0008,0x0008), ds)); std::string sLocalizerString = "LOCALIZER"; if (std::string(sTmpString).find(sLocalizerString) != std::string::npos) { - std::cout << "Loading Localizer: "<GetOutput()->GetDirection(); - std::cout<<"LocalizerImagDirectionDCM " <SetInputImage(caster2DInput->GetOutput() ); m_Duplicator->Update(); +// InternalImageType::IndexType pIndex; +// pIndex[0]=200; +// pIndex[1]=300; +// InternalImageType::PixelType pPixel = +// m_Duplicator->GetOutput()->GetPixel(pIndex); + //std::cout<<"processro PIXEL: "<GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<SetComputeZYX(true); - InputTransform->SetIdentity(); - - TransformType::OutputVectorType translation; - translation[0] = m_Translations[0]; - translation[1] = m_Translations[1]; - translation[2] = m_Translations[2]; - - InputTransform->SetTranslation(translation); - - const double dtr = (atan(1.0) * 4.0) / 180.0; - InputTransform->SetRotation( - dtr * m_Rotations[0], - dtr * m_Rotations[1], - dtr * m_Rotations[2]); - - ImageType3D::PointType m_TransformOrigin; - m_TransformOrigin.Fill(0.); - InputTransform->SetCenter( - m_TransformOrigin ); - - ImageType3D::PointType NewOriginTranslations = - InputTransform->TransformPoint(m_COR); - - ImageType3D::PointType DeltaNewOrigin = - NewOriginTranslations - m_COR; - TransformType::Pointer m_OutputTransform = - TransformType::New(); - m_OutputTransform ->SetComputeZYX(true); - m_OutputTransform ->SetIdentity(); - translation[0] = DeltaNewOrigin[0]; - translation[1] = DeltaNewOrigin[1]; - translation[2] = DeltaNewOrigin[2]; - m_OutputTransform->SetTranslation(translation); - m_OutputTransform->SetRotation( - dtr * m_Rotations[0], - dtr * m_Rotations[1], - dtr * m_Rotations[2]); +Optimizer::ParametersType +itkImageProcessor::GetFinalR23Parameters(){ - m_OutputTransform->SetCenter(m_COR); - InputTransform = NULL; + if (!m_TransformMetaInfo) { + itkExceptionMacro(<< "m_TransformMetaInfo not present"); + } - // m_OutputTransform.Print(std::cout); - return m_OutputTransform; -} + Optimizer::ParametersType pPars(6); + pPars.fill(0.); -itkImageProcessor::TransformType::Pointer -itkImageProcessor::CalculateInternalTransform( - ImageType3D::PointType m_TranslationOffset, //IEC - ImageType3D::PointType m_RotationOffset, //IEC - ImageType3D::PointType m_TranslationUser, //IEC - ImageType3D::PointType m_RotationUser, //IEC - ImageType3D::PointType m_ProjectionTransformCenter, //LPS - ImageType3D::PointType m_UserTransformCenter, //LPS - ImageType3D::PointType m_OffsetTransformCenter, //LPS - InternalImageType::DirectionType m_IECtoLPSDirections - ) -{ + itk::Optimizer::ParametersType currentPosition = + m_R23->GetCurrentPosition(); - //Convert all inputs into LPS + switch (m_r23MetaInfo->GetDegreeOfFreedom()) { + case THREE_DOF: + pPars[3] = currentPosition[0] + - m_TransformMetaInfo->GetHiddenTranslations()[0]; + pPars[4] = currentPosition[1] + - m_TransformMetaInfo->GetHiddenTranslations()[1]; + pPars[5] = currentPosition[2] + - m_TransformMetaInfo->GetHiddenTranslations()[2]; + break; - ImageType3D::PointType m_TOffsetLPS = - m_IECtoLPSDirections * m_TranslationOffset; + case SIX_DOF: + pPars[0] = currentPosition[0] + - m_TransformMetaInfo->GetHiddenRotations()[0]; + pPars[1] = currentPosition[1] + - m_TransformMetaInfo->GetHiddenRotations()[1]; + pPars[2] = currentPosition[2] + - m_TransformMetaInfo->GetHiddenRotations()[2]; + pPars[3] = currentPosition[3] + - m_TransformMetaInfo->GetHiddenTranslations()[0]; + pPars[4] = currentPosition[4] + - m_TransformMetaInfo->GetHiddenTranslations()[1]; + pPars[5] = currentPosition[5] + - m_TransformMetaInfo->GetHiddenTranslations()[2]; + break; - ImageType3D::PointType m_ROffsetLPS = - m_IECtoLPSDirections * m_RotationOffset; + default: + pPars.fill(0.); + break; + } - ImageType3D::PointType m_TUserLPS = - m_IECtoLPSDirections * m_TranslationUser; - - ImageType3D::PointType m_RUserLPS = - m_IECtoLPSDirections * m_RotationUser; - - // Map user to the projection center - TransformType::Pointer m_UserTransform = - MapTransformToNewOrigin ( - m_ProjectionTransformCenter - m_UserTransformCenter, - m_TUserLPS, - m_RUserLPS - ); - - // Map offset to the projection center - TransformType::Pointer m_OffsetTransform = - MapTransformToNewOrigin ( - m_ProjectionTransformCenter - m_UserTransformCenter,//m_OffsetTransformCenter, - m_TOffsetLPS, - m_ROffsetLPS - ); - - TransformType::Pointer m_OutputTransform = - TransformType::New(); - m_OutputTransform ->SetComputeZYX(true); - m_OutputTransform ->SetIdentity(); - - m_OutputTransform->SetTranslation( - m_OffsetTransform->GetTranslation() + - m_UserTransform->GetTranslation() - ); - - m_OutputTransform->SetRotation( - m_OffsetTransform->GetAngleX() + - m_UserTransform->GetAngleX(), - m_OffsetTransform->GetAngleY() + - m_UserTransform->GetAngleY(), - m_OffsetTransform->GetAngleZ() + - m_UserTransform->GetAngleZ() - ); - - m_OutputTransform->SetCenter( - m_ProjectionTransformCenter); return - m_OutputTransform; + pPars; + } +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() { - std::cout<< "itkImageProcessor::InitializeProjector()" <GetLPS2IECDirections().GetTranspose(); + // IECtoLPS_Directions = + // m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); TransformType::Pointer CurrTransform; - //CurrTransform->SetComputeZYX(true); - //CurrTransform->SetIdentity(); - if(m_RTMetaInfo == NULL) - { - - //AAAAAA TODOOOOO CERCA - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage1MetaInfo->GetProjectionOriginLPS(); - - //ImageType3D::PointType pFakeIsoLPS = - // m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter1()); - - - CurrTransform = - CalculateInternalTransform( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - pFakeIsoLPS, - ZeroPoint, - IECtoLPS_Directions - ); - } else { - - CurrTransform = - CalculateInternalTransform( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetImportOffset(), - IECtoLPS_Directions - ); - - } + CurrTransform= CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf1->GetpCalProjCenter(), + m_InternalTransf1->GetpRTIsocenter(), + m_InternalTransf1->GetIECtoLPSDirs()); transform1->SetComputeZYX(true); @@ -1331,10 +1217,9 @@ void itkImageProcessor::InitializeProjector() CurrTransform->GetAngleZ() ); - transform1->SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); - transform1->Print(std::cout); + // 2D Image 1 interpolator1->SetProjectionAngle( @@ -1343,55 +1228,18 @@ void itkImageProcessor::InitializeProjector() m_DRTImage1MetaInfo->GetSCD()); interpolator1->SetThreshold( m_DRTGeometryMetaInfo->GetIntensityThreshold() ); + interpolator1->SetPanelOffset( + m_DRTGeometryMetaInfo->GetPanel1Offset()); + interpolator1->SetTransform(transform1); interpolator1->Initialize(); - - - - - - - //CurrTransform->SetComputeZYX(true); - //CurrTransform->SetIdentity(); - - if(m_RTMetaInfo == NULL) - { - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage2MetaInfo->GetProjectionOriginLPS(); - - //ImageType3D::PointType pFakeIsoLPS = - // m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter2()); - - - - CurrTransform = - CalculateInternalTransform( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - pFakeIsoLPS, - ZeroPoint, - IECtoLPS_Directions - ); - } else { - - CurrTransform = - CalculateInternalTransform( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetImportOffset(), - IECtoLPS_Directions - ); - - } + CurrTransform= CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf2->GetpCalProjCenter(), + m_InternalTransf2->GetpRTIsocenter(), + m_InternalTransf2->GetIECtoLPSDirs()); transform2->SetComputeZYX(true); @@ -1417,37 +1265,18 @@ void itkImageProcessor::InitializeProjector() m_DRTImage2MetaInfo->GetSCD() ); interpolator2->SetThreshold( m_DRTGeometryMetaInfo->GetIntensityThreshold() ); + interpolator2->SetPanelOffset( + m_DRTGeometryMetaInfo->GetPanel2Offset()); + interpolator2->SetTransform(transform2); interpolator2->Initialize(); - - - - - - // // 2D Image 2 - // interpolator2->SetProjectionAngle( - // dtr * - // m_DRTImage2MetaInfo->GetProjectionAngleLPS() ); - // interpolator2->SetFocalPointToIsocenterDistance( - // m_DRTImage2MetaInfo->GetSCD()); - // interpolator2->SetThreshold( - // m_DRTGeometryMetaInfo->GetIntensityThreshold() - // ); - // interpolator2->SetTransform(transform); - // interpolator2->Initialize(); - - - - - resampleFilter1 = ResampleFilterType::New(); resampleFilter1->SetInput( image3DIn); resampleFilter1->SetDefaultPixelValue(0); - - resampleFilter1->SetNumberOfWorkUnits(8); + resampleFilter1->SetNumberOfWorkUnits(iNWUnits); // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance // have been set before registration. Here we only need to replace the initial @@ -1466,7 +1295,7 @@ void itkImageProcessor::InitializeProjector() resampleFilter2->SetInput( image3DIn); resampleFilter2->SetDefaultPixelValue(0); - resampleFilter2->SetNumberOfWorkUnits(8); + resampleFilter2->SetNumberOfWorkUnits(iNWUnits); // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance @@ -1496,31 +1325,23 @@ void itkImageProcessor::InitializeProjector() filter1->ChangeOriginOn(); filter1->UpdateOutputInformation(); - std::cout<< "itkImageProcessor::InitializeProjector() " <Update(); filter2->SetInput( resampleFilter2 ->GetOutput() ); filter2->SetOutputDirection(m_DRTImage2MetaInfo->GetImageDirectionsLPS() ); filter2->ChangeDirectionOn(); - filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() ); + filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() ); filter2->ChangeOriginOn(); filter2->UpdateOutputInformation(); - filter1->Update(); filter2->Update(); - - imageDRT1In = filter1->GetOutput(); imageDRT2In = filter2->GetOutput(); - - std::cout<< "itkImageProcessor::InitializeProjector() END" <SetImportOffset(pZero); + m_TransformMetaInfo->SetHiddenRotations(pZero); + + return 0; +} + + void itkImageProcessor::loadRTPlanInfo( double dIsoX, double dIsoY, double dIsoZ, double dLAT, double dVRT ,double dLNG){ if(m_RTMetaInfo != NULL){ + std::cout << " NEVER HERE loadRTPlanInfo m_RTMetaInfo" << std::endl; //TODO } if(m_CTMetaInfo != NULL){ + std::cout << " ALWAYS HERE loadRTPlanInfo m_CTMetaInfo" << std::endl; //TODO } if(m_DRTGeometryMetaInfo != NULL){ + std::cout << " ALWAYS HERE loadRTPlanInfo m_DRTGeometryMetaInfo" << std::endl; //TODO } @@ -1574,23 +1413,32 @@ void itkImageProcessor::loadRTPlanInfo( m_RTMetaInfo->SetIsocenterLPS(Point); Point[0] = dLAT; - Point[1] = dVRT; - Point[2] = dLNG; + Point[1] = dLNG; + Point[2] = dVRT; m_RTMetaInfo->SetIsocenterIECS(Point); + + ImageType3D::PointType + pZero (3); + pZero.Fill(0.); + m_CTMetaInfo->SetImportOffset( CalcImportVolumeOffset( m_RTMetaInfo->GetIsocenterIECS(), m_CTMetaInfo->GetLPS2IECDirections(), m_RTMetaInfo->GetIsocenterLPS(), m_DRTGeometryMetaInfo->GetIECS2IECScannerT(), - m_DRTGeometryMetaInfo->GetIECS2IECScannerR() + (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ? + m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero) ) ); - std::cout<< "m_CTMetaInfo->GetImportOffset() " - <GetImportOffset() <SetHiddenRotations( + (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ? + m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero) + ); + this->UpdateProjectionGeometryMeta(); this->InitializeProjector(); @@ -1600,18 +1448,22 @@ void itkImageProcessor::loadRTPlanInfo( void itkImageProcessor::UpdateProjectionGeometryMeta(){ if(m_CTMetaInfo == NULL){ + return; //TODO. } if(m_DRTGeometryMetaInfo == NULL){ + return; //TODO. } if(m_DRTImage1MetaInfo == NULL){ + return; //TODO. } if(m_DRTImage2MetaInfo == NULL){ + return; //TODO. } @@ -1632,20 +1484,17 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTGeometryMetaInfo->GetProjectionCenter() ); - std::cout<< "///////////////// PGEOM META BEG ///////////////" <GetProjectionCenter() <ConvertIECPointToLPSPoint( m_DRTGeometryMetaInfo->GetProjectionCenterOffset1()); - + // std::cout<< "///////////////// PGEOM META BEG ///////////////" <GetProjectionCenter() <SetProjectionAngleLPS( + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() + + m_DRTGeometryMetaInfo->GetProjectionAngle1OffsetIEC()) + ); - + m_DRTImage1MetaInfo->SetSCD( + m_DRTGeometryMetaInfo->GetSCD1() + + m_DRTGeometryMetaInfo->GetSCD1Offset() + ); m_DRTImage1MetaInfo->SetProjectionOriginLPS( CalibratedIsocenterLPS); @@ -1675,26 +1535,19 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTImage1MetaInfo->SetProjectionOriginLPSZero( CalibratedIsocenterZeroLPS); - // This is based on the calibrated iso to be sure... - std::vector vBounds = m_CTMetaInfo->CalculateRegions( - NominalIsocenter_wZcorrectionLPS//m_DRTImage1MetaInfo->GetProjectionOriginLPS() - ); - - std::cout<<"Bounds1: " - <SetSizeWithBounds( - vBounds.data(), + NULL,//vBounds.data(), m_DRTGeometryMetaInfo->GetDRT1Size(), m_DRTGeometryMetaInfo->GetDRT1Spacing() ); + ImageType3D::PointType PanelOffsetIEC1; + + PanelOffsetIEC1[0] = -m_DRTGeometryMetaInfo->GetPanel1Offset(); + PanelOffsetIEC1[1] = 0.0; + PanelOffsetIEC1[2] = 0.0; + // This HAS TO be calculated from the nominal isocenter. m_DRTImage1MetaInfo->SetOriginLPS( CalcDRTImageOrigin( @@ -1704,32 +1557,41 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),// m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) + ); - - m_DRTImage1MetaInfo->SetImageDirectionsLPS( CalcDRTImageDirections( this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),//m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); - std::cout<<"IsocenterOffsetLPS"<< IsocenterOffsetLPS <GetSize()"<< m_DRTImage1MetaInfo->GetSize() <GetSpacing()"<< m_DRTImage1MetaInfo->GetSpacing() <GetOrigin()"<< m_DRTImage1MetaInfo->GetOrigin() <GetCOV()"<< m_DRTImage1MetaInfo->GetCOV() <GetOriginLPS()"<< m_DRTImage1MetaInfo->GetOriginLPS() <GetProjectionAngleLPS()"<GetProjectionAngleLPS()<GetImageDirectionsLPS()"<< m_DRTImage1MetaInfo->GetImageDirectionsLPS() <SetpCalProjCenter( + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); + m_InternalTransf1->SetpRTIsocenter( + m_DRTImage1MetaInfo->GetProjectionOriginLPS()); + InternalImageType::DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions); + + } else { + m_InternalTransf1->SetpCalProjCenter( + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); + m_InternalTransf1->SetpRTIsocenter( + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS()); + InternalImageType::DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions); + + } // SECOND @@ -1753,39 +1615,52 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ IsocenterOffsetLPS[2]; + m_DRTImage2MetaInfo->SetProjectionAngleLPS( + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC() + + m_DRTGeometryMetaInfo->GetProjectionAngle2OffsetIEC()) + ); + + m_DRTImage2MetaInfo->SetSCD( + m_DRTGeometryMetaInfo->GetSCD2()+ + m_DRTGeometryMetaInfo->GetSCD2Offset() + ); m_DRTImage2MetaInfo->SetProjectionOriginLPS( CalibratedIsocenterLPS); - // m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter2()) - // ); m_DRTImage2MetaInfo->SetProjectionOriginLPSZero( CalibratedIsocenterZeroLPS); - // m_CTMetaInfo->GetProjectionOriginLPSZero( - // m_DRTGeometryMetaInfo->GetProjectionCenter2()) - // ); - vBounds.clear(); - vBounds = m_CTMetaInfo->CalculateRegions( - NominalIsocenter_wZcorrectionLPS//m_DRTImage2MetaInfo->GetProjectionOriginLPS() - ); - - std::cout<<"Bounds2: " - <SetSizeWithBounds( - vBounds.data(), + NULL,//vBounds.data(), m_DRTGeometryMetaInfo->GetDRT2Size(), m_DRTGeometryMetaInfo->GetDRT2Spacing() ); + + ImageType3D::PointType PanelOffsetIEC2; + + PanelOffsetIEC2[0] = -m_DRTGeometryMetaInfo->GetPanel2Offset(); + PanelOffsetIEC2[1] = 0.0; + PanelOffsetIEC2[2] = 0.0; + + + m_DRTImage2MetaInfo->SetOriginLPS( + CalcDRTImageOrigin( + m_DRTImage2MetaInfo->GetOrigin(), + m_DRTImage2MetaInfo->GetCOV(), + NominalIsocenter_wZcorrectionLPS,//m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , + Std_DRT2LPS + )); + + m_DRTImage2MetaInfo->SetOriginLPS( CalcDRTImageOrigin( m_DRTImage2MetaInfo->GetOrigin(), @@ -1794,9 +1669,9 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) + ); m_DRTImage2MetaInfo->SetImageDirectionsLPS( @@ -1804,22 +1679,31 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); - std::cout<<"IsocenterOffsetLPS"<< IsocenterOffsetLPS <GetSize()"<< m_DRTImage2MetaInfo->GetSize() <GetSpacing()"<< m_DRTImage2MetaInfo->GetSpacing() <GetOrigin()"<< m_DRTImage2MetaInfo->GetOrigin() <GetCOV()"<< m_DRTImage2MetaInfo->GetCOV() <GetOriginLPS()"<< m_DRTImage2MetaInfo->GetOriginLPS() <GetProjectionAngleLPS()"<GetProjectionAngleLPS()<GetImageDirectionsLPS()"<< m_DRTImage2MetaInfo->GetImageDirectionsLPS() <SetpCalProjCenter( + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero()); + m_InternalTransf2->SetpRTIsocenter( + m_DRTImage2MetaInfo->GetProjectionOriginLPS()); + InternalImageType::DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions); + } else { + m_InternalTransf2->SetpCalProjCenter( + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero()); + m_InternalTransf2->SetpRTIsocenter( + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS()); + InternalImageType::DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions); + } } @@ -1850,17 +1734,20 @@ itkImageProcessor::CalcDRTImageOrigin( } + void itkImageProcessor::GetProjectionImages(){ - std::cout<< "itkImageProcessor::GetProjectionImages" <SetComputeZYX(true); transform1->SetIdentity(); - InternalImageType::DirectionType IECtoLPS_Directions; - - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); TransformType::Pointer CurrTransform; - //CurrTransform->SetComputeZYX(true); - //CurrTransform->SetIdentity(); - - if(m_RTMetaInfo == NULL) - { - - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage1MetaInfo->GetProjectionOriginLPS(); - //m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter()); - - std::cout<<"pFakeIsoLPS: "<GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - pFakeIsoLPS, - ZeroPoint, - IECtoLPS_Directions - ); - } else { - - CurrTransform = - CalculateInternalTransform( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetImportOffset(), - IECtoLPS_Directions - ); - - } - + CurrTransform= CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf1->GetpCalProjCenter(), + m_InternalTransf1->GetpRTIsocenter(), + m_InternalTransf1->GetIECtoLPSDirs()); transform1->SetComputeZYX(true); @@ -1926,21 +1775,18 @@ void itkImageProcessor::GetProjectionImages(){ CurrTransform->GetAngleY(), CurrTransform->GetAngleZ() ); - - std::cout<< "itkImageProcessor::GetProjectionImages" <SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); - std::cout<<"Angle1 "<GetProjectionAngleLPS() <GetProjectionOriginLPSZero() <GetProjectionOriginLPS()<GetProjectionAngleLPS() <GetProjectionOriginLPSZero() <GetProjectionOriginLPS()<GetAngleX() <<" "<< +// CurrTransform->GetAngleY() <<" "<< +// CurrTransform->GetAngleZ() <<" "<< std::endl; + +// std::cout<< "CurrTransform Translations: "<< +// CurrTransform->GetTranslation()<GetCenter()<Print(std::cout); @@ -1951,55 +1797,14 @@ void itkImageProcessor::GetProjectionImages(){ interpolator1->Initialize(); - - - transform2->SetComputeZYX(true); transform2->SetIdentity(); - - //CurrTransform->SetComputeZYX(true); - //CurrTransform->SetIdentity(); - - if(m_RTMetaInfo == NULL) - { - - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage2MetaInfo->GetProjectionOriginLPS(); - //m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter()); - - std::cout<<"pFakeIsoLPS: "<GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - pFakeIsoLPS, - ZeroPoint, - IECtoLPS_Directions - ); - - std::cout<<"pFakeIsoLPS: "<GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetImportOffset(), - IECtoLPS_Directions - ); - - } - + CurrTransform = CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf2->GetpCalProjCenter(), + m_InternalTransf2->GetpRTIsocenter(), + m_InternalTransf2->GetIECtoLPSDirs()); transform2->SetComputeZYX(true); @@ -2013,34 +1818,21 @@ void itkImageProcessor::GetProjectionImages(){ CurrTransform->GetAngleZ() ); - transform2 ->SetCenter( m_DRTImage2MetaInfo->GetProjectionOriginLPSZero()); - //transform2 ->Print(std::cout); - - - // // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance // // have been set before registration. Here we only need to replace the initial // // transform with the final transform. interpolator2->SetTransform(transform2); //finalTransform); interpolator2->Initialize(); - std::cout<<"pFakeIsoLPS: "<Update(); filter2->Update(); imageDRT1In = filter1->GetOutput(); imageDRT2In = filter2->GetOutput(); - std::cout<< "itkImageProcessor::GetProjectionImages END" <GetInverseTransform()->GetMatrix()<GetInverseTransform()->GetTranslation()<; RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New(); rescaler1->SetOutputMinimum(0); - rescaler1->SetOutputMaximum(255); + rescaler1->SetOutputMaximum(32768); rescaler1->SetInput(m_PASourceDupli->GetOutput()); rescaler1->Update(); @@ -2173,44 +1968,45 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK() toVTKLocalizer1->SetInput(rescaler1->GetOutput()); toVTKLocalizer1->Update(); - if(false) { - using ImageRegionType3D = ImageType3D::RegionType; - using SizeType3D = ImageRegionType3D::SizeType; - using ImageDirectionType3D = ImageType3D::DirectionType; +// if(false) { +//// using ImageRegionType3D = ImageType3D::RegionType; +//// using SizeType3D = ImageRegionType3D::SizeType; +//// using ImageDirectionType3D = ImageType3D::DirectionType; - ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion(); - SizeType3D size3D = region3D.GetSize(); - ImageType3D::PointType origin3D = rescaler1->GetOutput()->GetOrigin(); - const itk::Vector resolution3D = rescaler1->GetOutput()->GetSpacing(); - ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection(); +//// ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion(); +//// SizeType3D size3D = region3D.GetSize(); +//// //ImageType3D::PointType origin3D = rescaler1->GetOutput()->GetOrigin(); +//// const itk::Vector resolution3D = rescaler1->GetOutput()->GetSpacing(); +//// ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection(); - /* calculate image size in 3D Space */ - using VectorType = itk::Vector; - VectorType Dest; - Dest[0]=(size3D[0]-1) * resolution3D[0]; - Dest[1]=(size3D[1]-1) * resolution3D[1]; - Dest[2]=(size3D[2]-1) * resolution3D[2]; +//// /* calculate image size in 3D Space */ +//// using VectorType = itk::Vector; +//// VectorType Dest; +//// Dest[0]=(size3D[0]-1) * resolution3D[0]; +//// Dest[1]=(size3D[1]-1) * resolution3D[1]; +//// Dest[2]=(size3D[2]-1) * resolution3D[2]; - ImageType3D::PointType LastVoxelPosition = - origin3D + (imagDirection * Dest); - std::cout<<" -------- Localizer 1 ITK --------"<GetOutput()->GetBounds(); +// // double* dBounds = toVTKLocalizer1->GetOutput()->GetBounds(); - std::cout<< "-------- Localizer 1 VTK --------" <GetOutput(); @@ -2220,12 +2016,12 @@ vtkImageData* itkImageProcessor::GetLocalizer2VTK() { - // Rescale the intensity of the projection images to 0-255 for output. + // Rescale the intensity of the projection images to 0-32768 for output. using RescaleFilterType = itk::RescaleIntensityImageFilter; RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New(); rescaler2->SetOutputMinimum(0); - rescaler2->SetOutputMaximum(255); + rescaler2->SetOutputMaximum(32768); rescaler2->SetInput(m_LATSourceDupli ->GetOutput()); rescaler2->Update(); @@ -2235,42 +2031,42 @@ vtkImageData* itkImageProcessor::GetLocalizer2VTK() - if(true) { - using ImageRegionType3D = ImageType3D::RegionType; - using SizeType3D = ImageRegionType3D::SizeType; - using ImageDirectionType3D = ImageType3D::DirectionType; +// if(true) { +// using ImageRegionType3D = ImageType3D::RegionType; +// using SizeType3D = ImageRegionType3D::SizeType; +// using ImageDirectionType3D = ImageType3D::DirectionType; - ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion(); - SizeType3D size3D = region3D.GetSize(); - ImageType3D::PointType origin3D = rescaler2->GetOutput()->GetOrigin(); - const itk::Vector resolution3D = rescaler2->GetOutput()->GetSpacing(); - ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection(); +// ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion(); +// SizeType3D size3D = region3D.GetSize(); +// ImageType3D::PointType origin3D = rescaler2->GetOutput()->GetOrigin(); +// const itk::Vector resolution3D = rescaler2->GetOutput()->GetSpacing(); +// ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection(); - /* calculate image size in 3D Space */ - using VectorType = itk::Vector; - VectorType Dest; - Dest[0]=(size3D[0]-1) * resolution3D[0]; - Dest[1]=(size3D[1]-1) * resolution3D[1]; - Dest[2]=(size3D[2]-1) * resolution3D[2]; +// /* calculate image size in 3D Space */ +// using VectorType = itk::Vector; +// VectorType Dest; +// Dest[0]=(size3D[0]-1) * resolution3D[0]; +// Dest[1]=(size3D[1]-1) * resolution3D[1]; +// Dest[2]=(size3D[2]-1) * resolution3D[2]; - ImageType3D::PointType LastVoxelPosition = - origin3D + (imagDirection * Dest); +// ImageType3D::PointType LastVoxelPosition = +// origin3D + (imagDirection * Dest); - // std::cout<<" -------- Localizer 2 ITK --------"<GetOutput()->GetBounds(); - // std::cout<< "-------- Localizer 2 VTK --------" <GetOutput()->GetBounds(); +// // std::cout<< "-------- Localizer 2 VTK --------" <; RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New(); rescaler1->SetOutputMinimum(0); - rescaler1->SetOutputMaximum(255); + rescaler1->SetOutputMaximum(32768); rescaler1->SetInput( imageDRT1In ); - rescaler1->Update(); + // using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator; + // auto imageCalculatorFilter = ImageCalculatorFilterType::New(); + // imageCalculatorFilter->SetImage(imageDRT1In); + // imageCalculatorFilter->Compute(); + + // std::cout<< "itkImageProcessor::imageDRT1In() " << + // imageCalculatorFilter <; + // auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New(); + // imageCalculatorFilter2->SetImage(rescaler1->GetOutput()); + // imageCalculatorFilter2->Compute(); + // std::cout<< "itkImageProcessor::imageDRT2In() " << + // imageCalculatorFilter2 <SetInput(rescaler1->GetOutput()); toVTK2D1->Update(); @@ -2337,11 +2155,6 @@ vtkImageData* itkImageProcessor::GetProjection1VTK() vtkMatrix4x4 * itkImageProcessor::GetProjection1VTKTransform() { - // std::cout<< "Composed " <GetComposedTransform()->Print(std::cout); - // std::cout<< "Inverse " <GetInverseTransform()->Print(std::cout); - m_Projection1VTKTransform->Identity(); for(int iIdx = 0; iIdx<3 ; iIdx++){ @@ -2395,12 +2208,19 @@ vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform() vtkImageData* itkImageProcessor::GetProjection2VTK() { - // Rescale the intensity of the projection images to 0-255 for output. + if(m_DRTImage2MetaInfo == NULL || + m_DRTGeometryMetaInfo == NULL || + m_TransformMetaInfo == NULL ){ + return NULL; + //TODO + } + + // Rescale the intensity of the projection images to 0-32768 for output. using RescaleFilterType = itk::RescaleIntensityImageFilter; RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New(); rescaler2->SetOutputMinimum(0); - rescaler2->SetOutputMaximum(255); + rescaler2->SetOutputMaximum(32768); rescaler2->SetInput( imageDRT2In ); rescaler2->Update(); @@ -2439,13 +2259,12 @@ vtkImageData* itkImageProcessor::GetProjection2VTK() - // double* dBounds = toVTK2D2->GetOutput()->GetBounds(); - - // std::cout<< "-------- Proj 2 --------" <GetOutput()->GetBounds(); + // std::cout<< "-------- Proj 2 --------" <GetOutput(); @@ -2479,7 +2298,7 @@ void itkImageProcessor::WriteProjectionImages() try { - std::cout << "Writing image 1 " << std::endl; + // std::cout << "Writing image 1 " << std::endl; writer1->Update(); } catch (itk::ExceptionObject & err) @@ -2493,7 +2312,7 @@ void itkImageProcessor::WriteProjectionImages() try { - std::cout << "Writing image 2" << std::endl; + // std::cout << "Writing image 2" << std::endl; writer2->Update(); } catch (itk::ExceptionObject & err) @@ -2515,13 +2334,12 @@ void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ) //todo } - ImageType3D::PointType Translations; Translations[0]= dX; Translations[1]= dY; Translations[2]= dZ; - m_TransformMetaInfo->SetTranslations(Translations); + m_TransformMetaInfo->SetUserTranslations(Translations); } void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) @@ -2531,10 +2349,167 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) Rotations[1]= dY; Rotations[2]= dZ; - m_TransformMetaInfo->SetRotations(Rotations); + m_TransformMetaInfo->SetUserRotations(Rotations); + + +} + +Optimizer::ParametersType +itkImageProcessor::GetUserIsocentricTransform(){ + + Optimizer::ParametersType + m_Pars(6); + m_Pars.fill(0.); + + if(!m_TransformMetaInfo){ + return m_Pars; + } + + m_Pars[0] = m_TransformMetaInfo->GetUserRotations()[0]; + m_Pars[1] = m_TransformMetaInfo->GetUserRotations()[1]; + m_Pars[2] = m_TransformMetaInfo->GetUserRotations()[2]; + + m_Pars[3] = m_TransformMetaInfo->GetUserTranslations()[0]; + m_Pars[4] = m_TransformMetaInfo->GetUserTranslations()[1]; + m_Pars[5] = m_TransformMetaInfo->GetUserTranslations()[2]; + + return + m_Pars; } + +TransformType::Pointer +itkImageProcessor::GetCompleteIsocentricTransform +(/*ImageType3D::PointType TransformCenter*/){ + + TransformType::Pointer + regTransform = TransformType::New(); + regTransform->SetComputeZYX(true); + regTransform->SetIdentity(); + + if(m_TransformMetaInfo == nullptr || + m_CTMetaInfo == nullptr){ + return nullptr; + } + + ImageType3D::PointType mISORTIEC = + m_CTMetaInfo->GetLPS2IECDirections() * m_RTMetaInfo->GetIsocenterLPS(); + + mISORTIEC[0] *= -1; + mISORTIEC[1] *= -1; + mISORTIEC[2] *= -1; + + TransformType::Pointer mTransf= + MapTransformToNewOrigin( + mISORTIEC, + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR() + ); + + itk::TransformMetaInformation::PointType + pTM = m_CTMetaInfo->GetLPS2IECDirections() * m_CTMetaInfo->GetImportOffset(); + + TransformType::OutputVectorType pTransl; + + pTransl[0] = mTransf->GetTranslation()[0] - pTM[0]; + pTransl[1] = mTransf->GetTranslation()[1] - pTM[1]; + pTransl[2] = mTransf->GetTranslation()[2] - pTM[2]; + + mTransf->SetTranslation(pTransl); + itk::PointpZero; + pZero.Fill(0.); + mTransf->SetCenter(pZero); + + return mTransf; +} + + + +void itkImageProcessor::SetRegionFixed1( + int iIdx0, int iIdx1, int iSz0, int iSz1){ + + auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion(); + + auto index1 = region1temp.GetIndex(); + auto size1 = region1temp.GetSize(); + + index1[0] = iIdx0; + index1[1] = iIdx1; + + size1[0] = iSz0; + size1[1] = iSz1; + size1[2] = 1; + + roiAutoReg1.SetIndex(index1); + roiAutoReg1.SetSize(size1); + + //std::cout << "itkImageProcessor " << std::endl; + //std::cout << roiAutoReg1 << std::endl; + + this->m_R23->SetroiAutoReg1(roiAutoReg1); + +} + +void itkImageProcessor::SetRegionFixed2( + int iIdx0, int iIdx1, int iSz0, int iSz1){ + + auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion(); + + auto index2 = region2temp.GetIndex(); + auto size2 = region2temp.GetSize(); + + index2[0] = iIdx0; + index2[1] = iIdx1; + + size2[0] = iSz0; + size2[1] = iSz1; + size2[2] = 1; + + roiAutoReg2.SetIndex(index2); + roiAutoReg2.SetSize(size2); + + //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 2be4ca8..b4f2400 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -38,20 +38,18 @@ gfattori 08.11.2021 #include "gdcmGlobal.h" - #include "vtkImageData.h" #include "vtkTransform.h" #include "vtkSmartPointer.h" +#include "itkQtIterationUpdate.h" + +#include "itkImageProcessorHelpers.h" + +#include "itkReg23.h" namespace itk { -/* reference string required for comparison with tag values */ -static const char *ImageOrientationStrings[] = { - "NotDefined", - "HFS", - "FFS" -}; class ITK_EXPORT itkImageProcessor : public itk::Object { @@ -71,13 +69,50 @@ public: /** Run-time type information (and related methods). */ itkTypeMacro(itkImageProcessor, Object); + + 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); + /** Input data load methods*/ - int load3DSerieFromFolder(const char* ); int load3DSerieFromFiles( std::vector ); + int unload3DVolumeAndMeta(); int load2D(const char *); + int unload2DAndMeta(int); + //int query2DimageReconstructionDiameter(const char*); void loadRTPlanInfo(double, double, double, double, double ,double); + int unloadRTPlanAndMeta(); /** Projection angles - Gantry angle IEC */ void SetProjectionAngle1IEC(double); @@ -96,8 +131,14 @@ public: double GetSCD1(); double GetSCD2(); - void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double); + /** Panel offsets - panel offsets IEC */ + void SetPanelOffsets(double,double); + /** Get projection angles - Gantry angle LPS + * Only meaningful after a 3D Volume is loaded */ + double GetPanelOffset1(); + double GetPanelOffset2(); + void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double); void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); /** Intensity threshold used for image projection, @@ -107,14 +148,49 @@ public: /** Custom settings of the projection images */ void SetCustom_2Dres(double,double,double,double); void SetCustom_2Dsize(int,int,int,int); + + /** Custom transform to map IEC of imaging device into IEC-Support + * Tx Ty Tz [mm] Rx Ry Rz [deg] + */ void SetCustom_ImportTransform(double, double, double, double, double, double); + /** Should rotations be used when computing import offset and hidden transform? */ + void SetUseRotationsForImportOffset(bool bVal); + + void SetCustom_UpdateMetaInfo(); + /** Set transform parameters for 3D Volume */ void SetInitialTranslations(double, double, double); void SetInitialRotations(double, double, double); + ///** Get transform parameters for 3D Volume */ + /** Get the complete transform, likely for export to SRO. + * This includes user and hidden contributions. + * Transform center has to be specified, Zero would result + * in Isocentric transform with center at the RT Iso. + */ + TransformType::Pointer + GetCompleteIsocentricTransform(/*ImageType3D::PointType TransformCenter*/); + + /** Get only the User component of the transform + * as parameters' list. Import offset and hidden not included. + * This can be use to update UI transform display + * ParametersType :: Rx Ry Rz Tx Ty Tz + */ + Optimizer::ParametersType + GetUserIsocentricTransform(); + + /** Return useful meta info for external use + * Consumer should check for nullptr... + */ + const CTVolumeImageMetaInformation::Pointer + GetCTMetaInfo(); + const RTGeometryMetaInformation::Pointer + GetRTMetaInfo(); + const DRTProjectionGeometryImageMetaInformation::Pointer + GetDRTGeoMetaInfo(); /** Initialize projection geometry */ void InitializeProjector(); @@ -125,8 +201,6 @@ public: /** Get Projection origin in LPS coordinates*/ const std::vector GetRTImportOffset(); - - /** Get the Localizer intensity window and * center for rendering */ double GetLocalizerDisplayWindowLevel(int ); @@ -150,13 +224,12 @@ public: void Write2DImages(); - protected: /** Various pixel types */ using InternalPixelType = float; using PixelType2D = short; using PixelType3D = short; - using OutputPixelType = unsigned char; + using OutputPixelType = unsigned short; using ImageType3D = itk::Image; using InternalImageType = itk::Image; @@ -171,7 +244,7 @@ private: /** Fill Meta after 3D volume load */ int fill3DVolumeMeta(InternalImageType::Pointer, - tPatOrientation); + tPatOrientation); /** Image types */ using ImageType2D = itk::Image; @@ -181,11 +254,12 @@ private: /** The following lines define each of the components used in the registration: The transform, optimizer, metric, interpolator and the registration method itself. */ - using TransformType = itk::Euler3DTransform; using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; using ResampleFilterType = itk::ResampleImageFilter; + using RoiForRegistration = itk::ImageRegion; + /** Image reader types */ using ImageReaderType2D = itk::ImageFileReader; using ImageReaderType3D = itk::ImageFileReader; @@ -202,23 +276,10 @@ private: /** ITK to VTK filter */ using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; - - TransformType::Pointer - CalculateInternalTransform( - ImageType3D::PointType m_TranslationOffset, - ImageType3D::PointType m_RotationOffset, - ImageType3D::PointType m_TranslationUser, - ImageType3D::PointType m_RotationUser, - ImageType3D::PointType m_ProjectionTransformCenter, - ImageType3D::PointType m_UserTransformCenter, - ImageType3D::PointType m_OffsetTransformCenter, - InternalImageType::DirectionType m_IECtoLPSDirections - ); - - TransformType::Pointer transform1, - transform2; + transform2, + IsocTransf; InterpolatorType::Pointer interpolator1, @@ -235,6 +296,8 @@ private: imageDRT1In, imageDRT2In; + + DuplicatorType::Pointer m_LATSourceDupli, m_PASourceDupli, @@ -265,13 +328,6 @@ private: ); - TransformType::Pointer - MapTransformToNewOrigin( - ImageType3D::PointType m_COR, - ImageType3D::PointType m_Translations, - ImageType3D::PointType m_Rotations - ); - double CalcProjectionAngleLPS( tPatOrientation pOrient, @@ -305,8 +361,6 @@ private: * m_Projection1VTKTransform, * m_Projection2VTKTransform; - - /** * Many meta containers */ @@ -330,6 +384,26 @@ private: TransformMetaInformation::Pointer m_TransformMetaInfo; + R23MetaInformation::Pointer + m_r23MetaInfo; + + InternalTransformMetaInformation::Pointer + m_InternalTransf1, + m_InternalTransf2; + + + RoiForRegistration + roiAutoReg1, + roiAutoReg2; + + + int + iNWUnits; + + + itk::itkReg23::Pointer + m_R23; + }; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp new file mode 100644 index 0000000..bb93324 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp @@ -0,0 +1,165 @@ +#include "itkImageProcessorHelpers.h" +#include +#include + +namespace itk { + + +TransformType::Pointer +MapTransformToNewOrigin( + ImageType3D::PointType m_COR, // Center of rotation for the proj geometry. this is my new origin. + ImageType3D::PointType m_Translations, + ImageType3D::PointType m_Rotations + ){ + + TransformType::Pointer InputTransform = TransformType::New(); + InputTransform->SetComputeZYX(true); + InputTransform->SetIdentity(); + + TransformType::OutputVectorType translation; + translation[0] = m_Translations[0]; + translation[1] = m_Translations[1]; + translation[2] = m_Translations[2]; + + InputTransform->SetTranslation(translation); + + const double dtr = (atan(1.0) * 4.0) / 180.0; + InputTransform->SetRotation( + dtr * m_Rotations[0], + dtr * m_Rotations[1], + dtr * m_Rotations[2]); + + ImageType3D::PointType m_TransformOrigin; + m_TransformOrigin.Fill(0.); + InputTransform->SetCenter( + m_TransformOrigin ); + + ImageType3D::PointType NewOriginTranslations = + InputTransform->TransformPoint(m_COR); + + ImageType3D::PointType DeltaNewOrigin = + NewOriginTranslations - m_COR; + + + TransformType::Pointer m_OutputTransform = + TransformType::New(); + m_OutputTransform ->SetComputeZYX(true); + m_OutputTransform ->SetIdentity(); + + translation[0] = DeltaNewOrigin[0]; + translation[1] = DeltaNewOrigin[1]; + translation[2] = DeltaNewOrigin[2]; + + m_OutputTransform->SetTranslation(translation); + m_OutputTransform->SetRotation( + dtr * m_Rotations[0], + dtr * m_Rotations[1], + dtr * m_Rotations[2]); + + m_OutputTransform->SetCenter(m_COR); + + InputTransform = NULL; + + return m_OutputTransform; +} + + + +TransformType::Pointer +CalculateInternalTransformV3( + ImageType3D::PointType m_Translation, //IEC + ImageType3D::PointType m_Rotation, //IEC + ImageType3D::PointType m_CalibratedProjectionCenter, //LPS + ImageType3D::PointType m_RTIsocenter, //LPS + InternalImageType::DirectionType m_IECtoLPSDirections + ) +{ + + //Convert all inputs into LPS + + ImageType3D::PointType m_TLPS = + m_IECtoLPSDirections * m_Translation; + + ImageType3D::PointType m_RLPS = + m_IECtoLPSDirections * m_Rotation; + + // Map offset to the projection center + TransformType::Pointer m_outputTransform = + MapTransformToNewOrigin ( + m_CalibratedProjectionCenter - m_RTIsocenter, + m_TLPS, + m_RLPS + ); + + m_outputTransform->SetCenter(m_CalibratedProjectionCenter); + + return m_outputTransform; + +} + +const char* queryStationName(const char* pcFName){ + + static std::string buffer; + + /* Check if we can open the file */ + gdcm::Reader R; + R.SetFileName(pcFName); + if(!R.Read()) + {std::cerr<<"ERROR: cannot read file: "<GetPointer(), bv->GetLength() ); + // Will be padded with at least one \0 + } + } + // Since return is a const char* the very first \0 will be considered + return buffer.c_str(); +} + +} diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h new file mode 100644 index 0000000..ea16394 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h @@ -0,0 +1,55 @@ + +#ifndef ITKIMAGEPROCESSORHELPERS_H +#define ITKIMAGEPROCESSORHELPERS_H + +#include "itkEuler3DTransform.h" +#include "itkImage.h" +#include "itkEuler3DTransform.h" + +#include "itkCommand.h" +#include "itkObject.h" +#include "itkObjectFactory.h" + +#include +#include + +namespace itk +{ + + +using TransformType = itk::Euler3DTransform; +constexpr static unsigned int Dimension = 3; +using PixelType3D = short; +using InternalPixelType = float; + +using ImageType3D = itk::Image; +using InternalImageType = itk::Image; + + +TransformType::Pointer +MapTransformToNewOrigin( + ImageType3D::PointType m_COR, + ImageType3D::PointType m_Translations, + ImageType3D::PointType m_Rotations + ); + + +TransformType::Pointer + CalculateInternalTransformV3( + ImageType3D::PointType m_Translation, + ImageType3D::PointType m_Rotation, + ImageType3D::PointType m_CalibratedProjectionCenter, + ImageType3D::PointType m_RTIsocenter, + InternalImageType::DirectionType m_IECtoLPSDirections + ); + +const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds); + +int query2DimageReconstructionDiameter(const char*); + +const char* queryStationName(const char*); + +} + + +#endif diff --git a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h new file mode 100644 index 0000000..89b231d --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h @@ -0,0 +1,192 @@ +#ifndef ITKQTITERATIONUPDATE_H +#define ITKQTITERATIONUPDATE_H + +#include "itkCommand.h" +#include +#include "itkPowellOptimizer.h" +#include "itkAmoebaOptimizer.h" +#include "itkExhaustiveOptimizer.h" + +#include "itkTwoProjectionImageRegistrationMethod.h" + + +class QObjectIterationUpdate : public QObject{ + Q_OBJECT +public: + QObjectIterationUpdate(){ + bAbortProcessCommand = false; + }; + + bool getAbortFlag(){ + return bAbortProcessCommand; + }; + void setAbortFlag(bool bVal){ + bAbortProcessCommand = bVal; + }; + +public slots: + void onAbortProcessCommand(){ + bAbortProcessCommand = true; + }; + void onIteration(double dI,double dX,double dY,double dZ){ + emit + sendRegistrationProgress(dI,dX,dY,dZ); + + } +private: + bool + bAbortProcessCommand; + + + + + +signals: +void sendRegistrationProgress(double,double,double,double); + +}; + +namespace itk +{ + + +class CommandIterationUpdate : public itk::Command { + // TODO: Move to own files. + + constexpr static unsigned int Dimension = 3; + +public: + QObjectIterationUpdate *objIterUpdate; + + using Self = CommandIterationUpdate; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer; + itkNewMacro(CommandIterationUpdate); + + using InternalPixelType = float; + using InternalImageType = itk::Image; + using ProcesssType = typename itk::TwoProjectionImageRegistrationMethod; + using ProcesssPointer = typename ProcesssType::Pointer; + + /** Set/Get the Process. */ + itkSetObjectMacro(Process, ProcesssType); + itkGetConstObjectMacro(Process, ProcesssType); + +private: + +protected: + CommandIterationUpdate() { + objIterUpdate = new QObjectIterationUpdate; + } + ProcesssPointer m_Process; + +public: + using OptimizerType = itk::PowellOptimizer; + using OptimizerPointer = const OptimizerType*; + + using AmoebaOptimizerType = itk::AmoebaOptimizer; + using AmoebaOptimizerPointer = const OptimizerType*; + + void + Execute(itk::Object* caller, const itk::EventObject& event) override + { + Execute((const itk::Object*)caller, event); + } + + void + Execute(const itk::Object* object, const itk::EventObject& event) override + { + if(objIterUpdate->getAbortFlag()){ + std::cout << "AbortRequest" << std::endl; + objIterUpdate->setAbortFlag(false); + throw itk::ProcessAborted(); + } +// std::cout << "Progress: " << this->m_Process->GetAbortGenerateData() << std::endl; + + auto optimizer = dynamic_cast(object); + + if (typeid(event) == typeid(itk::IterationEvent)) { + + //Feedback from the optimizer executed at the end of every itteration + // currently just print the result into the cout. We might add + // functionality to register and emit signals to update the UI. + + std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl; + auto oldprecision = std::cout.precision(); + std::cout.precision(std::numeric_limits::digits10 + 2); + std::cout << "Similarity: " << optimizer->GetValue() << std::endl; + std::cout.precision(oldprecision); + std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl; + +// objIterUpdate->onIteration( +// optimizer->GetCurrentIteration()+1, +// optimizer->GetCurrentPosition()[0], +// optimizer->GetCurrentPosition()[2], +// -optimizer->GetCurrentPosition()[1] +// ); + objIterUpdate->onIteration( + optimizer->GetCurrentIteration()+1, + optimizer->GetCurrentPosition()[0], + optimizer->GetCurrentPosition()[1], + optimizer->GetCurrentPosition()[2] + ); + } + return; + } +}; + +class ExhaustiveCommandIterationUpdate : public itk::Command { + // TODO: Move to own files. + +public: + using Self = ExhaustiveCommandIterationUpdate; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer; + itkNewMacro(Self); + + std::ostream* os; + +protected: + ExhaustiveCommandIterationUpdate() = default; + +public: + using OptimizerType = itk::ExhaustiveOptimizer; + ; + using OptimizerPointer = const OptimizerType*; + + void set_stream(std::ostream& stream) + { + os = &stream; + } + + void + Execute(itk::Object* caller, const itk::EventObject& event) override + { + Execute((const itk::Object*)caller, event); + } + + void + Execute(const itk::Object* object, const itk::EventObject& event) override + { + auto optimizer = dynamic_cast(object); + if (typeid(event) == typeid(itk::IterationEvent)) { + + //crude LPS to IEC transform for HFS. + auto position = optimizer->GetCurrentPosition(); + + auto oldprecision = os->precision(); + os->precision(std::numeric_limits::digits10 + 2); + *os << optimizer->GetCurrentValue(); + os->precision(oldprecision); + //*os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl; + *os << "\t" << position[0] << "\t" << position[1] << "\t" << position[2] << std::endl; + + } + return; + } +}; + +} + +#endif + diff --git a/reg23Topograms/itkDTRrecon/itkReg23.cpp b/reg23Topograms/itkDTRrecon/itkReg23.cpp new file mode 100644 index 0000000..e599807 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkReg23.cpp @@ -0,0 +1,393 @@ +#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); + + 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 <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.; + +} + + + + +} 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 diff --git a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h index 86579e3..6c48c3b 100644 --- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h +++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h @@ -182,6 +182,10 @@ public: itkSetMacro(Threshold, double); itkGetMacro(Threshold, double); + /** Set and get the Panel Offset */ + itkSetMacro(PanelOffset, double); + itkGetMacro(PanelOffset, double); + /** Check if a point is inside the image buffer. * \warning For efficiency, no validity checking of * the input image pointer is done. */ @@ -232,6 +236,7 @@ protected: double m_Threshold; double m_FocalPointToIsocenterDistance; // Focal point to isocenter distance double m_ProjectionAngle; // Linac gantry rotation angle in radians + double m_PanelOffset; private: void diff --git a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx index e6eb7df..ee26045 100644 --- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx +++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx @@ -62,6 +62,7 @@ gSiddonJacobsRayCastInterpolateImageFunction::gSiddonJac m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm. m_ProjectionAngle = 0.; // Angle in radians betweeen projection central axis and reference axis m_Threshold = 0.; // Intensity threshold, below which is ignored. + m_PanelOffset = 0.; m_SourcePoint[0] = 0.; m_SourcePoint[1] = 0.; @@ -148,14 +149,15 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c // m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint); } + PointType PointReq = point; + PointReq[0] += m_PanelOffset; - drrPixelWorld = m_InverseTransform->TransformPoint(point); + drrPixelWorld = m_InverseTransform->TransformPoint(PointReq); PointType SlidingSourcePoint = m_SourcePoint; SlidingSourcePoint[0] = 0.; SlidingSourcePoint[1] = point[1]; SlidingSourcePoint[2] = 0.; - //PointType SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint); PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint); //std::cout<<"SourceWorld: "<::Evaluate(c regionCT = inputPtr->GetLargestPossibleRegion(); sizeCT = regionCT.GetSize(); + // If Pixel position (in mm) is outside bounds of CT (zero-based) + // assign 0 at the pixel and move on. Ensures regular spacing of resulting + // DRR + float xSizeCT = sizeCT[0] * ctPixelSpacing[0]; + float ySizeCT = sizeCT[1] * ctPixelSpacing[1]; + float zSizeCT = sizeCT[2] * ctPixelSpacing[2]; + float xDrrPix = drrPixelWorld[0];float yDrrPix = drrPixelWorld[1];float zDrrPix = drrPixelWorld[2]; +// if(zDrrPix < 0 /*|| yDrrPix < 0 || xDrrPix < 0*/) +// std::cout << drrPixelWorld[0]<<" " <xSizeCT || yDrrPix>ySizeCT || zDrrPix>zSizeCT){ + pixval = static_cast(0.0); + return pixval; + } // calculate the detector position for this pixel center by moving // 2*m_FocalPointToIsocenterDistance from the source in the pixel // directions diff --git a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx index c3bd7bb..602abf9 100644 --- a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx +++ b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx @@ -66,6 +66,22 @@ vtkContourTopogramProjectionFilter::~vtkContourTopogramProjectionFilter() } } +void vtkContourTopogramProjectionFilter::cleanUpFilter(){ + + this->dProjectionLPSOffset[0]=0.; + this->dProjectionLPSOffset[1]=0.; + this->dProjectionLPSOffset[2]=0.; + + this->dImportOffsetLPS[0] = 0.; + this->dImportOffsetLPS[1] = 0.; + this->dImportOffsetLPS[2] = 0.; + + this->dSCD = 0.; + + this->m_RefTransform = nullptr; + this->m_Transform = nullptr; + +} int vtkContourTopogramProjectionFilter::RequestData( vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector) @@ -94,6 +110,7 @@ int vtkContourTopogramProjectionFilter::RequestData( // If no points, then nothing to do. if (points == nullptr) { + std::cout << "Cannot Project; no input points" << std::endl; vtkDebugMacro("Cannot Project; no input points"); return 1; } @@ -101,6 +118,7 @@ int vtkContourTopogramProjectionFilter::RequestData( // If reference transform, then nothing to do. if (m_RefTransform == nullptr) { + std::cout << "Cannot Project; no input reference projection transform" << std::endl; vtkDebugMacro("Cannot Project; no input reference projection transform"); return 1; } @@ -108,6 +126,7 @@ int vtkContourTopogramProjectionFilter::RequestData( // If transform, then nothing to do. if (m_Transform == nullptr) { + std::cout << "Cannot Project; no input projection transform" << std::endl; vtkDebugMacro("Cannot Project; no input projection transform"); return 1; } diff --git a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h index 568d308..da0f84a 100644 --- a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h +++ b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h @@ -41,6 +41,8 @@ public: void SetSCD(const double dVal); void SetImportOffsetLPS(const double *); + void cleanUpFilter(); + protected: vtkContourTopogramProjectionFilter(); ~vtkContourTopogramProjectionFilter() override;