Get Projection inverse transform as vtk4x4Matrix and Projection center LPS for further processing outside of the itkImageProcessor. Required for implementing vtk contour projection based on the same transform.

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.

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 SetProjectionAngle1IEC(double);
-    void SetProjectionAngle2(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<PixelType3D, Dimension>;
 
@@ -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;
+
+
 };
 
 

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);