Isocentric transform

* Meta Transform container
* Mapping from User Isocentric IEC correction to Projection center LPS (6DOF)
* Cleanup continues

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

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,

reg23Topograms/itkDTRrecon/itkImageProcessor.cpp

@@ -86,11 +86,13 @@ itkImageProcessor::itkImageProcessor()
 
     //    std::cout<<"itkImageProcessor::itkImageProcessor contructor"<<std::endl;
     transform = TransformType::New();
+    transform->SetIdentity();
+
     interpolator1 = InterpolatorType::New();
     interpolator2 = InterpolatorType::New();
 
-    TZero[0]=TZero[1]=TZero[2]=0.;
+    //    TZero[0]=TZero[1]=TZero[2]=0.;
-    RZero[0]=RZero[1]=RZero[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() <<std::endl;
+    /* initialise empy plan info */
-    std::cout<< m_DRTGeometryMetaInfo->GetDRT1Spacing() <<std::endl;
+    m_RTMetaInfo = RTGeometryMetaInformation::New();
 
 
     this->UpdateProjectionGeometryMeta();
 
 
-    std::cout<< m_DRTGeometryMetaInfo->GetDRT1Size() <<std::endl;
-    std::cout<< m_DRTGeometryMetaInfo->GetDRT1Spacing() <<std::endl;
-
-    std::cout<< m_DRTImage1MetaInfo->GetSize() <<std::endl;
-    std::cout<< m_DRTImage1MetaInfo->GetSpacing() <<std::endl;
-    std::cout<< m_DRTImage1MetaInfo->GetOrigin() <<std::endl;
-
-    std::cout<< m_DRTImage2MetaInfo->GetSize() <<std::endl;
-    std::cout<< m_DRTImage2MetaInfo->GetSpacing() <<std::endl;
-    std::cout<< m_DRTImage2MetaInfo->GetOrigin() <<std::endl;
-
-
     // To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image
     // to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
     // geometry is fully defined. The origin of the CT image becomes irrelavent.
@@ -959,32 +990,212 @@ 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;
+
+    //    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 <<std::endl;
+    std::cout<< "m_UserTransformCenter "<< m_UserTransformCenter <<std::endl;
+    std::cout<< "m_TUserLPS "<< m_TUserLPS <<std::endl;
+
+    // Map user to the projection center
+    TransformType::Pointer m_UserTransform =
+            MapTransformToNewOrigin (
+                m_ProjectionTransformCenter - m_UserTransformCenter,
+                m_TUserLPS,
+                m_RUserLPS
+                );
+
+    std::cout<< "m_UserTransform->GetCenter() "<< m_UserTransform->GetCenter()  <<std::endl;
+    std::cout<< "m_UserTransform->GetTranslation() "<< m_UserTransform->GetTranslation() <<std::endl;
+    std::cout<< "m_UserTransform->GetAngleX() "<< m_UserTransform->GetAngleX() <<std::endl;
+
+    // Map offset to the projection center
+    TransformType::Pointer m_OffsetTransform =
+            MapTransformToNewOrigin (
+                m_ProjectionTransformCenter - m_OffsetTransformCenter,
+                m_TOffsetLPS,
+                m_ROffsetLPS
+                );
+
+    TransformType::Pointer m_OutputTransform =
+            TransformType::New();
+    m_OutputTransform ->SetComputeZYX(true);
+    m_OutputTransform ->SetIdentity();
+
+    m_OutputTransform->SetTranslation(
+                //m_OffsetTransform->GetTranslation() +
+                m_UserTransform->GetTranslation()
+                );
+
+    m_OutputTransform->SetRotation(
+                //m_OffsetTransform->GetAngleX() +
+                m_UserTransform->GetAngleX(),
+                //m_OffsetTransform->GetAngleY() +
+                m_UserTransform->GetAngleY(),
+                //m_OffsetTransform->GetAngleZ() +
+                m_UserTransform->GetAngleZ()
+                );
+
+    m_OutputTransform->SetCenter(
+                m_ProjectionTransformCenter);
+
+    return
+            m_OutputTransform;
+}
 
 
 void itkImageProcessor::InitializeProjector()
 {
 
-    transform->SetComputeZYX(true);
+    std::cout<< "itkImageProcessor::InitializeProjector()" <<std::endl;
-
+    if(m_DRTImage1MetaInfo == NULL ||
-    TransformType::OutputVectorType translation;
+            m_DRTImage2MetaInfo == NULL ||
-    translation[0] = TZero[0];
+            m_DRTGeometryMetaInfo == NULL ||
-    translation[1] = TZero[1];
+            m_TransformMetaInfo == NULL     ){
-    translation[2] = TZero[2];
+        //TODO
-
-    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
     }
 
+    ImageType3D::PointType ZeroPoint;
+    ZeroPoint.Fill(0.);
+
+    std::cout<< "itkImageProcessor::InitializeProjector()" <<std::endl;
+
+
+    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::InitializeProjector()" <<std::endl;
+
+//    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);
+    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() " <<std::endl;
+
+
+
+    /* 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();
+    filter2->Update();
+
+    std::cout<< "itkImageProcessor::InitializeProjector()  AFTER UPDATE" <<std::endl;
+
+    imageDRT1In = filter1->GetOutput();
+    imageDRT2In = filter2->GetOutput();
+
+
+    std::cout<< "itkImageProcessor::InitializeProjector() END" <<std::endl;
+
 }
 
 
@@ -1044,7 +1338,15 @@ void itkImageProcessor::loadRTPlanInfo(
         //TODO
     }
 
-    m_RTMetaInfo = RTGeometryMetaInformation::New();
+    if(m_CTMetaInfo != NULL){
+        //TODO
+    }
+
+    if(m_DRTGeometryMetaInfo != NULL){
+        //TODO
+    }
+
+    //m_RTMetaInfo = RTGeometryMetaInformation::New();
 
     ImageType3D::PointType
             Point;
@@ -1083,6 +1385,18 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
         //TODO.
     }
 
+    if(m_DRTGeometryMetaInfo == NULL){
+        //TODO.
+    }
+
+    if(m_DRTImage1MetaInfo == NULL){
+        //TODO.
+    }
+
+    if(m_DRTImage2MetaInfo == NULL){
+        //TODO.
+    }
+
     m_DRTImage1MetaInfo->SetProjectionOriginLPS(
                 CalculateProjectionCenterLPS(
                     m_CTMetaInfo->GetImportOffset(),
@@ -1181,6 +1495,9 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
     std::cout<< "Projection Origin LPS" << m_DRTImage2MetaInfo->GetProjectionOriginLPS() <<std::endl;
     std::cout<< "Origin LPS" << m_DRTImage2MetaInfo->GetOriginLPS() <<std::endl;
     std::cout<< "Directions LPS" << m_DRTImage2MetaInfo->GetImageDirectionsLPS() <<std::endl;
+
+
+
 }
 
 
@@ -1212,109 +1529,89 @@ itkImageProcessor::CalcDRTImageOrigin(
 
 void itkImageProcessor::GetProjectionImages(){
 
-    //    std::cout<< "itkImageProcessor::GetProjectionImages" <<std::endl;
+        std::cout<< "itkImageProcessor::GetProjectionImages" <<std::endl;
-
+    if(m_DRTImage1MetaInfo == NULL ||
-    const double dtr = (atan(1.0) * 4.0) / 180.0;
+            m_DRTImage2MetaInfo == NULL ||
-    TransformType::Pointer finalTransform = TransformType::New();
+            m_DRTGeometryMetaInfo == NULL ||
-
+            m_TransformMetaInfo == NULL     ){
-    finalTransform->SetComputeZYX(true);
+        //TODO
-    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
     }
 
-    finalTransform->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" <<std::endl;
+
+//    transform->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());
 
-
+    //transform ->Print(std::cout);
-    using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
-
-    // 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);
-
     // 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);
 
-
+    //    // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
-    resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize());
+    //    // have been set before registration. Here we only need to replace the initial
-    resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing());
+    //    // transform with the final transform.
-    resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin());
+    interpolator2->SetTransform(transform); //finalTransform);
-
-    // 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);
     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" <<std::endl;
 
 }
 
@@ -1768,31 +2066,52 @@ void itkImageProcessor::WriteProjectionImages()
 void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ)
 {
 
-    /* input is in IEC */
+    if(m_TransformMetaInfo == NULL){
-    ImageType3D::PointType IECTranslations;
+        //todo
-    IECTranslations[0] = dX;
+    }
-    IECTranslations[1] = dY;
-    IECTranslations[2] = dZ;
 
-    InternalImageType::DirectionType IECtoLPS_Directions;
+    //    /* input is in IEC */
+    //    ImageType3D::PointType IECTranslations;
+    //    IECTranslations[0] = dX;
+    //    IECTranslations[1] = dY;
+    //    IECTranslations[2] = dZ;
 
-    IECtoLPS_Directions =
+    //    InternalImageType::DirectionType IECtoLPS_Directions;
-            m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
 
-    ImageType3D::PointType LPSTranslations =
+    //    IECtoLPS_Directions =
-            IECtoLPS_Directions * IECTranslations;
+    //            m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
 
-    TZero[0]= LPSTranslations[0];
+    //    ImageType3D::PointType LPSTranslations =
-    TZero[1]= LPSTranslations[1];
+    //            IECtoLPS_Directions * IECTranslations;
-    TZero[2]= LPSTranslations[2];
 
+    //  //  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;
+    ImageType3D::PointType Rotations;
-    RZero[1]=dY;
+    Rotations[0]= dX;
-    RZero[2]=dZ;
+    Rotations[1]= dY;
+    Rotations[2]= dZ;
+
+    m_TransformMetaInfo->SetRotations(Rotations);
+
+    //    RZero[0]=dX;
+    //    RZero[1]=dY;
+    //    RZero[2]=dZ;
+
 }
 
 

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 <double> 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<PixelType3D, Dimension>;
-    using InternalImageType = itk::Image<InternalPixelType, Dimension>;
 
 
 protected:
+    /** Various pixel types */
+        using InternalPixelType = float;
+        using PixelType2D = short;
+        using PixelType3D = short;
+        using OutputPixelType = unsigned char;
+
+        using ImageType3D = itk::Image<PixelType3D, Dimension>;
+        using InternalImageType = itk::Image<InternalPixelType, Dimension>;
+
     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<double>;
     using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
+    using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
+
 
     /** Image reader types */
     using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
@@ -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