From b0beda81b86f9465215b60c4b15de3942be09de7 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Wed, 3 May 2023 16:43:52 +0200 Subject: [PATCH 01/12] REG working --- .../itkDTRrecon/itkImageProcessor.cpp | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 8af1bf8..63925cc 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -1495,11 +1495,6 @@ void itkImageProcessor::InitializeRegistration( TransformType::Pointer CurrTransform; if(m_RTMetaInfo == NULL) { - - //AAAAAA TODOOOOO CERCA - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage1MetaInfo->GetProjectionOriginLPS(); - CurrTransform = CalculateInternalTransformV2( ZeroPoint, @@ -1545,9 +1540,6 @@ void itkImageProcessor::InitializeRegistration( /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/ if(m_RTMetaInfo == NULL) { - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage2MetaInfo->GetProjectionOriginLPS(); - CurrTransform = CalculateInternalTransformV2( ZeroPoint, @@ -1743,10 +1735,20 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) m_OptmizerValue = optimizer->GetValue(); + auto finalParameters = transform1->GetParameters(); - optimizer->GetCurrentPosition(); + //optimizer->GetCurrentPosition(); finalParameters = finalParameters - startParameters; + std::cout<<"startParameters "<GetCenter()<GetTranslation()<GetCenter()<GetTranslation()<SetParameters(finalParameters); @@ -1784,6 +1786,7 @@ void itkImageProcessor::InitializeProjector() m_DRTGeometryMetaInfo == NULL || m_TransformMetaInfo == NULL ){ //TODO + return; } const double dtr = (atan(1.0) * 4.0) / 180.0; @@ -2181,7 +2184,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),// m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) @@ -2192,7 +2194,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),//m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); @@ -2260,7 +2261,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS )); @@ -2273,7 +2273,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) @@ -2284,7 +2283,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); @@ -2326,6 +2324,7 @@ void itkImageProcessor::GetProjectionImages(){ m_DRTGeometryMetaInfo == NULL || m_TransformMetaInfo == NULL ){ //TODO + return; } const double dtr = (atan(1.0) * 4.0) / 180.0; From 9f126e8a34e34dba7b6a2170333e3ef6fd92303b Mon Sep 17 00:00:00 2001 From: Proton local user Date: Thu, 11 May 2023 09:41:10 +0200 Subject: [PATCH 02/12] debug commented --- reg23Topograms/itkDTRrecon/itkImageProcessor.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 63925cc..92e9645 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -2384,6 +2384,13 @@ void itkImageProcessor::GetProjectionImages(){ transform1 ->SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); + std::cout<< "CurrTransform Rotations: "<< + CurrTransform->GetAngleX() <<" "<< + CurrTransform->GetAngleY() <<" "<< + CurrTransform->GetAngleZ() <<" "<< std::endl; + + std::cout<< "CurrTransform Translations: "<< + CurrTransform->GetTranslation()<Print(std::cout); // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance From 07db1193150445bf0c1804a07e691576ff606c26 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Fri, 12 May 2023 13:29:02 +0200 Subject: [PATCH 03/12] backup before changing meta reg --- .../itkTwoProjectionImageRegistrationMethod.hxx | 2 +- .../autoreg/itkgTwoImageToOneImageMetric.hxx | 16 ++++++++-------- reg23Topograms/itkDTRrecon/itkImageProcessor.cpp | 5 ++++- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx index b6f4b12..81a7e91 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -201,7 +201,7 @@ void TwoProjectionImageRegistrationMethod::OnInternal // if (m_CancelRequest) // filter->SetAbortGenerateData(true); // may be handled by filter - std::cout << "Porca Madonna " << std::endl; + std::cout << "OnInternalFilterProgressReceptor :: filter (TURE) " << std::endl; } } diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index 90a9548..0e9c600 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -97,14 +97,14 @@ bool gTwoImageToOneImageMetric::SetTransformParameter break; } -// std::cout << "New Transform Parameters = " << 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; + std::cout << "New Transform Parameters = " << 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; transformParameters[0] = RotationAlongX; transformParameters[1] = RotationAlongY; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 36699a1..0049101 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -1734,8 +1734,11 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) m_OptmizerValue = optimizer->GetValue(); - auto finalParameters = transform1->GetParameters(); + std::cout<<"-->startParameters "<finalParameters1 "<GetParameters()<finalParameters2 "<GetParameters()<GetCurrentPosition(); finalParameters = finalParameters - startParameters; From c67281a10a038e0892e8e825fa9692d963cf9c73 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Fri, 12 May 2023 14:00:32 +0200 Subject: [PATCH 04/12] new R23MetaInfo so that R23 does not need TransformMeta anymore --- .../itkDTRrecon/DRTMetaInformation.cpp | 26 +++++++- .../itkDTRrecon/DRTMetaInformation.h | 64 +++++++++++++++---- .../itkTwoProjectionImageRegistrationMethod.h | 9 ++- .../autoreg/itkgTwoImageToOneImageMetric.h | 6 +- .../autoreg/itkgTwoImageToOneImageMetric.hxx | 4 +- .../itkDTRrecon/itkImageProcessor.cpp | 17 +++-- .../itkDTRrecon/itkImageProcessor.h | 23 +------ 7 files changed, 103 insertions(+), 46 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 55ada69..9b0c1dc 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -518,7 +518,27 @@ RTGeometryMetaInformation } +R23MetaInformation:: +R23MetaInformation (){ + m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + +} + + +void +R23MetaInformation +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + + +R23MetaInformation +::~R23MetaInformation () +{ + +} TransformMetaInformation :: TransformMetaInformation (){ @@ -531,11 +551,11 @@ TransformMetaInformation (){ m_UserTranslations.Fill(0.); - m_ReferenceTransform.Fill(0.); + //m_ReferenceTransform.Fill(0.); - m_CurrentTransform.Fill(0.); + //m_CurrentTransform.Fill(0.); - m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; + // m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN; } diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 1f77f24..917a7f2 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -514,6 +514,48 @@ private: }; +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; + + itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + + +protected: + + tDegreeOfFreedomEnum + m_DegreeOfFreedom; + + /** 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{ @@ -547,14 +589,14 @@ public: itkSetMacro(UserRotations,PointType); itkGetMacro(UserRotations,PointType); - itkSetMacro(ReferenceTransform,TransformMatrixType); - itkGetMacro(ReferenceTransform,TransformMatrixType); +// itkSetMacro(ReferenceTransform,TransformMatrixType); +// itkGetMacro(ReferenceTransform,TransformMatrixType); - itkSetMacro(CurrentTransform,TransformMatrixType); - itkGetMacro(CurrentTransform,TransformMatrixType); +// itkSetMacro(CurrentTransform,TransformMatrixType); +// itkGetMacro(CurrentTransform,TransformMatrixType); - itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); - itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); +// itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); +// itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); protected: @@ -564,12 +606,12 @@ protected: m_UserTranslations, m_UserRotations; - TransformMatrixType - m_ReferenceTransform, - m_CurrentTransform; +// TransformMatrixType +// m_ReferenceTransform, +// m_CurrentTransform; - tDegreeOfFreedomEnum - m_DegreeOfFreedom; +// tDegreeOfFreedomEnum +// m_DegreeOfFreedom; /** Default Constructor **/ TransformMetaInformation (); diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h index 3dd6241..7770224 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -170,8 +170,8 @@ public: itkGetConstObjectMacro(Interpolator2, InterpolatorType); /** Set/Get the Meta informations. */ - itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation); - itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation); + itkSetObjectMacro(TransformMetaInfo, R23MetaInformation); + itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation); /** Set/Get the output filters. */ itkSetObjectMacro(Filter1, ChangeInformationFilterType); @@ -253,7 +253,10 @@ private: InterpolatorPointer m_Interpolator1; InterpolatorPointer m_Interpolator2; - TransformMetaInformation::Pointer + //TransformMetaInformation::Pointer + // m_TransformMetaInfo; + + R23MetaInformation::Pointer m_TransformMetaInfo; ChangeInformationFilterPointer diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h index 9cb93cb..c3d8712 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h @@ -185,10 +185,10 @@ public: itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long); /** Connect the DRTGeometryMetaInfo. */ - itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation); + itkSetObjectMacro(TransformMetaInfo, R23MetaInformation); /** Get a pointer to the DRTGeometryMetaInfo. */ - itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation); + itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation); /** Set the region over which the metric will be computed */ itkSetMacro(FixedImageRegion1, FixedImageRegionType); @@ -274,7 +274,7 @@ protected: mutable FixedImageMaskPointer m_FixedImageMask2; mutable MovingImageMaskPointer m_MovingImageMask; - TransformMetaInformation::Pointer + R23MetaInformation::Pointer m_TransformMetaInfo; ChangeInformationFilterPointer diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index 0e9c600..d4ea48a 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -113,7 +113,9 @@ bool gTwoImageToOneImageMetric::SetTransformParameter transformParameters[4] = TranslationAlongY; transformParameters[5] = TranslationAlongZ; - bool transformValid = (std::abs(TranslationAlongX) < m_MaxTranslation) && (std::abs(TranslationAlongY) < m_MaxTranslation) && (std::abs(TranslationAlongZ) < m_MaxTranslation); + bool transformValid = (std::abs(TranslationAlongX) < m_MaxTranslation) && + (std::abs(TranslationAlongY) < m_MaxTranslation) && + (std::abs(TranslationAlongZ) < m_MaxTranslation); if (transformValid) { m_Transform1->SetParameters(transformParameters); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 0049101..66326dd 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -118,6 +118,9 @@ itkImageProcessor::itkImageProcessor() 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 @@ -284,13 +287,16 @@ double itkImageProcessor::GetPanelOffset2(){ void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof) { - m_TransformMetaInfo->SetDegreeOfFreedom(dof); + // m_TransformMetaInfo->SetDegreeOfFreedom(dof); + m_r23MetaInfo->SetDegreeOfFreedom(dof); } tDegreeOfFreedomEnum itkImageProcessor::GetDegreeOfFreedom() { - return m_TransformMetaInfo->GetDegreeOfFreedom(); + return + m_r23MetaInfo->GetDegreeOfFreedom(); + //m_TransformMetaInfo->GetDegreeOfFreedom(); } void itkImageProcessor::SetCustom_ImportTransform(double dTx, @@ -428,6 +434,7 @@ int itkImageProcessor::unload3DVolumeAndMeta(){ m_TransformMetaInfo = NULL; m_TransformMetaInfo = TransformMetaInformation::New(); + return 1; } @@ -1436,7 +1443,8 @@ void itkImageProcessor::InitializeRegistration( metric->SetSubtractMean(true); metric->SetMaxTranslation(maxTranslation); - m_TransformMetaInfo->SetDegreeOfFreedom(dof); +// m_TransformMetaInfo->SetDegreeOfFreedom(dof); + m_r23MetaInfo->SetDegreeOfFreedom(dof); mimetric->ComputeGradientOff(); @@ -1479,7 +1487,8 @@ void itkImageProcessor::InitializeRegistration( } - registration->SetTransformMetaInfo(m_TransformMetaInfo); +// registration->SetTransformMetaInfo(m_TransformMetaInfo); + registration->SetTransformMetaInfo(m_r23MetaInfo); registration->SetFilter1(filter1); registration->SetFilter2(filter2); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 553212e..5a25620 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -269,18 +269,6 @@ private: TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo); -// 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 -// ); - /* Calculate the transform used in siddon. * The isocentric transform is mapped to the calibrated center of projection */ TransformType::Pointer @@ -294,8 +282,6 @@ private: InternalImageType::DirectionType m_IECtoLPSDirections ); - - TransformType::Pointer transform1, transform2; @@ -360,13 +346,6 @@ private: ); -// ImageType3D::PointType -// CalcDRTImageOffset( -// ImageType3D::PointType m_DRTOffset, -// double dAngle, -// InternalImageType::DirectionType stdDRT2LPS -// ); - TransformType::Pointer MapTransformToNewOrigin( ImageType3D::PointType m_COR, @@ -433,6 +412,8 @@ private: TransformMetaInformation::Pointer m_TransformMetaInfo; + R23MetaInformation::Pointer + m_r23MetaInfo; double m_OptmizerValue; int m_MaxNumberOfIterations; From 9a8faa0693d58c9f90501f9f4824631073ec4344 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Fri, 12 May 2023 14:55:17 +0200 Subject: [PATCH 05/12] Simplified CalcInternalTransformV3 which makes use of T and R calculated by the transform Meta --- .../itkDTRrecon/DRTMetaInformation.cpp | 34 +- .../itkDTRrecon/DRTMetaInformation.h | 28 +- .../itkDTRrecon/itkImageProcessor.cpp | 402 +++++++++--------- .../itkDTRrecon/itkImageProcessor.h | 9 + 4 files changed, 261 insertions(+), 212 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 9b0c1dc..8395221 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -543,9 +543,9 @@ R23MetaInformation TransformMetaInformation :: TransformMetaInformation (){ - m_Translations.Fill(0.); + m_HiddenTranslations.Fill(0.); - m_Rotations.Fill(0.); + m_HiddenRotations.Fill(0.); m_UserRotations.Fill(0.); @@ -559,6 +559,36 @@ TransformMetaInformation (){ } +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; + +} + void TransformMetaInformation ::PrintSelf(std::ostream& os, itk::Indent indent) const diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 917a7f2..f797b30 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -577,11 +577,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); @@ -589,29 +589,17 @@ public: itkSetMacro(UserRotations,PointType); itkGetMacro(UserRotations,PointType); -// itkSetMacro(ReferenceTransform,TransformMatrixType); -// itkGetMacro(ReferenceTransform,TransformMatrixType); - -// itkSetMacro(CurrentTransform,TransformMatrixType); -// itkGetMacro(CurrentTransform,TransformMatrixType); - -// itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); -// itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum); + PointType GetT(); + PointType GetR(); protected: PointType - m_Translations, - m_Rotations, + m_HiddenTranslations, + m_HiddenRotations, m_UserTranslations, m_UserRotations; -// TransformMatrixType -// m_ReferenceTransform, -// m_CurrentTransform; - -// tDegreeOfFreedomEnum -// m_DegreeOfFreedom; /** Default Constructor **/ TransformMetaInformation (); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 66326dd..c828274 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -136,7 +136,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; @@ -287,7 +287,7 @@ double itkImageProcessor::GetPanelOffset2(){ void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof) { - // m_TransformMetaInfo->SetDegreeOfFreedom(dof); + // m_TransformMetaInfo->SetDegreeOfFreedom(dof); m_r23MetaInfo->SetDegreeOfFreedom(dof); } @@ -296,7 +296,7 @@ itkImageProcessor::GetDegreeOfFreedom() { return m_r23MetaInfo->GetDegreeOfFreedom(); - //m_TransformMetaInfo->GetDegreeOfFreedom(); + //m_TransformMetaInfo->GetDegreeOfFreedom(); } void itkImageProcessor::SetCustom_ImportTransform(double dTx, @@ -523,11 +523,11 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ 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(); -// } + // 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(); @@ -538,7 +538,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ } catch (itk::ExceptionObject & ex) { -// std::cout << ex << std::endl; + // std::cout << ex << std::endl; return EXIT_FAILURE; } @@ -755,22 +755,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(); @@ -911,7 +911,7 @@ void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2) m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1); m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2); - std::cout << "SetProjectionAngleOffsetIEC " << dOff1 << " " << dOff2 << std::endl; + std::cout << "SetProjectionAngleOffsetIEC " << dOff1 << " " << dOff2 << std::endl; } @@ -986,7 +986,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 " <Update(); -// std::cout<< " ////////////// 2D Topo BEG " <GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<SetCenter(m_CalibratedProjectionCenter); + + return m_outputTransform; + +} + + void itkImageProcessor::InitializeRegistration( double stepLength, double maxTranslation, @@ -1443,7 +1476,7 @@ void itkImageProcessor::InitializeRegistration( metric->SetSubtractMean(true); metric->SetMaxTranslation(maxTranslation); -// m_TransformMetaInfo->SetDegreeOfFreedom(dof); + // m_TransformMetaInfo->SetDegreeOfFreedom(dof); m_r23MetaInfo->SetDegreeOfFreedom(dof); mimetric->ComputeGradientOff(); @@ -1487,7 +1520,7 @@ void itkImageProcessor::InitializeRegistration( } -// registration->SetTransformMetaInfo(m_TransformMetaInfo); + // registration->SetTransformMetaInfo(m_TransformMetaInfo); registration->SetTransformMetaInfo(m_r23MetaInfo); registration->SetFilter1(filter1); @@ -1504,28 +1537,27 @@ void itkImageProcessor::InitializeRegistration( if(m_RTMetaInfo == NULL) { CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + IECtoLPS_Directions); + } else { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), + CurrTransform = + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), IECtoLPS_Directions ); + + } @@ -1548,29 +1580,26 @@ void itkImageProcessor::InitializeRegistration( /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/ if(m_RTMetaInfo == NULL) { - CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + CurrTransform = + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + IECtoLPS_Directions); + } else { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint, - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); + CurrTransform = + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), + IECtoLPS_Directions + ); } @@ -1707,10 +1736,10 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) } try { -// timer.Start("Registration"); + // timer.Start("Registration"); // Start the registration. registration->StartRegistration(); -// timer.Stop("Registration"); + // timer.Stop("Registration"); } catch (itk::ExceptionObject& err) { registration->ResetPipeline(); @@ -1818,29 +1847,25 @@ void itkImageProcessor::InitializeProjector() if(m_RTMetaInfo == NULL) { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - ZeroPoint , - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + CurrTransform = CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + IECtoLPS_Directions); + } else { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); - + CurrTransform = + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), + IECtoLPS_Directions + ); } @@ -1877,28 +1902,25 @@ void itkImageProcessor::InitializeProjector() { CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + IECtoLPS_Directions); + } else { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint, - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), + CurrTransform = CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), IECtoLPS_Directions ); + } @@ -1985,7 +2007,7 @@ void itkImageProcessor::InitializeProjector() filter1->ChangeOriginOn(); filter1->UpdateOutputInformation(); -// std::cout<< "itkImageProcessor::InitializeProjector() " <Update(); @@ -2033,7 +2055,14 @@ int itkImageProcessor::unloadRTPlanAndMeta(){ std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; m_RTMetaInfo = NULL; - return 1; + ImageType3D::PointType + pZero; + pZero.Fill(0.); + + m_CTMetaInfo->SetImportOffset(pZero); + m_TransformMetaInfo->SetHiddenRotations(pZero); + + return 0; } @@ -2082,6 +2111,9 @@ void itkImageProcessor::loadRTPlanInfo( ) ); + m_TransformMetaInfo->SetHiddenRotations( + m_DRTGeometryMetaInfo->GetIECS2IECScannerR()); + this->UpdateProjectionGeometryMeta(); this->InitializeProjector(); @@ -2133,11 +2165,11 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint( m_DRTGeometryMetaInfo->GetProjectionCenterOffset1()); -// std::cout<< "///////////////// PGEOM META BEG ///////////////" <GetProjectionCenter() <GetProjectionCenter() <SetProjectionAngleLPS( this->CalcProjectionAngleLPS( @@ -2334,7 +2366,7 @@ itkImageProcessor::CalcDRTImageOrigin( void itkImageProcessor::GetProjectionImages(){ -// std::cout<< "itkImageProcessor::GetProjectionImages" <GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + IECtoLPS_Directions); } else { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint , - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); + CurrTransform = + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), + IECtoLPS_Directions + ); } @@ -2425,23 +2453,18 @@ void itkImageProcessor::GetProjectionImages(){ { CurrTransform = - CalculateInternalTransformV2( - ZeroPoint, - ZeroPoint, - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions - ); + CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + IECtoLPS_Directions); } else { - CurrTransform = CalculateInternalTransformV2( - ZeroPoint, - m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), + CurrTransform = CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), IECtoLPS_Directions @@ -2634,22 +2657,22 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK() ImageType3D::PointType LastVoxelPosition = origin3D + (imagDirection * Dest); -// std::cout<<" -------- Localizer 1 ITK --------"<GetOutput()->GetBounds(); -// std::cout<< "-------- Localizer 1 VTK --------" <SetInput( imageDRT1In ); rescaler1->Update(); -// using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator; -// auto imageCalculatorFilter = ImageCalculatorFilterType::New(); -// imageCalculatorFilter->SetImage(imageDRT1In); -// imageCalculatorFilter->Compute(); + // 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 <; + // auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New(); + // imageCalculatorFilter2->SetImage(rescaler1->GetOutput()); + // imageCalculatorFilter2->Compute(); + // std::cout<< "itkImageProcessor::imageDRT2In() " << + // imageCalculatorFilter2 <SetInput(rescaler1->GetOutput()); toVTK2D1->Update(); @@ -2908,12 +2931,12 @@ vtkImageData* itkImageProcessor::GetProjection2VTK() - //double* dBounds = toVTK2D2->GetOutput()->GetBounds(); -// std::cout<< "-------- Proj 2 --------" <GetOutput()->GetBounds(); + // std::cout<< "-------- Proj 2 --------" <GetOutput(); @@ -2947,7 +2970,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) @@ -2961,7 +2984,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) @@ -2983,13 +3006,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) @@ -2999,18 +3021,19 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) Rotations[1]= dY; Rotations[2]= dZ; - m_TransformMetaInfo->SetRotations(Rotations); + m_TransformMetaInfo->SetUserRotations(Rotations); } +// THIS IS READ FOR EXPORT. TO FIX. double* itkImageProcessor::GetTransformParameters(){ ImageType3D::PointType Translations; ImageType3D::PointType Rotations; - Translations = m_TransformMetaInfo->GetTranslations(); - Rotations = m_TransformMetaInfo->GetRotations(); + Translations = m_TransformMetaInfo->GetUserTranslations(); + Rotations = m_TransformMetaInfo->GetUserRotations(); dTransfParam[0] = Translations[0]; dTransfParam[1] = Translations[1]; @@ -3023,27 +3046,26 @@ double* itkImageProcessor::GetTransformParameters(){ return dTransfParam; } -// Doubts about this user thing +// THESE ARE CALLED AT THE END OF REG void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) { - ImageType3D::PointType translations = m_TransformMetaInfo->GetTranslations(); ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations(); - dX = translations[0] + userTranslations[0]; - dY = translations[1] + userTranslations[1]; - dZ = translations[2] + userTranslations[2]; + dX = userTranslations[0]; + dY = userTranslations[1]; + dZ = userTranslations[2]; } +// THESE ARE CALLED AT THE END OF REG void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ) { - ImageType3D::PointType rotations = m_TransformMetaInfo->GetRotations(); ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations(); - dRX = rotations[0] + userRotations[0]; - dRY = rotations[1] + userRotations[1]; - dRZ = rotations[2] + userRotations[2]; -} + dRX = userRotations[0]; + dRY = userRotations[1]; + dRZ = userRotations[2]; +} double itkImageProcessor::GetOptimizerValue() { diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 5a25620..72af105 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -282,6 +282,15 @@ private: InternalImageType::DirectionType m_IECtoLPSDirections ); + TransformType::Pointer + CalculateInternalTransformV3( + ImageType3D::PointType m_Translation, + ImageType3D::PointType m_Rotation, + ImageType3D::PointType m_CalibratedProjectionCenter, + ImageType3D::PointType m_RTIsocenter, + InternalImageType::DirectionType m_IECtoLPSDirections + ); + TransformType::Pointer transform1, transform2; From 6ad0344791e162433f70f9abb4bac787cb7fd01e Mon Sep 17 00:00:00 2001 From: Proton local user Date: Fri, 12 May 2023 16:13:12 +0200 Subject: [PATCH 06/12] MapTransformToNewOrigin and CalcInternalTransformV3 moved to another file, ouside of scout processor --- reg23Topograms/itkDTRrecon/CMakeLists.txt | 2 + .../itkTwoProjectionImageRegistrationMethod.h | 6 + ...tkTwoProjectionImageRegistrationMethod.hxx | 9 +- .../autoreg/itkgTwoImageToOneImageMetric.h | 6 + .../autoreg/itkgTwoImageToOneImageMetric.hxx | 13 +- .../itkDTRrecon/itkImageProcessor.cpp | 377 +++++++++--------- .../itkDTRrecon/itkImageProcessor.h | 45 +-- .../itkDTRrecon/itkImageProcessorHelpers.cpp | 93 +++++ .../itkDTRrecon/itkImageProcessorHelpers.h | 38 ++ 9 files changed, 377 insertions(+), 212 deletions(-) create mode 100644 reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp create mode 100644 reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index fd530ce..e97344c 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -16,6 +16,7 @@ INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS}) SET(SRCS itkImageProcessor.cpp + itkImageProcessorHelpers.cpp vtkContourTopogramProjectionFilter.cxx DRTMetaInformation.cpp @@ -23,6 +24,7 @@ SET(SRCS SET(HDR itkImageProcessor.h + itkImageProcessorHelpers.h itkQtIterationUpdate.h itkgSiddonJacobsRayCastInterpolateImageFunction.h itkgSiddonJacobsRayCastInterpolateImageFunction.hxx diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h index 7770224..10d638b 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -163,6 +163,10 @@ public: 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); @@ -248,6 +252,8 @@ private: FixedImageConstPointer m_FixedImage1; FixedImageConstPointer m_FixedImage2; + TransformPointer m_IsocIECTransform; + TransformPointer m_Transform1; TransformPointer m_Transform2; InterpolatorPointer m_Interpolator1; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx index 81a7e91..fa127aa 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -38,6 +38,7 @@ TwoProjectionImageRegistrationMethod::TwoProjectionIm 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_Interpolator1 = nullptr; // has to be provided by the user. m_Interpolator2 = nullptr; // has to be provided by the user. @@ -80,6 +81,7 @@ void TwoProjectionImageRegistrationMethod::SetFixedIm m_FixedImageRegionDefined2 = true; } + /* * Initialize by setting the interconnects between components. */ @@ -119,7 +121,8 @@ void TwoProjectionImageRegistrationMethod::Initialize // Connect the transform to the Decorator. auto* transformOutput = static_cast(this->ProcessObject::GetOutput(0)); - transformOutput->Set(m_Transform1.GetPointer()); + //transformOutput->Set(m_Transform1.GetPointer()); + transformOutput->Set(m_IsocIECTransform.GetPointer()); if (!m_Interpolator1) { itkExceptionMacro(<< "Interpolator1 is not present"); @@ -147,6 +150,7 @@ void TwoProjectionImageRegistrationMethod::Initialize 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->SetFilter1(m_Filter1); m_Metric->SetFilter2(m_Filter2); @@ -174,7 +178,8 @@ void TwoProjectionImageRegistrationMethod::Initialize m_Optimizer->SetScales(m_Metric->GetWeightings()); m_Optimizer->SetCostFunction(m_Metric); - m_InitialOptimizerParameters = m_Metric->GetParameters(); + m_InitialOptimizerParameters = + m_Metric->GetParameters(); m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters); } diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h index c3d8712..752f7c1 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h @@ -169,6 +169,10 @@ public: /** Get a pointer to the Transform2. */ itkGetConstObjectMacro(Transform2, TransformType); + itkSetObjectMacro(IsocTransf, TransformType); + itkGetObjectMacro(IsocTransf, TransformType); + + /** Connect the Interpolator. */ itkSetObjectMacro(Interpolator1, InterpolatorType); @@ -261,6 +265,8 @@ protected: mutable TransformPointer m_Transform1; mutable TransformPointer m_Transform2; + mutable TransformPointer m_IsocTransf; + InterpolatorPointer m_Interpolator1; InterpolatorPointer m_Interpolator2; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index d4ea48a..8cade2d 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -55,7 +55,8 @@ bool gTwoImageToOneImageMetric::SetTransformParameter itkExceptionMacro(<< "Transform 2 has not been assigned"); } - auto transformParameters = m_Transform1->GetParameters(); + //auto transformParameters = m_Transform1->GetParameters(); + auto transformParameters = m_IsocTransf->GetParameters(); double TranslationAlongX = transformParameters[3]; double TranslationAlongY = transformParameters[4]; double TranslationAlongZ = transformParameters[5]; @@ -106,6 +107,10 @@ bool gTwoImageToOneImageMetric::SetTransformParameter std::cout << " Rotation Along Z = " << RotationAlongZ / dtr << " deg" << std::endl; std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; + + //HERE WE HAVE TO CALCULATE INTERNAL TRANSFORMS!!! + // GIOVANNI + transformParameters[0] = RotationAlongX; transformParameters[1] = RotationAlongY; transformParameters[2] = RotationAlongZ; @@ -281,10 +286,12 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric: } template -typename gTwoImageToOneImageMetric::ParametersType gTwoImageToOneImageMetric::GetParameters() const +typename gTwoImageToOneImageMetric:: +ParametersType gTwoImageToOneImageMetric::GetParameters() const { - ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ +// 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: diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index c828274..0c2c739 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -1281,64 +1281,63 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg) } -itkImageProcessor::TransformType::Pointer -itkImageProcessor::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 - ){ +//itkImageProcessor::TransformType::Pointer +//itkImageProcessor::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::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]; +// TransformType::OutputVectorType translation; +// translation[0] = m_Translations[0]; +// translation[1] = m_Translations[1]; +// translation[2] = m_Translations[2]; - InputTransform->SetTranslation(translation); +// 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]); +// 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 m_TransformOrigin; +// m_TransformOrigin.Fill(0.); +// InputTransform->SetCenter( +// m_TransformOrigin ); - ImageType3D::PointType NewOriginTranslations = - InputTransform->TransformPoint(m_COR); +// ImageType3D::PointType NewOriginTranslations = +// InputTransform->TransformPoint(m_COR); - ImageType3D::PointType DeltaNewOrigin = - NewOriginTranslations - m_COR; +// ImageType3D::PointType DeltaNewOrigin = +// NewOriginTranslations - m_COR; - TransformType::Pointer m_OutputTransform = - TransformType::New(); - m_OutputTransform ->SetComputeZYX(true); - m_OutputTransform ->SetIdentity(); +// 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]; +// 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->SetTranslation(translation); +// m_OutputTransform->SetRotation( +// dtr * m_Rotations[0], +// dtr * m_Rotations[1], +// dtr * m_Rotations[2]); - m_OutputTransform->SetCenter(m_COR); +// m_OutputTransform->SetCenter(m_COR); - InputTransform = NULL; +// InputTransform = NULL; - // m_OutputTransform.Print(std::cout); - return m_OutputTransform; -} +// return m_OutputTransform; +//} // This External User Transform thing need to be checked out @@ -1369,88 +1368,88 @@ void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer tr } -itkImageProcessor::TransformType::Pointer -itkImageProcessor::CalculateInternalTransformV2( - ImageType3D::PointType m_TranslationOffset, //IEC - ImageType3D::PointType m_RotationOffset, //IEC - ImageType3D::PointType m_TranslationUser, //IEC - ImageType3D::PointType m_RotationUser, //IEC - ImageType3D::PointType m_CalibratedProjectionCenter, //LPS - ImageType3D::PointType m_RTIsocenter, //LPS - InternalImageType::DirectionType m_IECtoLPSDirections - ) -{ +//itkImageProcessor::TransformType::Pointer +//itkImageProcessor::CalculateInternalTransformV2( +// ImageType3D::PointType m_TranslationOffset, //IEC +// ImageType3D::PointType m_RotationOffset, //IEC +// ImageType3D::PointType m_TranslationUser, //IEC +// ImageType3D::PointType m_RotationUser, //IEC +// ImageType3D::PointType m_CalibratedProjectionCenter, //LPS +// ImageType3D::PointType m_RTIsocenter, //LPS +// InternalImageType::DirectionType m_IECtoLPSDirections +// ) +//{ - //Convert all inputs into LPS +// //Convert all inputs into LPS - ImageType3D::PointType m_TOffsetLPS = - m_IECtoLPSDirections * m_TranslationOffset; +// ImageType3D::PointType m_TOffsetLPS = +// m_IECtoLPSDirections * m_TranslationOffset; - ImageType3D::PointType m_ROffsetLPS = - m_IECtoLPSDirections * m_RotationOffset; +// ImageType3D::PointType m_ROffsetLPS = +// m_IECtoLPSDirections * m_RotationOffset; - ImageType3D::PointType m_TUserLPS = - m_IECtoLPSDirections * m_TranslationUser; +// ImageType3D::PointType m_TUserLPS = +// m_IECtoLPSDirections * m_TranslationUser; - ImageType3D::PointType m_RUserLPS = - m_IECtoLPSDirections * m_RotationUser; +// ImageType3D::PointType m_RUserLPS = +// m_IECtoLPSDirections * m_RotationUser; - TransformType::OutputVectorType translation; - translation[0] = m_TOffsetLPS[0] + m_TUserLPS[0]; - translation[1] = m_TOffsetLPS[1] + m_TUserLPS[1]; - translation[2] = m_TOffsetLPS[2] + m_TUserLPS[2]; +// TransformType::OutputVectorType translation; +// translation[0] = m_TOffsetLPS[0] + m_TUserLPS[0]; +// translation[1] = m_TOffsetLPS[1] + m_TUserLPS[1]; +// translation[2] = m_TOffsetLPS[2] + m_TUserLPS[2]; - TransformType::OutputVectorType rotations; - rotations[0] = m_ROffsetLPS[0] + m_RUserLPS[0]; - rotations[1] = m_ROffsetLPS[1] + m_RUserLPS[1]; - rotations[2] = m_ROffsetLPS[2] + m_RUserLPS[2]; +// TransformType::OutputVectorType rotations; +// rotations[0] = m_ROffsetLPS[0] + m_RUserLPS[0]; +// rotations[1] = m_ROffsetLPS[1] + m_RUserLPS[1]; +// rotations[2] = m_ROffsetLPS[2] + m_RUserLPS[2]; - // Map offset to the projection center - TransformType::Pointer m_outputTransform = - MapTransformToNewOrigin ( - m_CalibratedProjectionCenter - m_RTIsocenter, - translation, - rotations - ); +// // Map offset to the projection center +// TransformType::Pointer m_outputTransform = +// MapTransformToNewOrigin ( +// m_CalibratedProjectionCenter - m_RTIsocenter, +// translation, +// rotations +// ); - m_outputTransform->SetCenter(m_CalibratedProjectionCenter); +// m_outputTransform->SetCenter(m_CalibratedProjectionCenter); - return m_outputTransform; +// return m_outputTransform; -} +//} -itkImageProcessor::TransformType::Pointer -itkImageProcessor::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 - ) -{ +//itkImageProcessor::TransformType::Pointer +//itkImageProcessor::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 +// //Convert all inputs into LPS - ImageType3D::PointType m_TLPS = - m_IECtoLPSDirections * m_Translation; +// ImageType3D::PointType m_TLPS = +// m_IECtoLPSDirections * m_Translation; - ImageType3D::PointType m_RLPS = - m_IECtoLPSDirections * m_Rotation; +// 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 - ); +// // Map offset to the projection center +// TransformType::Pointer m_outputTransform = +// MapTransformToNewOrigin ( +// m_CalibratedProjectionCenter - m_RTIsocenter, +// m_TLPS, +// m_RLPS +// ); - m_outputTransform->SetCenter(m_CalibratedProjectionCenter); +// m_outputTransform->SetCenter(m_CalibratedProjectionCenter); - return m_outputTransform; +// return m_outputTransform; -} +//} void itkImageProcessor::InitializeRegistration( @@ -1528,94 +1527,109 @@ void itkImageProcessor::InitializeRegistration( //CalculateInternalTransform(transform1, m_DRTImage1MetaInfo); /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ1 *********/ - ImageType3D::PointType ZeroPoint; - ZeroPoint.Fill(0.); - InternalImageType::DirectionType IECtoLPS_Directions; - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); - TransformType::Pointer CurrTransform; - if(m_RTMetaInfo == NULL) - { - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions); +// ImageType3D::PointType ZeroPoint; +// ZeroPoint.Fill(0.); +// InternalImageType::DirectionType IECtoLPS_Directions; +// IECtoLPS_Directions = +// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); +// TransformType::Pointer CurrTransform; +// if(m_RTMetaInfo == NULL) +// { +// CurrTransform = +// CalculateInternalTransformV3( +// m_TransformMetaInfo->GetT(), +// m_TransformMetaInfo->GetR(), +// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), +// m_DRTImage1MetaInfo->GetProjectionOriginLPS(), +// IECtoLPS_Directions); - } else { +// } else { - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); +// CurrTransform = +// CalculateInternalTransformV3( +// m_TransformMetaInfo->GetT(), +// m_TransformMetaInfo->GetR(), +// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), +// m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), +// IECtoLPS_Directions +// ); - } +// } - transform1->SetComputeZYX(true); - transform1->SetIdentity(); - transform1->SetTranslation( - CurrTransform->GetTranslation()); - transform1->SetRotation( - CurrTransform->GetAngleX(), - CurrTransform->GetAngleY(), - CurrTransform->GetAngleZ() - ); - transform1->SetCenter( - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); +// transform1->SetComputeZYX(true); +// transform1->SetIdentity(); +// transform1->SetTranslation( +// CurrTransform->GetTranslation()); +// transform1->SetRotation( +// CurrTransform->GetAngleX(), +// CurrTransform->GetAngleY(), +// CurrTransform->GetAngleZ() +// ); +// transform1->SetCenter( +// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); +// /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ + + +// //CalculateInternalTransform(transform2, m_DRTImage2MetaInfo); +// /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/ +// if(m_RTMetaInfo == NULL) +// { +// CurrTransform = +// CalculateInternalTransformV3( +// m_TransformMetaInfo->GetT(), +// m_TransformMetaInfo->GetR(), +// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), +// m_DRTImage2MetaInfo->GetProjectionOriginLPS(), +// IECtoLPS_Directions); + + +// } else { + + +// CurrTransform = +// CalculateInternalTransformV3( +// m_TransformMetaInfo->GetT(), +// m_TransformMetaInfo->GetR(), +// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), +// m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), +// IECtoLPS_Directions +// ); + + +// } +// transform2->SetComputeZYX(true); +// transform2->SetIdentity(); +// transform2->SetTranslation( +// CurrTransform->GetTranslation()); +// transform2->SetRotation( +// CurrTransform->GetAngleX(), +// CurrTransform->GetAngleY(), +// CurrTransform->GetAngleZ() +// ); +// transform2->SetCenter( +// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ); /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ + IsocTransf = TransformType::New(); + IsocTransf->SetRotation( + m_TransformMetaInfo->GetR()[0], + m_TransformMetaInfo->GetR()[1], + m_TransformMetaInfo->GetR()[2] + ); - //CalculateInternalTransform(transform2, m_DRTImage2MetaInfo); - /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/ - if(m_RTMetaInfo == NULL) - { - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions); + TransformType::OutputVectorType TranslV; + TranslV[0] = m_TransformMetaInfo->GetT()[0]; + TranslV[1] = m_TransformMetaInfo->GetT()[1]; + TranslV[2] = m_TransformMetaInfo->GetT()[2]; + IsocTransf->SetTranslation(TranslV); - } else { - - - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); - - - } - transform2->SetComputeZYX(true); - transform2->SetIdentity(); - transform2->SetTranslation( - CurrTransform->GetTranslation()); - transform2->SetRotation( - CurrTransform->GetAngleX(), - CurrTransform->GetAngleY(), - CurrTransform->GetAngleZ() - ); - transform2->SetCenter( - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ); - /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ - + registration->SetIsocIECTransform(IsocTransf); if (verbose) { registration->DebugOn(); @@ -2055,8 +2069,7 @@ int itkImageProcessor::unloadRTPlanAndMeta(){ std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl; m_RTMetaInfo = NULL; - ImageType3D::PointType - pZero; + ImageType3D::PointType pZero; pZero.Fill(0.); m_CTMetaInfo->SetImportOffset(pZero); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 72af105..e83945c 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -51,6 +51,8 @@ gfattori 08.11.2021 #include "itkQtIterationUpdate.h" +#include "itkImageProcessorHelpers.h" + namespace itk { @@ -271,29 +273,22 @@ private: /* Calculate the transform used in siddon. * The isocentric transform is mapped to the calibrated center of projection */ - TransformType::Pointer - CalculateInternalTransformV2( - ImageType3D::PointType m_TranslationOffset, - ImageType3D::PointType m_RotationOffset, - ImageType3D::PointType m_TranslationUser, - ImageType3D::PointType m_RotationUser, - ImageType3D::PointType m_CalibratedProjectionCenter, - ImageType3D::PointType m_RTIsocenter, - InternalImageType::DirectionType m_IECtoLPSDirections - ); +// TransformType::Pointer +// CalculateInternalTransformV2( +// ImageType3D::PointType m_TranslationOffset, +// ImageType3D::PointType m_RotationOffset, +// ImageType3D::PointType m_TranslationUser, +// ImageType3D::PointType m_RotationUser, +// ImageType3D::PointType m_CalibratedProjectionCenter, +// ImageType3D::PointType m_RTIsocenter, +// InternalImageType::DirectionType m_IECtoLPSDirections +// ); - TransformType::Pointer - CalculateInternalTransformV3( - ImageType3D::PointType m_Translation, - ImageType3D::PointType m_Rotation, - ImageType3D::PointType m_CalibratedProjectionCenter, - ImageType3D::PointType m_RTIsocenter, - InternalImageType::DirectionType m_IECtoLPSDirections - ); TransformType::Pointer transform1, - transform2; + transform2, + IsocTransf; InterpolatorType::Pointer interpolator1, @@ -355,12 +350,12 @@ private: ); - TransformType::Pointer - MapTransformToNewOrigin( - ImageType3D::PointType m_COR, - ImageType3D::PointType m_Translations, - ImageType3D::PointType m_Rotations - ); +// TransformType::Pointer +// MapTransformToNewOrigin( +// ImageType3D::PointType m_COR, +// ImageType3D::PointType m_Translations, +// ImageType3D::PointType m_Rotations +// ); double CalcProjectionAngleLPS( diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp new file mode 100644 index 0000000..018ecf1 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp @@ -0,0 +1,93 @@ +#include "itkImageProcessorHelpers.h" + +itk::TransformType::Pointer +itk::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; +} + + + +itk::TransformType::Pointer +itk::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; + +} diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h new file mode 100644 index 0000000..2b8aec8 --- /dev/null +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h @@ -0,0 +1,38 @@ + +#ifndef ITKIMAGEPROCESSORHELPERS_H +#define ITKIMAGEPROCESSORHELPERS_H + +#include "itkEuler3DTransform.h" +#include "itkImage.h" + + +namespace itk +{ +constexpr static unsigned int Dimension = 3; +using PixelType3D = short; +using InternalPixelType = float; + +using ImageType3D = itk::Image; +using InternalImageType = itk::Image; + +using TransformType = itk::Euler3DTransform; + +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 + ); +} + +#endif From d4c800dfcd65e867376c553d7add749fcac1521a Mon Sep 17 00:00:00 2001 From: Proton local user Date: Fri, 12 May 2023 23:36:24 +0200 Subject: [PATCH 07/12] R23 optimizes the isocentric IEC transform, internally the two internal transform are computed at every iteration. As a result, the feedback to regDialog is now already correct, no need for flips or sign change. CalculateExternalTransform becomes obsolete... and has been commented out. Methods for gui to get latest paremeters have been rewritten. --- .../itkDTRrecon/DRTMetaInformation.cpp | 23 + .../itkDTRrecon/DRTMetaInformation.h | 59 ++ .../itkTwoProjectionImageRegistrationMethod.h | 10 + ...tkTwoProjectionImageRegistrationMethod.hxx | 27 +- .../autoreg/itkgTwoImageToOneImageMetric.h | 13 + .../autoreg/itkgTwoImageToOneImageMetric.hxx | 104 ++- .../itkDTRrecon/itkImageProcessor.cpp | 613 ++++++------------ .../itkDTRrecon/itkImageProcessor.h | 43 +- .../itkDTRrecon/itkImageProcessorHelpers.cpp | 81 ++- .../itkDTRrecon/itkImageProcessorHelpers.h | 73 ++- .../itkDTRrecon/itkQtIterationUpdate.h | 14 +- 11 files changed, 586 insertions(+), 474 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 8395221..4981437 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -517,6 +517,29 @@ 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 (){ diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index f797b30..5beb57a 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -514,6 +514,65 @@ 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{ diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h index 10d638b..d389042 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -177,6 +177,12 @@ public: 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); @@ -259,6 +265,10 @@ private: InterpolatorPointer m_Interpolator1; InterpolatorPointer m_Interpolator2; + InternalTransformMetaInformation::Pointer + m_internalMeta1, + m_internalMeta2; + //TransformMetaInformation::Pointer // m_TransformMetaInfo; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx index fa127aa..b871033 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -40,6 +40,8 @@ TwoProjectionImageRegistrationMethod::TwoProjectionIm 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. @@ -144,6 +146,14 @@ void TwoProjectionImageRegistrationMethod::Initialize 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); @@ -152,6 +162,9 @@ void TwoProjectionImageRegistrationMethod::Initialize 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); @@ -206,7 +219,7 @@ void TwoProjectionImageRegistrationMethod::OnInternal // if (m_CancelRequest) // filter->SetAbortGenerateData(true); // may be handled by filter - std::cout << "OnInternalFilterProgressReceptor :: filter (TURE) " << std::endl; + std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl; } } @@ -217,22 +230,10 @@ template void TwoProjectionImageRegistrationMethod::StartRegistration() { -// typedef itk::MemberCommand ITKCommandType; -// typedef ITKCommandType::Pointer MemberPointer ; - - -// ITKCommandType::Pointer cmd = ITKCommandType::New(); -// itk::Command cmd= itk::Command::Pointer(this); - -// cmd->SetCallbackFunction(this, -// (ITKCommandType::Pointer)&TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor); - itk::ProgressReporter progress( this, 1, 1000,100); -// this->SetAbortGenerateData(true); - ParametersType empty(1); empty.Fill(0.0); try { diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h index 752f7c1..2c08287 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h @@ -40,6 +40,8 @@ #include #include +#include + namespace itk { /** \class gTwoImageToOneImageMetric @@ -194,6 +196,13 @@ public: /** 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); @@ -283,6 +292,10 @@ protected: R23MetaInformation::Pointer m_TransformMetaInfo; + InternalTransformMetaInformation::Pointer + m_internalMeta1, + m_internalMeta2; + ChangeInformationFilterPointer m_Filter1, m_Filter2; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index 8cade2d..b2c2f2b 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -32,6 +32,8 @@ gTwoImageToOneImageMetric::gTwoImageToOneImageMetric( 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; @@ -55,6 +57,13 @@ bool gTwoImageToOneImageMetric::SetTransformParameter 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"); + } + //auto transformParameters = m_Transform1->GetParameters(); auto transformParameters = m_IsocTransf->GetParameters(); double TranslationAlongX = transformParameters[3]; @@ -64,6 +73,8 @@ bool gTwoImageToOneImageMetric::SetTransformParameter double RotationAlongY = transformParameters[1]; double RotationAlongZ = transformParameters[2]; + //std::cout << "m_IsocTransf PARS: "<< transformParameters <GetDegreeOfFreedom()) { case X_ONLY: TranslationAlongX = parameters[0]; @@ -102,33 +113,73 @@ bool gTwoImageToOneImageMetric::SetTransformParameter 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 << " 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; //HERE WE HAVE TO CALCULATE INTERNAL TRANSFORMS!!! // GIOVANNI - transformParameters[0] = RotationAlongX; - transformParameters[1] = RotationAlongY; - transformParameters[2] = RotationAlongZ; - transformParameters[3] = TranslationAlongX; - transformParameters[4] = TranslationAlongY; - transformParameters[5] = TranslationAlongZ; + ImageType3D::PointType + pT; + pT[0] = TranslationAlongX; + pT[1] = TranslationAlongY; + pT[2] = TranslationAlongZ; - bool transformValid = (std::abs(TranslationAlongX) < m_MaxTranslation) && - (std::abs(TranslationAlongY) < m_MaxTranslation) && - (std::abs(TranslationAlongZ) < m_MaxTranslation); + ImageType3D::PointType + pR; + pR[0] = RotationAlongX; + pR[1] = RotationAlongY; + pR[2] = RotationAlongZ; - if (transformValid) { - m_Transform1->SetParameters(transformParameters); - m_Transform2->SetParameters(transformParameters); - } else { - std::cout << " Transform invalid, out of range!" << std::endl; - std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; - } + 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 << "m_Transform1 PARS: "<< m_Transform1->GetParameters()<GetpCalProjCenter(), + m_internalMeta2->GetpRTIsocenter(), + m_internalMeta2->GetIECtoLPSDirs()); + + m_Transform2->SetIdentity(); + m_Transform2->SetParameters( + CurrTransform->GetParameters()); + m_Transform2->SetCenter( + m_internalMeta2->GetpCalProjCenter()); + // std::cout << "m_Transform2 PARS: "<< m_Transform2->GetParameters()<SetParameters(transformParameters); +// m_Transform2->SetParameters(transformParameters); +// } else { +// std::cout << " Transform invalid, out of range!" << std::endl; +// std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; +// } return transformValid; } @@ -136,18 +187,25 @@ template void gTwoImageToOneImageMetric::Initialize() { - if (!m_Transform1) { + if (!m_IsocTransf) { itkExceptionMacro(<< "Transform is not present"); } if (!m_Transform1) { - itkExceptionMacro(<< "Transform is not present"); + itkExceptionMacro(<< "Transform1 is not present"); } if (!m_Transform2) { - itkExceptionMacro(<< "Transform is not present"); + 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"); } @@ -244,7 +302,7 @@ unsigned int gTwoImageToOneImageMetric::GetNumberOfPa break; } - return m_Transform1->GetNumberOfParameters(); + return m_IsocTransf->GetNumberOfParameters(); } template diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 0c2c739..494217f 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -68,6 +68,10 @@ itkImageProcessor::itkImageProcessor() transform2 = TransformType::New(); transform2->SetIdentity(); + m_InternalTransf1 = InternalTransformMetaInformation::New(); + m_InternalTransf2 = InternalTransformMetaInformation::New(); + + interpolator1 = InterpolatorType::New(); interpolator2 = InterpolatorType::New(); @@ -1281,174 +1285,33 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg) } -//itkImageProcessor::TransformType::Pointer -//itkImageProcessor::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; -//} // This External User Transform thing need to be checked out -void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo) -{ - //crude hack to get from internal transform in LPS to external transform in IEC. - //result is stored in m_TransformMetaInfo - - //TODO: support for projectionTransformCenter and userTransformCenter which are not the same - //currentlz we support only pFakeIsoLPS situations. - - InternalImageType::DirectionType LPStoIEC_Directions; - LPStoIEC_Directions = m_CTMetaInfo->GetLPS2IECDirections(); - - auto parameters = transform->GetParameters(); - - ImageType3D::PointType translationUser; - translationUser[0] = parameters[3]; - translationUser[1] = parameters[4]; - translationUser[2] = parameters[5]; - ImageType3D::PointType rotationUser; - rotationUser[0] = parameters[0]; - rotationUser[1] = parameters[1]; - rotationUser[2] = parameters[2]; - - m_TransformMetaInfo->SetUserTranslations(LPStoIEC_Directions * translationUser); - m_TransformMetaInfo->SetUserRotations(LPStoIEC_Directions * rotationUser); -} - - -//itkImageProcessor::TransformType::Pointer -//itkImageProcessor::CalculateInternalTransformV2( -// ImageType3D::PointType m_TranslationOffset, //IEC -// ImageType3D::PointType m_RotationOffset, //IEC -// ImageType3D::PointType m_TranslationUser, //IEC -// ImageType3D::PointType m_RotationUser, //IEC -// ImageType3D::PointType m_CalibratedProjectionCenter, //LPS -// ImageType3D::PointType m_RTIsocenter, //LPS -// InternalImageType::DirectionType m_IECtoLPSDirections -// ) +//void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo) //{ +// //crude hack to get from internal transform in LPS to external transform in IEC. +// //result is stored in m_TransformMetaInfo -// //Convert all inputs into LPS +// //TODO: support for projectionTransformCenter and userTransformCenter which are not the same +// //currentlz we support only pFakeIsoLPS situations. -// ImageType3D::PointType m_TOffsetLPS = -// m_IECtoLPSDirections * m_TranslationOffset; +// InternalImageType::DirectionType LPStoIEC_Directions; +// LPStoIEC_Directions = m_CTMetaInfo->GetLPS2IECDirections(); -// ImageType3D::PointType m_ROffsetLPS = -// m_IECtoLPSDirections * m_RotationOffset; +// auto parameters = transform->GetParameters(); -// ImageType3D::PointType m_TUserLPS = -// m_IECtoLPSDirections * m_TranslationUser; - -// ImageType3D::PointType m_RUserLPS = -// m_IECtoLPSDirections * m_RotationUser; - -// TransformType::OutputVectorType translation; -// translation[0] = m_TOffsetLPS[0] + m_TUserLPS[0]; -// translation[1] = m_TOffsetLPS[1] + m_TUserLPS[1]; -// translation[2] = m_TOffsetLPS[2] + m_TUserLPS[2]; - -// TransformType::OutputVectorType rotations; -// rotations[0] = m_ROffsetLPS[0] + m_RUserLPS[0]; -// rotations[1] = m_ROffsetLPS[1] + m_RUserLPS[1]; -// rotations[2] = m_ROffsetLPS[2] + m_RUserLPS[2]; - -// // Map offset to the projection center -// TransformType::Pointer m_outputTransform = -// MapTransformToNewOrigin ( -// m_CalibratedProjectionCenter - m_RTIsocenter, -// translation, -// rotations -// ); - -// m_outputTransform->SetCenter(m_CalibratedProjectionCenter); - -// return m_outputTransform; - -//} - - -//itkImageProcessor::TransformType::Pointer -//itkImageProcessor::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; +// ImageType3D::PointType translationUser; +// translationUser[0] = parameters[3]; +// translationUser[1] = parameters[4]; +// translationUser[2] = parameters[5]; +// ImageType3D::PointType rotationUser; +// rotationUser[0] = parameters[0]; +// rotationUser[1] = parameters[1]; +// rotationUser[2] = parameters[2]; +// m_TransformMetaInfo->SetUserTranslations(LPStoIEC_Directions * translationUser); +// m_TransformMetaInfo->SetUserRotations(LPStoIEC_Directions * rotationUser); //} @@ -1521,103 +1384,19 @@ void itkImageProcessor::InitializeRegistration( // registration->SetTransformMetaInfo(m_TransformMetaInfo); registration->SetTransformMetaInfo(m_r23MetaInfo); + registration->SetinternalMeta1(m_InternalTransf1); + registration->SetinternalMeta2(m_InternalTransf2); registration->SetFilter1(filter1); registration->SetFilter2(filter2); - //CalculateInternalTransform(transform1, m_DRTImage1MetaInfo); - /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ1 *********/ -// ImageType3D::PointType ZeroPoint; -// ZeroPoint.Fill(0.); -// InternalImageType::DirectionType IECtoLPS_Directions; -// IECtoLPS_Directions = -// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); -// TransformType::Pointer CurrTransform; -// if(m_RTMetaInfo == NULL) -// { -// CurrTransform = -// CalculateInternalTransformV3( -// m_TransformMetaInfo->GetT(), -// m_TransformMetaInfo->GetR(), -// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), -// m_DRTImage1MetaInfo->GetProjectionOriginLPS(), -// IECtoLPS_Directions); - - -// } else { - -// CurrTransform = -// CalculateInternalTransformV3( -// m_TransformMetaInfo->GetT(), -// m_TransformMetaInfo->GetR(), -// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), -// m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), -// IECtoLPS_Directions -// ); - - - -// } - - -// transform1->SetComputeZYX(true); -// transform1->SetIdentity(); -// transform1->SetTranslation( -// CurrTransform->GetTranslation()); -// transform1->SetRotation( -// CurrTransform->GetAngleX(), -// CurrTransform->GetAngleY(), -// CurrTransform->GetAngleZ() -// ); -// transform1->SetCenter( -// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); - -// /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ - - -// //CalculateInternalTransform(transform2, m_DRTImage2MetaInfo); -// /********* BEGIN CALCULATION OF INTERNAL TRANSFORM FOR PROJ2 *********/ -// if(m_RTMetaInfo == NULL) -// { -// CurrTransform = -// CalculateInternalTransformV3( -// m_TransformMetaInfo->GetT(), -// m_TransformMetaInfo->GetR(), -// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), -// m_DRTImage2MetaInfo->GetProjectionOriginLPS(), -// IECtoLPS_Directions); - - -// } else { - - -// CurrTransform = -// CalculateInternalTransformV3( -// m_TransformMetaInfo->GetT(), -// m_TransformMetaInfo->GetR(), -// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), -// m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), -// IECtoLPS_Directions -// ); - - -// } -// transform2->SetComputeZYX(true); -// transform2->SetIdentity(); -// transform2->SetTranslation( -// CurrTransform->GetTranslation()); -// transform2->SetRotation( -// CurrTransform->GetAngleX(), -// CurrTransform->GetAngleY(), -// CurrTransform->GetAngleZ() -// ); -// transform2->SetCenter( -// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ); - /********* END OF CALCULATE INTERNAL TRANSFORM FOR PROJ1 *********/ IsocTransf = TransformType::New(); + IsocTransf->SetComputeZYX(true); + IsocTransf->SetIdentity(); + IsocTransf->SetRotation( - m_TransformMetaInfo->GetR()[0], + m_TransformMetaInfo->GetR()[0], m_TransformMetaInfo->GetR()[1], m_TransformMetaInfo->GetR()[2] ); @@ -1690,8 +1469,8 @@ void itkImageProcessor::InitializeRegistration( const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651. const double stepLength = 0.1; - auto parameters = transform1->GetParameters(); - transform1->SetParameters(parameters); + //auto parameters = transform1->GetParameters(); + //transform1->SetParameters(parameters); ExhaustiveOptimizerType::StepsType steps(3); steps[0] = numberOfSteps; @@ -1717,7 +1496,7 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) // TODO: Check if the registartion pipeline has been initialized using ParametersType = RegistrationType::ParametersType; - auto startParameters = transform1->GetParameters(); + //auto startParameters = transform1->GetParameters(); time_t t = time(0); // get time now struct tm* now = localtime(&t); @@ -1784,54 +1563,106 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) fs.close(); + + m_OptmizerValue = optimizer->GetValue(); + auto finalParameters = optimizer->GetCurrentPosition(); - auto finalParameters = transform1->GetParameters(); - std::cout<<"-->startParameters "<finalParameters1 "<GetParameters()<finalParameters2 "<GetParameters()<GetCurrentPosition(); - finalParameters = finalParameters - startParameters; - - std::cout<<"startParameters "<GetCenter()<GetTranslation()<GetCenter()<GetTranslation()<SetParameters(finalParameters); + // auto finalParameters = transform1->GetParameters(); + // std::cout<<"-->startParameters "<finalParameters1 "<GetParameters()<finalParameters2 "<GetParameters()<GetCurrentPosition(); + // finalParameters = finalParameters - startParameters; - const double translationAlongX = finalParameters[3]; - const double translationAlongY = finalParameters[4]; - const double translationAlongZ = finalParameters[5]; - const double rotationAlongX = finalParameters[0]; - const double rotationAlongY = finalParameters[1]; - const double rotationAlongZ = finalParameters[2]; + // std::cout<<"startParameters "<GetCenter()<GetTranslation()<GetCenter()<GetTranslation()<GetCurrentIteration(); - std::cout << "Result = " << std::endl; - std::cout << " Rotation Along X = " << rotationAlongX / dtr << " deg" << std::endl; - std::cout << " Rotation Along Y = " << rotationAlongY / dtr << " deg" << std::endl; - std::cout << " Rotation Along Z = " << rotationAlongZ / dtr << " deg" << std::endl; - std::cout << " Translation X = " << translationAlongX << " mm" << std::endl; - std::cout << " Translation Y = " << translationAlongY << " mm" << std::endl; - std::cout << " Translation Z = " << translationAlongZ << " mm" << std::endl; - std::cout << " Number Of Iterations = " << numberOfIterations << std::endl; - std::cout << " Metric value = " << m_OptmizerValue << std::endl; + // TransformType::Pointer offsetTransform = TransformType::New(); + // offsetTransform->SetParameters(finalParameters); + + // CalculateExternalUserTransform(offsetTransform, m_DRTImage1MetaInfo); + + // const double translationAlongX = finalParameters[3]; + // const double translationAlongY = finalParameters[4]; + // const double translationAlongZ = finalParameters[5]; + // const double rotationAlongX = finalParameters[0]; + // const double rotationAlongY = finalParameters[1]; + // const double rotationAlongZ = finalParameters[2]; + + // const int numberOfIterations = optimizer->GetCurrentIteration(); + + // std::cout << "Result = " << std::endl; + // std::cout << " Rotation Along X = " << rotationAlongX / dtr << " deg" << std::endl; + // std::cout << " Rotation Along Y = " << rotationAlongY / dtr << " deg" << std::endl; + // std::cout << " Rotation Along Z = " << rotationAlongZ / dtr << " deg" << std::endl; + // std::cout << " Translation X = " << translationAlongX << " mm" << std::endl; + // std::cout << " Translation Y = " << translationAlongY << " mm" << std::endl; + // std::cout << " Translation Z = " << translationAlongZ << " mm" << std::endl; + // std::cout << " Number Of Iterations = " << numberOfIterations << std::endl; + // std::cout << " Metric value = " << m_OptmizerValue << std::endl; std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl; return 0; } +Optimizer::ParametersType +itkImageProcessor::GetFinalR23Parameters(){ + + Optimizer::ParametersType pPars(6); + pPars.fill(0.); + + if(optimizer == nullptr){ + return pPars; + } + + switch (m_r23MetaInfo->GetDegreeOfFreedom()) { + case THREE_DOF: + pPars[3] = optimizer->GetCurrentPosition()[0] + - m_TransformMetaInfo->GetHiddenTranslations()[0]; + pPars[4] = optimizer->GetCurrentPosition()[1] + - m_TransformMetaInfo->GetHiddenTranslations()[1]; + pPars[5] = optimizer->GetCurrentPosition()[2] + - m_TransformMetaInfo->GetHiddenTranslations()[2]; + break; + + case SIX_DOF: + pPars[0] = optimizer->GetCurrentPosition()[0] + - m_TransformMetaInfo->GetHiddenRotations()[0]; + pPars[1] = optimizer->GetCurrentPosition()[1] + - m_TransformMetaInfo->GetHiddenRotations()[1]; + pPars[2] = optimizer->GetCurrentPosition()[2] + - m_TransformMetaInfo->GetHiddenRotations()[2]; + pPars[3] = optimizer->GetCurrentPosition()[3] + - m_TransformMetaInfo->GetHiddenTranslations()[0]; + pPars[4] = optimizer->GetCurrentPosition()[4] + - m_TransformMetaInfo->GetHiddenTranslations()[1]; + pPars[5] = optimizer->GetCurrentPosition()[5] + - m_TransformMetaInfo->GetHiddenTranslations()[2]; + break; + + default: + pPars.fill(0.); + break; + } + + + return + pPars; + +} + void itkImageProcessor::InitializeProjector() { @@ -1851,37 +1682,19 @@ void itkImageProcessor::InitializeProjector() ZeroPoint.Fill(0.); - InternalImageType::DirectionType IECtoLPS_Directions; + // InternalImageType::DirectionType IECtoLPS_Directions; - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + // IECtoLPS_Directions = + // m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); TransformType::Pointer CurrTransform; - if(m_RTMetaInfo == NULL) - { - - CurrTransform = CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions); - - - } else { - - - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); - - } + CurrTransform= CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf1->GetpCalProjCenter(), + m_InternalTransf1->GetpRTIsocenter(), + m_InternalTransf1->GetIECtoLPSDirs()); transform1->SetComputeZYX(true); @@ -1912,30 +1725,12 @@ void itkImageProcessor::InitializeProjector() interpolator1->SetTransform(transform1); interpolator1->Initialize(); - if(m_RTMetaInfo == NULL) - { - - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions); - - - - } else { - - CurrTransform = CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); - - } + CurrTransform= CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf2->GetpCalProjCenter(), + m_InternalTransf2->GetpRTIsocenter(), + m_InternalTransf2->GetIECtoLPSDirs()); transform2->SetComputeZYX(true); @@ -2258,6 +2053,29 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ Std_DRT2LPS) ); + if(m_RTMetaInfo == NULL) + { + m_InternalTransf1->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 @@ -2347,6 +2165,29 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ Std_DRT2LPS) ); + if(m_RTMetaInfo == NULL) + { + m_InternalTransf2->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); + + } + } @@ -2398,35 +2239,13 @@ void itkImageProcessor::GetProjectionImages(){ transform1->SetComputeZYX(true); transform1->SetIdentity(); - InternalImageType::DirectionType IECtoLPS_Directions; - - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); TransformType::Pointer CurrTransform; - - - if(m_RTMetaInfo == NULL) - { - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions); - - } else { - - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); - - } + CurrTransform= CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf1->GetpCalProjCenter(), + m_InternalTransf1->GetpRTIsocenter(), + m_InternalTransf1->GetIECtoLPSDirs()); transform1->SetComputeZYX(true); @@ -2460,32 +2279,12 @@ void itkImageProcessor::GetProjectionImages(){ transform2->SetComputeZYX(true); transform2->SetIdentity(); - - - if(m_RTMetaInfo == NULL) - { - - CurrTransform = - CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - IECtoLPS_Directions); - - } else { - - CurrTransform = CalculateInternalTransformV3( - m_TransformMetaInfo->GetT(), - m_TransformMetaInfo->GetR(), - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS(), - IECtoLPS_Directions - ); - - - } - + CurrTransform = CalculateInternalTransformV3( + m_TransformMetaInfo->GetT(), + m_TransformMetaInfo->GetR(), + m_InternalTransf2->GetpCalProjCenter(), + m_InternalTransf2->GetpRTIsocenter(), + m_InternalTransf2->GetIECtoLPSDirs()); transform2->SetComputeZYX(true); @@ -3059,26 +2858,26 @@ double* itkImageProcessor::GetTransformParameters(){ return dTransfParam; } -// THESE ARE CALLED AT THE END OF REG -void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) -{ - ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations(); +//// THESE ARE CALLED AT THE END OF REG +//void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) +//{ +// ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations(); - dX = userTranslations[0]; - dY = userTranslations[1]; - dZ = userTranslations[2]; -} +// dX = userTranslations[0]; +// dY = userTranslations[1]; +// dZ = userTranslations[2]; +//} -// THESE ARE CALLED AT THE END OF REG -void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ) -{ - ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations(); +//// THESE ARE CALLED AT THE END OF REG +//void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ) +//{ +// ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations(); - dRX = userRotations[0]; - dRY = userRotations[1]; - dRZ = userRotations[2]; +// dRX = userRotations[0]; +// dRY = userRotations[1]; +// dRZ = userRotations[2]; -} +//} double itkImageProcessor::GetOptimizerValue() { diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index e83945c..285e756 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -139,13 +139,17 @@ public: void SetInitialTranslations(double, double, double); void SetInitialRotations(double, double, double); - /** Get transform parameters for 3D Volume */ - void GetFinalTranslations(double&, double&, double&); - void GetFinalRotations(double&, double&, double&); - /** Set user transform parameters for 3D Volume in LPS*/ - void SetUserTranslationsLPS(double, double, double); - void SetUserRotationsLPS(double, double, double); + Optimizer::ParametersType + GetFinalR23Parameters(); + + /** Get transform parameters for 3D Volume */ +// void GetFinalTranslations(double&, double&, double&); +// void GetFinalRotations(double&, double&, double&); + +// /** Set user transform parameters for 3D Volume in LPS*/ +// void SetUserTranslationsLPS(double, double, double); +// void SetUserRotationsLPS(double, double, double); /** Initialize the registration pipeline*/ void InitializeRegistration(double, double, eDegreeOfFreedomType); @@ -238,6 +242,8 @@ private: using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; using ResampleFilterType = itk::ResampleImageFilter; + //using gTransformType = itk::gEuler3DTransform; + /** Optimizer which tries to find the minimn (Powell Optimizer)*/ using AmoebaOptimizerType = itk::AmoebaOptimizer; /** Optimizer which samples the whol space */ @@ -266,23 +272,10 @@ private: /** ITK to VTK filter */ using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; - void - CalculateExternalUserTransform( - TransformType::Pointer transform, - DRTImageMetaInformation::Pointer imageMetaInfo); - - /* Calculate the transform used in siddon. - * The isocentric transform is mapped to the calibrated center of projection */ -// TransformType::Pointer -// CalculateInternalTransformV2( -// ImageType3D::PointType m_TranslationOffset, -// ImageType3D::PointType m_RotationOffset, -// ImageType3D::PointType m_TranslationUser, -// ImageType3D::PointType m_RotationUser, -// ImageType3D::PointType m_CalibratedProjectionCenter, -// ImageType3D::PointType m_RTIsocenter, -// InternalImageType::DirectionType m_IECtoLPSDirections -// ); +// void +// CalculateExternalUserTransform( +// TransformType::Pointer transform, +// DRTImageMetaInformation::Pointer imageMetaInfo); TransformType::Pointer @@ -419,6 +412,10 @@ private: R23MetaInformation::Pointer m_r23MetaInfo; + InternalTransformMetaInformation::Pointer + m_InternalTransf1, + m_InternalTransf2; + double m_OptmizerValue; int m_MaxNumberOfIterations; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp index 018ecf1..730ceb6 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp @@ -1,7 +1,78 @@ #include "itkImageProcessorHelpers.h" -itk::TransformType::Pointer -itk::MapTransformToNewOrigin( +namespace itk { + +//template +//gEuler3DTransform ::gEuler3DTransform() +//{ +// //m_CalibPrjCenter = nullptr; // has to be provided by the user. +// //m_RTIso = nullptr; // has to be provided by the user. +// // m_IECtoLPSDirs =; // has to be provided by the user. +//} + + +//template +//gEuler3DTransform::gEuler3DTransform +//(const MatrixType & matrix, const OutputPointType & offset) +//{ + + +//} + +//template +//gEuler3DTransform::gEuler3DTransform +//(unsigned int parametersDimension) +// : Superclass(parametersDimension) + +//{ + +//} + +//template +//void gEuler3DTransform ::SetParameters(const ParametersType & parameters) +//{ + +//// if(m_CalibPrjCenter != nullptr && +//// m_RTIso != nullptr //&& +//// //m_IECtoLPSDirs != nullptr +//// ){ + +//// itk::Point pT; +//// pT[0] = parameters[3]; +//// pT[1] = parameters[4]; +//// pT[2] = parameters[5]; + +//// itk::Point pR; +//// pR[0] = parameters[0]; +//// pR[1] = parameters[1]; +//// pR[2] = parameters[2]; + +//// TransformType::Pointer m_OutputTransform = +//// CalculateInternalTransformV3( +//// pT, +//// pR, +//// m_CalibPrjCenter, +//// m_RTIso, +//// m_IECtoLPSDirs); +//// Superclass::SetParameters(m_OutputTransform->GetParameters()); +//// }else{ +// Superclass::SetParameters(parameters); +//// } + + +//} + +//template +//void gEuler3DTransform ::PrintSelf(std::ostream& os, Indent indent) const +//{ +// Superclass::PrintSelf(os, indent); +// //os << indent << "ComputeGradient: " << static_cast::PrintType>(m_ComputeGradient) +// // << std::endl; +//} + + +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 @@ -60,8 +131,8 @@ itk::MapTransformToNewOrigin( -itk::TransformType::Pointer -itk::CalculateInternalTransformV3( +TransformType::Pointer +CalculateInternalTransformV3( ImageType3D::PointType m_Translation, //IEC ImageType3D::PointType m_Rotation, //IEC ImageType3D::PointType m_CalibratedProjectionCenter, //LPS @@ -91,3 +162,5 @@ itk::CalculateInternalTransformV3( return m_outputTransform; } + +} diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h index 2b8aec8..aa8a355 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h @@ -4,10 +4,80 @@ #include "itkEuler3DTransform.h" #include "itkImage.h" +#include "itkEuler3DTransform.h" +#include "itkCommand.h" +#include "itkObject.h" +#include "itkObjectFactory.h" namespace itk { + + + + + +//template +//class ITK_TEMPLATE_EXPORT gEuler3DTransform : +// public Euler3DTransform +//{ + +//public: +// ITK_DISALLOW_COPY_AND_ASSIGN(gEuler3DTransform); +// /** Standard class type aliases. */ +// using Self = gEuler3DTransform; +// using Superclass = Euler3DTransform; +// using Pointer = SmartPointer; +// using ConstPointer = SmartPointer; + +// /** New macro for creation of through a Smart Pointer. */ +// itkNewMacro(Self); + +// /** Run-time type information (and related methods). */ +// itkTypeMacro(gEuler3DTransform, Euler3DTransform); + +// // static constexpr unsigned int ParametersDimension = 6; + +// using MatrixType = typename Superclass::MatrixType; +// using OutputPointType = typename Superclass::OutputPointType; +// using ParametersType = typename Superclass::ParametersType; +// using PointType = typename itk::Point; + +// /** Set/Get the transformation from a container of parameters +// * This is typically used by optimizers. There are 6 parameters. The first +// * three represent the angles to rotate around the coordinate axis, and the +// * last three represents the offset. */ +// void +// SetParameters(const ParametersType & parameters) override; + +// itkSetMacro (CalibPrjCenter, PointType); +// itkSetMacro (RTIso, PointType); +// itkSetMacro(IECtoLPSDirs,MatrixType); + + +//protected: +// gEuler3DTransform(const MatrixType & matrix, const OutputPointType & offset); +// gEuler3DTransform(unsigned int parametersDimension); +// gEuler3DTransform(); + +// ~gEuler3DTransform() override = default; + +// PointType +// m_CalibPrjCenter, +// m_RTIso; + +// MatrixType +// m_IECtoLPSDirs; + + +// void +// PrintSelf(std::ostream & os, Indent indent) const override; + +//}; + + + +using TransformType = itk::Euler3DTransform; constexpr static unsigned int Dimension = 3; using PixelType3D = short; using InternalPixelType = float; @@ -15,7 +85,6 @@ using InternalPixelType = float; using ImageType3D = itk::Image; using InternalImageType = itk::Image; -using TransformType = itk::Euler3DTransform; TransformType::Pointer MapTransformToNewOrigin( @@ -33,6 +102,8 @@ TransformType::Pointer ImageType3D::PointType m_RTIsocenter, InternalImageType::DirectionType m_IECtoLPSDirections ); + + } #endif diff --git a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h index 91e982a..89b231d 100644 --- a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h +++ b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h @@ -118,11 +118,17 @@ public: 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()[2], - -optimizer->GetCurrentPosition()[1] + optimizer->GetCurrentPosition()[1], + optimizer->GetCurrentPosition()[2] ); } return; @@ -172,7 +178,9 @@ public: 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[2] << "\t" << -position[1] << std::endl; + *os << "\t" << position[0] << "\t" << position[1] << "\t" << position[2] << std::endl; + } return; } From 9b616b3b550f4001abe330f155d6e622c8d1424d Mon Sep 17 00:00:00 2001 From: Proton local user Date: Mon, 15 May 2023 17:02:18 +0200 Subject: [PATCH 08/12] a number of autoREG and geometry fixes - automatic registration actually makes use of ROI. Before were ignored. - cleanup of user and hidden transform and pass to UI after registration. - change of internal image interpolator used in autoREG. Now autoREG performs well. - general number of working units variable. - GetFinalRotations, GetFinalTranslations and GetTransformParameters become obsolete. - GetFinalR23Parameters and GetCompleteIsocentricTransform and GetUserIsocentricTransform are provided. - only GetCompleteIsocentricTransform is missing implementation. --- ...ualInformationTwoImageToOneImageMetric.hxx | 58 ++++-- .../itkTwoProjectionImageRegistrationMethod.h | 7 +- ...tkTwoProjectionImageRegistrationMethod.hxx | 43 +++-- .../autoreg/itkgTwoImageToOneImageMetric.hxx | 114 ++++++----- .../itkDTRrecon/itkImageProcessor.cpp | 181 ++++++++++-------- .../itkDTRrecon/itkImageProcessor.h | 34 ++-- .../itkDTRrecon/itkImageProcessorHelpers.cpp | 68 ------- .../itkDTRrecon/itkImageProcessorHelpers.h | 63 ------ 8 files changed, 250 insertions(+), 318 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx index 089c7f9..561225f 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx @@ -93,33 +93,61 @@ MutualInformationTwoImageToOneImageMetric::GetValue() 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; - - using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction; + /* + * 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 = NearestNeighborType::New(); + auto dummyInterpolator1 = LinearInterpType::New(); auto metric1 = InternalMetricType::New(); metric1->SetTransform(dummyTransform); metric1->SetTransformParameters(dummyParameters); metric1->SetInterpolator(dummyInterpolator1); - auto fixedImageRegion1 = fixedImage1->GetBufferedRegion(); - metric1->SetFixedImageRegion(fixedImageRegion1); + //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% - metric1->UseAllPixelsOn(); - metric1->SetNumberOfHistogramBins(30); + // since we have a ROI, then we should not set allPixels to TRUE. + //metric1->UseAllPixelsOn(); + metric1->SetNumberOfHistogramBins(50); metric1->SetFixedImage(fixedImage1); metric1->SetMovingImage(this->m_Filter1->GetOutput()); @@ -132,26 +160,26 @@ MutualInformationTwoImageToOneImageMetric::GetValue() this->m_Interpolator2->Initialize(); this->m_Filter2->Update(); - auto dummyInterpolator2 = NearestNeighborType::New(); + auto dummyInterpolator2 = LinearInterpType::New(); auto metric2 = InternalMetricType::New(); metric2->SetTransform(dummyTransform); metric2->SetTransformParameters(dummyParameters); metric2->SetInterpolator(dummyInterpolator2); - auto fixedImageRegion2 = fixedImage1->GetBufferedRegion(); - metric2->SetFixedImageRegion(fixedImageRegion2); +// 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(30); - + //metric2->SetNumberOfSpatialSamples(numberOfSamples); + //metric2->UseAllPixelsOn(); + metric2->SetNumberOfHistogramBins(50); metric2->SetFixedImage(fixedImage2); metric2->SetMovingImage(this->m_Filter2->GetOutput()); - metric2->Initialize(); auto measure2 = metric2->GetValue(dummyParameters); diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h index d389042..9ddf934 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h @@ -246,8 +246,8 @@ protected: itkSetMacro(LastOptimizerParameters, ParametersType); /** Entry-point for internal ITK filter observer. **/ - void OnInternalFilterProgressReceptor(itk::Object *caller, - const itk::EventObject &event); +// void OnInternalFilterProgressReceptor(itk::Object *caller, +// const itk::EventObject &event); private: @@ -269,9 +269,6 @@ private: m_internalMeta1, m_internalMeta2; - //TransformMetaInformation::Pointer - // m_TransformMetaInfo; - R23MetaInformation::Pointer m_TransformMetaInfo; diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx index b871033..2c9b6e0 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx @@ -160,6 +160,7 @@ void TwoProjectionImageRegistrationMethod::Initialize 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); @@ -197,31 +198,31 @@ void TwoProjectionImageRegistrationMethod::Initialize m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters); } -template -void TwoProjectionImageRegistrationMethod::OnInternalFilterProgressReceptor( - itk::Object *caller, const itk::EventObject &event) -{ - itk::ProcessObject *filter = dynamic_cast(caller); +//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(!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 (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 +//// if (m_CancelRequest) +//// filter->SetAbortGenerateData(true); // may be handled by filter - std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl; - } -} +// std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl; +// } +//} /* * Starts the Registration Process diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index b2c2f2b..d8b0cea 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -56,13 +56,15 @@ bool gTwoImageToOneImageMetric::SetTransformParameter 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"); - } + 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(); @@ -109,7 +111,7 @@ bool gTwoImageToOneImageMetric::SetTransformParameter break; } - std::cout << "New Transform Parameters = " << std::endl; + 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; @@ -119,67 +121,72 @@ bool gTwoImageToOneImageMetric::SetTransformParameter std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; - //HERE WE HAVE TO CALCULATE INTERNAL TRANSFORMS!!! - // GIOVANNI - + // individual R and T components of the IsoC transform ImageType3D::PointType - pT; + pT; pT[0] = TranslationAlongX; pT[1] = TranslationAlongY; pT[2] = TranslationAlongZ; ImageType3D::PointType - pR; + pR; pR[0] = RotationAlongX; pR[1] = RotationAlongY; pR[2] = RotationAlongZ; - TransformType::Pointer CurrTransform; - CurrTransform = CalculateInternalTransformV3( - pT, - pR, - m_internalMeta1->GetpCalProjCenter(), - m_internalMeta1->GetpRTIsocenter(), - m_internalMeta1->GetIECtoLPSDirs()); + // 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); - m_Transform1->SetIdentity(); - m_Transform1->SetParameters( - CurrTransform->GetParameters()); - m_Transform1->SetCenter( + 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 << "m_Transform1 PARS: "<< m_Transform1->GetParameters()<GetpCalProjCenter(), - m_internalMeta2->GetpRTIsocenter(), - m_internalMeta2->GetIECtoLPSDirs()); - m_Transform2->SetIdentity(); - m_Transform2->SetParameters( - CurrTransform->GetParameters()); - m_Transform2->SetCenter( - m_internalMeta2->GetpCalProjCenter()); - // std::cout << "m_Transform2 PARS: "<< m_Transform2->GetParameters()<GetCenter()<GetParameters()<SetParameters(transformParameters); -// m_Transform2->SetParameters(transformParameters); -// } else { -// std::cout << " Transform invalid, out of range!" << std::endl; -// std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl; -// } + // Calculate internal transform2 + CurrTransform = CalculateInternalTransformV3( + pT, + pR, + m_internalMeta2->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; } @@ -203,8 +210,8 @@ void gTwoImageToOneImageMetric::Initialize() itkExceptionMacro(<< "m_internalMeta1 is not present"); } if(!m_internalMeta2) { - itkExceptionMacro(<< "m_internalMeta2 is not present"); - } + itkExceptionMacro(<< "m_internalMeta2 is not present"); + } if (!m_Interpolator1) { itkExceptionMacro(<< "Interpolator1 is not present"); @@ -328,6 +335,7 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric: weightings[2] = 1. / dtr; break; case SIX_DOF: + //G: ITHINKTHOSEARE FLIPPED!! weightings[0] = 1.; weightings[1] = 1.; weightings[2] = 1.; @@ -348,7 +356,7 @@ typename gTwoImageToOneImageMetric:: ParametersType gTwoImageToOneImageMetric::GetParameters() const { -// ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ + // 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()) { diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 494217f..6e98f98 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -61,6 +61,7 @@ const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds) itkImageProcessor::itkImageProcessor() { + iNWUnits = 1; // std::cout<<"itkImageProcessor::itkImageProcessor contructor"<SetProjectionCenterOffset1(Point3D); m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D); - - // // TEST MAPPING - // { - - // ImageType3D::PointType m_COR; - // ImageType3D::PointType m_TransformOrigin; // Origin of the user defined trasnform. likely isocenter. - // ImageType3D::PointType m_Translations; - // ImageType3D::PointType m_Rotations; - - // 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); - - // }; - // setup the Automatic Registration metric = MetricType::New(); mimetric = MIMetricType::New(); @@ -239,6 +203,13 @@ itkImageProcessor::~itkImageProcessor() } +void itkImageProcessor::SetNumberOfWorkingUnits(int iN){ + if(iN>0){ + iNWUnits = iN; + } +} + + void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os, indent); @@ -482,7 +453,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(); @@ -621,7 +592,7 @@ int itkImageProcessor::load3DSerieFromFolder(const char * pcDirName) imageSeriesReader3D->SetImageIO(dicomIO); imageSeriesReader3D->SetFileNames(fileNames); - imageSeriesReader3D->SetNumberOfWorkUnits(12); + imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits); imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt imageSeriesReader3D->Update(); @@ -1287,33 +1258,6 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg) -// This External User Transform thing need to be checked out -//void itkImageProcessor::CalculateExternalUserTransform(TransformType::Pointer transform, DRTImageMetaInformation::Pointer imageMetaInfo) -//{ -// //crude hack to get from internal transform in LPS to external transform in IEC. -// //result is stored in m_TransformMetaInfo - -// //TODO: support for projectionTransformCenter and userTransformCenter which are not the same -// //currentlz we support only pFakeIsoLPS situations. - -// InternalImageType::DirectionType LPStoIEC_Directions; -// LPStoIEC_Directions = m_CTMetaInfo->GetLPS2IECDirections(); - -// auto parameters = transform->GetParameters(); - -// ImageType3D::PointType translationUser; -// translationUser[0] = parameters[3]; -// translationUser[1] = parameters[4]; -// translationUser[2] = parameters[5]; -// ImageType3D::PointType rotationUser; -// rotationUser[0] = parameters[0]; -// rotationUser[1] = parameters[1]; -// rotationUser[2] = parameters[2]; - -// m_TransformMetaInfo->SetUserTranslations(LPStoIEC_Directions * translationUser); -// m_TransformMetaInfo->SetUserRotations(LPStoIEC_Directions * rotationUser); -//} - void itkImageProcessor::InitializeRegistration( double stepLength, @@ -1338,6 +1282,7 @@ void itkImageProcessor::InitializeRegistration( metric->SetSubtractMean(true); metric->SetMaxTranslation(maxTranslation); + // m_TransformMetaInfo->SetDegreeOfFreedom(dof); m_r23MetaInfo->SetDegreeOfFreedom(dof); @@ -1364,6 +1309,8 @@ void itkImageProcessor::InitializeRegistration( registration->SetFixedImageRegion1(roiAutoReg1); registration->SetFixedImageRegion2(roiAutoReg2); + mimetric->SetFixedImageRegion1(roiAutoReg1); + mimetric->SetFixedImageRegion2(roiAutoReg2); std::cout << "using ROIs for Auto Reg " << std::endl; @@ -1379,6 +1326,8 @@ void itkImageProcessor::InitializeRegistration( registration->SetFixedImageRegion1(region1); registration->SetFixedImageRegion2(region2); + mimetric->SetFixedImageRegion1(roiAutoReg1); + mimetric->SetFixedImageRegion2(roiAutoReg2); } @@ -1620,6 +1569,11 @@ int itkImageProcessor::StartRegistration(std::string extraInfo) Optimizer::ParametersType itkImageProcessor::GetFinalR23Parameters(){ + + if (!m_TransformMetaInfo) { + itkExceptionMacro(<< "m_TransformMetaInfo not present"); + } + Optimizer::ParametersType pPars(6); pPars.fill(0.); @@ -1767,7 +1721,7 @@ void itkImageProcessor::InitializeProjector() 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 @@ -1786,7 +1740,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 @@ -2268,6 +2222,9 @@ void itkImageProcessor::GetProjectionImages(){ std::cout<< "CurrTransform Translations: "<< CurrTransform->GetTranslation()<GetCenter()<Print(std::cout); // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance @@ -2838,26 +2795,87 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) } -// THIS IS READ FOR EXPORT. TO FIX. -double* itkImageProcessor::GetTransformParameters(){ +Optimizer::ParametersType +itkImageProcessor::GetUserIsocentricTransform(){ - ImageType3D::PointType Translations; - ImageType3D::PointType Rotations; + Optimizer::ParametersType + m_Pars(6); + m_Pars.fill(0.); - Translations = m_TransformMetaInfo->GetUserTranslations(); - Rotations = m_TransformMetaInfo->GetUserRotations(); + if(!m_TransformMetaInfo){ + return m_Pars; + } - dTransfParam[0] = Translations[0]; - dTransfParam[1] = Translations[1]; - dTransfParam[2] = Translations[2]; + m_Pars[0] = m_TransformMetaInfo->GetUserRotations()[0]; + m_Pars[1] = m_TransformMetaInfo->GetUserRotations()[1]; + m_Pars[2] = m_TransformMetaInfo->GetUserRotations()[2]; - dTransfParam[3] = Rotations[0]; - dTransfParam[4] = Rotations[1]; - dTransfParam[5] = Rotations[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; - return dTransfParam; } + + +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 regTransform; + } + + std::cout<<" --- GetCompleteIsocentricTransform --- "<GetR()[0] <<" "<< + m_TransformMetaInfo->GetR()[1] <<" "<< + m_TransformMetaInfo->GetR()[2] <<" "<GetT()[0] <<" "<< + m_TransformMetaInfo->GetT()[1] <<" "<< + m_TransformMetaInfo->GetT()[2] <<" "<GetImportOffset()[0]<<" "<< + m_CTMetaInfo->GetImportOffset()[1]<<" "<< + m_CTMetaInfo->GetImportOffset()[2]<<" "<GetUserTranslations(); +// Rotations = m_TransformMetaInfo->GetUserRotations(); + +// dTransfParam[0] = Translations[0]; +// dTransfParam[1] = Translations[1]; +// dTransfParam[2] = Translations[2]; + +// dTransfParam[3] = Rotations[0]; +// dTransfParam[4] = Rotations[1]; +// dTransfParam[5] = Rotations[2]; + +// return dTransfParam; +//} + //// THESE ARE CALLED AT THE END OF REG //void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) //{ @@ -2963,6 +2981,7 @@ void itkImageProcessor::SetRegionFixed1( std::cout << "itkImageProcessor " << std::endl; std::cout << roiAutoReg1 << std::endl; + } void itkImageProcessor::SetRegionFixed2( diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 285e756..825ba79 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -78,6 +78,7 @@ public: CommandIterationUpdate::Pointer optimizerObserver; + void SetNumberOfWorkingUnits(int iN); /** Input data load methods*/ int load3DSerieFromFolder(const char* ); @@ -156,8 +157,21 @@ public: /** Start the registration process*/ int StartRegistration(std::string extraInfo); - /** Get transform parameters for 3D Volume */ - double* GetTransformParameters(); + ///** Get transform parameters for 3D Volume */ + //double* GetTransformParameters(); + /** 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 */ + Optimizer::ParametersType + GetUserIsocentricTransform(); + /** Initialize projection geometry */ void InitializeProjector(); @@ -238,7 +252,7 @@ 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 TransformType = itk::Euler3DTransform; using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; using ResampleFilterType = itk::ResampleImageFilter; @@ -343,13 +357,6 @@ private: ); -// TransformType::Pointer -// MapTransformToNewOrigin( -// ImageType3D::PointType m_COR, -// ImageType3D::PointType m_Translations, -// ImageType3D::PointType m_Rotations -// ); - double CalcProjectionAngleLPS( tPatOrientation pOrient, @@ -383,8 +390,8 @@ private: * m_Projection1VTKTransform, * m_Projection2VTKTransform; - /*Transformation Parameters */ - double dTransfParam[6]; + ///*Transformation Parameters */ + // double dTransfParam[6]; /** * Many meta containers @@ -430,6 +437,9 @@ private: bool m_UseMutualInformation; bool m_UseDumptoFile; + int + iNWUnits; + }; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp index 730ceb6..eb9ac1a 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp @@ -2,74 +2,6 @@ namespace itk { -//template -//gEuler3DTransform ::gEuler3DTransform() -//{ -// //m_CalibPrjCenter = nullptr; // has to be provided by the user. -// //m_RTIso = nullptr; // has to be provided by the user. -// // m_IECtoLPSDirs =; // has to be provided by the user. -//} - - -//template -//gEuler3DTransform::gEuler3DTransform -//(const MatrixType & matrix, const OutputPointType & offset) -//{ - - -//} - -//template -//gEuler3DTransform::gEuler3DTransform -//(unsigned int parametersDimension) -// : Superclass(parametersDimension) - -//{ - -//} - -//template -//void gEuler3DTransform ::SetParameters(const ParametersType & parameters) -//{ - -//// if(m_CalibPrjCenter != nullptr && -//// m_RTIso != nullptr //&& -//// //m_IECtoLPSDirs != nullptr -//// ){ - -//// itk::Point pT; -//// pT[0] = parameters[3]; -//// pT[1] = parameters[4]; -//// pT[2] = parameters[5]; - -//// itk::Point pR; -//// pR[0] = parameters[0]; -//// pR[1] = parameters[1]; -//// pR[2] = parameters[2]; - -//// TransformType::Pointer m_OutputTransform = -//// CalculateInternalTransformV3( -//// pT, -//// pR, -//// m_CalibPrjCenter, -//// m_RTIso, -//// m_IECtoLPSDirs); -//// Superclass::SetParameters(m_OutputTransform->GetParameters()); -//// }else{ -// Superclass::SetParameters(parameters); -//// } - - -//} - -//template -//void gEuler3DTransform ::PrintSelf(std::ostream& os, Indent indent) const -//{ -// Superclass::PrintSelf(os, indent); -// //os << indent << "ComputeGradient: " << static_cast::PrintType>(m_ComputeGradient) -// // << std::endl; -//} - TransformType::Pointer MapTransformToNewOrigin( diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h index aa8a355..04f504e 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h @@ -14,69 +14,6 @@ namespace itk { - - - -//template -//class ITK_TEMPLATE_EXPORT gEuler3DTransform : -// public Euler3DTransform -//{ - -//public: -// ITK_DISALLOW_COPY_AND_ASSIGN(gEuler3DTransform); -// /** Standard class type aliases. */ -// using Self = gEuler3DTransform; -// using Superclass = Euler3DTransform; -// using Pointer = SmartPointer; -// using ConstPointer = SmartPointer; - -// /** New macro for creation of through a Smart Pointer. */ -// itkNewMacro(Self); - -// /** Run-time type information (and related methods). */ -// itkTypeMacro(gEuler3DTransform, Euler3DTransform); - -// // static constexpr unsigned int ParametersDimension = 6; - -// using MatrixType = typename Superclass::MatrixType; -// using OutputPointType = typename Superclass::OutputPointType; -// using ParametersType = typename Superclass::ParametersType; -// using PointType = typename itk::Point; - -// /** Set/Get the transformation from a container of parameters -// * This is typically used by optimizers. There are 6 parameters. The first -// * three represent the angles to rotate around the coordinate axis, and the -// * last three represents the offset. */ -// void -// SetParameters(const ParametersType & parameters) override; - -// itkSetMacro (CalibPrjCenter, PointType); -// itkSetMacro (RTIso, PointType); -// itkSetMacro(IECtoLPSDirs,MatrixType); - - -//protected: -// gEuler3DTransform(const MatrixType & matrix, const OutputPointType & offset); -// gEuler3DTransform(unsigned int parametersDimension); -// gEuler3DTransform(); - -// ~gEuler3DTransform() override = default; - -// PointType -// m_CalibPrjCenter, -// m_RTIso; - -// MatrixType -// m_IECtoLPSDirs; - - -// void -// PrintSelf(std::ostream & os, Indent indent) const override; - -//}; - - - using TransformType = itk::Euler3DTransform; constexpr static unsigned int Dimension = 3; using PixelType3D = short; From 646ba2e0af41e9c4fb363511ad0e35ed082f6dd3 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Mon, 15 May 2023 18:09:04 +0200 Subject: [PATCH 09/12] Calculate Transfrom for SRO isocetric mapped to pCT Zero. --- ...ualInformationTwoImageToOneImageMetric.hxx | 30 +++++------ .../autoreg/itkgTwoImageToOneImageMetric.hxx | 24 ++++----- .../itkDTRrecon/itkImageProcessor.cpp | 52 +++++++++++++++---- .../itkDTRrecon/itkImageProcessor.h | 8 +-- 4 files changed, 70 insertions(+), 44 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx index 561225f..2e18525 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx @@ -102,12 +102,12 @@ MutualInformationTwoImageToOneImageMetric::GetValue() } - std::cout<< "----- itkMutualInformationTwoImage GetValue -----" <m_Transform1.GetPointer()->GetCenter() <m_Transform1.GetPointer()->GetParameters() <m_Transform1.GetPointer()->GetCenter() <m_Transform1.GetPointer()->GetParameters() <m_Interpolator1->SetTransform(this->m_Transform1); this->m_Interpolator1->Initialize(); @@ -139,7 +139,7 @@ MutualInformationTwoImageToOneImageMetric::GetValue() //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; + //std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl; auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion(); unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels(); @@ -169,7 +169,7 @@ MutualInformationTwoImageToOneImageMetric::GetValue() // auto fixedImageRegion2 = fixedImage2->GetBufferedRegion(); // metric2->SetFixedImageRegion(fixedImageRegion2); metric2->SetFixedImageRegion(Superclass::GetFixedImageRegion2()); - std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl; + //std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl; movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion(); numberOfPixels = movingImageRegion.GetNumberOfPixels(); @@ -191,13 +191,13 @@ MutualInformationTwoImageToOneImageMetric::GetValue() 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); +// 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; } diff --git a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx index d8b0cea..770ba54 100644 --- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx +++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx @@ -111,14 +111,14 @@ bool gTwoImageToOneImageMetric::SetTransformParameter 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; +// 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 @@ -162,10 +162,10 @@ bool gTwoImageToOneImageMetric::SetTransformParameter m_Transform1->SetCenter( m_internalMeta1->GetpCalProjCenter()); - std::cout<<"----- itkgTwoImageToOne SetTransformParameters -----"<GetCenter()<GetParameters()<GetCenter()<GetParameters()<SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); - std::cout<< "CurrTransform Rotations: "<< - CurrTransform->GetAngleX() <<" "<< - CurrTransform->GetAngleY() <<" "<< - CurrTransform->GetAngleZ() <<" "<< std::endl; +// std::cout<< "CurrTransform Rotations: "<< +// CurrTransform->GetAngleX() <<" "<< +// CurrTransform->GetAngleY() <<" "<< +// CurrTransform->GetAngleZ() <<" "<< std::endl; - std::cout<< "CurrTransform Translations: "<< - CurrTransform->GetTranslation()<GetCenter()<GetTranslation()<GetCenter()<Print(std::cout); @@ -2823,7 +2823,7 @@ itkImageProcessor::GetUserIsocentricTransform(){ TransformType::Pointer itkImageProcessor::GetCompleteIsocentricTransform -(ImageType3D::PointType TransformCenter){ +(/*ImageType3D::PointType TransformCenter*/){ TransformType::Pointer regTransform = TransformType::New(); @@ -2851,6 +2851,38 @@ itkImageProcessor::GetCompleteIsocentricTransform m_CTMetaInfo->GetImportOffset()[2]<<" "<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); + + std::cout<<"mTransf R"<< mTransf->GetMatrix().GetVnlMatrix() << std::endl; + std::cout<<"mTransf T"<< mTransf->GetTranslation() << std::endl; + + return mTransf; } diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 825ba79..88ba2fc 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -165,7 +165,7 @@ public: * in Isocentric transform with center at the RT Iso. */ TransformType::Pointer - GetCompleteIsocentricTransform(ImageType3D::PointType TransformCenter); + GetCompleteIsocentricTransform(/*ImageType3D::PointType TransformCenter*/); /** Get only the User component of the transform * as parameters' list */ @@ -286,12 +286,6 @@ private: /** ITK to VTK filter */ using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; -// void -// CalculateExternalUserTransform( -// TransformType::Pointer transform, -// DRTImageMetaInformation::Pointer imageMetaInfo); - - TransformType::Pointer transform1, transform2, From e97a24642720ebc9d00791c03458d683de925e03 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Tue, 16 May 2023 13:57:58 +0200 Subject: [PATCH 10/12] Crash on close patient and removal of uncessary load from dir - flipped unload sequence, Scouts, RTP, RTS, CT - removed unused load functions --- .../itkDTRrecon/itkImageProcessor.cpp | 139 ------------------ .../itkDTRrecon/itkImageProcessor.h | 25 ++-- 2 files changed, 12 insertions(+), 152 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index d39df10..f74f9db 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -524,141 +524,6 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ } -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(iNWUnits); - 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(); - - if (m_VolumeSourceDupli) { - m_VolumeSourceDupli = NULL; - m_VolumeSourceDupli = DuplicatorType::New(); - } - - m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput()); - m_VolumeSourceDupli->Update(); - - caster3D = NULL; - - } - - } - 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); - -} - /** Fill Meta after 3D volume load */ int itkImageProcessor::fill3DVolumeMeta( InternalImageType::Pointer m_InputImage, @@ -2951,10 +2816,6 @@ double itkImageProcessor::GetCurrentMetricValue() } } -void itkImageProcessor::SetROI(double, double, double, double) -{ - //TODO: Not implemanted yet. ROI are hard coded and can be enabled using SetFullROI. -} void itkImageProcessor::SetOptimizer(std::string optimizer) { diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 88ba2fc..c451a42 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -78,10 +78,10 @@ public: CommandIterationUpdate::Pointer optimizerObserver; + /** 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(); @@ -144,21 +144,12 @@ public: Optimizer::ParametersType GetFinalR23Parameters(); - /** Get transform parameters for 3D Volume */ -// void GetFinalTranslations(double&, double&, double&); -// void GetFinalRotations(double&, double&, double&); - -// /** Set user transform parameters for 3D Volume in LPS*/ -// void SetUserTranslationsLPS(double, double, double); -// void SetUserRotationsLPS(double, double, double); - /** Initialize the registration pipeline*/ void InitializeRegistration(double, double, eDegreeOfFreedomType); /** Start the registration process*/ int StartRegistration(std::string extraInfo); ///** Get transform parameters for 3D Volume */ - //double* GetTransformParameters(); /** Get the complete transform, likely for export to SRO. * This includes user and hidden contributions. * Transform center has to be specified, Zero would result @@ -168,7 +159,10 @@ public: GetCompleteIsocentricTransform(/*ImageType3D::PointType TransformCenter*/); /** Get only the User component of the transform - * as parameters' list */ + * 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(); @@ -204,18 +198,23 @@ public: void WriteProjectionImages(); void Write2DImages(); + /** Auto Reg23 methods */ + /** Get the current cost function value from the optimizer*/ double GetOptimizerValue(); /** Get the cost function value for the current transform*/ double GetCurrentMetricValue(); - void SetROI(double, double, double, double); - + /** Set Optimizer type */ void SetOptimizer(std::string); + /** Set Metric type */ void SetMetric(std::string); + /** Use full image as ROI */ void SetFullROI(bool); + /** Maximum number of iterations for the optimizer */ void SetMaxNumberOfIterations(int); + /** Define the region to be used as ROI on fixed images */ void SetRegionFixed1(int,int,int,int); void SetRegionFixed2(int,int,int,int); From 21751524bb20b1143493ed35191bed715fdb1991 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Tue, 16 May 2023 17:03:26 +0200 Subject: [PATCH 11/12] Calculation of offset in PAT and rework of signals to approval. - offset calc is verified against outside script. To be careful now in case we want to exclude the rotations from import offset. For PAT calc we need the complete matrix, therefore we may need to filter the rotations somewhere else. - initial work on signal of current transform. was based on readout of spin box. i would use internal storage and check it against spin. - i leave it broken so that we know where to continue --- .../itkDTRrecon/itkImageProcessor.cpp | 104 +++++++----------- .../itkDTRrecon/itkImageProcessor.h | 14 ++- 2 files changed, 53 insertions(+), 65 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index f74f9db..f69b637 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -204,12 +204,36 @@ itkImageProcessor::~itkImageProcessor() } void itkImageProcessor::SetNumberOfWorkingUnits(int iN){ - if(iN>0){ + 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); @@ -262,7 +286,6 @@ double itkImageProcessor::GetPanelOffset2(){ void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof) { - // m_TransformMetaInfo->SetDegreeOfFreedom(dof); m_r23MetaInfo->SetDegreeOfFreedom(dof); } @@ -271,7 +294,6 @@ itkImageProcessor::GetDegreeOfFreedom() { return m_r23MetaInfo->GetDegreeOfFreedom(); - //m_TransformMetaInfo->GetDegreeOfFreedom(); } void itkImageProcessor::SetCustom_ImportTransform(double dTx, @@ -2697,24 +2719,24 @@ itkImageProcessor::GetCompleteIsocentricTransform if(m_TransformMetaInfo == nullptr || m_CTMetaInfo == nullptr){ - return regTransform; + return nullptr; } - std::cout<<" --- GetCompleteIsocentricTransform --- "<GetR()[0] <<" "<< - m_TransformMetaInfo->GetR()[1] <<" "<< - m_TransformMetaInfo->GetR()[2] <<" "<GetT()[0] <<" "<< - m_TransformMetaInfo->GetT()[1] <<" "<< - m_TransformMetaInfo->GetT()[2] <<" "<GetImportOffset()[0]<<" "<< - m_CTMetaInfo->GetImportOffset()[1]<<" "<< - m_CTMetaInfo->GetImportOffset()[2]<<" "<GetR()[0] <<" "<< +// m_TransformMetaInfo->GetR()[1] <<" "<< +// m_TransformMetaInfo->GetR()[2] <<" "<GetT()[0] <<" "<< +// m_TransformMetaInfo->GetT()[1] <<" "<< +// m_TransformMetaInfo->GetT()[2] <<" "<GetImportOffset()[0]<<" "<< +// m_CTMetaInfo->GetImportOffset()[1]<<" "<< +// m_CTMetaInfo->GetImportOffset()[2]<<" "<GetLPS2IECDirections() * m_RTMetaInfo->GetIsocenterLPS(); @@ -2744,56 +2766,10 @@ itkImageProcessor::GetCompleteIsocentricTransform pZero.Fill(0.); mTransf->SetCenter(pZero); - std::cout<<"mTransf R"<< mTransf->GetMatrix().GetVnlMatrix() << std::endl; - std::cout<<"mTransf T"<< mTransf->GetTranslation() << std::endl; - return mTransf; } -// THIS IS READ FOR EXPORT. TO FIX. -//double* itkImageProcessor::GetTransformParameters(){ - -// ImageType3D::PointType Translations; -// ImageType3D::PointType Rotations; - - - -// Translations = m_TransformMetaInfo->GetUserTranslations(); -// Rotations = m_TransformMetaInfo->GetUserRotations(); - -// dTransfParam[0] = Translations[0]; -// dTransfParam[1] = Translations[1]; -// dTransfParam[2] = Translations[2]; - -// dTransfParam[3] = Rotations[0]; -// dTransfParam[4] = Rotations[1]; -// dTransfParam[5] = Rotations[2]; - -// return dTransfParam; -//} - -//// THESE ARE CALLED AT THE END OF REG -//void itkImageProcessor::GetFinalTranslations(double& dX, double& dY, double& dZ) -//{ -// ImageType3D::PointType userTranslations = m_TransformMetaInfo->GetUserTranslations(); - -// dX = userTranslations[0]; -// dY = userTranslations[1]; -// dZ = userTranslations[2]; -//} - -//// THESE ARE CALLED AT THE END OF REG -//void itkImageProcessor::GetFinalRotations(double& dRX, double& dRY, double& dRZ) -//{ -// ImageType3D::PointType userRotations = m_TransformMetaInfo->GetUserRotations(); - -// dRX = userRotations[0]; -// dRY = userRotations[1]; -// dRZ = userRotations[2]; - -//} - double itkImageProcessor::GetOptimizerValue() { return m_OptmizerValue; diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index c451a42..43f6108 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -140,7 +140,8 @@ public: void SetInitialTranslations(double, double, double); void SetInitialRotations(double, double, double); - + /** Returns user components only of the autoReg + * */ Optimizer::ParametersType GetFinalR23Parameters(); @@ -166,6 +167,17 @@ public: 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(); From ab29d81e5e81658a7d3677112ba796f360e8e4e1 Mon Sep 17 00:00:00 2001 From: Proton local user Date: Wed, 17 May 2023 10:49:22 +0200 Subject: [PATCH 12/12] UseOfRotations for hidden transform and station name - itkDRT has method to set whether rotations should be used for import offset and hidden transform. Results checked against external code for final PAT conversion. - helpers for query of station name and Fov moved to helpers functions file - preparation of station name check in scout processor load of 2D and 3D images --- .../itkDTRrecon/DRTMetaInformation.cpp | 2 + .../itkDTRrecon/DRTMetaInformation.h | 7 ++ .../itkDTRrecon/itkImageProcessor.cpp | 76 +++++++++---------- .../itkDTRrecon/itkImageProcessor.h | 24 +++--- .../itkDTRrecon/itkImageProcessorHelpers.cpp | 67 ++++++++++++++++ .../itkDTRrecon/itkImageProcessorHelpers.h | 9 +++ 6 files changed, 132 insertions(+), 53 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 4981437..ff6ae1c 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -474,6 +474,8 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_ProjectionCenter.Fill(0.); + this->m_UseRotationsForHiddenTransform = true; + } void diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index 5beb57a..0e6e9cd 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -408,6 +408,10 @@ public: itkSetMacro(ProjectionCenter, PointType); itkGetMacro(ProjectionCenter, PointType); + itkSetMacro(UseRotationsForHiddenTransform, bool); + itkGetMacro(UseRotationsForHiddenTransform, bool); + + protected: double @@ -447,6 +451,9 @@ protected: m_ProjectionCenterOffset2; + bool + m_UseRotationsForHiddenTransform; + /** Default Constructor **/ DRTProjectionGeometryImageMetaInformation (); /** Default Destructor **/ diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index f69b637..ecaf8ed 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -36,29 +36,6 @@ gfattori 08.11.2021 namespace itk { - -//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; @@ -296,6 +273,17 @@ itkImageProcessor::GetDegreeOfFreedom() 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, @@ -1060,25 +1048,25 @@ int itkImageProcessor::load2D(const char * pcFName){ -int itkImageProcessor::query2DimageReconstructionDiameter(const char *pcFName){ +//int itkImageProcessor::query2DimageReconstructionDiameter(const char *pcFName){ - /* Check if we can open the file */ - gdcm::Reader R; - R.SetFileName(pcFName); - if(!R.Read()) - { cerr<<"ERROR: cannot read file: "<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) ) ); m_TransformMetaInfo->SetHiddenRotations( - m_DRTGeometryMetaInfo->GetIECS2IECScannerR()); + (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ? + m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero) + ); this->UpdateProjectionGeometryMeta(); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 43f6108..df9f45b 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -87,7 +87,7 @@ public: int load2D(const char *); int unload2DAndMeta(int); - int query2DimageReconstructionDiameter(const char*); + //int query2DimageReconstructionDiameter(const char*); void loadRTPlanInfo(double, double, double, double, double ,double); int unloadRTPlanAndMeta(); @@ -121,7 +121,6 @@ public: tDegreeOfFreedomEnum GetDegreeOfFreedom(); void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double); - void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); /** Intensity threshold used for image projection, @@ -131,8 +130,16 @@ 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(); @@ -172,10 +179,8 @@ public: */ const CTVolumeImageMetaInformation::Pointer GetCTMetaInfo(); - const RTGeometryMetaInformation::Pointer GetRTMetaInfo(); - const DRTProjectionGeometryImageMetaInformation::Pointer GetDRTGeoMetaInfo(); @@ -230,9 +235,6 @@ public: void SetRegionFixed1(int,int,int,int); void SetRegionFixed2(int,int,int,int); - /** Optimizer which tries to find the minimn (Powell Optimizer)*/ - using OptimizerType = itk::PowellOptimizer; - protected: /** Various pixel types */ using InternalPixelType = float; @@ -263,12 +265,11 @@ 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 gTransformType = itk::gEuler3DTransform; - + /** Optimizer which tries to find the minimun (Powell Optimizer)*/ + using OptimizerType = itk::PowellOptimizer; /** Optimizer which tries to find the minimn (Powell Optimizer)*/ using AmoebaOptimizerType = itk::AmoebaOptimizer; /** Optimizer which samples the whol space */ @@ -395,9 +396,6 @@ private: * m_Projection1VTKTransform, * m_Projection2VTKTransform; - ///*Transformation Parameters */ - // double dTransfParam[6]; - /** * Many meta containers */ diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp index eb9ac1a..bb93324 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp @@ -1,4 +1,6 @@ #include "itkImageProcessorHelpers.h" +#include +#include namespace itk { @@ -95,4 +97,69 @@ CalculateInternalTransformV3( } +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 index 04f504e..ea16394 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h @@ -10,6 +10,9 @@ #include "itkObject.h" #include "itkObjectFactory.h" +#include +#include + namespace itk { @@ -40,7 +43,13 @@ TransformType::Pointer 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