diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 3a4b34f..5d79ebc 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -149,6 +149,11 @@ itkImageProcessor::itkImageProcessor() rtCouchOffset.Fill(0.); ProjectionCenterFixedAxes.Fill(0.); + + m_Projection1VTKTransform = vtkMatrix4x4::New(); + m_Projection1VTKTransform->Identity(); + m_Projection2VTKTransform = vtkMatrix4x4::New(); + m_Projection2VTKTransform->Identity(); } itkImageProcessor::~itkImageProcessor() @@ -171,6 +176,10 @@ void itkImageProcessor::SetSCD(double dDist) d_scd = dDist; } +double itkImageProcessor::GetSCD(){ + return d_scd; +} + void itkImageProcessor::SetCustom_ImportTransform(double dTx, double dTy, double dTz, @@ -530,6 +539,12 @@ void itkImageProcessor::ApplyVolumeImportTransform(){ } +const double* itkImageProcessor::GetProjectionOriginLPS() +{ + return + m_3DProjectionOriginLPS.GetDataPointer(); +} + double itkImageProcessor::CalcProjectionAngleLPS( tPatOrientation pOrient, @@ -597,16 +612,25 @@ itkImageProcessor::CalcProjectionAngleLPS( } -void itkImageProcessor::SetProjectionAngle1(double dGantryAngle) +void itkImageProcessor::SetProjectionAngle1IEC(double dGantryAngle) { projAngle1_IEC = dGantryAngle; } -void itkImageProcessor::SetProjectionAngle2(double dGantryAngle) +double itkImageProcessor::GetProjectionAngle1LPS(){ + + return projAngle1; +} + +void itkImageProcessor::SetProjectionAngle2IEC(double dGantryAngle) { projAngle2_IEC= dGantryAngle; } +double itkImageProcessor::GetProjectionAngle2LPS(){ + return projAngle2; +} + int itkImageProcessor::load2D(const char * pcFName){ /* Check if we can open the file */ @@ -899,7 +923,8 @@ itkImageProcessor::CalcDRTImageDirections( itkImageProcessor::ImageType3D::PointType itkImageProcessor::CalcDRTImageOrigin(itkImageProcessor::InternalImageType::Pointer m_Image, itkImageProcessor::ImageType3D::PointType m_VolCOOV, - double dAngle + double dAngle, + InternalImageType::DirectionType stdDRT2LPS ){ using ImageRegionType3D = itkImageProcessor::ImageType3D::RegionType; @@ -931,7 +956,7 @@ itkImageProcessor::CalcDRTImageOrigin(itkImageProcessor::InternalImageType::Poin DRTCOV[2] = (( m_Image ->GetOrigin())[2]+LastVoxelPosition[2])/2; itkImageProcessor::InternalImageType::DirectionType DRT2LPS - = this->CalcDRTImageDirections(dAngle, Std_DRT2LPS); + = this->CalcDRTImageDirections(dAngle, stdDRT2LPS); ImageType3D::PointType NewOrigin = m_VolCOOV + DRT2LPS * (m_Image ->GetOrigin() - DRTCOV); @@ -1004,9 +1029,6 @@ void itkImageProcessor::GetProjectionImages(){ resampleFilter1->SetOutputOrigin(origin); - - - // Do the same thing for the output image 2. ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New(); resampleFilter2->SetInput( @@ -1073,7 +1095,8 @@ void itkImageProcessor::GetProjectionImages(){ ImageType3D::PointType NewOrigin = this->CalcDRTImageOrigin(resampleFilter1->GetOutput(), m_3DProjectionOriginLPS, - projAngle1); + projAngle1, + Std_DRT2LPS); InternalImageType::DirectionType DRT2LPS1 = this->CalcDRTImageDirections(projAngle1, Std_DRT2LPS); @@ -1093,7 +1116,8 @@ void itkImageProcessor::GetProjectionImages(){ NewOrigin = this->CalcDRTImageOrigin(resampleFilter2->GetOutput(), m_3DProjectionOriginLPS, - projAngle2); + projAngle2, + Std_DRT2LPS); DRT2LPS1 = this->CalcDRTImageDirections(projAngle2, @@ -1382,6 +1406,65 @@ vtkImageData* itkImageProcessor::GetProjection1VTK() toVTK2D1->GetOutput(); } +vtkMatrix4x4 * itkImageProcessor::GetProjection1VTKTransform() +{ + + m_Projection1VTKTransform->Identity(); + + for(int iIdx = 0; iIdx<3 ; iIdx++){ + for(int jIdx = 0; jIdx<3 ; jIdx++){ + m_Projection1VTKTransform->SetElement(iIdx,jIdx, + interpolator1->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]); + } + } + + m_Projection1VTKTransform->SetElement(0,3, + interpolator1->GetInverseTransform()->GetTranslation()[0]); + + m_Projection1VTKTransform->SetElement(1,3, + interpolator1->GetInverseTransform()->GetTranslation()[1]); + + m_Projection1VTKTransform->SetElement(2,3, + interpolator1->GetInverseTransform()->GetTranslation()[2]); + + + //m_Projection1VTKTransform->Print(std::cout); + + return + m_Projection1VTKTransform; + + +} + +vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform() +{ + + m_Projection2VTKTransform->Identity(); + + for(int iIdx = 0; iIdx<3 ; iIdx++){ + for(int jIdx = 0; jIdx<3 ; jIdx++){ + m_Projection2VTKTransform->SetElement(iIdx,jIdx, + interpolator2->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]); + } + } + + m_Projection2VTKTransform->SetElement(0,3, + interpolator2->GetInverseTransform()->GetTranslation()[0]); + + m_Projection2VTKTransform->SetElement(1,3, + interpolator2->GetInverseTransform()->GetTranslation()[1]); + + m_Projection2VTKTransform->SetElement(2,3, + interpolator2->GetInverseTransform()->GetTranslation()[2]); + + + + return + m_Projection2VTKTransform; + +} + + vtkImageData* itkImageProcessor::GetProjection2VTK() { // Rescale the intensity of the projection images to 0-255 for output. diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 3f32352..a654f56 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -39,6 +39,9 @@ gfattori 08.11.2021 #include "itkImageSeriesReader.h" +#include "vtkImageData.h" +#include "vtkTransform.h" + namespace itk { @@ -75,19 +78,23 @@ public: /** Projection angles - Gantry angle IEC */ - void SetProjectionAngle1(double); - void SetProjectionAngle2(double); + void SetProjectionAngle1IEC(double); + void SetProjectionAngle2IEC(double); + /** Get projection angles - Gantry angle LPS + * Only meaningful after a 3D Volume is loaded */ + double GetProjectionAngle1LPS(); + double GetProjectionAngle2LPS(); /** Source to axis distance - SAD * Often called source to object (SOD), (0018,1111) */ void SetSCD(double); + double GetSCD(); /** Intensity threshold used for image projection, * any voxel below threshold is ignored */ void SetIntensityThreshold(double); /** Custom settings of the projection images */ -// void SetCustom_Isocenter(double,double,double); void SetCustom_2Dres(double,double,double,double); void SetCustom_2Dsize(int,int,int,int); void SetCustom_RTIsocenter_LPS(double, double, double); @@ -103,6 +110,9 @@ public: /** Initialize projection geometry */ void InitializeProjector(); + /** Get Projection origin in LPS coordinates*/ + const double* GetProjectionOriginLPS(); + /** Apply transform to CT volume. * This moves the volume based on RT Plan info * if those are available. */ @@ -122,6 +132,9 @@ public: vtkImageData* GetLocalizer1VTK(); vtkImageData* GetLocalizer2VTK(); + vtkMatrix4x4 * GetProjection1VTKTransform(); + vtkMatrix4x4 * GetProjection2VTKTransform(); + /** Debug writers */ void WriteProjectionImages(); void Write2DImages(); @@ -156,10 +169,6 @@ private: void operator=(const Self&); //purposely not implemented - - //void Finalize(); - - /** Image types */ using ImageType2D = itk::Image; @@ -233,6 +242,9 @@ private: InternalImageType::DirectionType LPStoIEC_Directions; + /* The following functions do not rely on class variables, + * Only input variables are used... */ + InternalImageType::DirectionType CalcDRTImageDirections(double angle, InternalImageType::DirectionType DRT2LPS); @@ -240,7 +252,8 @@ private: ImageType3D::PointType CalcDRTImageOrigin(InternalImageType::Pointer m_Image, ImageType3D::PointType m_VolCOOV, - double dAngle + double dAngle, + InternalImageType::DirectionType stdDRT2LPS ); double @@ -308,6 +321,13 @@ private: InternalImageType::DirectionType Std_DRT2LPS; + + + vtkMatrix4x4 + * m_Projection1VTKTransform, + * m_Projection2VTKTransform; + + }; diff --git a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h index 666c0ea..f3da9ba 100644 --- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h +++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h @@ -166,6 +166,9 @@ public: /** Get a pointer to the Transform. */ itkGetConstObjectMacro(Transform, TransformType); + /** Get a pointer to the Inverse Transform. used to calculate the rays */ + itkGetConstObjectMacro(InverseTransform, TransformType); + /** Set and get the focal point to isocenter distance in mm */ itkSetMacro(FocalPointToIsocenterDistance, double); itkGetMacro(FocalPointToIsocenterDistance, double);