From c08dba7cf02785d9996a87dffeab05195c1acd1b Mon Sep 17 00:00:00 2001 From: fattori_g Date: Wed, 23 Feb 2022 22:11:21 +0100 Subject: [PATCH] Isocentric transform * Meta Transform container * Mapping from User Isocentric IEC correction to Projection center LPS (6DOF) * Cleanup continues --- .../itkDTRrecon/DRTMetaInformation.cpp | 4 + .../itkDTRrecon/DRTMetaInformation.h | 24 +- .../itkDTRrecon/itkImageProcessor.cpp | 613 +++++++++++++----- .../itkDTRrecon/itkImageProcessor.h | 66 +- 4 files changed, 541 insertions(+), 166 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index c8ab0be..6f8a6d2 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -321,6 +321,10 @@ TransformMetaInformation (){ m_Rotations.Fill(0.); + m_UserRotations.Fill(0.); + + m_UserTranslations.Fill(0.); + m_ReferenceTransform.Fill(0.); m_CurrentTransform.Fill(0.); diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index c98ea86..b270e7b 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -428,12 +428,32 @@ public: void PrintSelf(std::ostream& os, itk::Indent indent) const; -protected: + itkSetMacro(Translations,PointType); + itkGetMacro(Translations,PointType); + itkSetMacro(Rotations,PointType); + itkGetMacro(Rotations,PointType); + + itkSetMacro(UserTranslations,PointType); + itkGetMacro(UserTranslations,PointType); + + itkSetMacro(UserRotations,PointType); + itkGetMacro(UserRotations,PointType); + + itkSetMacro(ReferenceTransform,TransformMatrixType); + itkGetMacro(ReferenceTransform,TransformMatrixType); + + itkSetMacro(CurrentTransform,TransformMatrixType); + itkGetMacro(CurrentTransform,TransformMatrixType); + + +protected: PointType m_Translations, - m_Rotations; + m_Rotations, + m_UserTranslations, + m_UserRotations; TransformMatrixType m_ReferenceTransform, diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 919aa82..d8b18c2 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -86,11 +86,13 @@ itkImageProcessor::itkImageProcessor() // std::cout<<"itkImageProcessor::itkImageProcessor contructor"<SetIdentity(); + interpolator1 = InterpolatorType::New(); interpolator2 = InterpolatorType::New(); - TZero[0]=TZero[1]=TZero[2]=0.; - RZero[0]=RZero[1]=RZero[2]=0.; + // TZero[0]=TZero[1]=TZero[2]=0.; + // RZero[0]=RZero[1]=RZero[2]=0.; toVTK2D1 = ITKtoVTKFilterType::New(); toVTK2D2 = ITKtoVTKFilterType::New(); @@ -120,7 +122,7 @@ itkImageProcessor::itkImageProcessor() - /* Set to NULL the metainfo that are filled in on load */ + /* Set to NULL the metainfo */ m_CTMetaInfo = NULL; m_TImage1MetaInfo = NULL; @@ -136,7 +138,9 @@ itkImageProcessor::itkImageProcessor() m_TransformMetaInfo = NULL; - /* Initialiser the projection geoemtry with defaults */ + m_TransformMetaInfo = TransformMetaInformation::New(); + + /* Initialise the projection geoemtry with defaults */ m_DRTGeometryMetaInfo = DRTProjectionGeometryImageMetaInformation::New(); m_DRTGeometryMetaInfo->SetSCD(570.); @@ -169,6 +173,45 @@ itkImageProcessor::itkImageProcessor() m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Point3D); m_DRTGeometryMetaInfo->SetIECS2IECScannerR(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); + + // }; + + } itkImageProcessor::~itkImageProcessor() @@ -440,7 +483,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) /* initialise DRT meta */ m_DRTImage1MetaInfo = DRTImageMetaInformation::New(); m_DRTImage1MetaInfo->SetProjectionAngleLPS( - CalcProjectionAngleLPS( + this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()) ); @@ -452,7 +495,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) m_DRTImage2MetaInfo = DRTImageMetaInformation::New(); m_DRTImage2MetaInfo->SetProjectionAngleLPS( - CalcProjectionAngleLPS( + this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ); @@ -460,25 +503,13 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) m_DRTGeometryMetaInfo->GetSCD()); - std::cout<< m_DRTGeometryMetaInfo->GetDRT1Size() <GetDRT1Spacing() <UpdateProjectionGeometryMeta(); - std::cout<< m_DRTGeometryMetaInfo->GetDRT1Size() <GetDRT1Spacing() <GetSize() <GetSpacing() <GetOrigin() <GetSize() <GetSpacing() <GetOrigin() <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; + + // m_OutputTransform.Print(std::cout); + return m_OutputTransform; +} + +itkImageProcessor::TransformType::Pointer +itkImageProcessor::CalculateInternalTransform( + ImageType3D::PointType m_TranslationOffset, //IEC + ImageType3D::PointType m_RotationOffset, //IEC + ImageType3D::PointType m_TranslationUser, //IEC + ImageType3D::PointType m_RotationUser, //IEC + ImageType3D::PointType m_ProjectionTransformCenter, //LPS + ImageType3D::PointType m_UserTransformCenter, //LPS + ImageType3D::PointType m_OffsetTransformCenter, //LPS + InternalImageType::DirectionType m_IECtoLPSDirections + ) +{ + + //Convert all inputs into LPS + + ImageType3D::PointType m_TOffsetLPS = + m_IECtoLPSDirections * m_TranslationOffset; + + ImageType3D::PointType m_ROffsetLPS = + m_IECtoLPSDirections * m_RotationOffset; + + ImageType3D::PointType m_TUserLPS = + m_IECtoLPSDirections * m_TranslationUser; + + ImageType3D::PointType m_RUserLPS = + m_IECtoLPSDirections * m_RotationUser; + + std::cout<< "m_ProjectionTransformCenter "<< m_ProjectionTransformCenter <GetCenter() "<< m_UserTransform->GetCenter() <GetTranslation() "<< m_UserTransform->GetTranslation() <GetAngleX() "<< m_UserTransform->GetAngleX() <SetComputeZYX(true); + m_OutputTransform ->SetIdentity(); + + m_OutputTransform->SetTranslation( + //m_OffsetTransform->GetTranslation() + + m_UserTransform->GetTranslation() + ); + + m_OutputTransform->SetRotation( + //m_OffsetTransform->GetAngleX() + + m_UserTransform->GetAngleX(), + //m_OffsetTransform->GetAngleY() + + m_UserTransform->GetAngleY(), + //m_OffsetTransform->GetAngleZ() + + m_UserTransform->GetAngleZ() + ); + + m_OutputTransform->SetCenter( + m_ProjectionTransformCenter); + + return + m_OutputTransform; +} void itkImageProcessor::InitializeProjector() { - transform->SetComputeZYX(true); - - TransformType::OutputVectorType translation; - translation[0] = TZero[0]; - translation[1] = TZero[1]; - translation[2] = TZero[2]; - - transform->SetTranslation(translation); - const double dtr = (atan(1.0) * 4.0) / 180.0; - transform->SetRotation( - dtr * RZero[0], - dtr * RZero[1], - dtr * RZero[2]); - - if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { - // TODO + std::cout<< "itkImageProcessor::InitializeProjector()" <GetLPS2IECDirections().GetTranspose(); + + TransformType::Pointer CurrTransform = + CalculateInternalTransform( + ZeroPoint , + ZeroPoint , + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + ZeroPoint, + IECtoLPS_Directions + ); + + transform->SetComputeZYX(true); + transform->SetIdentity(); + + transform->SetTranslation( + CurrTransform->GetTranslation()); + transform->SetRotation( + CurrTransform->GetAngleX(), + CurrTransform->GetAngleY(), + CurrTransform->GetAngleZ() + ); +// transform->SetCenter( +// CurrTransform->GetCenter() +// ); + std::cout<< "itkImageProcessor::InitializeProjector()" <GetTranslations()[0]; //TZero[0]; +// translation[1] = m_TransformMetaInfo->GetTranslations()[1]; //TZero[1]; +// translation[2] = m_TransformMetaInfo->GetTranslations()[2]; //TZero[2]; + +// transform->SetTranslation(translation); + const double dtr = (atan(1.0) * 4.0) / 180.0; +// transform->SetRotation( +// dtr * m_TransformMetaInfo->GetRotations()[0],// RZero[0], +// dtr * m_TransformMetaInfo->GetRotations()[1], // RZero[1], +// dtr * m_TransformMetaInfo->GetRotations()[2]); //RZero[2]); + +// if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != +// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { +// // TODO +// } + transform->SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); + transform->Print(std::cout); // 2D Image 1 interpolator1->SetProjectionAngle( @@ -1010,6 +1221,89 @@ void itkImageProcessor::InitializeProjector() interpolator2->SetTransform(transform); interpolator2->Initialize(); + + + + + resampleFilter1 = ResampleFilterType::New(); + resampleFilter1->SetInput( + image3DIn); + resampleFilter1->SetDefaultPixelValue(0); + + resampleFilter1->SetNumberOfWorkUnits(8); + + // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance + // have been set before registration. Here we only need to replace the initial + // transform with the final transform. + interpolator1->SetTransform(transform); + interpolator1->Initialize(); + resampleFilter1->SetInterpolator(interpolator1); + + + resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize()); + resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing()); + resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin()); + + + resampleFilter2 = ResampleFilterType::New(); + resampleFilter2->SetInput( + image3DIn); + resampleFilter2->SetDefaultPixelValue(0); + resampleFilter2->SetNumberOfWorkUnits(8); + + + // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance + // have been set before registration. Here we only need to replace the initial + // transform with the final transform. + interpolator2->SetTransform(transform); + interpolator2->Initialize(); + resampleFilter2->SetInterpolator(interpolator2); + + resampleFilter2->SetSize(m_DRTImage2MetaInfo->GetSize()); + resampleFilter2->SetOutputSpacing(m_DRTImage2MetaInfo->GetSpacing()); + resampleFilter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOrigin()); + + filter1 = + ChangeInformationFilterType::New(); + + filter2 = + ChangeInformationFilterType::New(); + + + /* First of all we need the center of the Projection images in its reference frame */ + resampleFilter1->Update(); + filter1->SetInput( resampleFilter1->GetOutput() ); + filter1->SetOutputDirection(m_DRTImage1MetaInfo->GetImageDirectionsLPS() );//IECtoLPS_Directions); + filter1->ChangeDirectionOn(); + filter1->SetOutputOrigin( m_DRTImage1MetaInfo->GetOriginLPS() ); + filter1->ChangeOriginOn(); + filter1->UpdateOutputInformation(); + + std::cout<< "itkImageProcessor::InitializeProjector() " <Update(); + filter2->SetInput( resampleFilter2 ->GetOutput() ); + filter2->SetOutputDirection(m_DRTImage2MetaInfo->GetImageDirectionsLPS() ); + filter2->ChangeDirectionOn(); + filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() ); + filter2->ChangeOriginOn(); + filter2->UpdateOutputInformation(); + + + filter1->Update(); + filter2->Update(); + + std::cout<< "itkImageProcessor::InitializeProjector() AFTER UPDATE" <GetOutput(); + imageDRT2In = filter2->GetOutput(); + + + std::cout<< "itkImageProcessor::InitializeProjector() END" <SetProjectionOriginLPS( CalculateProjectionCenterLPS( m_CTMetaInfo->GetImportOffset(), @@ -1181,6 +1495,9 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ std::cout<< "Projection Origin LPS" << m_DRTImage2MetaInfo->GetProjectionOriginLPS() <GetOriginLPS() <GetImageDirectionsLPS() <SetComputeZYX(true); - TransformType::OutputVectorType translation; - translation[0] = TZero[0]; - translation[1] = TZero[1]; - translation[2] = TZero[2]; - - finalTransform->SetTranslation(translation); - finalTransform->SetRotation(dtr * RZero[0], dtr * RZero[1], dtr * RZero[2]); - - // The transform is determined by the parameters and the rotation center. - - if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { - // TODO + std::cout<< "itkImageProcessor::GetProjectionImages" <SetCenter( + const double dtr = (atan(1.0) * 4.0) / 180.0; + + ImageType3D::PointType ZeroPoint; + ZeroPoint.Fill(0.); + + + transform->SetComputeZYX(true); + transform->SetIdentity(); + + InternalImageType::DirectionType IECtoLPS_Directions; + + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + + TransformType::Pointer CurrTransform = + CalculateInternalTransform( + ZeroPoint, + ZeroPoint, + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + ZeroPoint, + IECtoLPS_Directions + ); + + transform->SetComputeZYX(true); + transform->SetIdentity(); + + transform->SetTranslation( + CurrTransform->GetTranslation()); + transform->SetRotation( + CurrTransform->GetAngleX(), + CurrTransform->GetAngleY(), + CurrTransform->GetAngleZ() + ); +// transform->SetCenter( +// CurrTransform->GetCenter() +// ); + std::cout<< "itkImageProcessor::GetProjectionImages" <SetComputeZYX(true); +// transform->SetIdentity(); +// TransformType::OutputVectorType translation; +// translation[0] = m_TransformMetaInfo->GetTranslations()[0]; //TZero[0]; +// translation[1] = m_TransformMetaInfo->GetTranslations()[1]; //TZero[1]; +// translation[2] = m_TransformMetaInfo->GetTranslations()[2]; //TZero[2]; + +// transform->SetTranslation(translation); +// transform->SetRotation( +// dtr * m_TransformMetaInfo->GetRotations()[0],// RZero[0], +// dtr * m_TransformMetaInfo->GetRotations()[1], // RZero[1], +// dtr * m_TransformMetaInfo->GetRotations()[2]); //RZero[2]); + + +// if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != +// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { +// // TODO +// } + + transform ->SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); - - using ResampleFilterType = itk::ResampleImageFilter; - - // The ResampleImageFilter is the driving force for the projection image generation. - ResampleFilterType::Pointer resampleFilter1 = ResampleFilterType::New(); - resampleFilter1->SetInput( - image3DIn); - resampleFilter1->SetDefaultPixelValue(0); - - resampleFilter1->SetNumberOfWorkUnits(8); - + //transform ->Print(std::cout); // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance // have been set before registration. Here we only need to replace the initial // transform with the final transform. - interpolator1->SetTransform(finalTransform); + interpolator1->SetTransform(transform); interpolator1->Initialize(); - resampleFilter1->SetInterpolator(interpolator1); - - resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize()); - resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing()); - resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin()); - - // Do the same thing for the output image 2. - ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New(); - resampleFilter2->SetInput( - image3DIn); - resampleFilter2->SetDefaultPixelValue(0); - resampleFilter2->SetNumberOfWorkUnits(8); - - - // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance - // have been set before registration. Here we only need to replace the initial - // transform with the final transform. - interpolator2->SetTransform(finalTransform); + // // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance + // // have been set before registration. Here we only need to replace the initial + // // transform with the final transform. + interpolator2->SetTransform(transform); //finalTransform); interpolator2->Initialize(); - resampleFilter2->SetInterpolator(interpolator2); - - resampleFilter2->SetSize(m_DRTImage2MetaInfo->GetSize()); - resampleFilter2->SetOutputSpacing(m_DRTImage2MetaInfo->GetSpacing()); - resampleFilter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOrigin()); - - resampleFilter1->Update(); - resampleFilter2->Update(); - - /* here probably we should move everything back to dicom */ - - /* Projection images are in a fancy reference system... - * With some knowledge (Projection angle, 3DCOV, etc) they can be - * brought back into Patient reference so that can directly - * overlay with the loaded Localizer images. - */ - - ChangeInformationFilterType::Pointer filter1 = - ChangeInformationFilterType::New(); - - ChangeInformationFilterType::Pointer filter2 = - ChangeInformationFilterType::New(); - - - /* First of all we need the center of the Projection images in its reference frame */ - resampleFilter1->Update(); - filter1->SetInput( resampleFilter1->GetOutput() ); - filter1->SetOutputDirection(m_DRTImage1MetaInfo->GetImageDirectionsLPS() );//IECtoLPS_Directions); - filter1->ChangeDirectionOn(); - filter1->SetOutputOrigin( m_DRTImage1MetaInfo->GetOriginLPS() ); - filter1->ChangeOriginOn(); - filter1->UpdateOutputInformation(); - - - - /* Again.. */ - resampleFilter2->Update(); - filter2->SetInput( resampleFilter2 ->GetOutput() ); - filter2->SetOutputDirection(m_DRTImage2MetaInfo->GetImageDirectionsLPS() ); - filter2->ChangeDirectionOn(); - filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() ); - filter2->ChangeOriginOn(); - filter2->UpdateOutputInformation(); filter1->Update(); @@ -1322,6 +1619,7 @@ void itkImageProcessor::GetProjectionImages(){ imageDRT1In = filter1->GetOutput(); imageDRT2In = filter2->GetOutput(); + std::cout<< "itkImageProcessor::GetProjectionImages" <GetLPS2IECDirections().GetTranspose(); + // InternalImageType::DirectionType IECtoLPS_Directions; - ImageType3D::PointType LPSTranslations = - IECtoLPS_Directions * IECTranslations; + // IECtoLPS_Directions = + // m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); - TZero[0]= LPSTranslations[0]; - TZero[1]= LPSTranslations[1]; - TZero[2]= LPSTranslations[2]; + // ImageType3D::PointType LPSTranslations = + // IECtoLPS_Directions * IECTranslations; + // // TZero[0]= LPSTranslations[0]; + // // TZero[1]= LPSTranslations[1]; + // // TZero[2]= LPSTranslations[2]; + + // m_TransformMetaInfo->SetTranslations(LPSTranslations); + + + ImageType3D::PointType Translations; + Translations[0]= dX; + Translations[1]= dY; + Translations[2]= dZ; + + m_TransformMetaInfo->SetTranslations(Translations); } void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) { - RZero[0]=dX; - RZero[1]=dY; - RZero[2]=dZ; + ImageType3D::PointType Rotations; + Rotations[0]= dX; + Rotations[1]= dY; + Rotations[2]= dZ; + + m_TransformMetaInfo->SetRotations(Rotations); + + // RZero[0]=dX; + // RZero[1]=dY; + // RZero[2]=dZ; + } diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index c141abd..ee1dfa0 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -105,6 +105,7 @@ public: void SetInitialTranslations(double, double, double); void SetInitialRotations(double, double, double); + /** Initialize projection geometry */ void InitializeProjector(); @@ -115,13 +116,6 @@ public: const std::vector GetRTImportOffset(); - /** Apply transform to CT volume. - * This moves the volume based on RT Plan info - * if those are available. */ - //void ApplyVolumeImportTransform(); - - - /** Get the Localizer intensity window and * center for rendering */ @@ -145,17 +139,18 @@ public: void WriteProjectionImages(); void Write2DImages(); - /** Various pixel types */ - using InternalPixelType = float; - using PixelType2D = short; - using PixelType3D = short; - using OutputPixelType = unsigned char; - - using ImageType3D = itk::Image; - using InternalImageType = itk::Image; protected: + /** Various pixel types */ + using InternalPixelType = float; + using PixelType2D = short; + using PixelType3D = short; + using OutputPixelType = unsigned char; + + using ImageType3D = itk::Image; + using InternalImageType = itk::Image; + itkImageProcessor(); virtual ~itkImageProcessor(); void PrintSelf(std::ostream& os, itk::Indent indent) const; @@ -175,6 +170,8 @@ private: the registration method itself. */ using TransformType = itk::Euler3DTransform; using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction; + using ResampleFilterType = itk::ResampleImageFilter; + /** Image reader types */ using ImageReaderType2D = itk::ImageFileReader; @@ -209,6 +206,19 @@ private: ImageType3D::PointType m_ProjectionCenter, bool bZero); + TransformType::Pointer + CalculateInternalTransform( + ImageType3D::PointType m_TranslationOffset, + ImageType3D::PointType m_RotationOffset, + ImageType3D::PointType m_TranslationUser, + ImageType3D::PointType m_RotationUser, + ImageType3D::PointType m_ProjectionTransformCenter, + ImageType3D::PointType m_UserTransformCenter, + ImageType3D::PointType m_OffsetTransformCenter, + InternalImageType::DirectionType m_IECtoLPSDirections + ); + + TransformType::Pointer transform; @@ -235,7 +245,9 @@ private: ChangeInformationFilterType::Pointer m_3DInputChangeInformationToZero; - + /** Apply transform to CT volume. + * This moves the volume based on RT Plan info + * if those are available. */ void UpdateProjectionGeometryMeta(); /* The following functions do not rely on class variables, @@ -255,11 +267,21 @@ private: ); + TransformType::Pointer + MapTransformToNewOrigin( + ImageType3D::PointType m_COR, + ImageType3D::PointType m_Translations, + ImageType3D::PointType m_Rotations + ); + double CalcProjectionAngleLPS( tPatOrientation pOrient, double pAngleIEC); + /** Apply transform to CT volume. + * This moves the volume based on RT Plan info + * if those are available. */ ImageType3D::PointType CalcImportVolumeOffset( ImageType3D::PointType rtCouchOffset, @@ -267,7 +289,17 @@ private: ImageType3D::PointType rtIsocenterLPS, ImageType3D::PointType IEC2DCMMapT); - double TZero[3], RZero[3]; + + // The ResampleImageFilter is the driving force for the projection image generation. + ResampleFilterType::Pointer + resampleFilter1, + resampleFilter2; + + ChangeInformationFilterType::Pointer + filter1, + filter2; + + //double TZero[3], RZero[3]; InternalImageType::DirectionType