Merge branch 'scoutRTUIDevelREG' into 'ScoutRTUIDevel'

MapTransformToNewOrigin and CalcInternalTransformV3 moved to another file,...

See merge request cpt_bioeng/drt!28

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

reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp

@@ -474,6 +474,8 @@ DRTProjectionGeometryImageMetaInformation(){
 
     this->m_ProjectionCenter.Fill(0.);
 
+    this->m_UseRotationsForHiddenTransform = true;
+
 }
 
 void
@@ -517,25 +519,98 @@ 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 (){
+
+    m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
+
+}
+
+
+void
+R23MetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+R23MetaInformation
+::~R23MetaInformation  ()
+{
+
+}
 
 TransformMetaInformation ::
 TransformMetaInformation (){
 
-    m_Translations.Fill(0.);
+    m_HiddenTranslations.Fill(0.);
 
-    m_Rotations.Fill(0.);
+    m_HiddenRotations.Fill(0.);
 
     m_UserRotations.Fill(0.);
 
     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;
+
+}
+
+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;
 
 }
 

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 **/
@@ -514,6 +521,107 @@ private:
 };
 
 
+
+
+
+class InternalTransformMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+     typedef InternalTransformMetaInformation Self;
+     typedef itk::Object Superclass;
+     typedef itk::SmartPointer<Self> Pointer;
+
+     typedef itk::Point<double, 3> PointType;
+     typedef itk::Matrix<double,3,3> 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{
+
+public:
+    /** standard typedefs **/
+     typedef R23MetaInformation Self;
+     typedef itk::Object Superclass;
+     typedef itk::SmartPointer<Self> 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{
 
@@ -535,11 +643,11 @@ public:
     void PrintSelf(std::ostream& os, itk::Indent indent) const;
 
 
-    itkSetMacro(Translations,PointType);
+    itkSetMacro(HiddenTranslations,PointType);
-    itkGetMacro(Translations,PointType);
+    itkGetMacro(HiddenTranslations,PointType);
 
-    itkSetMacro(Rotations,PointType);
+    itkSetMacro(HiddenRotations,PointType);
-    itkGetMacro(Rotations,PointType);
+    itkGetMacro(HiddenRotations,PointType);
 
     itkSetMacro(UserTranslations,PointType);
     itkGetMacro(UserTranslations,PointType);
@@ -547,29 +655,17 @@ public:
     itkSetMacro(UserRotations,PointType);
     itkGetMacro(UserRotations,PointType);
 
-    itkSetMacro(ReferenceTransform,TransformMatrixType);
+    PointType GetT();
-    itkGetMacro(ReferenceTransform,TransformMatrixType);
+    PointType GetR();
-
-    itkSetMacro(CurrentTransform,TransformMatrixType);
-    itkGetMacro(CurrentTransform,TransformMatrixType);
-
-    itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
-    itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
 
 protected:
 
     PointType
-    m_Translations,
+    m_HiddenTranslations,
-    m_Rotations,
+    m_HiddenRotations,
     m_UserTranslations,
     m_UserRotations;
 
-    TransformMatrixType
-    m_ReferenceTransform,
-    m_CurrentTransform;
-
-    tDegreeOfFreedomEnum
-        m_DegreeOfFreedom;
 
     /** Default Constructor **/
     TransformMetaInformation ();

reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx

@@ -93,33 +93,61 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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 -----" <<std::endl;
+//    std::cout<< "Mutual GetValue T1 :: center : " <<
+//                this->m_Transform1.GetPointer()->GetCenter() <<std::endl;
+//    std::cout<< "Mutual GetValue T1 :: Pars : " <<
+//                this->m_Transform1.GetPointer()->GetParameters() <<std::endl;
+//    std::cout<< "--------------------------------" <<std::endl;
+
     this->m_Interpolator1->SetTransform(this->m_Transform1);
     this->m_Interpolator1->Initialize();
     this->m_Filter1->Update();
 
     using InternalMetricType = itk::MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>;
-
+    /*
-    using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction<TMovingImage, double>;
+     * 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<TMovingImage, double>;
+    /*
+     * LinearInterpolateImageFunction performs better for us, performance degradation not noticeable.
+     * */
+    using LinearInterpType = itk::LinearInterpolateImageFunction<TMovingImage, double>;
 
     //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();
+    //auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
-    metric1->SetFixedImageRegion(fixedImageRegion1);
+    //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<unsigned int>(numberOfPixels * 0.50); //100%
-    metric1->UseAllPixelsOn();
+    // since we have a ROI, then we should not set allPixels to TRUE.
-    metric1->SetNumberOfHistogramBins(30);
+    //metric1->UseAllPixelsOn();
+    metric1->SetNumberOfHistogramBins(50);
 
     metric1->SetFixedImage(fixedImage1);
     metric1->SetMovingImage(this->m_Filter1->GetOutput());
@@ -132,26 +160,26 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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();
+//    auto fixedImageRegion2 = fixedImage2->GetBufferedRegion();
-    metric2->SetFixedImageRegion(fixedImageRegion2);
+//    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<unsigned int>(numberOfPixels * 0.50); //100%
-    //%metric2->SetNumberOfSpatialSamples(numberOfSamples);
+    //metric2->SetNumberOfSpatialSamples(numberOfSamples);
-    metric2->UseAllPixelsOn();
+    //metric2->UseAllPixelsOn();
-    metric2->SetNumberOfHistogramBins(30);
+    metric2->SetNumberOfHistogramBins(50);
-
     metric2->SetFixedImage(fixedImage2);
     metric2->SetMovingImage(this->m_Filter2->GetOutput());
-
     metric2->Initialize();
 
     auto measure2 = metric2->GetValue(dummyParameters);
@@ -163,13 +191,13 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
 
     auto measure = (measure1 + measure2) / 2.0;
 
-    auto oldprecision = std::cout.precision();
+//    auto oldprecision = std::cout.precision();
-    std::cout.precision(std::numeric_limits<double>::digits10 + 2);
+//    std::cout.precision(std::numeric_limits<double>::digits10 + 2);
-    std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
+//    std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
-    std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
+//    std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
-    std::cout << "Measure: " << measure << std::endl;
+//    std::cout << "Measure: " << measure << std::endl;
-    std::cout << "=======================================================" << std::endl;
+//    std::cout << "=======================================================" << std::endl;
-    std::cout.precision(oldprecision);
+//    std::cout.precision(oldprecision);
     return measure;
 }
 

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);
@@ -170,8 +174,14 @@ public:
     itkGetConstObjectMacro(Interpolator2, InterpolatorType);
 
     /** Set/Get the Meta informations. */
-    itkSetObjectMacro(TransformMetaInfo, TransformMetaInformation);
+    itkSetObjectMacro(TransformMetaInfo, R23MetaInformation);
-    itkGetConstObjectMacro(TransformMetaInfo, TransformMetaInformation);
+    itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation);
+
+    itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation);
+    itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation);
+
+    itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation);
+    itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation);
 
     /** Set/Get the output filters. */
     itkSetObjectMacro(Filter1, ChangeInformationFilterType);
@@ -236,8 +246,8 @@ protected:
     itkSetMacro(LastOptimizerParameters, ParametersType);
 
     /** Entry-point for internal ITK filter observer. **/
-    void OnInternalFilterProgressReceptor(itk::Object *caller,
+//    void OnInternalFilterProgressReceptor(itk::Object *caller,
-        const itk::EventObject &event);
+//        const itk::EventObject &event);
 
 
 private:
@@ -248,12 +258,18 @@ private:
     FixedImageConstPointer m_FixedImage1;
     FixedImageConstPointer m_FixedImage2;
 
+    TransformPointer m_IsocIECTransform;
+
     TransformPointer m_Transform1;
     TransformPointer m_Transform2;
     InterpolatorPointer m_Interpolator1;
     InterpolatorPointer m_Interpolator2;
 
-    TransformMetaInformation::Pointer
+    InternalTransformMetaInformation::Pointer
+    m_internalMeta1,
+    m_internalMeta2;
+
+    R23MetaInformation::Pointer
         m_TransformMetaInfo;
 
     ChangeInformationFilterPointer

reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx

@@ -38,7 +38,10 @@ TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::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_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.
@@ -80,6 +83,7 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedIm
     m_FixedImageRegionDefined2 = true;
 }
 
+
 /*
  * Initialize by setting the interconnects between components.
  */
@@ -119,7 +123,8 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize
     // Connect the transform to the Decorator.
     auto* transformOutput = static_cast<TransformOutputType*>(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");
@@ -141,13 +146,26 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::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);
     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);
+    m_Metric->SetinternalMeta2(m_internalMeta2);
+
     m_Metric->SetFilter1(m_Filter1);
     m_Metric->SetFilter2(m_Filter2);
 
@@ -174,36 +192,37 @@ void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::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);
 }
 
-template <typename TFixedImage, typename TMovingImage>
+//template <typename TFixedImage, typename TMovingImage>
-void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
+//void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
-    itk::Object *caller, const itk::EventObject &event)
+//    itk::Object *caller, const itk::EventObject &event)
-{
+//{
-  itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(caller);
+//  itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(caller);
 
-  if(!itk::ProgressEvent().CheckEvent(&event) || !filter)
+//  if(!itk::ProgressEvent().CheckEvent(&event) || !filter)
-    return;
+//    return;
 
-  if (filter)
+//  if (filter)
-  {
+//  {
-//    double p = m_CurrentProgressStart;
+////    double p = m_CurrentProgressStart;
-//    // filter->GetProgress() in [0;1]
+////    // filter->GetProgress() in [0;1]
-//    p += m_CurrentProgressSpan * filter->GetProgress();
+////    p += m_CurrentProgressSpan * filter->GetProgress();
-//    // NOTE: filter is buggy - not in [0;1] if multi-threading is activated!
+////    // NOTE: filter is buggy - not in [0;1] if multi-threading is activated!
-//    if (p > (m_CurrentProgressStart + m_CurrentProgressSpan))
+////    if (p > (m_CurrentProgressStart + m_CurrentProgressSpan))
-//      p = m_CurrentProgressStart + m_CurrentProgressSpan;
+////      p = m_CurrentProgressStart + m_CurrentProgressSpan;
-//    TaskProgressInfo(m_CurrentProgressDirection, p);
+////    TaskProgressInfo(m_CurrentProgressDirection, p);
 
-//    if (m_CancelRequest)
+////    if (m_CancelRequest)
-//      filter->SetAbortGenerateData(true); // may be handled by filter
+////      filter->SetAbortGenerateData(true); // may be handled by filter
 
-      std::cout << "Porca Madonna " << std::endl;
+//      std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
-  }
+//  }
-}
+//}
 
 /*
  * Starts the Registration Process
@@ -212,22 +231,10 @@ template <typename TFixedImage, typename TMovingImage>
 void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
 {
 
-//    typedef itk::MemberCommand<TwoProjectionImageRegistrationMethod> 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 {

reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h

@@ -40,6 +40,8 @@
 #include <linux/limits.h>
 #include <unistd.h>
 
+#include <itkImageProcessorHelpers.h>
+
 namespace itk {
 
 /** \class gTwoImageToOneImageMetric
@@ -169,6 +171,10 @@ public:
     /** Get a pointer to the Transform2.  */
     itkGetConstObjectMacro(Transform2, TransformType);
 
+    itkSetObjectMacro(IsocTransf, TransformType);
+    itkGetObjectMacro(IsocTransf, TransformType);
+
+
     /** Connect the Interpolator. */
     itkSetObjectMacro(Interpolator1, InterpolatorType);
 
@@ -185,10 +191,17 @@ 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);
+
+
+    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);
@@ -261,6 +274,8 @@ protected:
 
     mutable TransformPointer m_Transform1;
     mutable TransformPointer m_Transform2;
+    mutable TransformPointer m_IsocTransf;
+
     InterpolatorPointer m_Interpolator1;
     InterpolatorPointer m_Interpolator2;
 
@@ -274,9 +289,13 @@ protected:
     mutable FixedImageMaskPointer m_FixedImageMask2;
     mutable MovingImageMaskPointer m_MovingImageMask;
 
-    TransformMetaInformation::Pointer
+    R23MetaInformation::Pointer
         m_TransformMetaInfo;
 
+    InternalTransformMetaInformation::Pointer
+    m_internalMeta1,
+    m_internalMeta2;
+
     ChangeInformationFilterPointer
         m_Filter1,
         m_Filter2;

reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx

@@ -32,6 +32,8 @@ gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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;
@@ -54,8 +56,18 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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");
+    }
+    if (!m_IsocTransf){
+        itkExceptionMacro(<< "m_IsocTransf  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];
@@ -63,6 +75,8 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameter
     double RotationAlongY = transformParameters[1];
     double RotationAlongZ = transformParameters[2];
 
+    //std::cout << "m_IsocTransf PARS: "<< transformParameters <<std::endl;
+
     switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
     case X_ONLY:
         TranslationAlongX = parameters[0];
@@ -97,31 +111,82 @@ bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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;
-//    std::cout << " Rotation Along X = " << RotationAlongX / 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 Y = " << RotationAlongY /*/ dtr*/ << " deg" << std::endl;
-//    std::cout << " Rotation Along Z = " << RotationAlongZ / dtr << " deg" << std::endl;
+//    std::cout << " Rotation Along Z = " << RotationAlongZ /*/ dtr*/ << " deg" << std::endl;
 //    std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
 
-    transformParameters[0] = RotationAlongX;
-    transformParameters[1] = RotationAlongY;
-    transformParameters[2] = RotationAlongZ;
-    transformParameters[3] = TranslationAlongX;
-    transformParameters[4] = TranslationAlongY;
-    transformParameters[5] = TranslationAlongZ;
 
-    bool transformValid = (std::abs(TranslationAlongX) < m_MaxTranslation) && (std::abs(TranslationAlongY) < m_MaxTranslation) && (std::abs(TranslationAlongZ) < m_MaxTranslation);
+    // individual R and T components of the IsoC transform
+    ImageType3D::PointType
+            pT;
+    pT[0]  = TranslationAlongX;
+    pT[1]  = TranslationAlongY;
+    pT[2]  = TranslationAlongZ;
+
+    ImageType3D::PointType
+            pR;
+    pR[0]  = RotationAlongX;
+    pR[1]  = RotationAlongY;
+    pR[2]  = RotationAlongZ;
+
+    // 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);
 
     if (transformValid) {
-        m_Transform1->SetParameters(transformParameters);
+        /*
-        m_Transform2->SetParameters(transformParameters);
+            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<<"----- itkgTwoImageToOne SetTransformParameters -----"<<std::endl;
+//        std::cout << "m_Transform1 :: Center: "<< m_Transform1->GetCenter()<<std::endl;
+//        std::cout << "m_Transform1 :: Pars: "<< m_Transform1->GetParameters()<<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;
 }
 
@@ -129,16 +194,23 @@ template <typename TFixedImage, typename TMovingImage>
 void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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) {
@@ -237,7 +309,7 @@ unsigned int gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetNumberOfPa
         break;
     }
 
-    return m_Transform1->GetNumberOfParameters();
+    return m_IsocTransf->GetNumberOfParameters();
 }
 
 template <typename TFixedImage, typename TMovingImage>
@@ -263,6 +335,7 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>:
         weightings[2] = 1. / dtr;
         break;
     case SIX_DOF:
+        //G: ITHINKTHOSEARE FLIPPED!!
         weightings[0] = 1.;
         weightings[1] = 1.;
         weightings[2] = 1.;
@@ -279,10 +352,12 @@ itk::Optimizer::ScalesType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>:
 }
 
 template <typename TFixedImage, typename TMovingImage>
-typename gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::ParametersType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetParameters() const
+typename gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::
+ParametersType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::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:

reg23Topograms/itkDTRrecon/itkImageProcessor.h

@@ -51,6 +51,8 @@ gfattori 08.11.2021
 
 #include "itkQtIterationUpdate.h"
 
+#include "itkImageProcessorHelpers.h"
+
 namespace itk
 {
 
@@ -76,15 +78,16 @@ 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<std::string> );
     int unload3DVolumeAndMeta();
 
     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();
@@ -118,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,
@@ -128,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();
 
 
@@ -137,21 +147,42 @@ public:
     void SetInitialTranslations(double, double, double);
     void SetInitialRotations(double, double, double);
 
-    /** Get transform parameters for 3D Volume */
+    /** Returns user components only of the autoReg
-    void GetFinalTranslations(double&, double&, double&);
+    * */
-    void GetFinalRotations(double&, double&, double&);
+    Optimizer::ParametersType
-
+    GetFinalR23Parameters();
-    /** 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 */
+    ///** 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. 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();
+
+    /** 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();
@@ -184,24 +215,26 @@ 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);
 
-    /** Optimizer which tries to find the minimn (Powell Optimizer)*/
-    using OptimizerType = itk::PowellOptimizer;
-
 protected:
     /** Various pixel types */
     using InternalPixelType = float;
@@ -232,10 +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<double>;
     using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
     using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
 
+    /** 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 */
@@ -264,41 +298,10 @@ private:
     /** ITK to VTK filter */
     using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
 
-    void
-    CalculateExternalUserTransform(
-            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
-      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
     transform1,
-    transform2;
+    transform2,
+    IsocTransf;
 
     InterpolatorType::Pointer
     interpolator1,
@@ -360,20 +363,6 @@ private:
             );
 
 
-//    ImageType3D::PointType
-//    CalcDRTImageOffset(
-//            ImageType3D::PointType m_DRTOffset,
-//            double dAngle,
-//            InternalImageType::DirectionType stdDRT2LPS
-//            );
-
-    TransformType::Pointer
-    MapTransformToNewOrigin(
-            ImageType3D::PointType m_COR,
-            ImageType3D::PointType m_Translations,
-            ImageType3D::PointType m_Rotations
-            );
-
     double
     CalcProjectionAngleLPS(
             tPatOrientation pOrient,
@@ -407,9 +396,6 @@ private:
     * m_Projection1VTKTransform,
     * m_Projection2VTKTransform;
 
-    /*Transformation Parameters */
-    double dTransfParam[6];
-
     /**
      * Many meta containers
      */
@@ -433,6 +419,12 @@ private:
     TransformMetaInformation::Pointer
     m_TransformMetaInfo;
 
+    R23MetaInformation::Pointer
+    m_r23MetaInfo;
+
+    InternalTransformMetaInformation::Pointer
+    m_InternalTransf1,
+    m_InternalTransf2;
 
     double m_OptmizerValue;
     int m_MaxNumberOfIterations;
@@ -448,6 +440,9 @@ private:
     bool m_UseMutualInformation;
     bool m_UseDumptoFile;
 
+    int
+    iNWUnits;
+
 };
 
 

reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp

@@ -0,0 +1,165 @@
+#include "itkImageProcessorHelpers.h"
+#include <gdcmReader.h>
+#include <stdio.h>
+
+namespace itk {
+
+
+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
+        ){
+
+    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;
+}
+
+
+
+TransformType::Pointer
+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;
+
+}
+
+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: "<<pcFName<<std::endl;
+        return "\n";
+    }
+
+    const gdcm::File &file = R.GetFile();
+    const gdcm::DataSet &ds = file.GetDataSet();
+
+    std::string s (gGetStringValueFromTag(gdcm::Tag(0x0008,0x1010), ds));
+    std::copy(s.rbegin(), s.rend(), std::back_inserter(buffer));
+
+    std::reverse(buffer.begin(),buffer.end());
+
+    return buffer.c_str();
+
+}
+
+int query2DimageReconstructionDiameter(const char *pcFName){
+
+    /* Check if we can open the file */
+    gdcm::Reader R;
+    R.SetFileName(pcFName);
+    if(!R.Read())
+    {std::cerr<<"ERROR: cannot read file: "<<pcFName<<std::endl;
+        return -1;
+    }
+
+    const gdcm::File &file = R.GetFile();
+    const gdcm::DataSet &ds = file.GetDataSet();
+
+    char* sTmpString = new char [255];
+    strcpy( sTmpString,gGetStringValueFromTag(gdcm::Tag(0x0018,0x1100), ds));
+
+    return std::atoi(sTmpString);
+}
+
+
+//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();
+}
+
+}

reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h

@@ -0,0 +1,55 @@
+
+#ifndef ITKIMAGEPROCESSORHELPERS_H
+#define ITKIMAGEPROCESSORHELPERS_H
+
+#include "itkEuler3DTransform.h"
+#include "itkImage.h"
+#include "itkEuler3DTransform.h"
+
+#include "itkCommand.h"
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+
+#include <gdcmTag.h>
+#include <gdcmDataSet.h>
+
+namespace itk
+{
+
+
+using TransformType = itk::Euler3DTransform<double>;
+constexpr static unsigned int Dimension = 3;
+using PixelType3D = short;
+using InternalPixelType = float;
+
+using ImageType3D = itk::Image<PixelType3D, Dimension>;
+using InternalImageType = itk::Image<InternalPixelType, Dimension>;
+
+
+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
+             );
+
+const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds);
+
+int query2DimageReconstructionDiameter(const char*);
+
+const char* queryStationName(const char*);
+
+}
+
+
+#endif

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<double>::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;
     }