clear folder structure

itkReg23DRT/CMakeLists.txt

@@ -0,0 +1,59 @@
+cmake_minimum_required(VERSION 2.8.12)
+
+SET(LIB_NAME itkDTRrecon)
+
+add_subdirectory(autoreg)
+
+find_package(ITK REQUIRED)
+include(${ITK_USE_FILE})
+INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
+
+find_package(VTK REQUIRED)
+INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS})
+
+FIND_PACKAGE(GDCM REQUIRED)
+INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS})
+
+SET(SRCS
+    itkImageProcessor.cpp  
+    itkImageProcessorHelpers.cpp
+    vtkContourTopogramProjectionFilter.cxx
+    DRTMetaInformation.cpp
+    itkReg23.cpp
+    itkReg23MetaInformation.cpp
+)
+
+SET(HDR
+    itkImageProcessor.h
+    itkImageProcessorHelpers.h
+    itkQtIterationUpdate.h
+    itkgSiddonJacobsRayCastInterpolateImageFunction.h
+    itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
+    vtkContourTopogramProjectionFilter.h
+    DRTMetaInformation.h
+    itkReg23.h
+    itkReg23MetaInformation.h
+)
+
+ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR})
+
+SET(LINK_LIBS
+    ${LIB_NAME}
+    ${VTK_LIBRARIES}
+    ${ITK_LIBRARIES}
+    gdcmMSFF
+    autoreg
+)
+
+TARGET_LINK_LIBRARIES(
+    ${LINK_LIBS}
+)
+
+
+target_include_directories (${LIB_NAME}
+    PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+
+set(EXTRA_LIBS ${EXTRA_LIBS} itkDTRrecon)
+set(EXTRA_LIBS ${EXTRA_LIBS} PARENT_SCOPE)

itkReg23DRT/DRTMetaInformation.cpp

@@ -0,0 +1,642 @@
+#include "DRTMetaInformation.h"
+
+
+static double HFS2IEC[9] = {
+    1, 0, 0,
+    0, 0, 1,
+    0, -1, 0};
+
+static double FFS2IEC[9] = {
+    -1, 0, 0,
+    0, 0, -1,
+    0, -1, 0};
+
+namespace itk
+{
+
+TopogramImageMetaInformation::
+TopogramImageMetaInformation(){
+
+    this->m_PatientOrientation =
+            tPatOrientation::NotDefined;
+
+    this->m_ProjectionOrientation =
+            tProjOrientationType::NA;
+
+    this->m_WLLevel = 0.;
+
+    this->m_WLWindow = 0.;
+
+    this->m_LPS2IECDirections.SetIdentity();
+
+}
+
+
+void
+TopogramImageMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+TopogramImageMetaInformation
+::~TopogramImageMetaInformation ()
+{
+
+}
+
+
+void TopogramImageMetaInformation ::
+SetPatientOrientation(tPatOrientation m_orient)
+{
+    this->m_PatientOrientation = m_orient;
+
+    for(int iIdx = 0 ; iIdx < 3; iIdx++){
+        for(int jIdx = 0 ; jIdx < 3; jIdx++){
+
+            switch (this->m_PatientOrientation) {
+            case tPatOrientation::HFS:
+                m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3];
+                break;
+            case tPatOrientation ::FFS:
+                m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3];
+                break;
+            default:
+                m_LPS2IECDirections.SetIdentity();
+                break;
+            }
+        }
+    }
+
+}
+
+
+DRTImageMetaInformation::
+DRTImageMetaInformation(){
+
+
+    this->m_Size.Fill(0.);
+
+    this->m_Spacing.Fill(0.);
+
+    //this->m_Origin.Fill(0.);
+
+    this->m_ProjectionOriginLPS.Fill(0.);
+
+    this->m_ProjectionOriginLPSZero.Fill(0.);
+
+    this->m_ProjectionAngleLPS = 0.;
+
+    this->m_SCD = 0.;
+
+    this->m_ImageDirectionsLPS.SetIdentity();
+
+    this->m_OriginLPS.Fill(0.);
+
+    this->m_PanelOffsetPGeo = 0.;
+
+}
+
+void
+DRTImageMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+DRTImageMetaInformation
+::~DRTImageMetaInformation  ()
+{
+
+}
+
+
+DRTImageMetaInformation::PointType
+DRTImageMetaInformation::GetOrigin()
+{
+
+    //PointType Origin;
+    PointType Origin;
+    double o2Dx, o2Dy;
+
+    o2Dx = ((double)this->m_Size[0] - 1.) / 2.;
+    o2Dy = ((double)this->m_Size[1] - 1.) / 2.;
+    // Compute the origin (in mm) of the 2D Image
+    Origin[0] = -(this->m_Spacing[0]) * o2Dx;
+    Origin[1] = -(this->m_Spacing[1]) * o2Dy;
+    Origin[2] = -(this->m_SCD) ;
+
+    return
+            Origin;
+}
+
+
+void DRTImageMetaInformation::SetSizeWithBounds(
+        double* dBounds,
+        SizeType MaxSize,
+        SpacingType Spacing)
+{
+
+//    double dXmin = dBounds[0];
+//    double dYmin = dBounds[1];
+//    double dZmin = dBounds[2];
+//    double dXmax = dBounds[3];
+//    double dYmax = dBounds[4];
+//    double dZmax = dBounds[5];
+
+//    double
+//            dNPixMin1 = 0.,
+//            dNPixMin2 = 0.,
+//            dNPixMax1 = 0.,
+//            dNPixMax2 = 0.;
+
+//    dNPixMin1 = abs(dXmin / Spacing[0]);
+//    dNPixMin2 = abs(dYmin / Spacing[0]);
+//    dNPixMax1 = abs(dXmax / Spacing[0]);
+//    dNPixMax2 = abs(dYmax / Spacing[0]);
+
+//    double dPixMaxOverlap =dNPixMin1;
+
+//    if(dPixMaxOverlap > dNPixMin2){
+//        dPixMaxOverlap = dNPixMin2;
+//    }
+
+//    if(dPixMaxOverlap > dNPixMax1){
+//        dPixMaxOverlap = dNPixMax1;
+//    }
+
+//    if(dPixMaxOverlap > dNPixMax2){
+//        dPixMaxOverlap = dNPixMax2;
+//    }
+
+//    if(floor(dPixMaxOverlap*2) > MaxSize[0]) {
+//        m_Size[0] = MaxSize[0];
+//    } else {
+//        m_Size[0] = floor(dPixMaxOverlap*2);
+//    }
+
+
+
+//    dNPixMin1 = abs(dZmin / Spacing[1]);
+//    dNPixMax1 = abs(dZmax / Spacing[1]);
+
+//    dPixMaxOverlap =dNPixMin1;
+//    if(dPixMaxOverlap > dNPixMax1){
+//        dPixMaxOverlap = dNPixMax1;
+//    }
+
+//    if(floor(dPixMaxOverlap*2) > MaxSize[1]) {
+//        m_Size[1] = MaxSize[1];
+//    } else {
+//        m_Size[1] = floor(dPixMaxOverlap*2);
+//    }
+
+    m_Size[0] = MaxSize[0];
+    m_Size[1] = MaxSize[1];
+    m_Size[2] = 1;
+
+    m_Spacing = Spacing;
+
+    return;
+}
+
+
+DRTImageMetaInformation::PointType
+DRTImageMetaInformation::GetLPStoProjectionGeoLPSOffset()
+{
+
+    PointType pOffset;
+
+    pOffset[0] = - m_ProjectionOriginLPS [0] + m_ProjectionOriginLPSZero [0];
+    pOffset[1] = - m_ProjectionOriginLPS [1] + m_ProjectionOriginLPSZero [1];
+    pOffset[2] = - m_ProjectionOriginLPS [2] + m_ProjectionOriginLPSZero [2];
+
+    return
+            pOffset;
+}
+
+DRTImageMetaInformation ::PointType
+DRTImageMetaInformation ::GetCOV()
+{
+
+    /* calculate image size in 3D Space by finding the last voxel position */
+    PointType Dest;
+    Dest[0]=(m_Size[0]-1) * m_Spacing [0];
+    Dest[1]=(m_Size[1]-1) * m_Spacing [1];
+    Dest[2]=(m_Size[2]-1) * m_Spacing [2];
+
+    PointType LastVoxelPosition = Dest;
+
+    PointType Origin = this->GetOrigin();
+    LastVoxelPosition [0] += Origin[0];
+    LastVoxelPosition [1] += Origin[1];
+    LastVoxelPosition [2] += Origin[2];
+
+    PointType image3DCOV;
+    image3DCOV[0] = (Origin [0] + LastVoxelPosition[0])/2.;
+    image3DCOV[1] = (Origin [1] + LastVoxelPosition[1])/2.;
+    image3DCOV[2] = (Origin [2] + LastVoxelPosition[2])/2.;
+
+    return
+            image3DCOV;
+}
+
+CTVolumeImageMetaInformation::PointType
+CTVolumeImageMetaInformation::ConvertIECPointToLPSPoint(
+        PointType IECPoint){
+
+    DirectionType IECtoLPS_Directions;
+    IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose();
+
+    PointType m_LPSPoint;
+    m_LPSPoint =
+            IECtoLPS_Directions * IECPoint;
+
+    return
+            m_LPSPoint;
+}
+
+CTVolumeImageMetaInformation::PointType
+CTVolumeImageMetaInformation::GetProjectionOriginLPS(
+        PointType ProjectionCenter)
+{
+
+    //    PointType Origin;
+    //    Origin = this->m_OriginLPS;
+    //    Origin = Origin - this->m_ImportOffset;
+
+    DirectionType IECtoLPS_Directions;
+    IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose();
+
+
+    PointType m_ProjectionOriginLPS;
+    m_ProjectionOriginLPS =
+            IECtoLPS_Directions * ProjectionCenter;
+
+    /* in longitudinal direction (Sup-Inf), we put it in the
+      * middle of the volume */
+    m_ProjectionOriginLPS[2] = this->GetCOV()[2] - this->m_ImportOffset[2];
+
+    return
+            m_ProjectionOriginLPS;
+
+}
+
+CTVolumeImageMetaInformation::PointType
+CTVolumeImageMetaInformation::GetProjectionOriginLPSZero(
+        PointType ProjectionCenter)
+{
+
+    PointType Origin;
+    Origin = this->m_OriginLPS;
+    Origin = Origin - this->m_ImportOffset;
+
+
+    DirectionType IECtoLPS_Directions;
+    IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose();
+
+    PointType m_ProjectionOriginLPS;
+    m_ProjectionOriginLPS =
+            IECtoLPS_Directions * ProjectionCenter;
+    /* in longitudinal direction (Sup-Inf), we put it in the
+      * middle of the volume */
+    m_ProjectionOriginLPS[2] = this->GetCOV()[2] - this->m_ImportOffset[2];
+
+
+    return
+            m_ProjectionOriginLPS -  Origin;
+
+}
+
+const std::vector<double>
+CTVolumeImageMetaInformation::CalculateRegions(PointType pProjectionOriginLPS)
+{
+
+    /* calculate image size in 3D Space by finding the last voxel position */
+    PointType Dest;
+    Dest[0]=(m_Size[0]-1) * m_Spacing [0];
+    Dest[1]=(m_Size[1]-1) * m_Spacing [1];
+    Dest[2]=(m_Size[2]-1) * m_Spacing [2];
+
+    PointType LastVoxelPosition =
+            (m_ImageDirections * Dest);
+
+    LastVoxelPosition [0] += m_OriginLPS[0];
+    LastVoxelPosition [1] += m_OriginLPS[1];
+    LastVoxelPosition [2] += m_OriginLPS[2];
+
+    std::vector <double>  vBounds;
+    vBounds.clear();
+    vBounds.push_back( pProjectionOriginLPS[0] - LastVoxelPosition [0] );
+    vBounds.push_back( pProjectionOriginLPS[1] - LastVoxelPosition [1] );
+    vBounds.push_back( pProjectionOriginLPS[2] - LastVoxelPosition [2] );
+    vBounds.push_back( pProjectionOriginLPS[0] - m_OriginLPS [0] );
+    vBounds.push_back( pProjectionOriginLPS[1] - m_OriginLPS [1] );
+    vBounds.push_back( pProjectionOriginLPS[2] - m_OriginLPS [2] );
+
+    return
+            vBounds;
+
+}
+
+CTVolumeImageMetaInformation::
+CTVolumeImageMetaInformation(){
+
+    this->m_PatientOrientation =
+            tPatOrientation::NotDefined;
+
+    this->m_Spacing.Fill(0.);
+
+    this->m_Size.Fill(0.);
+
+    this->m_OriginLPS.Fill(0.);
+
+    this->m_ImageDirections.SetIdentity();
+
+    this->m_LPS2IECDirections.SetIdentity();
+
+    this->m_ImportOffset.Fill(0.);
+}
+
+
+void
+CTVolumeImageMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+CTVolumeImageMetaInformation
+::~CTVolumeImageMetaInformation   ()
+{
+
+}
+
+void CTVolumeImageMetaInformation::
+SetPatientOrientation(tPatOrientation m_orient)
+{
+    this->m_PatientOrientation = m_orient;
+
+    for(int iIdx = 0 ; iIdx < 3; iIdx++){
+        for(int jIdx = 0 ; jIdx < 3; jIdx++){
+
+            switch (this->m_PatientOrientation) {
+            case tPatOrientation::HFS:
+                m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3];
+                break;
+            case tPatOrientation ::FFS:
+                m_LPS2IECDirections.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3];
+                break;
+            default:
+                m_LPS2IECDirections.SetIdentity();
+                break;
+            }
+        }
+    }
+
+
+}
+
+CTVolumeImageMetaInformation::PointType
+CTVolumeImageMetaInformation::GetCOV()
+{
+
+    /* calculate image size in 3D Space by finding the last voxel position */
+    PointType Dest;
+    Dest[0]=(m_Size[0]-1) * m_Spacing [0];
+    Dest[1]=(m_Size[1]-1) * m_Spacing [1];
+    Dest[2]=(m_Size[2]-1) * m_Spacing [2];
+
+    PointType LastVoxelPosition =
+            (m_ImageDirections * Dest);
+
+    LastVoxelPosition [0] += m_OriginLPS[0];
+    LastVoxelPosition [1] += m_OriginLPS[1];
+    LastVoxelPosition [2] += m_OriginLPS[2];
+
+    PointType image3DCOVLPS;
+    image3DCOVLPS[0] = (m_OriginLPS[0] + LastVoxelPosition[0])/2.;
+    image3DCOVLPS[1] = (m_OriginLPS[1] + LastVoxelPosition[1])/2.;
+    image3DCOVLPS[2] = (m_OriginLPS[2] + LastVoxelPosition[2])/2.;
+
+    return
+            image3DCOVLPS;
+}
+
+CTVolumeImageMetaInformation::PointType
+CTVolumeImageMetaInformation::GetOriginWOffset(){
+    return
+            this->m_OriginLPS - this->m_ImportOffset;
+}
+
+
+DRTProjectionGeometryImageMetaInformation::
+DRTProjectionGeometryImageMetaInformation(){
+
+    this->m_ProjectionAngle1IEC = 0.;
+
+    this->m_ProjectionAngle2IEC = 0.;
+
+    this->m_ProjectionAngle1OffsetIEC = 0.;
+
+    this->m_ProjectionAngle2OffsetIEC = 0.;
+
+    this->m_SCD1 = 0.;
+
+    this->m_SCD2 = 0.;
+
+    this->m_SCD1Offset = 0.;
+
+    this->m_SCD2Offset = 0.;
+
+    this->m_Panel1Offset = 0.;
+
+    this->m_Panel2Offset = 0.;
+
+    this->m_IntensityThreshold=0.;
+
+    this->m_DRT1Size.Fill(0.);
+
+    this->m_DRT2Size.Fill(0.);
+
+    this->m_DRT1Spacing.Fill(0.);
+
+    this->m_DRT2Spacing.Fill(0.);
+
+    this->m_IECS2IECScannerT.Fill(0.);
+
+    this->m_IECS2IECScannerR.Fill(0.);
+
+    this->m_ProjectionCenterOffset1.Fill(0.);
+
+    this->m_ProjectionCenterOffset2.Fill(0.);
+
+    this->m_ProjectionCenter.Fill(0.);
+
+    this->m_UseRotationsForHiddenTransform = true;
+
+}
+
+void
+DRTProjectionGeometryImageMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+DRTProjectionGeometryImageMetaInformation
+::~DRTProjectionGeometryImageMetaInformation    ()
+{
+
+}
+
+
+
+
+RTGeometryMetaInformation::
+RTGeometryMetaInformation(){
+
+    this->m_IsocenterLPS.Fill(0.);
+
+    this->m_IsocenterIECS.Fill(0.);
+
+
+}
+
+void
+RTGeometryMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+RTGeometryMetaInformation
+::~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;
+    m_OptimizerType = tOptimizerTypeEnum::POWELL;
+    m_MetricType = tMetricTypeEnum::NCC;
+}
+
+
+void
+R23MetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+R23MetaInformation
+::~R23MetaInformation  ()
+{
+
+}
+
+TransformMetaInformation ::
+TransformMetaInformation (){
+
+    m_HiddenTranslations.Fill(0.);
+
+    m_HiddenRotations.Fill(0.);
+
+    m_UserRotations.Fill(0.);
+
+    m_UserTranslations.Fill(0.);
+
+    //m_ReferenceTransform.Fill(0.);
+
+    //m_CurrentTransform.Fill(0.);
+
+ //   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;
+
+}
+
+void
+TransformMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+TransformMetaInformation
+::~TransformMetaInformation ()
+{
+
+}
+
+
+
+
+
+
+
+
+
+}

itkReg23DRT/DRTMetaInformation.h

@@ -0,0 +1,736 @@
+#ifndef DRTMETAINFORMATION_H
+#define DRTMETAINFORMATION_H
+
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "itkSmartPointer.h"
+#include "itkMacro.h"
+
+#include "itkImage.h"
+
+namespace itk
+{
+
+typedef enum eProjectionOrientationType{
+    NA = 0,
+    LAT=2,
+    PA=1
+} tProjOrientationType;
+
+/** Enum type for Orientation */
+typedef enum eImageOrientationType{
+    NotDefined = 0,
+    HFS = 1,
+    FFS = 2
+} tPatOrientation;
+
+
+typedef enum eHandlingRotImpTransform {
+    ALWAYS_USE = 0,
+    NEVER_USE =1,
+    STATION_DEPENDENT=2
+} tHandlingRotImpTransform;
+
+typedef enum eDegreeOfFreedomType {
+    UNKNOWN = 0,
+    THREE_DOF =1,
+    SIX_DOF =2 ,
+    X_ONLY =3 ,
+    Y_ONLY =4,
+    Z_ONLY =5 ,
+    ROTATION_ONLY =6,
+    OTHER =7
+} tDegreeOfFreedomEnum;
+
+
+typedef enum eOptimizerType{
+    POWELL = 0,
+    AMOEBA = 1,
+    EXHAUSTIVE =2
+} tOptimizerTypeEnum;
+
+typedef enum eMetricType{
+    NCC = 0,
+    MI = 1
+} tMetricTypeEnum;
+
+
+inline int GetNumberOfDOF(eDegreeOfFreedomType dof)
+{
+    switch (dof) {
+    case SIX_DOF:
+        return 6;
+        break;
+    case THREE_DOF:
+    case ROTATION_ONLY:
+        return 3;
+        break;
+    case X_ONLY:
+    case Y_ONLY:
+    case Z_ONLY:
+        return 1;
+        break;
+    default:
+        return 0;
+    }
+    return 0;
+};
+
+
+class DegreeOfFreedom : public itk::Object {
+
+public:
+    typedef DegreeOfFreedom Self;
+    typedef itk::Object Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+
+    itkNewMacro(Self);
+
+protected:
+    bool
+        m_TranslationX,
+        m_TranslationY,
+        m_TranslationZ,
+        m_RotationX,
+        m_RotationY,
+        m_RotationZ;
+    tDegreeOfFreedomEnum m_dof;
+};
+
+class TopogramImageMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef TopogramImageMetaInformation  Self;
+    typedef itk::Object Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::Matrix<double, 3, 3> DirectionType;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(TopogramImageMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    void SetPatientOrientation(tPatOrientation);
+    itkGetEnumMacro(PatientOrientation, tPatOrientation);
+
+    itkSetEnumMacro(ProjectionOrientation, tProjOrientationType);
+    itkGetEnumMacro(ProjectionOrientation, tProjOrientationType);
+
+    itkSetMacro(WLLevel,double);
+    itkGetMacro(WLLevel,double);
+
+    itkSetMacro(WLWindow,double);
+    itkGetMacro(WLWindow,double);
+
+    itkSetMacro(LPS2IECDirections,DirectionType);
+    itkGetMacro(LPS2IECDirections,DirectionType);
+
+protected:
+
+
+    tPatOrientation
+    m_PatientOrientation;
+
+    tProjOrientationType
+    m_ProjectionOrientation;
+
+    /** preferred level (windowing) **/
+    double m_WLLevel;
+    /** preferred window (windowing) **/
+    double m_WLWindow;
+
+    DirectionType
+    m_LPS2IECDirections;
+
+
+    /** Default Constructor **/
+    TopogramImageMetaInformation  ();
+    /** Default Destructor **/
+    virtual ~TopogramImageMetaInformation  ();
+
+private:
+    /** purposely not implemented **/
+    TopogramImageMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+class DRTImageMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef DRTImageMetaInformation Self;
+    typedef itk::Object Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::Point<double, 3> PointType;
+    typedef itk::Vector<double, 3> SpacingType;
+    //typedef itk::FixedArray<unsigned int, 3> SizeType;
+    typedef itk::Size<3> SizeType;
+    typedef itk::Matrix<double, 3, 3> DirectionType;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(DRTImageMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+    itkSetMacro(Size,SizeType);
+    itkGetMacro(Size,SizeType);
+
+    itkSetMacro(Spacing,SpacingType);
+    itkGetMacro(Spacing,SpacingType);
+
+    //itkSetMacro(Origin,PointType);
+    //itkGetMacro(Origin,PointType);
+
+    itkSetMacro(OriginLPS,PointType);
+    itkGetMacro(OriginLPS,PointType);
+
+    itkSetMacro(ProjectionOriginLPS,PointType);
+    itkGetMacro(ProjectionOriginLPS,PointType);
+
+    itkSetMacro(ProjectionOriginLPSZero,PointType);
+    itkGetMacro(ProjectionOriginLPSZero,PointType);
+
+    itkSetMacro(PanelOffsetPGeo,double);
+    itkGetMacro(PanelOffsetPGeo,double);
+
+    itkSetMacro(ProjectionAngleLPS,double);
+    itkGetMacro(ProjectionAngleLPS,double);
+
+    itkSetMacro(SCD,double);
+    itkGetMacro(SCD,double);
+
+    itkSetMacro(ImageDirectionsLPS,DirectionType);
+    itkGetMacro(ImageDirectionsLPS,DirectionType);
+
+    PointType GetCOV();
+    PointType GetLPStoProjectionGeoLPSOffset();
+
+    void SetSizeWithBounds(double*, SizeType, SpacingType);
+
+    PointType GetOrigin();
+
+protected:
+
+    SizeType
+    m_Size;
+
+    SpacingType
+    m_Spacing;
+
+    PointType
+    //m_Origin,
+    m_OriginLPS,
+    m_ProjectionOriginLPS,
+    m_ProjectionOriginLPSZero;
+
+    double
+    m_ProjectionAngleLPS,
+    m_SCD,
+    m_PanelOffsetPGeo;
+
+    DirectionType
+        m_ImageDirectionsLPS;
+
+    /** Default Constructor **/
+    DRTImageMetaInformation ();
+    /** Default Destructor **/
+    virtual ~DRTImageMetaInformation ();
+
+private:
+    /** purposely not implemented **/
+    DRTImageMetaInformation(const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+class CTVolumeImageMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef CTVolumeImageMetaInformation  Self;
+    typedef itk::Object Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::Matrix<double, 3, 3> DirectionType;
+    typedef itk::Point<double, 3> PointType;
+    typedef itk::Vector<double, 3> SpacingType;
+    typedef itk::Size<3> SizeType;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(CTVolumeImageMetaInformation , itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    void SetPatientOrientation(tPatOrientation);
+    itkGetEnumMacro(PatientOrientation, tPatOrientation);
+
+    itkSetMacro(Spacing,SpacingType);
+    itkGetMacro(Spacing,SpacingType);
+
+    itkSetMacro(Size,SizeType);
+    itkGetMacro(Size,SizeType);
+
+    itkSetMacro(OriginLPS,PointType);
+    itkGetMacro(OriginLPS,PointType);
+
+    itkSetMacro(ImportOffset,PointType);
+    itkGetMacro(ImportOffset,PointType);
+
+    itkSetMacro(ImageDirections,DirectionType);
+    itkGetMacro(ImageDirections,DirectionType);
+
+    itkGetMacro(LPS2IECDirections,DirectionType);
+
+    PointType GetCOV();
+    PointType GetOriginWOffset();
+
+    const std::vector <double> CalculateRegions(PointType pProjectionOriginLPS);
+
+
+    PointType ConvertIECPointToLPSPoint(
+            PointType IECPoint);
+
+    PointType GetProjectionOriginLPS(
+            PointType ProjectionCenter);
+
+    PointType GetProjectionOriginLPSZero(
+            PointType ProjectionCenter);
+
+protected:
+
+    tPatOrientation
+    m_PatientOrientation;
+
+    SpacingType
+    m_Spacing;
+
+    SizeType
+    m_Size;
+
+    PointType
+    m_OriginLPS,
+    m_ImportOffset;
+
+    DirectionType
+    m_ImageDirections,
+    m_LPS2IECDirections;
+
+    /** Default Constructor **/
+    CTVolumeImageMetaInformation  ();
+    /** Default Destructor **/
+    virtual ~CTVolumeImageMetaInformation  ();
+
+private:
+    /** purposely not implemented **/
+    CTVolumeImageMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+class DRTProjectionGeometryImageMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef DRTProjectionGeometryImageMetaInformation   Self;
+    typedef itk::Object Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::Point<double, 3> PointType;
+    typedef itk::Vector<double, 3> SpacingType;
+    typedef itk::Size<3> SizeType;
+
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(DRTProjectionGeometryImageMetaInformation  , itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+    itkSetMacro(ProjectionAngle1IEC, double);
+    itkGetMacro(ProjectionAngle1IEC, double);
+
+    itkSetMacro(ProjectionAngle2IEC, double);
+    itkGetMacro(ProjectionAngle2IEC, double);
+
+    itkSetMacro(ProjectionAngle1OffsetIEC, double);
+    itkGetMacro(ProjectionAngle1OffsetIEC, double);
+
+    itkSetMacro(ProjectionAngle2OffsetIEC, double);
+    itkGetMacro(ProjectionAngle2OffsetIEC, double);
+
+    itkSetMacro(SCD1, double);
+    itkGetMacro(SCD1, double);
+
+    itkSetMacro(SCD2, double);
+    itkGetMacro(SCD2, double);
+
+    itkSetMacro(SCD1Offset, double);
+    itkGetMacro(SCD1Offset, double);
+
+    itkSetMacro(SCD2Offset, double);
+    itkGetMacro(SCD2Offset, double);
+
+    itkSetMacro(Panel1Offset, double);
+    itkGetMacro(Panel1Offset, double);
+
+    itkSetMacro(Panel2Offset, double);
+    itkGetMacro(Panel2Offset, double);
+
+    itkSetMacro(IntensityThreshold, double);
+    itkGetMacro(IntensityThreshold, double);
+
+    itkSetMacro(DRT1Size, SizeType);
+    itkGetMacro(DRT1Size, SizeType);
+
+    itkSetMacro(DRT2Size, SizeType);
+    itkGetMacro(DRT2Size, SizeType);
+
+    itkSetMacro(DRT1Spacing, SpacingType);
+    itkGetMacro(DRT1Spacing, SpacingType);
+
+    itkSetMacro(DRT2Spacing, SpacingType);
+    itkGetMacro(DRT2Spacing, SpacingType);
+
+    itkSetMacro(IECS2IECScannerT, PointType);
+    itkGetMacro(IECS2IECScannerT, PointType);
+
+    itkSetMacro(IECS2IECScannerR, PointType);
+    itkGetMacro(IECS2IECScannerR, PointType);
+
+    itkSetMacro(ProjectionCenterOffset1, PointType);
+    itkGetMacro(ProjectionCenterOffset1, PointType);
+
+    itkSetMacro(ProjectionCenterOffset2, PointType);
+    itkGetMacro(ProjectionCenterOffset2, PointType);
+
+    itkSetMacro(ProjectionCenter, PointType);
+    itkGetMacro(ProjectionCenter, PointType);
+
+    itkSetMacro(UseRotationsForHiddenTransform, bool);
+    itkGetMacro(UseRotationsForHiddenTransform, bool);
+
+    itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
+    itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
+
+    itkSetEnumMacro(HandleRotationImpOffset, tHandlingRotImpTransform);
+    itkGetEnumMacro(HandleRotationImpOffset, tHandlingRotImpTransform);
+
+protected:
+
+    double
+    m_ProjectionAngle1IEC,
+    m_ProjectionAngle2IEC,
+    m_ProjectionAngle1OffsetIEC,
+    m_ProjectionAngle2OffsetIEC,
+    m_SCD1,
+    m_SCD2,
+    m_SCD1Offset,
+    m_SCD2Offset,
+    m_Panel1Offset,
+    m_Panel2Offset,
+    m_IntensityThreshold;
+
+
+    SizeType
+    m_DRT1Size,
+    m_DRT2Size;
+
+    SpacingType
+    m_DRT1Spacing,
+    m_DRT2Spacing;
+
+
+    PointType
+    /* Transform between IEC Support and
+     * IEC Scanner frame of reference */
+    m_IECS2IECScannerT,
+    m_IECS2IECScannerR;
+
+    PointType
+    /* center of projection in an IEC reference at
+     * Patient Origin of fixed images. Positionin scanner */
+    m_ProjectionCenter,
+    m_ProjectionCenterOffset1,
+    m_ProjectionCenterOffset2;
+
+
+    bool
+    m_UseRotationsForHiddenTransform;
+
+    tDegreeOfFreedomEnum
+    m_DegreeOfFreedom;
+
+    tHandlingRotImpTransform
+    m_HandleRotationImpOffset;
+
+    /** Default Constructor **/
+    DRTProjectionGeometryImageMetaInformation   ();
+    /** Default Destructor **/
+    virtual ~DRTProjectionGeometryImageMetaInformation   ();
+
+private:
+    /** purposely not implemented **/
+    DRTProjectionGeometryImageMetaInformation  (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+
+class RTGeometryMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef RTGeometryMetaInformation Self;
+    typedef itk::Object Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::Point<double, 3> PointType;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(RTGeometryMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    itkSetMacro(IsocenterLPS,PointType);
+    itkGetMacro(IsocenterLPS,PointType);
+
+    itkSetMacro(IsocenterIECS,PointType);
+    itkGetMacro(IsocenterIECS,PointType);
+
+
+
+protected:
+
+
+    PointType
+    /* RT Plan isocenter in LPS coordinates */
+    m_IsocenterLPS,
+    /* RT Plan couch SETUP offset corresponding to
+     * the RT Plan isocenter */
+    m_IsocenterIECS;
+
+    /** Default Constructor **/
+    RTGeometryMetaInformation ();
+    /** Default Destructor **/
+    virtual ~RTGeometryMetaInformation ();
+
+private:
+    /** purposely not implemented **/
+    RTGeometryMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+
+
+
+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 override;
+
+     itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
+     itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
+
+     itkSetEnumMacro(OptimizerType, tOptimizerTypeEnum);
+     itkGetEnumMacro(OptimizerType, tOptimizerTypeEnum);
+
+
+     itkSetEnumMacro(MetricType, tMetricTypeEnum);
+     itkGetEnumMacro(MetricType, tMetricTypeEnum);
+
+
+
+
+protected:
+
+    tDegreeOfFreedomEnum
+    m_DegreeOfFreedom;
+
+    tOptimizerTypeEnum
+    m_OptimizerType;
+
+    tMetricTypeEnum
+    m_MetricType;
+
+    /** 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{
+
+public:
+    /** standard typedefs **/
+    typedef TransformMetaInformation 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(TransformMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    itkSetMacro(HiddenTranslations,PointType);
+    itkGetMacro(HiddenTranslations,PointType);
+
+    itkSetMacro(HiddenRotations,PointType);
+    itkGetMacro(HiddenRotations,PointType);
+
+    itkSetMacro(UserTranslations,PointType);
+    itkGetMacro(UserTranslations,PointType);
+
+    itkSetMacro(UserRotations,PointType);
+    itkGetMacro(UserRotations,PointType);
+
+    PointType GetT();
+    PointType GetR();
+
+protected:
+
+    PointType
+    m_HiddenTranslations,
+    m_HiddenRotations,
+    m_UserTranslations,
+    m_UserRotations;
+
+
+    /** Default Constructor **/
+    TransformMetaInformation ();
+    /** Default Destructor **/
+    virtual ~TransformMetaInformation ();
+
+private:
+    /** purposely not implemented **/
+    TransformMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+};
+
+#endif

itkReg23DRT/autoreg/CMakeLists.txt

@@ -0,0 +1,31 @@
+SET(LIB_NAME autoreg)
+#SET(CMAKE_CXX_FLAGS "-std=c++11 -fPIC")
+
+find_package(ITK REQUIRED)
+include(${ITK_USE_FILE})
+INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
+
+
+SET(HDR
+itkDRTHelpers.h
+itkgTwoImageToOneImageMetric.h
+itkgTwoImageToOneImageMetric.hxx
+itkMutualInformationTwoImageToOneImageMetric.h
+itkMutualInformationTwoImageToOneImageMetric.hxx
+itkNormalizedCorrelationTwoImageToOneImageMetric.h
+itkNormalizedCorrelationTwoImageToOneImageMetric.hxx
+itkTwoProjectionImageRegistrationMethod.h
+itkTwoProjectionImageRegistrationMethod.hxx
+)
+
+ADD_LIBRARY(${LIB_NAME} ${HDR})
+
+target_link_libraries(
+    ${LIB_NAME}
+    ${ITK_LIBRARIES}
+)
+
+target_include_directories (${LIB_NAME}
+    PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+
+

itkReg23DRT/autoreg/itkDRTHelpers.h

@@ -0,0 +1,99 @@
+#ifndef ITKDRTHELPERS_H
+#define ITKDRTHELPERS_H
+
+#include "itkMath.h"
+#include <cstddef>
+#include <cstring>
+#include <iterator>
+
+#ifdef __GNUC__
+#define VARIABLE_IS_NOT_USED __attribute__ ((unused))
+#else
+#define VARIABLE_IS_NOT_USED
+#endif
+
+//Function to get method and class name for logging purposes.
+template <size_t FL, size_t PFL>
+const char* computeMethodName(const char (&function)[FL], const char (&prettyFunction)[PFL])
+{
+    using reverse_ptr = std::reverse_iterator<const char*>;
+    thread_local static char result[PFL];
+    const char* locFuncName = std::search(prettyFunction, prettyFunction + PFL - 1, function, function + FL - 1);
+    const char* locClassName = std::find(reverse_ptr(locFuncName), reverse_ptr(prettyFunction), ' ').base();
+    const char* endFuncName = std::find(locFuncName, prettyFunction + PFL - 1, '(');
+    result[0] = '\0';
+    std::strncat(result, locClassName, endFuncName - locClassName + 1);
+    std::strcat(result, ")");
+    return result;
+}
+#define __COMPACT_PRETTY_FUNCTION__ computeMethodName(__FUNCTION__, __PRETTY_FUNCTION__)
+
+namespace itk {
+/* reference string required for comparison with tag values */
+static const char VARIABLE_IS_NOT_USED *ImageOrientationStrings[] = {
+    "NotDefined",
+    "HFS",
+    "FFS",
+    "HFP",
+    "FFP",
+};
+
+// constant for converting degrees into radians
+const float dtr = itk::Math::pi_over_180;
+
+static const bool verbose = false;
+
+/* this is in the end a IEC to HFS...
+ * but we keep this for ourselves...
+ */
+static double VARIABLE_IS_NOT_USED Standard_DRT2LPS[9] = {
+    1, 0, 0,
+    0, 0, -1,
+    0, 1, 0
+};
+
+static double VARIABLE_IS_NOT_USED PAElementsIEC[9] = {
+    1, 0, 0,
+    0, -1, 0,
+    0, 0, -1
+};
+
+static double VARIABLE_IS_NOT_USED LATElementsIEC[9] = {
+    0, 0, -1,
+    0, -1, 0,
+    -1, 0, 0
+};
+
+static double VARIABLE_IS_NOT_USED HFS2IEC[9] = {
+    1, 0, 0,
+    0, 0, 1,
+    0, -1, 0
+};
+
+static double VARIABLE_IS_NOT_USED FFS2IEC[9] = {
+    -1, 0, 0,
+    0, 0, -1,
+    0, -1, 0
+};
+
+static double VARIABLE_IS_NOT_USED HFP2IEC[9] = {
+    -1, 0, 0,
+    0, 0, 1,
+    0, 1, 0
+};
+
+static double VARIABLE_IS_NOT_USED FFP2IEC[9] = {
+    1, 0, 0,
+    0, 0, -1,
+    0, 1, 0
+};
+
+static double VARIABLE_IS_NOT_USED PAT2IEC[9] = {
+    1, 0, 0,
+    0, 0, 1,
+    0, -1, 0
+};
+
+}
+
+#endif // ITKDRTHELPERS_H

itkReg23DRT/autoreg/itkMutualInformationTwoImageToOneImageMetric.h

@@ -0,0 +1,119 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkMutualInformationTwoImageToOneImageMetric_h
+#define itkMutualInformationTwoImageToOneImageMetric_h
+
+#include "itkCovariantVector.h"
+#include "itkMacro.h"
+#include "itkPoint.h"
+#include "itkgTwoImageToOneImageMetric.h"
+
+namespace itk {
+/** \class NormalizedCorrelationTwoImageToOneImageMetric
+ * \brief Computes similarity between two fixed images and one moving image
+ *
+ * This metric computes the correlation between pixels in the two fixed images
+ * and pixels in the moving image. The spatial correspondance between
+ * two fixed images and the moving image is established through a Transform. Pixel values are
+ * taken from the fixed images, their positions are mapped to the moving
+ * image and result in general in non-grid position on it. Values at these
+ * non-grid position of the moving image are interpolated using user-selected
+ * Interpolators. The correlation is normalized by the autocorrelations of both
+ * the fixed and moving images.
+ *
+ * \ingroup RegistrationMetrics
+ * \ingroup TwoProjectionRegistration
+ */
+template <typename TFixedImage, typename TMovingImage>
+class MutualInformationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric<TFixedImage, TMovingImage> {
+public:
+    ITK_DISALLOW_COPY_AND_ASSIGN(MutualInformationTwoImageToOneImageMetric);
+
+    /** Standard class type alias. */
+    using Self = MutualInformationTwoImageToOneImageMetric;
+    using Superclass = gTwoImageToOneImageMetric<TFixedImage, TMovingImage>;
+
+    using Pointer = SmartPointer<Self>;
+    using ConstPointer = SmartPointer<const Self>;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(MutualInformationTwoImageToOneImageMetric, Object);
+
+    /** Types transferred from the base class */
+    using RealType = typename Superclass::RealType;
+    using TransformType = typename Superclass::TransformType;
+    using TransformPointer = typename Superclass::TransformPointer;
+    using TransformParametersType = typename Superclass::TransformParametersType;
+    using TransformJacobianType = typename Superclass::TransformJacobianType;
+    using GradientPixelType = typename Superclass::GradientPixelType;
+
+    using MeasureType = typename Superclass::MeasureType;
+    using DerivativeType = typename Superclass::DerivativeType;
+    using FixedImageType = typename Superclass::FixedImageType;
+    using MovingImageType = typename Superclass::MovingImageType;
+    using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
+    using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
+
+    /** Get the derivatives of the match measure. */
+    void
+    GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override;
+
+    /**  Get the value for single valued optimizers. */
+    MeasureType
+    GetValue(const TransformParametersType& parameters) const override;
+
+    /**  Get the value using the current transforms. */
+    MeasureType
+    GetValue() const;
+
+    /**  Get value and derivatives for multiple valued optimizers. */
+    void
+    GetValueAndDerivative(const TransformParametersType& parameters,
+        MeasureType& Value,
+        DerivativeType& Derivative) const override;
+
+    /** Set/Get SubtractMean boolean. If true, the sample mean is subtracted
+   * from the sample values in the cross-correlation formula and
+   * typically results in narrower valleys in the cost fucntion.
+   * Default value is false. */
+    itkSetMacro(SubtractMean, bool);
+    itkSetMacro(NumberOfHistogramBins, int);
+    itkGetConstReferenceMacro(SubtractMean, bool);
+    itkBooleanMacro(SubtractMean);
+
+protected:
+    MutualInformationTwoImageToOneImageMetric();
+    ~MutualInformationTwoImageToOneImageMetric() override = default;
+    void
+    PrintSelf(std::ostream& os, Indent indent) const override;
+
+private:
+    bool m_SubtractMean;
+    int m_NumberOfHistogramBins;
+};
+
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkMutualInformationTwoImageToOneImageMetric.hxx"
+#endif
+
+#endif

itkReg23DRT/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx

@@ -0,0 +1,239 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkMutualInformationTwoImageToOneImageMetric_hxx
+#define itkMutualInformationTwoImageToOneImageMetric_hxx
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkMattesMutualInformationImageToImageMetric.h"
+#include "itkMutualInformationTwoImageToOneImageMetric.h"
+#include "itkNearestNeighborInterpolateImageFunction.h"
+#include "itkNormalizeImageFilter.h"
+#include "itkTranslationTransform.h"
+
+#include "itkProgressReporter.h"
+
+#include <limits>
+
+namespace itk {
+
+template <typename TFixedImage, typename TMovingImage>
+MutualInformationTwoImageToOneImageMetric<TFixedImage,
+    TMovingImage>::MutualInformationTwoImageToOneImageMetric()
+{
+    m_SubtractMean = false;
+    m_NumberOfHistogramBins = 50;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+typename MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
+MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue(
+    const TransformParametersType& parameters) const
+{
+    FixedImageConstPointer fixedImage1 = this->m_FixedImage1;
+
+    if (!fixedImage1) {
+        itkExceptionMacro(<< "Fixed image1 has not been assigned");
+    }
+
+    FixedImageConstPointer fixedImage2 = this->m_FixedImage2;
+
+    if (!fixedImage2) {
+        itkExceptionMacro(<< "Fixed image2 has not been assigned");
+    }
+
+    if (!this->m_Filter1) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    if (!this->m_Filter2) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    this->SetTransformParameters(parameters);
+
+    return GetValue();
+}
+
+template <typename TFixedImage, typename TMovingImage>
+typename MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
+MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue() const
+{
+
+    FixedImageConstPointer fixedImage1 = this->m_FixedImage1;
+
+    if (!fixedImage1) {
+        itkExceptionMacro(<< "Fixed image1 has not been assigned");
+    }
+
+    FixedImageConstPointer fixedImage2 = this->m_FixedImage2;
+
+    if (!fixedImage2) {
+        itkExceptionMacro(<< "Fixed image2 has not been assigned");
+    }
+
+    if (!this->m_Filter1) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    if (!this->m_Filter2) {
+        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>;
+    /*
+     * 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 = LinearInterpType::New();
+
+    auto metric1 = InternalMetricType::New();
+    metric1->SetTransform(dummyTransform);
+    metric1->SetTransformParameters(dummyParameters);
+    metric1->SetInterpolator(dummyInterpolator1);
+
+    //auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
+    //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%
+    // since we have a ROI, then we should not set allPixels to TRUE.
+    //metric1->UseAllPixelsOn();
+    //std::cout << "m_NumberOfHistogramBins " << m_NumberOfHistogramBins << std::endl;
+    metric1->SetNumberOfHistogramBins(m_NumberOfHistogramBins);
+
+//    InternalImageType::IndexType pIndex;
+//    pIndex[0]=200;
+//    pIndex[1]=300;
+//    InternalImageType::PixelType pPixel = fixedImage1->GetPixel(pIndex);
+//    std::cout<<"MATTES PIXEL: "<<pPixel <<std::endl;
+
+    metric1->SetFixedImage(fixedImage1);
+    metric1->SetMovingImage(this->m_Filter1->GetOutput());
+
+    metric1->Initialize();
+
+    auto measure1 = metric1->GetValue(dummyParameters);
+
+    this->m_Interpolator2->SetTransform(this->m_Transform2);
+    this->m_Interpolator2->Initialize();
+    this->m_Filter2->Update();
+
+    auto dummyInterpolator2 = LinearInterpType::New();
+    auto metric2 = InternalMetricType::New();
+    metric2->SetTransform(dummyTransform);
+    metric2->SetTransformParameters(dummyParameters);
+    metric2->SetInterpolator(dummyInterpolator2);
+
+//    auto fixedImageRegion2 = fixedImage2->GetBufferedRegion();
+//    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->UseAllPixelsOn();
+    metric2->SetNumberOfHistogramBins(m_NumberOfHistogramBins);
+    metric2->SetFixedImage(fixedImage2);
+    metric2->SetMovingImage(this->m_Filter2->GetOutput());
+    metric2->Initialize();
+
+    auto measure2 = metric2->GetValue(dummyParameters);
+
+    //std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl;
+    //std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl;
+
+    // Calculate the measure value between fixed image 2 and the moving image
+
+    auto measure = (measure1 + measure2) / 2.0;
+
+//    auto oldprecision = std::cout.precision();
+//    std::cout.precision(std::numeric_limits<double>::digits10 + 2);
+//    std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
+//    std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
+//    std::cout << "Measure: " << measure << std::endl;
+//    std::cout << "=======================================================" << std::endl;
+//    std::cout.precision(oldprecision);
+    return measure;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetDerivative(
+    const TransformParametersType& itkNotUsed(parameters),
+    DerivativeType& itkNotUsed(derivative)) const
+{
+    // under construction
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValueAndDerivative(
+    const TransformParametersType& itkNotUsed(parameters),
+    MeasureType& itkNotUsed(value),
+    DerivativeType& itkNotUsed(derivative)) const
+{
+    // under construction
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os,
+    Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+    os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
+}
+
+} // end namespace itk
+
+#endif

itkReg23DRT/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h

@@ -0,0 +1,127 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkNormalizedCorrelationTwoImageToOneImageMetric_h
+#define itkNormalizedCorrelationTwoImageToOneImageMetric_h
+
+#include "itkCovariantVector.h"
+#include "itkMacro.h"
+#include "itkNormalizedCorrelationImageToImageMetric.hxx"
+#include "itkPoint.h"
+#include "itkgTwoImageToOneImageMetric.h"
+
+namespace itk {
+/** \class NormalizedCorrelationTwoImageToOneImageMetric
+ * \brief Computes similarity between two fixed images and one moving image
+ *
+ * This metric computes the correlation between pixels in the two fixed images
+ * and pixels in the moving image. The spatial correspondance between
+ * two fixed images and the moving image is established through a Transform. Pixel values are
+ * taken from the fixed images, their positions are mapped to the moving
+ * image and result in general in non-grid position on it. Values at these
+ * non-grid position of the moving image are interpolated using user-selected
+ * Interpolators. The correlation is normalized by the autocorrelations of both
+ * the fixed and moving images.
+ *
+ * \ingroup RegistrationMetrics
+ * \ingroup TwoProjectionRegistration
+ */
+template <typename TFixedImage, typename TMovingImage>
+class NormalizedCorrelationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric<TFixedImage, TMovingImage> {
+public:
+    ITK_DISALLOW_COPY_AND_ASSIGN(NormalizedCorrelationTwoImageToOneImageMetric);
+
+    /** Standard class type alias. */
+    using Self = NormalizedCorrelationTwoImageToOneImageMetric;
+    using Superclass = gTwoImageToOneImageMetric<TFixedImage, TMovingImage>;
+
+    using Pointer = SmartPointer<Self>;
+    using ConstPointer = SmartPointer<const Self>;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(NormalizedCorrelationTwoImageToOneImageMetric, Object);
+
+    /** Types transferred from the base class */
+    using RealType = typename Superclass::RealType;
+    using TransformType = typename Superclass::TransformType;
+    using TransformPointer = typename Superclass::TransformPointer;
+    using TransformParametersType = typename Superclass::TransformParametersType;
+    using TransformJacobianType = typename Superclass::TransformJacobianType;
+    using GradientPixelType = typename Superclass::GradientPixelType;
+
+    using MeasureType = typename Superclass::MeasureType;
+    using DerivativeType = typename Superclass::DerivativeType;
+
+    using FixedImageType = typename Superclass::FixedImageType;
+    using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
+    using FixedImageMaskPointer = typename Superclass::FixedImageMaskPointer;
+    using FixedImageRegionType = typename Superclass::FixedImageRegionType;
+
+    using MovingImageType = typename Superclass::MovingImageType;
+    using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
+    using MovingImageMaskPointer = typename Superclass::MovingImageMaskPointer;
+
+    using ChangeInformationFilterPointer = typename Superclass::ChangeInformationFilterPointer;
+
+    /** Get the derivatives of the match measure. */
+    void
+    GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override;
+
+    /**  Get the value for single valued optimizers, after applying paramters to the transforms. */
+    MeasureType
+    GetValue(const TransformParametersType& parameters) const override;
+
+    /**  Get the value using the current transforms. */
+    MeasureType
+    GetValue() const;
+
+    /**  Get value and derivatives for multiple valued optimizers. */
+    void
+    GetValueAndDerivative(const TransformParametersType& parameters,
+        MeasureType& Value,
+        DerivativeType& Derivative) const override;
+
+    /** Set/Get SubtractMean boolean. If true, the sample mean is subtracted
+   * from the sample values in the cross-correlation formula and
+   * typically results in narrower valleys in the cost fucntion.
+   * Default value is false. */
+    itkSetMacro(SubtractMean, bool);
+    itkGetConstReferenceMacro(SubtractMean, bool);
+    itkBooleanMacro(SubtractMean);
+
+protected:
+    NormalizedCorrelationTwoImageToOneImageMetric();
+    ~NormalizedCorrelationTwoImageToOneImageMetric() override = default;
+    void
+    PrintSelf(std::ostream& os, Indent indent) const override;
+
+private:
+    bool m_SubtractMean;
+
+    MeasureType CalculateMeasure(int imageId) const;
+};
+
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkNormalizedCorrelationTwoImageToOneImageMetric.hxx"
+#endif
+
+#endif

itkReg23DRT/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx

@@ -0,0 +1,317 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkNormalizedCorrelationTwoImageToOneImageMetric_hxx
+#define itkNormalizedCorrelationTwoImageToOneImageMetric_hxx
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
+
+#
+
+namespace itk {
+
+template <typename TFixedImage, typename TMovingImage>
+NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage,
+    TMovingImage>::NormalizedCorrelationTwoImageToOneImageMetric()
+{
+    m_SubtractMean = false;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
+NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::CalculateMeasure(int imageId) const
+{
+
+    using FixedIteratorType = itk::ImageRegionConstIteratorWithIndex<FixedImageType>;
+    using MovingIteratorType = itk::ImageRegionConstIteratorWithIndex<MovingImageType>;
+    using AccumulateType = typename NumericTraits<MeasureType>::AccumulateType;
+
+    FixedImageConstPointer fixedImage;
+    ChangeInformationFilterPointer filter;
+    FixedImageMaskPointer fixedImageMask;
+    FixedImageRegionType fixedImageRegion;
+    MeasureType measure;
+
+    typename FixedImageType::IndexType indexFixed;
+    typename MovingImageType::IndexType indexMoving;
+    typename Superclass::InputPointType inputPoint;
+
+    if (imageId == 1) {
+        this->m_Interpolator1->SetTransform(this->m_Transform1);
+        this->m_Interpolator1->Initialize();
+
+        fixedImage = this->m_FixedImage1;
+        fixedImageMask = this->m_FixedImageMask1;
+        fixedImageRegion = this->GetFixedImageRegion1();
+        filter = this->m_Filter1;
+    } else if (imageId == 2) {
+        this->m_Interpolator2->SetTransform(this->m_Transform2);
+        this->m_Interpolator2->Initialize();
+
+        fixedImage = this->m_FixedImage2;
+        fixedImageMask = this->m_FixedImageMask2;
+        fixedImageRegion = this->GetFixedImageRegion2();
+        filter = this->m_Filter2;
+    } else {
+        itkExceptionMacro(<< "Non supported image id.");
+    }
+
+    FixedIteratorType ti(fixedImage, fixedImageRegion);
+
+    this->m_NumberOfPixelsCounted = 0;
+    int numberOfPixelsVisited = 0;
+
+    filter->Update();
+    MovingImageConstPointer imageDRTIn = filter->GetOutput();
+    //this->WriteImage(imageDRTIn, "", "");
+
+    AccumulateType sff = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType smm = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType sfm = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType sf = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType sm = NumericTraits<AccumulateType>::ZeroValue();
+
+    RealType movingValueMin = NumericTraits<RealType>::max();
+    RealType fixedValueMin = NumericTraits<RealType>::max();
+    RealType movingValueMax = NumericTraits<RealType>::NonpositiveMin();
+    RealType fixedValueMax = NumericTraits<RealType>::NonpositiveMin();
+
+    RealType maxX = NumericTraits<RealType>::NonpositiveMin();
+    RealType maxY = NumericTraits<RealType>::NonpositiveMin();
+    RealType maxZ = NumericTraits<RealType>::NonpositiveMin();
+    RealType minX = NumericTraits<RealType>::max();
+    RealType minY = NumericTraits<RealType>::max();
+    RealType minZ = NumericTraits<RealType>::max();
+
+    AccumulateType m_meanCurr_f = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_meanCurr_m = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_meanPrev_f = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_meanPrev_m = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_sPrev_f = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_sPrev_m = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_sCurr_f = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_sCurr_m = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_varianceCurr_f = NumericTraits<AccumulateType>::ZeroValue();
+    AccumulateType m_varianceCurr_m = NumericTraits<AccumulateType>::ZeroValue();
+
+    while (!ti.IsAtEnd()) {
+
+        indexFixed = ti.GetIndex();
+        ++numberOfPixelsVisited;
+
+        fixedImage->TransformIndexToPhysicalPoint(indexFixed, inputPoint);
+
+        if (fixedImageMask && !fixedImageMask->IsInsideInWorldSpace(inputPoint)) {
+            ++ti;
+            continue;
+        }
+
+        if (this->m_MovingImageMask && !this->m_MovingImageMask->IsInsideInWorldSpace(inputPoint)) {
+            ++ti;
+            continue;
+        }
+
+        const RealType fixedValue = ti.Get();
+
+        if (fixedValue == -1024) {
+            ++ti;
+            continue;
+        }
+
+        if (imageDRTIn->TransformPhysicalPointToIndex(inputPoint, indexMoving)) {
+
+            double x = inputPoint[0];
+            double y = inputPoint[1];
+            double z = inputPoint[2];
+            maxX = maxX > x ? maxX : x;
+            minX = minX < x ? minX : x;
+            maxY = maxY > y ? maxY : y;
+            minY = minY < y ? minY : y;
+            maxZ = maxZ > z ? maxZ : z;
+            minZ = minZ < z ? minZ : z;
+
+            const RealType movingValue = imageDRTIn->GetPixel(indexMoving);
+
+            if (NumericTraits<short>::NonpositiveMin() == movingValue) {
+                ++ti;
+                continue;
+            }
+
+            movingValueMin = movingValue < movingValueMin ? movingValue : movingValueMin;
+            movingValueMax = movingValue > movingValueMax ? movingValue : movingValueMax;
+            fixedValueMin = fixedValue < fixedValueMin ? fixedValue : fixedValueMin;
+            fixedValueMax = fixedValue > fixedValueMax ? fixedValue : fixedValueMax;
+
+            sff += fixedValue * fixedValue;
+            smm += movingValue * movingValue;
+            sfm += fixedValue * movingValue;
+            if (this->m_SubtractMean) {
+                sf += fixedValue;
+                sm += movingValue;
+            }
+            this->m_NumberOfPixelsCounted++;
+
+            if (this->m_NumberOfPixelsCounted == 1) {
+                m_meanCurr_f = fixedValue;
+                m_meanCurr_m = movingValue;
+            } else {
+                m_meanPrev_f = m_meanCurr_f;
+                m_meanPrev_m = m_meanCurr_m;
+                m_sPrev_f = m_sCurr_f;
+                m_sPrev_m = m_sCurr_m;
+
+                m_meanCurr_f = m_meanPrev_f + (fixedValue - m_meanPrev_f) / this->m_NumberOfPixelsCounted;
+                m_sCurr_f = m_sPrev_f + (fixedValue - m_meanPrev_f) * (fixedValue - m_meanCurr_f);
+                m_varianceCurr_f = m_sCurr_f / (this->m_NumberOfPixelsCounted - 1);
+
+                m_meanCurr_m = m_meanPrev_m + (movingValue - m_meanPrev_m) / this->m_NumberOfPixelsCounted;
+                m_sCurr_m = m_sPrev_m + (movingValue - m_meanPrev_m) * (movingValue - m_meanCurr_m);
+                m_varianceCurr_m = m_sCurr_m / (this->m_NumberOfPixelsCounted - 1);
+            }
+        }
+
+        ++ti;
+    }
+
+    if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0) {
+        sff -= (sf * sf / this->m_NumberOfPixelsCounted);
+        smm -= (sm * sm / this->m_NumberOfPixelsCounted);
+        sfm -= (sf * sm / this->m_NumberOfPixelsCounted);
+    }
+
+    RealType denom = -1.0 * sqrt(sff * smm);
+
+    if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
+        measure = sfm / denom;
+    } else {
+        measure = NumericTraits<MeasureType>::Zero;
+    }
+
+//    std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl;
+//    std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl;
+//    std::cout << "max coordinates: " << maxX << " " << maxX << " " << maxZ << std::endl;
+//    std::cout << "min coordinates: " << minX << " " << minY << " " << minZ << std::endl;
+
+    return measure;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
+NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue(
+    const TransformParametersType& parameters) const
+{
+    if (!this->m_FixedImage1) {
+        itkExceptionMacro(<< "Fixed image1 has not been assigned");
+    }
+
+    if (!this->m_FixedImage2) {
+        itkExceptionMacro(<< "Fixed image2 has not been assigned");
+    }
+
+    if (!this->m_Filter1) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    if (!this->m_Filter2) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    if (!this->SetTransformParameters(parameters)) {
+        std::cout << "Measure: 0.0, transform Invalid!" << std::endl;
+        std::cout << "=======================================================" << std::endl;
+
+        return 0.0;
+    }
+    return GetValue();
+}
+
+template <typename TFixedImage, typename TMovingImage>
+typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
+NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue() const
+{
+
+    if (!this->m_FixedImage1) {
+        itkExceptionMacro(<< "Fixed image1 has not been assigned");
+    }
+
+    if (!this->m_FixedImage2) {
+        itkExceptionMacro(<< "Fixed image2 has not been assigned");
+    }
+
+    if (!this->m_Filter1) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    if (!this->m_Filter2) {
+        itkExceptionMacro(<< "Fixed filter1 has not been assigned");
+    }
+
+    // Calculate the measure value between fixed image 1 and the moving image
+    //auto oldprecision = std::cout.precision();
+
+
+//    std::cout<<"Region "  <<this->GetFixedImageRegion1() <<std::endl;
+//    std::cout<<"Buffered Region " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
+//    std::cout<<"Buffered Region " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
+
+    MeasureType measure1 = CalculateMeasure(1);
+//    std::cout.precision(std::numeric_limits<double>::digits10 + 2);
+//    std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
+//    std::cout.precision(oldprecision);
+
+    MeasureType measure2 = CalculateMeasure(2);
+//    std::cout.precision(std::numeric_limits<double>::digits10 + 2);
+//    std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
+
+    auto measure = (measure1 + measure2) / 2.0;
+
+//    std::cout << "Measure: " << measure << std::endl;
+//    std::cout << "=======================================================" << std::endl;
+//    std::cout.precision(oldprecision);
+
+    return measure;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetDerivative(
+    const TransformParametersType& itkNotUsed(parameters),
+    DerivativeType& itkNotUsed(derivative)) const
+{
+    // under construction
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValueAndDerivative(
+    const TransformParametersType& itkNotUsed(parameters),
+    MeasureType& itkNotUsed(value),
+    DerivativeType& itkNotUsed(derivative)) const
+{
+    // under construction
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os,
+    Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+    os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
+}
+
+} // end namespace itk
+
+#endif

itkReg23DRT/autoreg/itkTwoProjectionImageRegistrationMethod.h

@@ -0,0 +1,294 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkTwoProjectionImageRegistrationMethod_h
+#define itkTwoProjectionImageRegistrationMethod_h
+
+#include "itkChangeInformationImageFilter.h"
+#include "itkDataObjectDecorator.h"
+#include "itkImage.h"
+#include "itkMacro.h"
+#include "itkProcessObject.h"
+#include "itkSingleValuedNonLinearOptimizer.h"
+#include "itkgTwoImageToOneImageMetric.h"
+
+#include "itkObject.h"
+#include "itkEventObject.h"
+
+#include <iostream>
+
+namespace itk {
+
+/** \class TwoProjectionImageRegistrationMethod
+ * \brief Base class for Image Registration Methods
+ *
+ * This Class define the generic interface for a registration method.
+ *
+ * This class is templated over the type of the two image to be
+ * registered. A generic Transform is used by this class. That allows
+ * to select at run time the particular type of transformation that
+ * is to be applied for registering the images.
+ *
+ * This method use a generic Metric in order to compare the two images.
+ * the final goal of the registration method is to find the set of
+ * parameters of the Transformation that optimizes the metric.
+ *
+ * The registration method also support a generic optimizer that can
+ * be selected at run-time. The only restriction for the optimizer is
+ * that it should be able to operate in single-valued cost functions
+ * given that the metrics used to compare images provide a single
+ * value as output.
+ *
+ * The terms : Fixed image and Moving image are used in this class
+ * to indicate what image is being mapped by the transform.
+ *
+ * This class uses the coordinate system of the Fixed image as a reference
+ * and searchs for a Transform that will map points from the space of the
+ * Fixed image to the space of the Moving image.
+ *
+ * For doing so, a Metric will be continously applied to compare the Fixed
+ * image with the Transformed Moving image. This process also requires to
+ * interpolate values from the Moving image.
+ *
+ * \ingroup RegistrationFilters
+ * \ingroup TwoProjectionRegistration
+ */
+template <typename TFixedImage, typename TMovingImage>
+class TwoProjectionImageRegistrationMethod : public ProcessObject {
+public:
+    ITK_DISALLOW_COPY_AND_ASSIGN(TwoProjectionImageRegistrationMethod);
+
+    /** Standard class type alias. */
+    using Self = TwoProjectionImageRegistrationMethod;
+    using Superclass = ProcessObject;
+    using Pointer = SmartPointer<Self>;
+    using ConstPointer = SmartPointer<const Self>;
+
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self);
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(TwoProjectionImageRegistrationMethod, ProcessObject);
+
+    /**  Type of the Fixed image. */
+    using FixedImageType = TFixedImage;
+    using FixedImageConstPointer = typename FixedImageType::ConstPointer;
+
+    /**  Type of the Moving image. */
+    using MovingImageType = TMovingImage;
+    using MovingImageConstPointer = typename MovingImageType::ConstPointer;
+    using ChangeInformationFilterType = ChangeInformationImageFilter<MovingImageType>;
+    using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer;
+
+    /**  Type of the metric. */
+    using MetricType = gTwoImageToOneImageMetric<FixedImageType, MovingImageType>;
+    using MetricPointer = typename MetricType::Pointer;
+    using FixedImageRegionType = typename MetricType::FixedImageRegionType;
+
+    /**  Type of the Transform . */
+    using TransformType = typename MetricType::TransformType;
+    using TransformPointer = typename TransformType::Pointer;
+
+    /** Type for the output: Using Decorator pattern for enabling
+   *  the Transform to be passed in the data pipeline */
+    using TransformOutputType = DataObjectDecorator<TransformType>;
+    using TransformOutputPointer = typename TransformOutputType::Pointer;
+    using TransformOutputConstPointer = typename TransformOutputType::ConstPointer;
+
+    /**  Type of the Interpolator. */
+    using InterpolatorType = typename MetricType::InterpolatorType;
+    using InterpolatorPointer = typename InterpolatorType::Pointer;
+
+    /**  Type of the optimizer. */
+    using OptimizerType = SingleValuedNonLinearOptimizer;
+
+    /** Type of the Transformation parameters This is the same type used to
+   *  represent the search space of the optimization algorithm */
+    using ParametersType = typename MetricType::TransformParametersType;
+
+    /** Smart Pointer type to a DataObject. */
+    using DataObjectPointer = typename DataObject::Pointer;
+
+    /** Method that initiates the registration. This will Initialize and ensure
+   * that all inputs the registration needs are in place, via a call to
+   * Initialize() will then start the optimization process via a call to
+   * StartOptimization()  */
+    void
+    StartRegistration();
+
+    /** Method that initiates the optimization process. */
+    void
+    StartOptimization();
+
+    /** Set/Get the Fixed images. */
+    void
+    SetFixedImage1(const FixedImageType* fixedImage1);
+    void
+    SetFixedImage2(const FixedImageType* fixedImage2);
+    itkGetConstObjectMacro(FixedImage1, FixedImageType);
+    itkGetConstObjectMacro(FixedImage2, FixedImageType);
+
+    /** Set/Get the Moving image. */
+    void
+    SetMovingImage(const MovingImageType* movingImage);
+    itkGetConstObjectMacro(MovingImage, MovingImageType);
+
+    /** Set/Get the Optimizer. */
+    itkSetObjectMacro(Optimizer, OptimizerType);
+    itkGetConstObjectMacro(Optimizer, OptimizerType);
+
+    /** Set/Get the Metric. */
+    itkSetObjectMacro(Metric, MetricType);
+    itkGetConstObjectMacro(Metric, MetricType);
+
+    /** Set/Get the Transfrom1. */
+    itkSetObjectMacro(Transform1, TransformType);
+    itkGetConstObjectMacro(Transform1, TransformType);
+
+    /** Set/Get the Transfrom2. */
+    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);
+    itkGetConstObjectMacro(Interpolator1, InterpolatorType);
+    itkGetConstObjectMacro(Interpolator2, InterpolatorType);
+
+    /** Set/Get the Meta informations. */
+    itkSetObjectMacro(TransformMetaInfo, R23MetaInformation);
+    itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation);
+
+    itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation);
+    itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation);
+
+    itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation);
+    itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation);
+
+    /** Set/Get the output filters. */
+    itkSetObjectMacro(Filter1, ChangeInformationFilterType);
+    itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
+    itkSetObjectMacro(Filter2, ChangeInformationFilterType);
+    itkGetConstObjectMacro(Filter2, ChangeInformationFilterType);
+
+    /** Set the region of the fixed image to be considered as region of
+   interest during the registration. This region will be passed to
+   the ImageMetric in order to restrict the metric computation to
+   consider only this region.
+   \warning The same region can also be set directly into the metric.
+   please avoid to set the region in both places since this can lead
+   to inconsistent configurations.  */
+    void
+    SetFixedImageRegion1(const FixedImageRegionType& region1);
+    void
+    SetFixedImageRegion2(const FixedImageRegionType& region2);
+    /** Get the region of the fixed image to be considered as region of
+   interest during the registration. This region will be passed to
+   the ImageMetric in order to restrict the metric computation to
+   consider only this region.  */
+    itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType);
+    itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType);
+    /** True if a region has been defined for the fixed image to which
+   the ImageMetric will limit its computation */
+    itkGetMacro(FixedImageRegionDefined1, bool);
+    itkGetMacro(FixedImageRegionDefined2, bool);
+    /** Turn on/off the use of a fixed image region to which
+   the ImageMetric will limit its computation.
+   \warning The region must have been previously defined using the
+   SetFixedImageRegion member function */
+    itkSetMacro(FixedImageRegionDefined1, bool);
+    itkSetMacro(FixedImageRegionDefined2, bool);
+
+    /** Initialize by setting the interconnects between the components. */
+    virtual void
+    Initialize();
+
+    /** Returns the transform resulting from the registration process  */
+    const TransformOutputType*
+    GetOutput() const;
+
+    /** Make a DataObject of the correct type to be used as the specified
+   * output. */
+    using Superclass::MakeOutput;
+    DataObjectPointer
+    MakeOutput(DataObjectPointerArraySizeType idx) override;
+
+protected:
+    TwoProjectionImageRegistrationMethod();
+    ~TwoProjectionImageRegistrationMethod() override = default;
+    void
+    PrintSelf(std::ostream& os, Indent indent) const override;
+
+    /** Method invoked by the pipeline in order to trigger the computation of
+   * the registration. */
+    void
+    GenerateData() override;
+
+    /** Provides derived classes with the ability to set this private var */
+    itkSetMacro(LastOptimizerParameters, ParametersType);
+
+    /** Entry-point for internal ITK filter observer. **/
+//    void OnInternalFilterProgressReceptor(itk::Object *caller,
+//        const itk::EventObject &event);
+
+
+private:
+    MetricPointer m_Metric;
+    OptimizerType::Pointer m_Optimizer;
+
+    MovingImageConstPointer m_MovingImage;
+    FixedImageConstPointer m_FixedImage1;
+    FixedImageConstPointer m_FixedImage2;
+
+    TransformPointer m_IsocIECTransform;
+
+    TransformPointer m_Transform1;
+    TransformPointer m_Transform2;
+    InterpolatorPointer m_Interpolator1;
+    InterpolatorPointer m_Interpolator2;
+
+    InternalTransformMetaInformation::Pointer
+    m_internalMeta1,
+    m_internalMeta2;
+
+    R23MetaInformation::Pointer
+        m_TransformMetaInfo;
+
+    ChangeInformationFilterPointer
+        m_Filter1,
+        m_Filter2;
+
+    ParametersType m_InitialOptimizerParameters;
+    ParametersType m_LastOptimizerParameters;
+
+    bool m_FixedImageRegionDefined1;
+    bool m_FixedImageRegionDefined2;
+    FixedImageRegionType m_FixedImageRegion1;
+    FixedImageRegionType m_FixedImageRegion2;
+};
+
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkTwoProjectionImageRegistrationMethod.hxx"
+#endif
+
+#endif

itkReg23DRT/autoreg/itkTwoProjectionImageRegistrationMethod.hxx

@@ -0,0 +1,381 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkTwoProjectionImageRegistrationMethod_hxx
+#define itkTwoProjectionImageRegistrationMethod_hxx
+
+#include "itkCommand.h"
+
+#include "itkTwoProjectionImageRegistrationMethod.h"
+
+namespace itk {
+
+// TwoProjectionImageRegistrationMethod from Jian Wu.
+// This part is not really needed and functionality should be moved to
+// to the itkImageProcessor.cpp. This is just abstraction with an additional step.
+
+template <typename TFixedImage, typename TMovingImage>
+TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TwoProjectionImageRegistrationMethod()
+{
+    this->SetNumberOfRequiredOutputs(1); // for the Transform
+
+    m_FixedImage1 = nullptr; // has to be provided by the user.
+    m_FixedImage2 = nullptr; // has to be provided by the user.
+    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.
+    m_Optimizer = nullptr; // has to be provided by the user.
+
+    m_InitialOptimizerParameters = ParametersType(1);
+    m_LastOptimizerParameters = ParametersType(1);
+
+    m_InitialOptimizerParameters.Fill(0.0f);
+    m_LastOptimizerParameters.Fill(0.0f);
+
+    m_FixedImageRegionDefined1 = false;
+    m_FixedImageRegionDefined2 = false;
+
+    TransformOutputPointer transformDecorator = static_cast<TransformOutputType*>(this->MakeOutput(0).GetPointer());
+
+    this->ProcessObject::SetNthOutput(0, transformDecorator.GetPointer());
+}
+
+/*
+ * Set the region of the fixed image 1 to be considered for registration
+ */
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImageRegion1(
+    const FixedImageRegionType& region1)
+{
+    m_FixedImageRegion1 = region1;
+    m_FixedImageRegionDefined1 = true;
+}
+
+/*
+ * Set the region of the fixed image 2 to be considered for registration
+ */
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImageRegion2(
+    const FixedImageRegionType& region2)
+{
+    m_FixedImageRegion2 = region2;
+    m_FixedImageRegionDefined2 = true;
+}
+
+
+/*
+ * Initialize by setting the interconnects between components.
+ */
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize()
+{
+
+    if (!m_FixedImage1) {
+        itkExceptionMacro(<< "FixedImage1 is not present");
+    }
+
+    if (!m_FixedImage2) {
+        itkExceptionMacro(<< "FixedImage2 is not present");
+    }
+
+    if (!m_MovingImage) {
+        itkExceptionMacro(<< "MovingImage is not present");
+    }
+
+    if (!m_Metric) {
+        itkExceptionMacro(<< "Metric is not present");
+    }
+
+    if (!m_Optimizer) {
+        itkExceptionMacro(<< "Optimizer is not present");
+    }
+
+    if (!m_Transform1) {
+        itkExceptionMacro(<< "Transform1 is not present");
+    }
+
+    if (!m_Transform2) {
+        itkExceptionMacro(<< "Transform2 is not present");
+    }
+
+    //
+    // Connect the transform to the Decorator.
+    auto* transformOutput = static_cast<TransformOutputType*>(this->ProcessObject::GetOutput(0));
+
+    //transformOutput->Set(m_Transform1.GetPointer());
+    transformOutput->Set(m_IsocIECTransform.GetPointer());
+
+    if (!m_Interpolator1) {
+        itkExceptionMacro(<< "Interpolator1 is not present");
+    }
+
+    if (!m_Interpolator2) {
+        itkExceptionMacro(<< "Interpolator2 is not present");
+    }
+
+    if (!m_Filter1) {
+        itkExceptionMacro(<< "Filter1 is not present");
+    }
+
+    if (!m_Filter1) {
+        itkExceptionMacro(<< "Filter2 is not present");
+    }
+
+    if (!m_TransformMetaInfo) {
+        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);
+
+    m_Metric->SetInterpolator1(m_Interpolator1);
+    m_Metric->SetInterpolator2(m_Interpolator2);
+
+    if (m_FixedImageRegionDefined1) {
+        m_Metric->SetFixedImageRegion1(m_FixedImageRegion1);
+        std::cout << "m_FixedImageRegionDefined1 is defined " << std::endl;
+    } else {
+        m_Metric->SetFixedImageRegion1(m_FixedImage1->GetBufferedRegion());
+    }
+
+    if (m_FixedImageRegionDefined2) {
+        m_Metric->SetFixedImageRegion2(m_FixedImageRegion2);
+        std::cout << "m_FixedImageRegionDefined2 is defined " << std::endl;
+    } else {
+        m_Metric->SetFixedImageRegion2(m_FixedImage2->GetBufferedRegion());
+    }
+
+    m_Metric->Initialize();
+
+    // Setup the optimizer
+    m_Optimizer->SetScales(m_Metric->GetWeightings());
+    m_Optimizer->SetCostFunction(m_Metric);
+
+    m_InitialOptimizerParameters =
+            m_Metric->GetParameters();
+
+    m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters);
+}
+
+//template <typename TFixedImage, typename TMovingImage>
+//void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
+//    itk::Object *caller, const itk::EventObject &event)
+//{
+//  itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(caller);
+
+//  if(!itk::ProgressEvent().CheckEvent(&event) || !filter)
+//    return;
+
+//  if (filter)
+//  {
+////    double p = m_CurrentProgressStart;
+////    // filter->GetProgress() in [0;1]
+////    p += m_CurrentProgressSpan * filter->GetProgress();
+////    // NOTE: filter is buggy - not in [0;1] if multi-threading is activated!
+////    if (p > (m_CurrentProgressStart + m_CurrentProgressSpan))
+////      p = m_CurrentProgressStart + m_CurrentProgressSpan;
+////    TaskProgressInfo(m_CurrentProgressDirection, p);
+
+////    if (m_CancelRequest)
+////      filter->SetAbortGenerateData(true); // may be handled by filter
+
+//      std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
+//  }
+//}
+
+/*
+ * Starts the Registration Process
+ */
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
+{
+
+    itk::ProgressReporter progress(
+                this,
+                1, 1000,100);
+
+    ParametersType empty(1);
+    empty.Fill(0.0);
+    try {
+        // initialize the interconnects between components
+        this->Initialize();
+        if (GetDebug()) {
+            m_Metric->Print(std::cout);
+            m_Optimizer->Print(std::cout);
+            this->Print(std::cout);
+        }
+    } catch (ExceptionObject& err) {
+        m_LastOptimizerParameters = empty;
+
+        // pass exception to caller
+        throw err;
+    }
+
+    this->StartOptimization();
+}
+
+/*
+ * Starts the Optimization process
+ */
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartOptimization()
+{
+    try {
+        // do the optimization
+        m_Optimizer->StartOptimization();
+
+    } catch (ExceptionObject& err) {
+        // An error has occurred in the optimization.
+        // Update the parameters
+        m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition();
+
+        // Pass exception to caller
+        throw err;
+    }
+
+    // get the results
+    m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition();
+
+    //Update transform1 and transform2 with the last transform parameters from the optimizer.
+    m_Metric->SetTransformParameters(m_LastOptimizerParameters);
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os, Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+    os << indent << "Metric: " << m_Metric.GetPointer() << std::endl;
+    os << indent << "Optimizer: " << m_Optimizer.GetPointer() << std::endl;
+    os << indent << "Transform1: " << m_Transform1.GetPointer() << std::endl;
+    os << indent << "Transform2: " << m_Transform2.GetPointer() << std::endl;
+    os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl;
+    os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl;
+    os << indent << "Fixed Image 1: " << m_FixedImage1.GetPointer() << std::endl;
+    os << indent << "Fixed Image 2: " << m_FixedImage2.GetPointer() << std::endl;
+    os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl;
+    os << indent << "Transform MetaInfo: " << m_TransformMetaInfo.GetPointer() << std::endl;
+    os << indent << "Fixed Image 1 Region Defined: " << m_FixedImageRegionDefined1 << std::endl;
+    os << indent << "Fixed Image 2 Region Defined: " << m_FixedImageRegionDefined2 << std::endl;
+    os << indent << "Fixed Image 1 Region: " << m_FixedImageRegion1 << std::endl;
+    os << indent << "Fixed Image 2 Region: " << m_FixedImageRegion2 << std::endl;
+    os << indent << "Initial Optimizer Parameters: " << m_InitialOptimizerParameters << std::endl;
+    os << indent << "Last Optimizer Parameters: " << m_LastOptimizerParameters << std::endl;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GenerateData()
+{
+
+    this->StartRegistration();
+}
+
+template <typename TFixedImage, typename TMovingImage>
+const typename TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TransformOutputType*
+TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GetOutput() const
+{
+    return static_cast<const TransformOutputType*>(this->ProcessObject::GetOutput(0));
+}
+
+template <typename TFixedImage, typename TMovingImage>
+DataObject::Pointer
+TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::MakeOutput(DataObjectPointerArraySizeType output)
+{
+    switch (output) {
+    case 0:
+        return static_cast<DataObject*>(TransformOutputType::New().GetPointer());
+        break;
+    default:
+        itkExceptionMacro("MakeOutput request for an output number larger than the expected number of outputs");
+        return nullptr;
+    }
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImage1(const FixedImageType* fixedImage1)
+{
+    itkDebugMacro("setting Fixed Image 1 to " << fixedImage1);
+
+    if (this->m_FixedImage1.GetPointer() != fixedImage1) {
+        this->m_FixedImage1 = fixedImage1;
+
+        // Process object is not const-correct so the const_cast is required here
+        this->ProcessObject::SetNthInput(0, const_cast<FixedImageType*>(fixedImage1));
+
+        this->Modified();
+    }
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImage2(const FixedImageType* fixedImage2)
+{
+    itkDebugMacro("setting Fixed Image 2 to " << fixedImage2);
+
+    if (this->m_FixedImage2.GetPointer() != fixedImage2) {
+        this->m_FixedImage2 = fixedImage2;
+
+        // Process object is not const-correct so the const_cast is required here
+        this->ProcessObject::SetNthInput(0, const_cast<FixedImageType*>(fixedImage2));
+
+        this->Modified();
+    }
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetMovingImage(const MovingImageType* movingImage)
+{
+    itkDebugMacro("setting Moving Image to " << movingImage);
+
+    if (this->m_MovingImage.GetPointer() != movingImage) {
+        this->m_MovingImage = movingImage;
+
+        // Process object is not const-correct so the const_cast is required here
+        this->ProcessObject::SetNthInput(1, const_cast<MovingImageType*>(movingImage));
+
+        this->Modified();
+    }
+}
+
+} // end namespace itk
+
+#endif

itkReg23DRT/autoreg/itkgTwoImageToOneImageMetric.h

@@ -0,0 +1,319 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkgTwoImageToOneImageMetric_h
+#define itkgTwoImageToOneImageMetric_h
+
+#include "DRTMetaInformation.h"
+#include "itkChangeInformationImageFilter.h"
+#include "itkDRTHelpers.h"
+#include "itkGradientRecursiveGaussianImageFilter.h"
+#include "itkImageBase.h"
+#include "itkInterpolateImageFunction.h"
+#include "itkOptimizer.h"
+#include "itkResampleImageFilter.h"
+#include "itkSingleValuedCostFunction.h"
+#include "itkSpatialObject.h"
+#include "itkTransform.h"
+#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
+
+#include "itkImageFileWriter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+#include <cstring>
+#include <iostream>
+#include <libgen.h>
+#include <linux/limits.h>
+#include <unistd.h>
+
+#include <itkImageProcessorHelpers.h>
+
+namespace itk {
+
+/** \class gTwoImageToOneImageMetric
+ * \brief Computes similarity between two fixed images and one fixed image.
+ *
+ * This Class is templated over the type of the two input images.
+ * It expects a Transform and two Interpolators to be plugged in.
+ * This particular class is the base class for a hierarchy of
+ * similarity metrics.
+ *
+ * This class computes a value that measures the similarity
+ * between two Fixed image and the transformed Moving images.
+ * The Interpolators are used to compute intensity values on
+ * non-grid positions resulting from mapping points through
+ * the Transform.
+ *
+ *
+ * \ingroup RegistrationMetrics
+ * \ingroup TwoProjectionRegistration
+ *
+ */
+
+template <typename TFixedImage, typename TMovingImage>
+class gTwoImageToOneImageMetric : public SingleValuedCostFunction {
+public:
+    ITK_DISALLOW_COPY_AND_ASSIGN(gTwoImageToOneImageMetric);
+
+    /** Standard class type alias. */
+    using Self = gTwoImageToOneImageMetric;
+    using Superclass = SingleValuedCostFunction;
+    using Pointer = SmartPointer<Self>;
+    using ConstPointer = SmartPointer<const Self>;
+
+    /** Type used for representing point components  */
+    using CoordinateRepresentationType = Superclass::ParametersValueType;
+
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(gTwoImageToOneImageMetric, SingleValuedCostFunction);
+
+    /**  Type of the moving Image. */
+    using MovingImageType = TMovingImage;
+    using MovingImagePixelType = typename TMovingImage::PixelType;
+    using MovingImageConstPointer = typename MovingImageType::ConstPointer;
+
+    /**  Type of the fixed Image. */
+    using FixedImageType = TFixedImage;
+    using FixedImageConstPointer = typename FixedImageType::ConstPointer;
+    using FixedImageRegionType = typename FixedImageType::RegionType;
+
+    /** Constants for the image dimensions */
+    static constexpr unsigned int MovingImageDimension = TMovingImage::ImageDimension;
+    static constexpr unsigned int FixedImageDimension = TFixedImage::ImageDimension;
+
+    /**  Type of the Transform Base class */
+    using TransformType = itk::Euler3DTransform<double>;
+
+    using TransformPointer = typename TransformType::Pointer;
+    using InputPointType = typename TransformType::InputPointType;
+    using OutputPointType = typename TransformType::OutputPointType;
+    using TransformParametersType = typename TransformType::ParametersType;
+    using TransformJacobianType = typename TransformType::JacobianType;
+
+    /**  Type of the Interpolator Base class */
+    using InterpolatorType = gSiddonJacobsRayCastInterpolateImageFunction<MovingImageType, CoordinateRepresentationType>;
+
+    /** Gaussian filter to compute the gradient of the Moving Image */
+    using RealType = typename NumericTraits<MovingImagePixelType>::RealType;
+    using GradientPixelType = CovariantVector<RealType, itkGetStaticConstMacro(MovingImageDimension)>;
+    using GradientImageType = Image<GradientPixelType, itkGetStaticConstMacro(MovingImageDimension)>;
+    using GradientImagePointer = SmartPointer<GradientImageType>;
+    using GradientImageFilterType = GradientRecursiveGaussianImageFilter<MovingImageType, GradientImageType>;
+    using GradientImageFilterPointer = typename GradientImageFilterType::Pointer;
+    using InterpolatorPointer = typename InterpolatorType::Pointer;
+
+    using ResampleImageFilterType = ResampleImageFilter<MovingImageType, MovingImageType>;
+    using ResampleImageFilterPointer = typename ResampleImageFilterType::Pointer;
+    using ChangeInformationFilterType = ChangeInformationImageFilter<MovingImageType>;
+    using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer;
+
+    /**  Type for the mask of the fixed image. Only pixels that are "inside"
+       this mask will be considered for the computation of the metric */
+    typedef SpatialObject<itkGetStaticConstMacro(FixedImageDimension)> FixedImageMaskType;
+    using FixedImageMaskPointer = typename FixedImageMaskType::Pointer;
+
+    /**  Type for the mask of the moving image. Only pixels that are "inside"
+       this mask will be considered for the computation of the metric */
+    typedef SpatialObject<itkGetStaticConstMacro(MovingImageDimension)> MovingImageMaskType;
+    using MovingImageMaskPointer = typename MovingImageMaskType::Pointer;
+
+    /**  Type of the measure. */
+    using MeasureType = Superclass::MeasureType;
+
+    /**  Type of the derivative. */
+    using DerivativeType = Superclass::DerivativeType;
+
+    /** Type of the Transformation parameters This is the same type used to
+   *  represent the search space of the optimization algorithm */
+    using ParametersType = Superclass::ParametersType;
+
+    /** Connect the Fixed Image.  */
+    itkSetConstObjectMacro(FixedImage1, FixedImageType);
+
+    /** Connect the Fixed Image.  */
+    itkSetConstObjectMacro(FixedImage2, FixedImageType);
+
+    /** Get the Fixed Image. */
+    itkGetConstObjectMacro(FixedImage1, FixedImageType);
+
+    /** Get the Fixed Image. */
+    itkGetConstObjectMacro(FixedImage2, FixedImageType);
+
+    /** Connect the Moving Image.  */
+    itkSetConstObjectMacro(MovingImage, MovingImageType);
+
+    /** Get the Moving Image. */
+    itkGetConstObjectMacro(MovingImage, MovingImageType);
+
+    /** Connect the Transform1. */
+    itkSetObjectMacro(Transform1, TransformType);
+
+    /** Get a pointer to the Transform1.  */
+    itkGetConstObjectMacro(Transform1, TransformType);
+
+    /** Connect the Transform2. */
+    itkSetObjectMacro(Transform2, TransformType);
+
+    /** Get a pointer to the Transform2.  */
+    itkGetConstObjectMacro(Transform2, TransformType);
+
+    itkSetObjectMacro(IsocTransf, TransformType);
+    itkGetObjectMacro(IsocTransf, TransformType);
+
+
+    /** Connect the Interpolator. */
+    itkSetObjectMacro(Interpolator1, InterpolatorType);
+
+    /** Get a pointer to the Interpolator.  */
+    itkGetConstObjectMacro(Interpolator1, InterpolatorType);
+
+    /** Connect the Interpolator. */
+    itkSetObjectMacro(Interpolator2, InterpolatorType);
+
+    /** Get a pointer to the Interpolator.  */
+    itkGetConstObjectMacro(Interpolator2, InterpolatorType);
+
+    /** Get the number of pixels considered in the computation. */
+    itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long);
+
+    /** Connect the DRTGeometryMetaInfo. */
+    itkSetObjectMacro(TransformMetaInfo, R23MetaInformation);
+
+    /** Get a pointer to the DRTGeometryMetaInfo.  */
+    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);
+
+    /** Set the region over which the metric will be computed */
+    itkSetMacro(FixedImageRegion2, FixedImageRegionType);
+
+    /** Get the region over which the metric will be computed */
+    itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType);
+
+    /** Get the region over which the metric will be computed */
+    itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType);
+
+    /** Set/Get the output filters. */
+    itkSetObjectMacro(Filter1, ChangeInformationFilterType);
+    itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
+    itkSetObjectMacro(Filter2, ChangeInformationFilterType);
+    itkGetConstObjectMacro(Filter2, ChangeInformationFilterType);
+
+    /** Set/Get the moving image mask. */
+    itkSetObjectMacro(MovingImageMask, MovingImageMaskType);
+    itkGetConstObjectMacro(MovingImageMask, MovingImageMaskType);
+
+    /** Set/Get the fixed image mask. */
+    itkSetObjectMacro(FixedImageMask1, FixedImageMaskType);
+    itkSetObjectMacro(FixedImageMask2, FixedImageMaskType);
+    itkGetConstObjectMacro(FixedImageMask1, FixedImageMaskType);
+    itkGetConstObjectMacro(FixedImageMask2, FixedImageMaskType);
+
+    /** Set/Get gradient computation. */
+    itkSetMacro(ComputeGradient, bool);
+    itkGetConstReferenceMacro(ComputeGradient, bool);
+    itkBooleanMacro(ComputeGradient);
+
+    itkSetMacro(MaxTranslation, double);
+    itkGetConstReferenceMacro(MaxTranslation, double);
+
+    /** Get Gradient Image. */
+    itkGetConstObjectMacro(GradientImage, GradientImageType);
+
+    /** Set the parameters defining the Transform. */
+    bool
+    SetTransformParameters(const ParametersType& parameters) const;
+
+    /** Return the number of parameters required by the Transform */
+
+    unsigned int GetNumberOfParameters() const;
+    itk::Optimizer::ScalesType GetWeightings() const;
+    ParametersType GetParameters() const;
+
+    /** Initialize the Metric by making sure that all the components
+   *  are present and plugged together correctly     */
+    virtual void
+    Initialize();
+
+protected:
+    gTwoImageToOneImageMetric();
+    ~gTwoImageToOneImageMetric() override = default;
+    void
+    PrintSelf(std::ostream& os, Indent indent) const override;
+
+    constexpr static unsigned int Dimension = 3;
+
+    mutable unsigned long m_NumberOfPixelsCounted;
+    mutable double m_bestMeassure;
+
+    FixedImageConstPointer m_FixedImage1;
+    FixedImageConstPointer m_FixedImage2;
+    MovingImageConstPointer m_MovingImage;
+
+    mutable TransformPointer m_Transform1;
+    mutable TransformPointer m_Transform2;
+    mutable TransformPointer m_IsocTransf;
+
+    InterpolatorPointer m_Interpolator1;
+    InterpolatorPointer m_Interpolator2;
+
+    bool m_ComputeGradient;
+    GradientImagePointer m_GradientImage;
+
+    using OutputPixelType = unsigned char;
+    using OutputImageType = itk::Image<OutputPixelType, Dimension>;
+
+    mutable FixedImageMaskPointer m_FixedImageMask1;
+    mutable FixedImageMaskPointer m_FixedImageMask2;
+    mutable MovingImageMaskPointer m_MovingImageMask;
+
+    R23MetaInformation::Pointer
+        m_TransformMetaInfo;
+
+    InternalTransformMetaInformation::Pointer
+    m_internalMeta1,
+    m_internalMeta2;
+
+    ChangeInformationFilterPointer
+        m_Filter1,
+        m_Filter2;
+
+    double m_MaxTranslation;
+
+    /* Writes the image to the disk. If path is empty it will be saved to /img/out/*/
+    void WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const;
+
+private:
+    FixedImageRegionType m_FixedImageRegion1;
+    FixedImageRegionType m_FixedImageRegion2;
+};
+
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkgTwoImageToOneImageMetric.hxx"
+#endif
+
+#endif

itkReg23DRT/autoreg/itkgTwoImageToOneImageMetric.hxx

@@ -0,0 +1,476 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkgTwoImageToOneImageMetric_hxx
+#define itkgTwoImageToOneImageMetric_hxx
+
+#include "itkgTwoImageToOneImageMetric.h"
+#include <limits>
+
+namespace itk {
+
+template <typename TFixedImage, typename TMovingImage>
+gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::gTwoImageToOneImageMetric()
+{
+    m_FixedImage1 = nullptr; // has to be provided by the user.
+    m_FixedImage2 = nullptr; // has to be provided by the user.
+    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_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;
+    m_ComputeGradient = true; // metric computes gradient by default
+    m_NumberOfPixelsCounted = 0; // initialize to zero
+    m_bestMeassure = std::numeric_limits<short>::max();
+    m_MaxTranslation = std::numeric_limits<short>::max();
+}
+
+/*
+ * Set the parameters that define a unique transform
+ */
+template <typename TFixedImage, typename TMovingImage>
+bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameters(const ParametersType& parameters) const
+{
+
+    if (!m_Transform1) {
+        itkExceptionMacro(<< "Transform 1  has not been assigned");
+    }
+    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_IsocTransf->GetParameters();
+    double TranslationAlongX = transformParameters[3];
+    double TranslationAlongY = transformParameters[4];
+    double TranslationAlongZ = transformParameters[5];
+    double RotationAlongX = transformParameters[0];
+    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];
+        break;
+    case Y_ONLY:
+        TranslationAlongY = parameters[0];
+        break;
+    case Z_ONLY:
+        TranslationAlongZ = parameters[0];
+        break;
+    case THREE_DOF:
+        TranslationAlongX = parameters[0];
+        TranslationAlongY = parameters[1];
+        TranslationAlongZ = parameters[2];
+        break;
+    case ROTATION_ONLY:
+        RotationAlongX = parameters[0];
+        RotationAlongY = parameters[1];
+        RotationAlongZ = parameters[2];
+        break;
+    case SIX_DOF:
+        TranslationAlongX = parameters[0];
+        TranslationAlongY = parameters[1];
+        TranslationAlongZ = parameters[2];
+        RotationAlongX = parameters[3];
+        RotationAlongY = parameters[4];
+        RotationAlongZ = parameters[5];
+        break;
+
+    default:
+        itkExceptionMacro(<< "Unsupported degree of freedeom");
+        break;
+    }
+
+//    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 Y = " << RotationAlongY /*/ dtr*/ << " deg" << std::endl;
+//    std::cout << " Rotation Along Z = " << RotationAlongZ /*/ dtr*/ << " deg" << std::endl;
+//    std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
+
+
+    // 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) {
+        /*
+            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;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
+{
+
+    if (!m_IsocTransf) {
+        itkExceptionMacro(<< "Transform is not present");
+    }
+
+    if (!m_Transform1) {
+        itkExceptionMacro(<< "Transform1 is not present");
+    }
+
+    if (!m_Transform2) {
+        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) {
+        itkExceptionMacro(<< "Interpolator1 is not present");
+    }
+    if (!m_Interpolator2) {
+        itkExceptionMacro(<< "Interpolator2 is not present");
+    }
+
+    if (!m_FixedImage1) {
+        itkExceptionMacro(<< "FixedImage1 is not present");
+    }
+
+    if (!m_FixedImage2) {
+        itkExceptionMacro(<< "FixedImage2 is not present");
+    }
+
+    if (m_FixedImageRegion1.GetNumberOfPixels() == 0) {
+        itkExceptionMacro(<< "FixedImageRegion1 is empty");
+    }
+
+    if (m_FixedImageRegion2.GetNumberOfPixels() == 0) {
+        itkExceptionMacro(<< "FixedImageRegion2 is empty");
+    }
+
+    // If the image is provided by a source, update the source.
+    if (m_MovingImage->GetSource()) {
+        m_MovingImage->GetSource()->Update();
+    }
+
+    // If the image is provided by a source, update the source.
+    if (m_FixedImage1->GetSource()) {
+        m_FixedImage1->GetSource()->Update();
+    }
+
+    if (m_FixedImage2->GetSource()) {
+        m_FixedImage2->GetSource()->Update();
+    }
+
+    //std::cout<<"FixedImageRegion1 "  <<this->GetFixedImageRegion1() <<std::endl;
+    //std::cout<<"m_FixedImage1 buffered " <<this->m_FixedImage1->GetBufferedRegion() <<std::endl;
+    //std::cout<<"Moving buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
+    //std::cout<<"GetFixedImage1 buffered " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
+    //std::cout<<"Filter1 buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
+
+    // Make sure the FixedImageRegion is within the FixedImage buffered region
+    if (!m_FixedImageRegion1.Crop(m_FixedImage1->GetBufferedRegion())) {
+        itkExceptionMacro(<< "FixedImageRegion1 does not overlap the fixed image buffered region");
+    }
+
+    if (!m_FixedImageRegion2.Crop(m_FixedImage2->GetBufferedRegion())) {
+        itkExceptionMacro(<< "FixedImageRegion2 does not overlap the fixed image buffered region");
+    }
+
+    this->SetTransformParameters(this->GetParameters());
+
+    if (m_ComputeGradient) {
+
+        GradientImageFilterPointer gradientFilter = GradientImageFilterType::New();
+
+        gradientFilter->SetInput(m_MovingImage);
+
+        const typename MovingImageType::SpacingType& spacing = m_MovingImage->GetSpacing();
+        double maximumSpacing = 0.0;
+        for (unsigned int i = 0; i < MovingImageDimension; i++) {
+            if (spacing[i] > maximumSpacing) {
+                maximumSpacing = spacing[i];
+            }
+        }
+        gradientFilter->SetSigma(maximumSpacing);
+        gradientFilter->SetNormalizeAcrossScale(true);
+
+        gradientFilter->Update();
+
+        m_GradientImage = gradientFilter->GetOutput();
+    }
+
+    // If there are any observers on the metric, call them to give the
+    // user code a chance to set parameters on the metric
+    this->InvokeEvent(InitializeEvent());
+}
+
+template <typename TFixedImage, typename TMovingImage>
+unsigned int gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetNumberOfParameters() const
+{
+    if (!m_TransformMetaInfo) {
+        itkExceptionMacro(<< "TransformMetaInfo is missing");
+    }
+
+    switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
+    case X_ONLY:
+    case Y_ONLY:
+    case Z_ONLY:
+        return 1;
+    case THREE_DOF:
+    case ROTATION_ONLY:
+        return 3;
+    case SIX_DOF:
+        return 6;
+    default:
+        itkExceptionMacro(<< "Unsupported degree of freedeom");
+        break;
+    }
+
+    return m_IsocTransf->GetNumberOfParameters();
+}
+
+template <typename TFixedImage, typename TMovingImage>
+itk::Optimizer::ScalesType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetWeightings() const
+{
+    itk::Optimizer::ScalesType weightings(this->GetNumberOfParameters());
+
+    switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
+    case X_ONLY:
+    case Y_ONLY:
+    case Z_ONLY:
+        weightings[0] = 1.0;
+        break;
+        ;
+    case THREE_DOF:
+        weightings[0] = 1.;
+        weightings[1] = 1.;
+        weightings[2] = 1.;
+        break;
+    case ROTATION_ONLY:
+        weightings[0] = 1. / dtr;
+        weightings[1] = 1. / dtr;
+        weightings[2] = 1. / dtr;
+        break;
+    case SIX_DOF:
+        //G: ITHINKTHOSEARE FLIPPED!!
+        weightings[0] = 1.;
+        weightings[1] = 1.;
+        weightings[2] = 1.;
+        weightings[3] = 1. / dtr;
+        weightings[4] = 1. / dtr;
+        weightings[5] = 1. / dtr;
+        break;
+
+    default:
+        itkExceptionMacro(<< "Unsupported degree of freedeom");
+        break;
+    }
+    return weightings;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+typename gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::
+ParametersType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetParameters() const
+{
+
+    //    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:
+        parameters[0] = parametersTransform[3];
+        break;
+    case Y_ONLY:
+        parameters[0] = parametersTransform[4];
+        break;
+    case Z_ONLY:
+        parameters[0] = parametersTransform[5];
+        break;
+    case THREE_DOF:
+        parameters[0] = parametersTransform[3];
+        parameters[1] = parametersTransform[4];
+        parameters[2] = parametersTransform[5];
+        break;
+    case ROTATION_ONLY:
+        parameters[0] = parametersTransform[0];
+        parameters[1] = parametersTransform[1];
+        parameters[2] = parametersTransform[2];
+        break;
+    case SIX_DOF:
+        parameters[0] = parametersTransform[3];
+        parameters[1] = parametersTransform[4];
+        parameters[2] = parametersTransform[5];
+        parameters[3] = parametersTransform[0];
+        parameters[4] = parametersTransform[1];
+        parameters[5] = parametersTransform[2];
+        break;
+
+    default:
+        itkExceptionMacro(<< "Unsupported degree of freedeom");
+        break;
+    }
+    return parameters;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os, Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+    os << indent << "ComputeGradient: " << static_cast<typename NumericTraits<bool>::PrintType>(m_ComputeGradient)
+       << std::endl;
+    os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl;
+    os << indent << "Fixed  Image 1: " << m_FixedImage1.GetPointer() << std::endl;
+    os << indent << "Fixed  Image 2: " << m_FixedImage2.GetPointer() << std::endl;
+    os << indent << "Gradient Image: " << m_GradientImage.GetPointer() << std::endl;
+    os << indent << "Transform1:    " << m_Transform1.GetPointer() << std::endl;
+    os << indent << "Transform2:    " << m_Transform2.GetPointer() << std::endl;
+    os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl;
+    os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl;
+    os << indent << "Filter 1: " << m_Filter1.GetPointer() << std::endl;
+    os << indent << "Filter 2: " << m_Filter2.GetPointer() << std::endl;
+    os << indent << "Transform MetaInfoe: " << m_TransformMetaInfo.GetPointer() << std::endl;
+    os << indent << "FixedImageRegion 1: " << m_FixedImageRegion1 << std::endl;
+    os << indent << "FixedImageRegion 2: " << m_FixedImageRegion2 << std::endl;
+    os << indent << "Moving Image Mask: " << m_MovingImageMask.GetPointer() << std::endl;
+    os << indent << "Fixed Image Mask 1: " << m_FixedImageMask1.GetPointer() << std::endl;
+    os << indent << "Fixed Image Mask 2: " << m_FixedImageMask2.GetPointer() << std::endl;
+    os << indent << "Number of Pixels Counted: " << m_NumberOfPixelsCounted << std::endl;
+}
+
+template <typename TFixedImage, typename TMovingImage>
+void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const
+{
+
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<MovingImageType, OutputImageType>;
+    using WriterType = itk::ImageFileWriter<OutputImageType>;
+
+    if (img) {
+
+        // Rescale the intensity of the projection images to 0-255 for output.
+        typename RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+        rescaler1->SetOutputMinimum(0);
+        rescaler1->SetOutputMaximum(255);
+        rescaler1->SetInput(img);
+
+        WriterType::Pointer writer1 = WriterType::New();
+
+        if (fileName.empty()) {
+
+            char result[PATH_MAX];
+            ssize_t count = readlink("/proc/self/exe", result, PATH_MAX);
+            if (count != -1) {
+                path = dirname(result);
+                path += "/img/out/";
+            }
+        }
+
+        if (fileName.empty()) {
+
+            fileName = path;
+            fileName += "2D1.tif";
+        }
+        writer1->SetFileName(fileName);
+        writer1->SetInput(rescaler1->GetOutput());
+
+        try {
+            std::cout << "Writing image: " << fileName << std::endl;
+            writer1->Update();
+        } catch (itk::ExceptionObject& err) {
+            std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+            std::cerr << err << std::endl;
+        }
+    }
+}
+
+} // end namespace itk
+
+#endif

itkReg23DRT/itkImageProcessor.h

@@ -0,0 +1,471 @@
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+#ifndef ITKIMAGEPROCESSOR_H
+#define ITKIMAGEPROCESSOR_H
+
+
+#include "DRTMetaInformation.h"
+
+#include "itkCommand.h"
+
+#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
+#include "itkEuler3DTransform.h"
+#include "itkResampleImageFilter.h"
+#include "itkCastImageFilter.h"
+#include "itkRescaleIntensityImageFilter.h"
+#include "itkIntensityWindowingImageFilter.h"
+#include "itkChangeInformationImageFilter.h"
+#include "itkGDCMImageIO.h"
+#include "itkMetaDataObject.h"
+#include "itkImageDuplicator.h"
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "itkSmartPointer.h"
+
+#include "itkImageToVTKImageFilter.h"
+
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkImageSeriesReader.h"
+
+#include "gdcmGlobal.h"
+
+#include "vtkImageData.h"
+#include "vtkTransform.h"
+#include "vtkSmartPointer.h"
+
+#include "itkQtIterationUpdate.h"
+
+#include "itkImageProcessorHelpers.h"
+
+#include "itkReg23.h"
+#include "itkReg23MetaInformation.h"
+
+#define IMG_VIS_MAXIMUM_RANGE 2048
+
+namespace itk
+{
+
+class ITK_EXPORT itkImageProcessor : public itk::Object
+{
+
+    constexpr static unsigned int Dimension = 3;
+
+public:
+    /** Standard "Self" typedef. */
+    typedef itkImageProcessor  Self;
+    /** Standard "Superclass" typedef. */
+    typedef Object Superclass;
+    /** Smart pointer typedef support */
+    typedef itk::SmartPointer<Self>  Pointer;
+    typedef itk::SmartPointer<const Self>  ConstPointer;
+    /** Method of creation through the object factory. */
+    itkNewMacro(Self);
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(itkImageProcessor, Object);
+
+
+    itkSetMacro(ResamplepCT, bool);
+    itkGetMacro(ResamplepCT, bool);
+
+    itkSetMacro(SamplingLNG, double);
+    itkGetMacro(SamplingLNG, double);
+
+    CommandIterationUpdate::Pointer
+    optimizerObserver;
+
+    ExhaustiveCommandIterationUpdate::Pointer
+    ExhaustiveOptimizerObserver;
+
+    /** Sets the degree of freedom for optimzation*/
+    tDegreeOfFreedomEnum GetDegreeOfFreedom();
+    /** Returns user components only of the autoReg
+    * */
+    Optimizer::ParametersType
+    GetFinalR23Parameters();
+    /** Define the region to be used as ROI on fixed images */
+    void SetRegionFixed1(int,int,int,int);
+    void SetRegionFixed2(int,int,int,int);
+    void ResetROIRegions();
+    void InizializeRegistration();
+
+    /** Maximum number of iterations for the optimizer */
+    itkGetMacro(R23, itkReg23::Pointer);
+
+
+    /** Set Optimizer type */
+    void SetOptimizer(std::string);
+    void SetOptimizer(tOptimizerTypeEnum);
+
+    /** Set Metric type */
+    void SetMetric(std::string);
+    void SetMetric(tMetricTypeEnum);
+
+    /** Set DOF */
+    void SetDegreeOfFreedom(tDegreeOfFreedomEnum);
+
+    /** Set Handle Of Rotation Import Offset */
+    void SetHandleRotationImportOffset(tHandlingRotImpTransform);
+
+    /** Set PowellOptimizer */
+    void SetPowellOptimParameters(double,double,double,int,int);
+
+    /** Set AmoebaOptimizer */
+    void SetAmoebaOptimParameters(double,double,double,int);
+
+    /** Set MI */
+    void SetMIMetricParameters(double,int);
+
+    /** Set NCC */
+    void SetNCCMetricParameters(double,bool);
+
+    /** Set number of logic CPU to be made available to interpolators*/
+    void SetNumberOfWorkingUnits(int iN);
+
+    /** Input data load methods*/
+    int load3DSerieFromFiles( std::vector<std::string> );
+    int unload3DVolumeAndMeta();
+
+    int load2D(const char *);
+    int unload2DAndMeta(int);
+    //int query2DimageReconstructionDiameter(const char*);
+
+    void loadRTPlanInfo(double, double, double, double, double ,double);
+    int unloadRTPlanAndMeta();
+
+    /** Projection angles - Gantry angle IEC */
+    void SetProjectionAngle1IEC(double);
+    void SetProjectionAngle2IEC(double);
+    void SetProjectionAngleOffsetIEC(double,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);
+    void SetSCDOffset(double,double);
+
+    double GetSCD1();
+    double GetSCD2();
+
+    /** Panel offsets - panel offsets IEC */
+    void SetPanelOffsets(double,double);
+
+    /** Panel offsets - takes into account patient orientation */
+    double GetPanelOffsetPGeo(tProjOrientationType ImgPrj);
+
+    /** Get projection angles - Gantry angle LPS
+     *  Only meaningful after a 3D Volume is loaded */
+    double GetPanelOffset1();
+    double GetPanelOffset2();
+
+    void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double);
+    void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double);
+
+    /** Intensity threshold used for image projection,
+     *  any voxel below threshold is ignored */
+    void SetIntensityThreshold(double);
+
+    /** 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();
+
+
+    /** Set transform parameters for 3D Volume */
+    void SetInitialTranslations(double, double, double);
+    void SetInitialRotations(double, double, double);
+
+    ///** Get transform parameters for 3D Volume */
+    /** 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();
+
+    /** Get Projection origin in LPS coordinates*/
+    const std::vector <double> GetLPStoProjectionGeoLPSOffset();
+
+    /** Get Projection origin in LPS coordinates*/
+    const std::vector <double> GetRTImportOffset();
+
+    /** Get the Localizer intensity window and
+     *  center for rendering */
+    double GetLocalizerDisplayWindowLevel(int );
+    double GetLocalizerDisplayWindowWidth(int );
+
+    /** Compute Digital Localizers */
+    void GetProjectionImages();
+
+    /** Conveniency methods to get VTK images for rendering */
+    vtkImageData* GetProjection1VTK();
+    vtkImageData* GetProjection2VTK();
+
+    /** Conveniency methods to get VTK images for rendering */
+    vtkImageData* GetProjection1VTKToWrite();
+    vtkImageData* GetProjection2VTKToWrite();
+
+    vtkImageData* GetLocalizer1VTK();
+    vtkImageData* GetLocalizer2VTK();
+
+    vtkMatrix4x4 * GetProjection1VTKTransform();
+    vtkMatrix4x4 * GetProjection2VTKTransform();
+
+    /** Debug writers */
+    void WriteProjectionImages();
+    void Write2DImages();
+
+
+protected:
+    /** Various pixel types */
+    using InternalPixelType = float;
+    using PixelType2D = short;
+    using PixelType3D = short;
+    using OutputPixelType = unsigned short;
+
+    using ImageType3D = itk::Image<PixelType3D, Dimension>;
+    using InternalImageType = itk::Image<InternalPixelType, Dimension>;
+
+    itkImageProcessor();
+    virtual ~itkImageProcessor();
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+    bool
+    m_ResamplepCT;
+    double
+    m_SamplingLNG;
+
+private:
+    itkImageProcessor(const Self&); //purposely not implemented
+    void operator=(const Self&); //purposely not implemented
+
+    /** Fill Meta after 3D volume load */
+    int fill3DVolumeMeta(InternalImageType::Pointer,
+                         tPatOrientation);
+
+    /** Image types */
+    using ImageType2D = itk::Image<PixelType3D, Dimension>;
+
+    using OutputImageType = itk::Image<OutputPixelType, Dimension>;
+
+    /** The following lines define each of the components used in the
+      registration: The transform, optimizer, metric, interpolator and
+      the registration method itself. */
+    using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
+    using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
+
+
+    using RoiForRegistration = itk::ImageRegion<Dimension>;
+
+    /** Image reader types */
+    using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
+    using ImageReaderType3D = itk::ImageFileReader<ImageType3D>;
+    using ImageSeriesReaderType3D = itk::ImageSeriesReader<ImageType3D>;
+    using ImageIOType = itk::GDCMImageIO;
+    using DuplicatorType = itk::ImageDuplicator<InternalImageType>;
+
+    /** Image input resampling Types */
+    using ResampleInputFilterType = itk::ResampleImageFilter<ImageType3D, ImageType3D>;
+    using LinearInterpolatorType = itk::LinearInterpolateImageFunction<ImageType3D, double>;
+    using SizeValueType = ImageType3D::SizeType::SizeValueType;
+
+    /** widely used filters: flip, cast, rescale, ... */
+    using Input2DRescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
+    using CastFilterType3D = itk::CastImageFilter<ImageType3D, InternalImageType>;
+    using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
+    using ChangeInformationFilterInput2DType = itk::ChangeInformationImageFilter<ImageType2D>;
+
+    /** ITK to VTK filter */
+    using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
+
+    TransformType::Pointer
+    transform1,
+    transform2,
+    IsocTransf;
+
+    InterpolatorType::Pointer
+    interpolator1,
+    interpolator2;
+
+    ITKtoVTKFilterType::Pointer
+    toVTK2D1,
+    toVTK2D2,
+    toVTKLocalizer1,
+    toVTKLocalizer2;
+
+    InternalImageType::Pointer
+    imageDRT1In,
+    imageDRT2In;
+
+
+
+    DuplicatorType::Pointer
+    m_LATSourceDupli,
+    m_PASourceDupli,
+    m_VolumeSourceDupli;
+
+    /** 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,
+     * Only input variables are used... */
+
+    InternalImageType::DirectionType
+    CalcDRTImageDirections(double angle,
+                           InternalImageType::DirectionType DRT2LPS);
+
+    ImageType3D::PointType
+    CalcDRTImageOrigin(
+            ImageType3D::PointType m_DRTOrigin,
+            ImageType3D::PointType m_DRTCOV,
+            ImageType3D::PointType m_ProjectionCenterLPS,
+            double dAngle,
+            InternalImageType::DirectionType stdDRT2LPS
+            );
+
+
+    double
+    CalcProjectionAngleLPS(
+            tPatOrientation pOrient,
+            double pAngleIEC);
+
+    double
+        CalcPanelOffsetPGeo(
+                tPatOrientation pOrient,
+                double dPOffset);
+    /** Apply transform to CT volume.
+     *  This moves the volume based on RT Plan info
+     *  if those are available. */
+    ImageType3D::PointType
+    CalcImportVolumeOffset(
+            ImageType3D::PointType rtCouchOffset,
+            InternalImageType::DirectionType VolumeLPStoIEC,
+            ImageType3D::PointType rtIsocenterLPS,
+            ImageType3D::PointType IEC2DCMMapT,
+            ImageType3D::PointType IEC2DCMMapR);
+
+
+    // The ResampleImageFilter is the driving force for the projection image generation.
+    ResampleFilterType::Pointer
+    resampleFilter1,
+    resampleFilter2;
+
+    ChangeInformationFilterType::Pointer
+    filter1,
+    filter2;
+
+    InternalImageType::DirectionType
+    Std_DRT2LPS;
+
+    vtkMatrix4x4
+    * m_Projection1VTKTransform,
+    * m_Projection2VTKTransform;
+
+    /**
+     * Many meta containers
+     */
+    CTVolumeImageMetaInformation::Pointer
+    m_CTMetaInfo;
+
+    DRTProjectionGeometryImageMetaInformation::Pointer
+    m_DRTGeometryMetaInfo;
+
+    DRTImageMetaInformation::Pointer
+    m_DRTImage1MetaInfo,
+    m_DRTImage2MetaInfo;
+
+    TopogramImageMetaInformation::Pointer
+    m_TImage1MetaInfo,
+    m_TImage2MetaInfo;
+
+    RTGeometryMetaInformation::Pointer
+    m_RTMetaInfo;
+
+    TransformMetaInformation::Pointer
+    m_TransformMetaInfo;
+
+    R23MetaInformation::Pointer
+    m_r23MetaInfo;
+
+    PowellOptimizerMetaInformation::Pointer
+    m_PowellOptimizerMetaInfo;
+
+    AmoebaOptimizerMetaInformation::Pointer
+    m_AmoebaOptimizerMetaInfo;
+
+    MIMetricMetaInformation::Pointer
+    m_MIMetricMetaInfo;
+
+    NCCMetricMetaInformation::Pointer
+    m_NCCMetricMetaInfo;
+
+    InternalTransformMetaInformation::Pointer
+    m_InternalTransf1,
+    m_InternalTransf2;
+
+
+    RoiForRegistration
+    roiAutoReg1,
+    roiAutoReg2;
+
+
+    int
+    iNWUnits;
+
+
+    itk::itkReg23::Pointer
+    m_R23;
+
+};
+
+
+}
+
+
+#endif

itkReg23DRT/itkImageProcessorHelpers.cpp

@@ -0,0 +1,166 @@
+#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;
+    buffer.clear();
+
+    /* 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.begin(), s.end(), buffer.begin());
+
+    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();
+
+    std::string buffer;
+
+    std::string s(gGetStringValueFromTag(gdcm::Tag(0x0018,0x1100), ds));
+    std::copy(s.begin(), s.end(), buffer.begin());
+
+    return std::stoi(buffer);
+}
+
+
+//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();
+}
+
+}

itkReg23DRT/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

itkReg23DRT/itkQtIterationUpdate.h

@@ -0,0 +1,224 @@
+#ifndef ITKQTITERATIONUPDATE_H
+#define ITKQTITERATIONUPDATE_H
+
+#include "itkCommand.h"
+#include <QObject>
+#include "itkPowellOptimizer.h"
+#include "itkAmoebaOptimizer.h"
+#include "itkExhaustiveOptimizer.h"
+
+#include "itkTwoProjectionImageRegistrationMethod.h"
+
+
+class QObjectIterationUpdate : public QObject{
+    Q_OBJECT
+public:
+    QObjectIterationUpdate(){
+        bAbortProcessCommand = false;
+    };
+
+    bool getAbortFlag(){
+        return bAbortProcessCommand;
+    };
+    void setAbortFlag(bool bVal){
+        bAbortProcessCommand = bVal;
+    };
+
+public slots:
+    void onAbortProcessCommand(){
+        bAbortProcessCommand = true;
+    };
+    void onIteration(itk::tOptimizerTypeEnum eType,double dI,double dX,double dY,double dZ){
+
+        switch (eType) {
+        case itk::eOptimizerType::POWELL:
+
+            emit
+            sendRegistrationProgress(0,dI,dX,dY,dZ);
+
+            break;
+
+        case itk::eOptimizerType::AMOEBA:
+
+            emit
+            sendRegistrationProgress(1,dI,dX,dY,dZ);
+
+             break;
+        default:
+            break;
+        }
+
+    }
+private:
+    bool
+    bAbortProcessCommand;
+
+
+
+
+
+signals:
+void sendRegistrationProgress(int,double,double,double,double);
+
+};
+
+namespace itk
+{
+
+
+class CommandIterationUpdate : public itk::Command {
+    // TODO: Move to own files.
+
+    constexpr static unsigned int Dimension = 3;
+
+public:
+    QObjectIterationUpdate *objIterUpdate;
+    tOptimizerTypeEnum
+    eOptimType;
+
+    int
+    iAmoebaIterations;
+
+    using Self = CommandIterationUpdate;
+    using Superclass = itk::Command;
+    using Pointer = itk::SmartPointer<Self>;
+    itkNewMacro(CommandIterationUpdate);
+
+    using InternalPixelType = float;
+    using InternalImageType = itk::Image<InternalPixelType, Dimension>;
+    using ProcesssType = typename itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
+    using ProcesssPointer = typename ProcesssType::Pointer;
+
+    /** Set/Get the Process. */
+    itkSetObjectMacro(Process, ProcesssType);
+    itkGetConstObjectMacro(Process, ProcesssType);
+
+private:
+
+protected:
+    CommandIterationUpdate() {
+        objIterUpdate = new QObjectIterationUpdate;
+    }
+    ProcesssPointer m_Process;
+
+public:
+    using OptimizerType = itk::PowellOptimizer;
+    using OptimizerPointer = const OptimizerType*;
+
+    using AmoebaOptimizerType = itk::AmoebaOptimizer;
+    using AmoebaOptimizerPointer = const AmoebaOptimizerType*;
+
+    void setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum eOptType){
+        eOptimType = eOptType;
+        iAmoebaIterations = 0;
+    }
+
+    void
+    Execute(itk::Object* caller, const itk::EventObject& event) override
+    {
+        Execute((const itk::Object*)caller, event);
+    }
+
+    void
+    Execute(const itk::Object* object, const itk::EventObject& event) override
+    {
+        if(objIterUpdate->getAbortFlag()){
+            std::cout << "AbortRequest" << std::endl;
+            objIterUpdate->setAbortFlag(false);
+            throw itk::ProcessAborted();
+        }
+//        std::cout << "Progress: " << this->m_Process->GetAbortGenerateData() << std::endl;
+
+        OptimizerPointer optPow;
+        AmoebaOptimizerPointer optAm;
+
+        switch (eOptimType) {
+        case eOptimizerType::POWELL:
+            optPow = dynamic_cast<OptimizerPointer>(object);
+            if (typeid(event) == typeid(itk::IterationEvent)) {
+
+                objIterUpdate->onIteration(
+                            eOptimType,
+                        optPow->GetCurrentIteration()+1,
+                        optPow->GetCurrentPosition()[0],
+                        optPow->GetCurrentPosition()[1],
+                        optPow->GetCurrentPosition()[2]
+                        );
+            }
+
+            break;
+        case eOptimizerType::AMOEBA:
+            optAm = dynamic_cast<AmoebaOptimizerPointer>(object);
+            iAmoebaIterations ++;
+
+            objIterUpdate->onIteration(
+                        eOptimType,
+                    iAmoebaIterations, //optAm->GetCachedValue(),
+                    optAm->GetCachedCurrentPosition()[0],
+                    optAm->GetCachedCurrentPosition()[1],
+                    optAm->GetCachedCurrentPosition()[2]
+                    );
+        break;
+        default:
+            break;
+        }
+
+        return;
+    }
+};
+
+class ExhaustiveCommandIterationUpdate : public itk::Command {
+    // TODO: Move to own files.
+
+public:
+    using Self = ExhaustiveCommandIterationUpdate;
+    using Superclass = itk::Command;
+    using Pointer = itk::SmartPointer<Self>;
+    itkNewMacro(Self);
+
+    std::ostream* os;
+
+protected:
+    ExhaustiveCommandIterationUpdate() = default;
+
+public:
+    using OptimizerType = itk::ExhaustiveOptimizer;
+    ;
+    using OptimizerPointer = const OptimizerType*;
+
+    void set_stream(std::ostream& stream)
+    {
+        os = &stream;
+    }
+
+    void
+    Execute(itk::Object* caller, const itk::EventObject& event) override
+    {
+        Execute((const itk::Object*)caller, event);
+    }
+
+    void
+    Execute(const itk::Object* object, const itk::EventObject& event) override
+    {
+        auto optimizer = dynamic_cast<OptimizerPointer>(object);
+        if (typeid(event) == typeid(itk::IterationEvent)) {
+
+            //crude LPS to IEC transform for HFS.
+            auto position = optimizer->GetCurrentPosition();
+
+            auto oldprecision = os->precision();
+            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[1] << "\t" << position[2] << std::endl;
+
+        }
+        return;
+    }
+};
+
+}
+
+#endif
+

itkReg23DRT/itkReg23.cpp

@@ -0,0 +1,439 @@
+#include "itkReg23.h"
+
+
+
+#include "itkTimeProbesCollectorBase.h"
+
+namespace itk
+{
+
+itkReg23::itkReg23()
+{
+    // Metrics
+    NCCmetric = MetricType::New();
+    MImetric = MIMetricType::New();
+
+
+    // Optimizers
+    PowellOptimizer = PowellOptimizerType::New();
+    AmoebaOptimizer = AmoebaOptimizerType::New();
+    ExhaustiveOptimizer = ExhaustiveOptimizerType::New();
+
+    // Registration
+    registration = RegistrationType::New();
+
+    //internalTransforms
+    Transform1 = TransformType::New();
+    Transform2 = TransformType::New();
+    IsocTransform = TransformType::New();
+
+    registration->SetTransform1(Transform1);
+    registration->SetTransform2(Transform2);
+
+
+
+    // to be provided by the user
+    m_r23Meta = nullptr;
+    m_Volume = nullptr;
+    m_PA = nullptr;
+    m_LAT = nullptr;
+    m_InternalTransf1 = nullptr;
+    m_InternalTransf2 = nullptr;
+    m_filter1 = nullptr;
+    m_filter2 = nullptr;
+    m_TransformMetaInfo = nullptr;
+    m_OptimizerObserver = nullptr;
+
+    //    m_UseFullROI = true; // if the full image ROI shall be used
+    m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file.
+
+
+    this->setROISizeToZero();
+}
+
+
+itkReg23::~itkReg23()
+{
+
+}
+
+void itkReg23::setROISizeToZero(){
+
+    ImageType3D::SizeType zeroSize;
+    zeroSize.Fill(0);
+
+    m_roiAutoReg1.SetSize(zeroSize);
+    m_roiAutoReg2.SetSize(zeroSize);
+
+}
+
+void itkReg23::PrintSelf(std::ostream& os, Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+void itkReg23::InitializeRegistration()
+{
+    std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
+
+    if (!m_PA) {
+        itkExceptionMacro(<< "PA topogram not present");
+    }
+
+    if (!m_LAT) {
+        itkExceptionMacro(<< "LAT topogram not present");
+    }
+
+    if (!m_Volume) {
+        itkExceptionMacro(<< "CT data not present");
+    }
+
+    if (!m_r23Meta) {
+        itkExceptionMacro(<< "r23Meta data not present");
+    }
+
+    if (!m_OptimizerObserver) {
+        itkExceptionMacro(<< "m_OptimizerObserver data not present");
+    }
+
+    if (!m_InternalTransf1) {
+        itkExceptionMacro(<< "m_InternalTransf1 data not present");
+    }
+
+    if (!m_InternalTransf2) {
+        itkExceptionMacro(<< "m_InternalTransf2 data not present");
+    }
+
+    if (!m_filter1) {
+        itkExceptionMacro(<< "m_filter1 data not present");
+    }
+
+    if (!m_filter2) {
+        itkExceptionMacro(<< "m_filter2 data not present");
+    }
+
+    if (!m_TransformMetaInfo) {
+        itkExceptionMacro(<< "m_TransformMetaInfo data not present");
+    }
+
+    if(!m_interpolator1){
+        itkExceptionMacro(<< "m_interpolator1 data not present");
+    }
+
+    if(!m_interpolator2){
+        itkExceptionMacro(<< "m_interpolator2 data not present");
+    }
+
+//    std::cout << "Print r23" <<
+//                 m_r23Meta->GetMetricType() << " " << m_r23Meta->GetDegreeOfFreedom() <<  std::endl;
+//    std::cout << "Print Powell" <<
+//                 m_PowellMeta->GetStepLength() << " " << m_PowellMeta->GetStepTolerance() <<  std::endl;
+//    std::cout << "Print Amoeba" <<
+//                m_AmoebaMeta->GetFunctionConvergenceTolerance() << " " << m_AmoebaMeta->GetSimplexDelta() <<  std::endl;
+//    std::cout << "Print NCC" <<
+//                m_NCCMeta->GetSubtractMean() << " " << m_NCCMeta->GetMaxTranslation() <<  std::endl;
+//    std::cout << "Print MI" <<
+//                 m_MIMeta->GetMaxTranslation() <<  std::endl;
+
+//    ExhaustiveOptimizerType::StepsType steps(3);
+//    const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
+    //const double stepLength = 0.1;
+
+    switch (m_r23Meta->GetOptimizerType()) {
+
+    case tOptimizerTypeEnum::POWELL:
+        std::cout<< "Using POWELL Optimizer" <<std::endl;
+
+        PowellOptimizer->SetMaximize(false); // for NCC and MI
+        PowellOptimizer->SetMaximumIteration(m_PowellMeta->GetMaxIterations());
+        PowellOptimizer->SetMaximumLineIteration(m_PowellMeta->GetMaximumLineInteration()); // for Powell's method
+        PowellOptimizer->SetStepLength(m_PowellMeta->GetStepLength());
+        PowellOptimizer->SetStepTolerance(m_PowellMeta->GetStepTolerance());
+        PowellOptimizer->SetValueTolerance(m_PowellMeta->GetValueTolerance());
+        PowellOptimizer->RemoveAllObservers();
+        PowellOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
+        m_OptimizerObserver->setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum::POWELL);
+
+        registration->SetOptimizer(PowellOptimizer);
+
+        break;
+
+    case tOptimizerTypeEnum::AMOEBA:
+        std::cout<< "Using AMOEBA Optimizer" <<std::endl;
+
+        AmoebaOptimizer->SetMinimize(true);
+        AmoebaOptimizer->SetMaximumNumberOfIterations(m_AmoebaMeta->GetMaxIterations());
+        AmoebaOptimizer->SetParametersConvergenceTolerance(m_AmoebaMeta->GetParametersConvergenceTolerance());
+        AmoebaOptimizer->SetFunctionConvergenceTolerance(m_AmoebaMeta->GetFunctionConvergenceTolerance());
+        AmoebaOptimizer->SetAutomaticInitialSimplex(false);
+        //Initial size of the simplex (otherwise it is a very small number and optimizer stops immeditaly)
+        // 2 mm / 2 degrees seems to be a good value which performs well.
+        if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 0) {
+            itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
+            return;
+        }else if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 3){
+            AmoebaOptimizerType::ParametersType simplexDelta3dof(3);
+            for (int i = 0; i < GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()); i++) {
+                simplexDelta3dof[i] = m_AmoebaMeta->GetSimplexDelta();
+            }
+            AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta3dof);
+
+        }else if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 6){
+            AmoebaOptimizerType::ParametersType simplexDelta6dof(6);
+            for (int i = 0; i < GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()); i++) {
+                simplexDelta6dof[i] = m_AmoebaMeta->GetSimplexDelta();
+            }
+            AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta6dof);
+        }
+
+        AmoebaOptimizer->RemoveAllObservers();
+        AmoebaOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
+        m_OptimizerObserver->setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum::AMOEBA);
+
+        registration->SetOptimizer(AmoebaOptimizer);
+
+        break;
+
+    case tOptimizerTypeEnum::EXHAUSTIVE:
+        std::cout<< "Using Extensive Optimizer" <<std::endl;
+
+//        steps[0] = numberOfSteps;
+//        steps[1] = numberOfSteps;
+//        steps[2] = numberOfSteps;
+//        ExhaustiveOptimizer->SetNumberOfSteps(steps);
+//        ExhaustiveOptimizer->SetStepLength(stepLength);
+
+//        registration->SetOptimizer(ExhaustiveOptimizer);
+
+        break;
+
+    default:
+        break;
+
+    }
+
+    switch (m_r23Meta->GetMetricType()) {
+
+    case tMetricTypeEnum::MI:
+        std::cout<< "Using MI Metric" <<std::endl;
+
+        registration->SetMetric(MImetric);
+        MImetric->SetNumberOfHistogramBins(m_MIMeta->GetNumberOfHistogramBins());
+        MImetric->ComputeGradientOff();
+        MImetric->SetMaxTranslation(m_MIMeta->GetMaxTranslation());
+
+        break;
+
+    case tMetricTypeEnum::NCC:
+        std::cout<< "Using NCC Metric" <<std::endl;
+        registration->SetMetric(NCCmetric);
+        NCCmetric->ComputeGradientOff();
+        NCCmetric->SetSubtractMean(true);
+        NCCmetric->SetMaxTranslation(m_NCCMeta->GetMaxTranslation());
+
+        break;
+
+    default:
+        break;
+
+    }
+
+    registration->SetFixedImage1(m_PA);
+    registration->SetFixedImage2(m_LAT);
+    registration->SetMovingImage(m_Volume);
+
+
+    if( this->GetroiAutoReg1().GetNumberOfPixels() == 0 ){
+        auto defaultRegion = m_PA->GetBufferedRegion();
+        this->SetroiAutoReg1(defaultRegion);
+        std::cout<< "Image1 region not defined. Using the whole image buffered region" <<std::endl;
+    }
+    if( this->GetroiAutoReg2().GetNumberOfPixels() == 0 ){
+        auto defaultRegion = m_LAT->GetBufferedRegion();
+        this->SetroiAutoReg2(defaultRegion);
+        std::cout<< "Image2 region not defined. Using the whole image buffered region" <<std::endl;
+    }
+
+    registration->SetFixedImageRegion1(m_roiAutoReg1);
+    registration->SetFixedImageRegion2(m_roiAutoReg2);
+
+    registration->SetTransformMetaInfo(m_r23Meta);
+
+    registration->SetinternalMeta1(m_InternalTransf1);
+    registration->SetinternalMeta2(m_InternalTransf2);
+
+    registration->SetInterpolator1(m_interpolator1);
+    registration->SetInterpolator2(m_interpolator2);
+
+    registration->SetFilter1(m_filter1);
+    registration->SetFilter2(m_filter2);
+
+
+    IsocTransform = TransformType::New();
+    IsocTransform->SetComputeZYX(true);
+    IsocTransform->SetIdentity();
+
+    IsocTransform->SetRotation(
+                m_TransformMetaInfo->GetR()[0],
+                m_TransformMetaInfo->GetR()[1],
+                m_TransformMetaInfo->GetR()[2]
+            );
+
+    TransformType::OutputVectorType TranslV;
+    TranslV[0] = m_TransformMetaInfo->GetT()[0];
+    TranslV[1] = m_TransformMetaInfo->GetT()[1];
+    TranslV[2] = m_TransformMetaInfo->GetT()[2];
+
+    IsocTransform->SetTranslation(TranslV);
+
+    registration->SetIsocIECTransform(IsocTransform);
+
+    //    if (verbose) {
+    //        registration->DebugOn();
+    //        registration->Print(std::cout);
+    //    }
+
+
+    std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
+}
+
+
+
+int itkReg23::StartRegistration(std::string extraInfo)
+{
+
+    std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
+    // TODO: Check if the registartion pipeline has been initialized
+    using ParametersType = RegistrationType::ParametersType;
+
+    //auto startParameters = transform1->GetParameters();
+
+    time_t t = time(0); // get time now
+    struct tm* now = localtime(&t);
+
+    bool useDumpToFile = false;
+
+    char buffer[255];
+    strftime(buffer, 255, "test-%F-%H%M%S.txt", now);
+
+    std::fstream fs;
+
+    if (useDumpToFile) {
+        fs.open(buffer, std::fstream::out);
+        fs << extraInfo;
+        fs << "Value\tX\tY\tZ " << std::endl;
+
+        //        if (r23Meta->GetOptimizerType() == tOptimizerTypeEnum::EXHAUSTIVE) {
+        //            ExhaustiveOptimizerObserver->set_stream(fs);
+        //        }
+    }
+
+    // Start the registration
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Create a timer to record calculation time.
+    itk::TimeProbesCollectorBase timer;
+
+    if (verbose) {
+        std::cout << "Starting the registration now..." << std::endl;
+    }
+
+    try {
+        //        timer.Start("Registration");
+        // Start the registration.
+        registration->StartRegistration();
+        //        timer.Stop("Registration");
+    } catch (itk::ExceptionObject& err) {
+
+        registration->ResetPipeline();
+        std::cout << "ExceptionObject caught !" << std::endl;
+        std::cout << err << std::endl;
+        return -1;
+    }
+
+    fs.close();
+
+
+    itk::Optimizer::ParametersType finalParameters =
+        this->GetCurrentPosition();
+    std::cout<<" FinalPars: "<< finalParameters <<std::endl;
+
+
+    std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
+
+    return 0;
+}
+
+double itkReg23::GetOptimizerValue()
+{
+    return m_OptmizerValue;
+}
+
+
+/** Get the cost function value for the current transform*/
+Optimizer::ParametersType itkReg23::GetCurrentPosition(){
+
+
+    itk::Optimizer::ParametersType finalParameters;
+
+    switch (m_r23Meta->GetOptimizerType()) {
+
+    case tOptimizerTypeEnum::POWELL:
+        finalParameters = PowellOptimizer->GetCurrentPosition();
+        m_OptmizerValue = PowellOptimizer->GetValue();
+        break;
+
+    case tOptimizerTypeEnum::AMOEBA:
+        finalParameters = AmoebaOptimizer->GetCurrentPosition();
+        m_OptmizerValue = AmoebaOptimizer->GetValue();
+        break;
+
+    case tOptimizerTypeEnum::EXHAUSTIVE:
+        finalParameters = ExhaustiveOptimizer->GetCurrentPosition();
+        break;
+
+    default:
+        break;
+
+    }
+
+    return finalParameters;
+}
+
+double itkReg23::GetCurrentMetricValue()
+{
+
+    switch (m_r23Meta->GetMetricType()) {
+
+    case tMetricTypeEnum::MI:
+        if(!MImetric){
+            itkExceptionMacro(<< "MImetric not present");
+            return -1.;
+        } else {
+            return MImetric->GetValue();
+        }
+
+        break;
+
+    case tMetricTypeEnum::NCC:
+        if(!NCCmetric){
+            itkExceptionMacro(<< "NCCmetric not present");
+            return -1.;
+        } else {
+            return NCCmetric->GetValue();
+        }
+        break;
+
+    default:
+        break;
+
+    }
+
+    return -1.;
+
+}
+
+
+
+
+}

itkReg23DRT/itkReg23.h

@@ -0,0 +1,224 @@
+#ifndef ITKREG23_H
+#define ITKREG23_h
+
+#include "itkCommand.h"
+
+#include "itkExhaustiveOptimizer.h"
+#include "itkMutualInformationTwoImageToOneImageMetric.h"
+#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
+#include "itkPowellOptimizer.h"
+#include "itkTwoProjectionImageRegistrationMethod.h"
+#include "itkAmoebaOptimizer.h"
+
+#include "itkCommand.h"
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "itkSmartPointer.h"
+
+#include "itkImageProcessorHelpers.h"
+#include "itkQtIterationUpdate.h"
+
+#include "itkReg23MetaInformation.h"
+
+
+namespace itk{
+
+class ITK_EXPORT itkReg23 : public itk::Object
+{
+
+public:
+
+    /** Standard "Self" typedef. */
+    typedef itkReg23  Self;
+    /** Standard "Superclass" typedef. */
+    typedef Object Superclass;
+    /** Smart pointer typedef support */
+    typedef itk::SmartPointer<Self>  Pointer;
+    typedef itk::SmartPointer<const Self>  ConstPointer;
+    /** Method of creation through the object factory. */
+    itkNewMacro(Self);
+    /** Run-time type information (and related methods). */
+    itkTypeMacro(itkReg23, Object);
+
+    using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
+    using RoiForRegistration = itk::ImageRegion<Dimension>;
+    using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
+
+    /* ---- User provided */
+
+    itkSetMacro(OptimizerObserver, CommandIterationUpdate::Pointer);
+    itkGetMacro(OptimizerObserver, CommandIterationUpdate::Pointer);
+
+    itkSetMacro(r23Meta, R23MetaInformation::Pointer);
+    itkGetMacro(r23Meta, R23MetaInformation::Pointer);
+
+    itkSetMacro(PowellMeta, PowellOptimizerMetaInformation::Pointer);
+    itkGetMacro(PowellMeta, PowellOptimizerMetaInformation::Pointer);
+
+    itkSetMacro(AmoebaMeta, AmoebaOptimizerMetaInformation::Pointer);
+    itkGetMacro(AmoebaMeta, AmoebaOptimizerMetaInformation::Pointer);
+
+    itkSetMacro(MIMeta, MIMetricMetaInformation::Pointer);
+    itkGetMacro(MIMeta, MIMetricMetaInformation::Pointer);
+
+    itkSetMacro(NCCMeta, NCCMetricMetaInformation::Pointer);
+    itkGetMacro(NCCMeta, NCCMetricMetaInformation::Pointer);
+
+
+    itkSetMacro(Volume, InternalImageType::Pointer);
+    itkGetMacro(Volume, InternalImageType::Pointer);
+
+    itkSetMacro(PA, InternalImageType::Pointer);
+    itkGetMacro(PA, InternalImageType::Pointer);
+
+    itkSetMacro(LAT, InternalImageType::Pointer);
+    itkGetMacro(LAT, InternalImageType::Pointer);
+
+    itkSetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer);
+    itkGetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer);
+
+    itkSetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer);
+    itkGetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer);
+
+    itkSetMacro(filter1, ChangeInformationFilterType::Pointer);
+    itkGetMacro(filter1, ChangeInformationFilterType::Pointer);
+
+    itkSetMacro(filter2, ChangeInformationFilterType::Pointer);
+    itkGetMacro(filter2, ChangeInformationFilterType::Pointer);
+
+    itkSetMacro(interpolator1, InterpolatorType::Pointer);
+    itkGetMacro(interpolator1, InterpolatorType::Pointer);
+
+    itkSetMacro(interpolator2, InterpolatorType::Pointer);
+    itkGetMacro(interpolator2, InterpolatorType::Pointer);
+
+    itkSetMacro(TransformMetaInfo, TransformMetaInformation::Pointer);
+    itkGetMacro(TransformMetaInfo, TransformMetaInformation::Pointer);
+
+    /* ---- User provided  END */
+    itkSetMacro(roiAutoReg1, RoiForRegistration);
+    itkGetMacro(roiAutoReg1, RoiForRegistration);
+
+    itkSetMacro(roiAutoReg2, RoiForRegistration);
+    itkGetMacro(roiAutoReg2, RoiForRegistration);
+
+
+
+    /** Auto Reg23 methods */
+
+    /** Initialize the registration pipeline*/
+    void InitializeRegistration();
+    /** Start the registration process*/
+    int StartRegistration(std::string extraInfo);
+    /** Get the current cost function value from the optimizer*/
+    double GetOptimizerValue();
+
+    /** Get the cost function value for the current transform*/
+    double GetCurrentMetricValue();
+
+    /** Get the cost function value for the current transform*/
+    Optimizer::ParametersType GetCurrentPosition();
+
+    /** Auto Reg23 methods END */
+
+    void setROISizeToZero();
+
+protected:
+    itkReg23();
+    virtual ~itkReg23();
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+private:
+    itkReg23(const Self&); //purposely not implemented
+    void operator=(const Self&); //purposely not implemented
+
+    /** Optimizer which tries to find the minimun (Powell Optimizer)*/
+    using PowellOptimizerType = itk::PowellOptimizer;
+    /** Optimizer which tries to find the minimn (Powell Optimizer)*/
+    using AmoebaOptimizerType = itk::AmoebaOptimizer;
+    /** Optimizer which samples the whol space */
+    using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer;
+    /** Metric to calculate similarites between images*/
+    using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
+    using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
+    /** The thing which actuall does the image registration*/
+    using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
+
+
+
+    RegistrationType::Pointer
+    registration;
+    MetricType::Pointer
+    NCCmetric;
+    MIMetricType::Pointer
+    MImetric;
+    PowellOptimizerType::Pointer
+    PowellOptimizer;
+    AmoebaOptimizerType::Pointer
+    AmoebaOptimizer;
+    ExhaustiveOptimizerType::Pointer
+    ExhaustiveOptimizer;
+
+    /* ---- User provided */
+    R23MetaInformation::Pointer
+    m_r23Meta;
+
+    PowellOptimizerMetaInformation::Pointer
+    m_PowellMeta;
+
+    AmoebaOptimizerMetaInformation::Pointer
+    m_AmoebaMeta;
+
+    MIMetricMetaInformation::Pointer
+    m_MIMeta;
+
+    NCCMetricMetaInformation::Pointer
+    m_NCCMeta;
+
+    InternalImageType::Pointer
+    m_Volume,
+    m_PA,
+    m_LAT;
+
+    InternalTransformMetaInformation::Pointer
+    m_InternalTransf1,
+    m_InternalTransf2;
+
+    ChangeInformationFilterType::Pointer
+    m_filter1,
+    m_filter2;
+
+    InterpolatorType::Pointer
+    m_interpolator1,
+    m_interpolator2;
+
+    TransformMetaInformation::Pointer
+    m_TransformMetaInfo;
+
+    CommandIterationUpdate::Pointer
+    m_OptimizerObserver;
+
+    /* ----  User provided  END */
+
+
+    TransformType::Pointer
+    IsocTransform,
+    Transform1,
+    Transform2;
+
+    RoiForRegistration
+    m_roiAutoReg1,
+    m_roiAutoReg2;
+
+    bool m_UseDumptoFile;
+
+
+    double m_OptmizerValue;
+
+
+};
+
+}
+
+#endif

itkReg23DRT/itkReg23MetaInformation.cpp

@@ -0,0 +1,104 @@
+#include "itkReg23MetaInformation.h"
+
+
+namespace itk {
+
+PowellOptimizerMetaInformation::
+PowellOptimizerMetaInformation(){
+
+    this->m_MaximumLineInteration = 4;
+    this->m_StepLength = 2.0;
+    this->m_StepTolerance = 0.01;
+    this->m_ValueTolerance = 0.000002;
+
+    this->m_MaxIterations = 6;
+}
+
+void
+PowellOptimizerMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+PowellOptimizerMetaInformation
+::~PowellOptimizerMetaInformation ()
+{
+
+}
+
+
+AmoebaOptimizerMetaInformation::
+AmoebaOptimizerMetaInformation(){
+
+    this->m_ParametersConvergenceTolerance = 0.1;
+    this->m_FunctionConvergenceTolerance = 0.000002;
+    this->m_SimplexDelta = 2.0;
+
+    this->m_MaxIterations = 100;
+}
+
+
+void
+AmoebaOptimizerMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+AmoebaOptimizerMetaInformation
+::~AmoebaOptimizerMetaInformation ()
+{
+
+}
+
+
+MIMetricMetaInformation::
+MIMetricMetaInformation(){
+
+    this->m_MaxTranslation = 100.0;
+    this->m_NumberOfHistogramBins = 50;
+
+}
+
+
+void
+MIMetricMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+MIMetricMetaInformation
+::~MIMetricMetaInformation ()
+{
+
+}
+
+NCCMetricMetaInformation::
+NCCMetricMetaInformation(){
+
+    this->m_MaxTranslation = 100.0;
+    this->m_SubtractMean = true;
+
+}
+
+
+void
+NCCMetricMetaInformation
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+
+NCCMetricMetaInformation
+::~NCCMetricMetaInformation ()
+{
+
+}
+
+}

itkReg23DRT/itkReg23MetaInformation.h

@@ -0,0 +1,232 @@
+#ifndef ITKREG23METAINFORMATION_H
+#define ITKREG23METAINFORMATION_H
+
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "itkSmartPointer.h"
+#include "itkMacro.h"
+
+
+namespace itk {
+
+
+class PowellOptimizerMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef PowellOptimizerMetaInformation  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(PowellOptimizerMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    itkSetMacro(StepTolerance,double);
+    itkGetMacro(StepTolerance,double);
+
+    itkSetMacro(ValueTolerance,double);
+    itkGetMacro(ValueTolerance,double);
+
+    itkSetMacro(StepLength,double);
+    itkGetMacro(StepLength,double);
+
+    itkSetMacro(MaximumLineInteration,int);
+    itkGetMacro(MaximumLineInteration,int);
+
+    itkSetEnumMacro(MaxIterations, int);
+    itkGetEnumMacro(MaxIterations, int);
+
+
+protected:
+
+    double m_StepTolerance;
+
+    double m_ValueTolerance;
+
+    double m_StepLength;
+
+    int m_MaximumLineInteration;
+
+    int
+    m_MaxIterations;
+
+    /** Default Constructor **/
+    PowellOptimizerMetaInformation  ();
+    /** Default Destructor **/
+    virtual ~PowellOptimizerMetaInformation  ();
+
+private:
+    /** purposely not implemented **/
+    PowellOptimizerMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+class AmoebaOptimizerMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef AmoebaOptimizerMetaInformation  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(AmoebaOptimizerMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    itkSetMacro(ParametersConvergenceTolerance,double);
+    itkGetMacro(ParametersConvergenceTolerance,double);
+
+    itkSetMacro(FunctionConvergenceTolerance,double);
+    itkGetMacro(FunctionConvergenceTolerance,double);
+
+    itkSetMacro(SimplexDelta,double);
+    itkGetMacro(SimplexDelta,double);
+
+    itkSetEnumMacro(MaxIterations, int);
+    itkGetEnumMacro(MaxIterations, int);
+
+protected:
+
+    double m_ParametersConvergenceTolerance;
+
+    double m_FunctionConvergenceTolerance;
+
+    double m_SimplexDelta;
+
+    int
+    m_MaxIterations;
+
+    /** Default Constructor **/
+    AmoebaOptimizerMetaInformation  ();
+    /** Default Destructor **/
+    virtual ~AmoebaOptimizerMetaInformation  ();
+
+private:
+    /** purposely not implemented **/
+    AmoebaOptimizerMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+class MIMetricMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef MIMetricMetaInformation  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(MIMetricMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    itkSetMacro(MaxTranslation,double);
+    itkGetMacro(MaxTranslation,double);
+
+    itkSetMacro(NumberOfHistogramBins,int);
+    itkGetMacro(NumberOfHistogramBins,int);
+
+protected:
+
+    double m_MaxTranslation;
+
+    int m_NumberOfHistogramBins;
+
+
+    /** Default Constructor **/
+    MIMetricMetaInformation  ();
+    /** Default Destructor **/
+    virtual ~MIMetricMetaInformation  ();
+
+private:
+    /** purposely not implemented **/
+    MIMetricMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+
+class NCCMetricMetaInformation :
+        public itk::Object{
+
+public:
+    /** standard typedefs **/
+    typedef NCCMetricMetaInformation  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(NCCMetricMetaInformation, itk::Object);
+
+    /** object information streaming **/
+    void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+
+    itkSetMacro(MaxTranslation,double);
+    itkGetMacro(MaxTranslation,double);
+
+    itkSetMacro(SubtractMean,bool);
+    itkGetMacro(SubtractMean,bool);
+
+protected:
+
+    double m_MaxTranslation;
+
+    bool m_SubtractMean;
+
+
+
+    /** Default Constructor **/
+    NCCMetricMetaInformation  ();
+    /** Default Destructor **/
+    virtual ~NCCMetricMetaInformation  ();
+
+private:
+    /** purposely not implemented **/
+    NCCMetricMetaInformation (const Self&);
+
+    /** purposely not implemented **/
+    void operator=(const Self&);
+};
+
+}
+
+
+
+
+
+
+
+
+#endif

itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.h

@@ -0,0 +1,260 @@
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+/*=========================================================================
+Calculate DRR from a CT dataset using incremental ray-tracing algorithm
+The algorithm was initially proposed by Robert Siddon and improved by
+Filip Jacobs etc.
+
+-------------------------------------------------------------------------
+References:
+
+R. L. Siddon, "Fast calculation of the exact radiological path for a
+threedimensional CT array," Medical Physics 12, 252-55 (1985).
+
+F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
+"A fast algorithm to calculate the exact radiological path through a pixel
+or voxel space," Journal of Computing and Information Technology ?
+CIT 6, 89-94 (1998).
+
+=========================================================================*/
+#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_h
+#define itkgSiddonJacobsRayCastInterpolateImageFunction_h
+
+#include "itkInterpolateImageFunction.h"
+#include "itkTransform.h"
+#include "itkVector.h"
+#include "itkEuler3DTransform.h"
+
+namespace itk
+{
+
+/** \class gSiddonJacobsRayCastInterpolateImageFunction
+ * \brief Projective interpolation of an image at specified positions.
+ *
+ * SiddonJacobsRayCastInterpolateImageFunction casts rays through a 3-dimensional
+ * image
+ * \warning This interpolator works for 3-dimensional images only.
+ *
+ * \ingroup ImageFunctions
+ * \ingroup TwoProjectionRegistration
+ */
+template <typename TInputImage, typename TCoordRep = float>
+class gSiddonJacobsRayCastInterpolateImageFunction : public InterpolateImageFunction<TInputImage, TCoordRep>
+{
+public:
+  ITK_DISALLOW_COPY_AND_ASSIGN(gSiddonJacobsRayCastInterpolateImageFunction);
+
+  /** Standard class type alias. */
+  using Self = gSiddonJacobsRayCastInterpolateImageFunction;
+  using Superclass = InterpolateImageFunction<TInputImage, TCoordRep>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** Constants for the image dimensions */
+  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
+
+
+  using TransformType = Euler3DTransform<TCoordRep>;
+
+  using TransformPointer = typename TransformType::Pointer;
+  using InputPointType = typename TransformType::InputPointType;
+  using OutputPointType = typename TransformType::OutputPointType;
+  using TransformParametersType = typename TransformType::ParametersType;
+  using TransformJacobianType = typename TransformType::JacobianType;
+
+  using PixelType = typename Superclass::InputPixelType;
+
+  using SizeType = typename TInputImage::SizeType;
+
+  using DirectionType = Vector<TCoordRep, 3>;
+
+  /**  Type of the Interpolator Base class */
+  using InterpolatorType = InterpolateImageFunction<TInputImage, TCoordRep>;
+
+  using InterpolatorPointer = typename InterpolatorType::Pointer;
+
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(gSiddonJacobsRayCastInterpolateImageFunction, InterpolateImageFunction);
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** OutputType type alias support. */
+  using OutputType = typename Superclass::OutputType;
+
+  /** InputImageType type alias support. */
+  using InputImageType = typename Superclass::InputImageType;
+
+  /** InputImageConstPointer type alias support. */
+  using InputImageConstPointer = typename Superclass::InputImageConstPointer;
+
+  /** RealType type alias support. */
+  using RealType = typename Superclass::RealType;
+
+  /** Dimension underlying input image. */
+  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
+
+  /** Point type alias support. */
+  using PointType = typename Superclass::PointType;
+
+  /** Index type alias support. */
+  using IndexType = typename Superclass::IndexType;
+
+  /** ContinuousIndex type alias support. */
+  using ContinuousIndexType = typename Superclass::ContinuousIndexType;
+
+
+  /** Internal floating point comparison accuracy **/
+  static const double EPSILON;
+
+  /** \brief
+   * Interpolate the image at a point position.
+   *
+   * Returns the interpolated image intensity at a
+   * specified point position. No bounds checking is done.
+   * The point is assume to lie within the image buffer.
+   *
+   * ImageFunction::IsInsideBuffer() can be used to check bounds before
+   * calling the method.
+   */
+  OutputType
+  Evaluate(const PointType & point) const override;
+
+  /** Interpolate the image at a continuous index position
+   *
+   * Returns the interpolated image intensity at a
+   * specified index position. No bounds checking is done.
+   * The point is assume to lie within the image buffer.
+   *
+   * Subclasses must override this method.
+   *
+   * ImageFunction::IsInsideBuffer() can be used to check bounds before
+   * calling the method.
+   */
+  OutputType
+  EvaluateAtContinuousIndex(const ContinuousIndexType & index) const override;
+
+  virtual void
+  Initialize();
+
+  /** Connect the Transform. */
+  itkSetObjectMacro(Transform, TransformType);
+  /** Get a pointer to the Transform.  */
+  itkGetConstObjectMacro(Transform, TransformType);
+
+  /** Get a pointer to the Inverse Transform. used to calculate the rays  */
+  itkGetConstObjectMacro(InverseTransform, TransformType);
+  itkGetConstObjectMacro(ComposedTransform, TransformType);
+
+  /** Set and get the focal point to isocenter distance in mm */
+  itkSetMacro(FocalPointToIsocenterDistance, double);
+  itkGetMacro(FocalPointToIsocenterDistance, double);
+
+  /** Set and get the Lianc grantry rotation angle in radians */
+  itkSetMacro(ProjectionAngle, double);
+  itkGetMacro(ProjectionAngle, double);
+
+  /** Set and get the Threshold */
+  itkSetMacro(Threshold, double);
+  itkGetMacro(Threshold, double);
+
+  /** Set and get the Panel Offset */
+  itkSetMacro(PanelOffset, double);
+  itkGetMacro(PanelOffset, double);
+
+  /** Check if a point is inside the image buffer.
+   * \warning For efficiency, no validity checking of
+   * the input image pointer is done. */
+  inline bool
+  IsInsideBuffer(const PointType &) const override
+  {
+    return true;
+  }
+  bool
+  IsInsideBuffer(const ContinuousIndexType &) const override
+  {
+    return true;
+  }
+  bool
+  IsInsideBuffer(const IndexType &) const override
+  {
+    return true;
+  }
+
+#if !defined(ITKV4_COMPATIBILITY)
+  SizeType
+  GetRadius() const override
+  {
+    const InputImageType * input = this->GetInputImage();
+    if (!input)
+    {
+      itkExceptionMacro("Input image required!");
+    }
+    return input->GetLargestPossibleRegion().GetSize();
+  }
+#endif
+
+protected:
+  gSiddonJacobsRayCastInterpolateImageFunction();
+
+  ~gSiddonJacobsRayCastInterpolateImageFunction() override = default;
+
+  void
+  PrintSelf(std::ostream & os, Indent indent) const override;
+
+  /// Transformation used to calculate the new focal point position
+  TransformPointer m_Transform; // Displacement of the volume
+  // Overall inverse transform used to calculate the ray position in the input space
+  TransformPointer m_InverseTransform;
+  TransformPointer m_ComposedTransform;  // Composed transform
+
+  // The threshold above which voxels along the ray path are integrated
+  double m_Threshold;
+  double m_FocalPointToIsocenterDistance; // Focal point to isocenter distance
+  double m_ProjectionAngle;               // Linac gantry rotation angle in radians
+  double m_PanelOffset;
+
+private:
+  void
+                   ComputeInverseTransform() const;
+  TransformPointer m_GantryRotTransform; // Gantry rotation transform
+  TransformPointer m_CamShiftTransform;  // Camera shift transform camRotTransform
+  TransformPointer m_CamRotTransform;    // Camera rotation transform
+  PointType        m_SourcePoint;        // Coordinate of the source in the standard Z projection geometry
+  PointType        m_SourceWorld;        // Coordinate of the source in the world coordinate system
+};
+
+} // namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkgSiddonJacobsRayCastInterpolateImageFunction.hxx"
+#endif
+
+#endif

itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx

@@ -0,0 +1,565 @@
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+/*=========================================================================
+The algorithm was initially proposed by Robert Siddon and improved by
+Filip Jacobs etc.
+
+-------------------------------------------------------------------------
+References:
+
+R. L. Siddon, "Fast calculation of the exact radiological path for a
+threedimensional CT array," Medical Physics 12, 252-55 (1985).
+
+F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
+"A fast algorithm to calculate the exact radiological path through a pixel
+or voxel space," Journal of Computing and Information Technology ?
+CIT 6, 89-94 (1998).
+
+=========================================================================*/
+
+#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
+#define itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
+
+#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
+
+#include "itkMath.h"
+#include <cstdlib>
+#include <unistd.h>
+
+namespace itk
+{
+
+
+template <typename TInputImage, typename TCoordRep>
+const double gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EPSILON =
+        0.00001;
+
+template <typename TInputImage, typename TCoordRep>
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::gSiddonJacobsRayCastInterpolateImageFunction()
+{
+    m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm.
+    m_ProjectionAngle = 0.;                  // Angle in radians betweeen projection central axis and reference axis
+    m_Threshold = 0.;                        // Intensity threshold, below which is ignored.
+    m_PanelOffset = 0.;
+
+    m_SourcePoint[0] = 0.;
+    m_SourcePoint[1] = 0.;
+    m_SourcePoint[2] = 0.;
+
+    m_InverseTransform = TransformType::New();
+    m_InverseTransform->SetComputeZYX(true);
+
+    m_ComposedTransform = TransformType::New();
+    m_ComposedTransform->SetComputeZYX(true);
+
+    m_GantryRotTransform = TransformType::New();
+    m_GantryRotTransform->SetComputeZYX(true);
+    m_GantryRotTransform->SetIdentity();
+
+    m_CamShiftTransform = TransformType::New();
+    m_CamShiftTransform->SetComputeZYX(true);
+    m_CamShiftTransform->SetIdentity();
+
+    m_CamRotTransform = TransformType::New();
+    m_CamRotTransform->SetComputeZYX(true);
+    m_CamRotTransform->SetIdentity();
+    // constant for converting degrees into radians
+    const float dtr = (atan(1.0) * 4.0) / 180.0;
+    m_CamRotTransform->SetRotation(dtr * (-90.0), 0.0, 0.0);
+
+    m_Threshold = 0;
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+void
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
+{
+    this->Superclass::PrintSelf(os, indent);
+
+    os << indent << "Threshold: " << m_Threshold << std::endl;
+    os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
+{
+    float     rayVector[3];
+    IndexType cIndex;
+
+    PointType  drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
+    OutputType pixval;
+
+
+    float firstIntersection[3];
+    float alphaX1, alphaXN, alphaXmin, alphaXmax;
+    float alphaY1, alphaYN, alphaYmin, alphaYmax;
+    float alphaZ1, alphaZN, alphaZmin, alphaZmax;
+    float alphaMin, alphaMax;
+    float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
+    float alphaUx, alphaUy, alphaUz;
+    float alphaIntersectionUp[3], alphaIntersectionDown[3];
+    float d12, value;
+    float firstIntersectionIndex[3];
+    int   firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
+    int   iU, jU, kU;
+
+
+    // Min/max values of the output pixel type AND these values
+    // represented as the output type of the interpolator
+    const OutputType minOutputValue = itk::NumericTraits<OutputType>::NonpositiveMin();
+    const OutputType maxOutputValue = itk::NumericTraits<OutputType>::max();
+
+
+    // If the volume was shifted, recalculate the overall inverse transform
+    unsigned long int interpMTime = this->GetMTime();
+    unsigned long int vTransformMTime = m_Transform->GetMTime();
+
+    if (interpMTime < vTransformMTime)
+    {
+        this->ComputeInverseTransform();
+        // The m_SourceWorld should be computed here to avoid the repeatedly calculation
+        // for each projection ray. However, we are in a const function, which prohibits
+        // the modification of class member variables. So the world coordiate of the source
+        // point is calculated for each ray as below. Performance improvement may be made
+        // by using a static variable?
+        // m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
+    }
+
+    PointType PointReq = point;
+    PointReq[0] += m_PanelOffset;
+
+    drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
+
+    // Get ths input pointers
+    InputImageConstPointer inputPtr = this->GetInputImage();
+
+    typename InputImageType::SizeType    sizeCT;
+    typename InputImageType::RegionType  regionCT;
+    typename InputImageType::SpacingType ctPixelSpacing;
+    typename InputImageType::PointType   ctOrigin;
+
+    ctPixelSpacing = inputPtr->GetSpacing();
+    ctOrigin = inputPtr->GetOrigin();
+    regionCT = inputPtr->GetLargestPossibleRegion();
+    sizeCT = regionCT.GetSize();
+
+
+    PointType SlidingSourcePoint = m_SourcePoint;
+    SlidingSourcePoint[0] = 0.;
+    // Simple sliding source model, the source follows the y of the pixel
+    // requested, which sets the Z position to look in the CT volume.
+    SlidingSourcePoint[1] = point[1];
+    SlidingSourcePoint[2] = 0.;
+
+    PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
+
+    PointType O(3);
+    O[0] = -ctPixelSpacing[0]/2.;
+    O[1] = -ctPixelSpacing[1]/2.;
+    O[2] = -ctPixelSpacing[2]/2.;
+
+    // be sure that drr pixel position is inside the volume
+    float xSizeCT = O[0] + (sizeCT[0]) * ctPixelSpacing[0];
+    float ySizeCT = O[1] + (sizeCT[1]) * ctPixelSpacing[1];
+    float zSizeCT = O[2] + (sizeCT[2]) * ctPixelSpacing[2];
+    float xDrrPix = drrPixelWorld[0];
+    float yDrrPix = drrPixelWorld[1];
+    float zDrrPix = drrPixelWorld[2];
+    if(
+            xDrrPix<= O[0] || xDrrPix>=xSizeCT ||
+            yDrrPix<= O[1] || yDrrPix>=ySizeCT ||
+            zDrrPix<= O[2] || zDrrPix>=zSizeCT
+            ){
+        pixval = static_cast<OutputType>(0.0);
+        return pixval;
+    }
+
+    // calculate the detector position for this pixel center by moving
+    // 2*m_FocalPointToIsocenterDistance from the source in the pixel
+    // directions
+    double p1[3], p2[3];
+
+    p1[0] = SourceWorld[0];
+    p1[1] = SourceWorld[1];
+    p1[2] = SourceWorld[2];
+
+    p2[0] = drrPixelWorld[0];
+    p2[1] = drrPixelWorld[1];
+    p2[2] = drrPixelWorld[2];
+
+    double P12D = sqrt(
+                (p2[0]-p1[0]) * (p2[0]-p1[0]) +
+            (p2[1]-p1[1]) * (p2[1]-p1[1]) +
+            (p2[2]-p1[2]) * (p2[2]-p1[2])
+            );
+
+    double p12V[3];
+    p12V[0] = p2[0] - p1[0];
+    p12V[1] = p2[1] - p1[1];
+    p12V[2] = p2[2] - p1[2];
+
+    double p12v [3];
+    p12v[0] = p12V[0]/P12D;
+    p12v[1] = p12V[1]/P12D;
+    p12v[2] = p12V[2]/P12D;
+
+    double pDest [3];
+    pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
+    pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
+    pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2 * p12v[2];
+
+
+    // The following is the Siddon-Jacob fast ray-tracing algorithm
+    //rayVector[0] = drrPixelWorld[0] - SourceWorld[0];
+    //rayVector[1] = drrPixelWorld[1] - SourceWorld[1];
+    // correct ray to reach detector for topogram geometry
+    // Differs from conventional FPD geometry (above) that has
+    // the panel at the end of the ray.
+    rayVector[0] = pDest[0] - SourceWorld[0];
+    rayVector[1] = pDest[1] - SourceWorld[1];
+    rayVector[2] = pDest[2] - SourceWorld[2];
+
+
+    /* Calculate the parametric  values of the first  and  the  last
+  intersection points of  the  ray  with the X,  Y, and Z-planes  that
+  define  the  CT volume. */
+    if (fabs(rayVector[0]) > EPSILON/* != 0*/)
+    {
+        //    alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
+        //    alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+        alphaX1 = ( O[0] - SourceWorld[0]) / rayVector[0];
+        alphaXN = ( O[0] + sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+        alphaXmin = std::min(alphaX1, alphaXN);
+        alphaXmax = std::max(alphaX1, alphaXN);
+    }
+    else
+    {
+        alphaXmin = -2;
+        alphaXmax = 2;
+    }
+
+    if (fabs(rayVector[1]) > EPSILON/* != 0*/)
+    {
+        alphaY1 = ( O[1] - SourceWorld[1]) / rayVector[1];
+        alphaYN = ( O[1] + sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+        //    alphaY1 = (0.0 - SourceWorld[1]) / rayVector[1];
+        //    alphaYN = (sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+        alphaYmin = std::min(alphaY1, alphaYN);
+        alphaYmax = std::max(alphaY1, alphaYN);
+    }
+    else
+    {
+        alphaYmin = -2;
+        alphaYmax = 2;
+    }
+
+    if (fabs(rayVector[2]) > EPSILON/* != 0*/)
+    {
+        //    alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
+        //    alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+        alphaZ1 = ( O[2] - SourceWorld[2]) / rayVector[2];
+        alphaZN = ( O[2] + sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+        alphaZmin = std::min(alphaZ1, alphaZN);
+        alphaZmax = std::max(alphaZ1, alphaZN);
+    }
+    else
+    {
+        alphaZmin = -2;
+        alphaZmax = 2;
+    }
+
+    /* Get the very first and the last alpha values when the ray
+  intersects with the CT volume. */
+    alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin);
+    alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax);
+
+    /* Calculate the parametric values of the first intersection point
+  of the ray with the X, Y, and Z-planes after the ray entered the
+  CT volume. */
+
+    firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0] - O[0] ; //ctPixelSpacing[0]/2.;
+    firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - O[1] ; //ctPixelSpacing[1]/2.;
+    firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2] - O[2] ; //ctPixelSpacing[2]/2.;
+
+    /* Transform world coordinate to the continuous index of the CT volume*/
+    firstIntersectionIndex[0] = firstIntersection[0] / ctPixelSpacing[0];
+    firstIntersectionIndex[1] = firstIntersection[1] / ctPixelSpacing[1];
+    firstIntersectionIndex[2] = firstIntersection[2] / ctPixelSpacing[2];
+
+    // now "correct" the indices (there are special cases):
+    if (alphaMin == alphaYmin) // j is special
+    {
+        if (alphaYmin == alphaY1)
+            firstIntersectionIndex [1] = 0;
+        else
+            // alpha_y_Ny
+            firstIntersectionIndex [1] = sizeCT[1] - 1;
+    }
+    else if (alphaMin == alphaXmin) // i is special
+    {
+        if (alphaXmin == alphaX1)
+            firstIntersectionIndex[0] = 0;
+        else
+            // alpha_x_Nx
+            firstIntersectionIndex[0] = sizeCT[0] - 1;
+    }
+    else if (alphaMin == alphaZmin) // k is special
+    {
+        if (alphaZmin == alphaZ1)
+            firstIntersectionIndex [2] = 0;
+        else
+            // alpha_z_Nz
+            firstIntersectionIndex [2] = sizeCT[2] - 1;
+    }
+
+    firstIntersectionIndexUp[0] = (int)ceil(firstIntersectionIndex[0]);
+    firstIntersectionIndexUp[1] = (int)ceil(firstIntersectionIndex[1]);
+    firstIntersectionIndexUp[2] = (int)ceil(firstIntersectionIndex[2]);
+
+    firstIntersectionIndexDown[0] = (int)floor(firstIntersectionIndex[0]);
+    firstIntersectionIndexDown[1] = (int)floor(firstIntersectionIndex[1]);
+    firstIntersectionIndexDown[2] = (int)floor(firstIntersectionIndex[2]);
+
+
+    if (fabs(rayVector[0]) <= EPSILON/* == 0*/)
+    {
+        alphaX = 2;
+    }
+    else
+    {
+        alphaIntersectionUp[0] =
+                (O[0] + firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+        alphaIntersectionDown[0] =
+                (O[0] + firstIntersectionIndexDown[0] * ctPixelSpacing[0]  - SourceWorld[0]) / rayVector[0];
+        alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
+    }
+
+    if (fabs(rayVector[1] ) <= EPSILON/* == 0*/)
+    {
+        alphaY = 2;
+    }
+    else
+    {
+        alphaIntersectionUp[1] =
+                (O[1] + firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+        alphaIntersectionDown[1] =
+                (O[1] + firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+        alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
+    }
+
+    if (fabs(rayVector[2]) <= EPSILON/* == 0*/)
+    {
+        alphaZ = 2;
+    }
+    else
+    {
+        alphaIntersectionUp[2] =
+                (O[2] + firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+        alphaIntersectionDown[2] =
+                (O[2] + firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+        alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
+    }
+
+    /* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */
+    if (fabs(rayVector[0]) > EPSILON /*!= 0*/)
+    {
+        alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
+    }
+    else
+    {
+        alphaUx = 999;
+    }
+    if (fabs(rayVector[1]) > EPSILON /*!= 0*/)
+    {
+        alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]);
+    }
+    else
+    {
+        alphaUy = 999;
+    }
+    if (fabs(rayVector[2]) > EPSILON /*!= 0*/)
+    {
+        alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]);
+    }
+    else
+    {
+        alphaUz = 999;
+    }
+    /* Calculate voxel index incremental values along the ray path. */
+
+    iU = (SourceWorld[0] < drrPixelWorld[0]) ? 1 : -1;
+    jU = (SourceWorld[1] < drrPixelWorld[1]) ? 1 : -1;
+    kU = (SourceWorld[2] < drrPixelWorld[2]) ? 1 : -1;
+
+    d12 = 0.0; /* Initialize the sum of the voxel intensities along the ray path to zero. */
+
+    /* Initialize the current ray position. */
+    alphaCmin = std::min(std::min(alphaX, alphaY), alphaZ);
+
+    /* Initialize the current voxel index. */
+    cIndex[0] = firstIntersectionIndexDown[0];
+    cIndex[1] = firstIntersectionIndexDown[1];
+    cIndex[2] = firstIntersectionIndexDown[2];
+
+    while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
+    {
+        /* Store the current ray position */
+        alphaCminPrev = alphaCmin;
+
+        if ((alphaX <= alphaY) && (alphaX <= alphaZ))
+        {
+            /* Current ray front intercepts with x-plane. Update alphaX. */
+            alphaCmin = alphaX;
+            cIndex[0] = cIndex[0] + iU;
+            alphaX = alphaX + alphaUx;
+        }
+        else if ((alphaY <= alphaX) && (alphaY <= alphaZ))
+        {
+            /* Current ray front intercepts with y-plane. Update alphaY. */
+            alphaCmin = alphaY;
+            cIndex[1] = cIndex[1] + jU;
+            alphaY = alphaY + alphaUy;
+        }
+        else
+        {
+            /* Current ray front intercepts with z-plane. Update alphaZ. */
+            alphaCmin = alphaZ;
+            cIndex[2] = cIndex[2] + kU;
+            alphaZ = alphaZ + alphaUz;
+        }
+
+        if ((cIndex[0] >= 0) && (cIndex[0] < static_cast<IndexValueType>(sizeCT[0])) && (cIndex[1] >= 0) &&
+                (cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
+                (cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
+        {
+            /* If it is a valid index, get the voxel intensity. */
+            value = static_cast<float>(inputPtr->GetPixel(cIndex));
+            if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
+            {
+                d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
+            }
+        }
+    }
+
+    if (d12 < minOutputValue)
+    {
+        pixval = minOutputValue;
+    }
+    else if (d12 > maxOutputValue)
+    {
+        pixval = maxOutputValue;
+    }
+    else
+    {
+        pixval = static_cast<OutputType>(d12);
+    }
+    return (pixval);
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateAtContinuousIndex(
+        const ContinuousIndexType & index) const
+{
+    OutputPointType point;
+    this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
+
+    return this->Evaluate(point);
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+void
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInverseTransform() const
+{
+    m_ComposedTransform->SetIdentity();
+    m_ComposedTransform->Compose(m_Transform, 0);
+
+    typename TransformType::InputPointType isocenter;
+    isocenter = m_Transform->GetCenter();
+
+//   //This is problably not necessary since we correct the ray.... to be checked..
+//    InputImageConstPointer inputPtr = this->GetInputImage();
+
+//    std::cout<<inputPtr<<std::endl;
+//      if(inputPtr != NULL)
+//    {
+//      isocenter[0] -= inputPtr->GetSpacing()[0]/2.;
+//      isocenter[1] -= inputPtr->GetSpacing()[1]/2.;
+//      isocenter[2] -= inputPtr->GetSpacing()[2]/2.;
+//      std::cout<<inputPtr->GetSpacing()<<std::endl;
+//    } else {
+//          isocenter[0] += 0.5;
+//          isocenter[1] += 0.5;
+//          isocenter[2] += 1;
+//      }
+    // An Euler 3D transform is used to rotate the volume to simulate the roation of the linac gantry.
+    // The rotation is about z-axis. After the transform, a AP projection geometry (projecting
+    // towards positive y direction) is established.
+    m_GantryRotTransform->SetRotation(0.0, 0.0, -m_ProjectionAngle);
+    m_GantryRotTransform->SetCenter(isocenter);
+    m_ComposedTransform->Compose(m_GantryRotTransform, 0);
+
+    // An Euler 3D transfrom is used to shift the source to the origin.
+    typename TransformType::OutputVectorType focalpointtranslation;
+    focalpointtranslation[0] = -isocenter[0];
+    focalpointtranslation[1] = m_FocalPointToIsocenterDistance - isocenter[1];
+    focalpointtranslation[2] = -isocenter[2];
+    m_CamShiftTransform->SetTranslation(focalpointtranslation);
+    m_ComposedTransform->Compose(m_CamShiftTransform, 0);
+
+    // A Euler 3D transform is used to establish the standard negative z-axis projection geometry. (By
+    // default, the camera is situated at the origin, points down the negative z-axis, and has an up-
+    // vector of (0, 1, 0).)
+
+    m_ComposedTransform->Compose(m_CamRotTransform, 0);
+
+    // The overall inverse transform is computed. The inverse transform will be used by the interpolation
+    // procedure.
+    m_ComposedTransform->GetInverse(m_InverseTransform);
+    this->Modified();
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+void
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Initialize()
+{
+    this->ComputeInverseTransform();
+    m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
+}
+
+} // namespace itk
+
+#endif

itkReg23DRT/vtkContourTopogramProjectionFilter.cxx

@@ -0,0 +1,255 @@
+#include "vtkContourTopogramProjectionFilter.h"
+
+
+#include "vtkCallbackCommand.h"
+#include "vtkCell.h"
+#include "vtkCellArray.h"
+#include "vtkCellData.h"
+#include "vtkContourHelper.h"
+#include "vtkContourValues.h"
+#include "vtkGarbageCollector.h"
+#include "vtkGenericCell.h"
+#include "vtkInformation.h"
+#include "vtkInformationVector.h"
+#include "vtkNew.h"
+#include "vtkObjectFactory.h"
+#include "vtkPointData.h"
+#include "vtkPolyData.h"
+#include "vtkPolyDataNormals.h"
+#include "vtkStreamingDemandDrivenPipeline.h"
+#include "vtkTimerLog.h"
+#include "vtkUniformGrid.h"
+#include "vtkGraph.h"
+
+#include "vtkTransformPolyDataFilter.h"
+#include <cmath>
+
+#include "vtkSmartPointer.h"
+#define VTK_CREATE(type, name) vtkSmartPointer<type> name = vtkSmartPointer<type>::New()
+
+
+vtkObjectFactoryNewMacro(vtkContourTopogramProjectionFilter);
+
+
+vtkContourTopogramProjectionFilter::vtkContourTopogramProjectionFilter()
+{
+
+    this->dProjectionLPSOffset[0]=0.;
+    this->dProjectionLPSOffset[1]=0.;
+    this->dProjectionLPSOffset[2]=0.;
+
+    this->dImportOffsetLPS[0] = 0.;
+    this->dImportOffsetLPS[1] = 0.;
+    this->dImportOffsetLPS[2] = 0.;
+
+    this->dSCD = 0.;
+    this->dPOffset = 0.;
+
+    this->m_RefTransform = nullptr;
+    this->m_Transform = nullptr;
+
+    // by default process active point scalars
+    this->SetInputArrayToProcess(
+                0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
+}
+
+
+vtkContourTopogramProjectionFilter::~vtkContourTopogramProjectionFilter()
+{
+    if (this->m_Transform)
+     {
+       this->m_Transform->UnRegister(this);
+     }
+
+    if (this->m_RefTransform)
+     {
+       this->m_RefTransform->UnRegister(this);
+     }
+}
+
+void vtkContourTopogramProjectionFilter::cleanUpFilter(){
+
+    this->dProjectionLPSOffset[0]=0.;
+    this->dProjectionLPSOffset[1]=0.;
+    this->dProjectionLPSOffset[2]=0.;
+
+    this->dImportOffsetLPS[0] = 0.;
+    this->dImportOffsetLPS[1] = 0.;
+    this->dImportOffsetLPS[2] = 0.;
+
+    this->dSCD = 0.;
+
+    this->m_RefTransform = nullptr;
+    this->m_Transform = nullptr;
+
+}
+
+int vtkContourTopogramProjectionFilter::RequestData(
+        vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
+{
+
+    // Get the info objects.
+    vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
+    vtkInformation* outInfo = outputVector->GetInformationObject(0);
+
+    // Get the input and output.
+    vtkPointSet* psInput = vtkPointSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
+    vtkGraph* graphInput = vtkGraph::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
+    vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
+
+    vtkPoints* points = nullptr;
+    if (psInput)
+    {
+        points = psInput->GetPoints();
+    }
+    else
+    {
+        points = graphInput->GetPoints();
+    }
+
+    // check input
+    // If no points, then nothing to do.
+    if (points == nullptr)
+    {
+//        std::cout << "Cannot Project; no input points" << std::endl;
+        vtkDebugMacro("Cannot Project; no input points");
+        return 1;
+    }
+
+    // If reference transform, then nothing to do.
+    if (m_RefTransform == nullptr)
+    {
+//        std::cout << "Cannot Project; no input reference projection transform" << std::endl;
+        vtkDebugMacro("Cannot Project; no input reference projection transform");
+        return 1;
+    }
+
+    // If transform, then nothing to do.
+    if (m_Transform == nullptr)
+    {
+//        std::cout << "Cannot Project; no input projection transform" << std::endl;
+        vtkDebugMacro("Cannot Project; no input projection transform");
+        return 1;
+    }
+
+    vtkPoints* PrjPoints = vtkPoints::New();
+    PrjPoints->SetNumberOfPoints( points->GetNumberOfPoints() );
+
+    double pp0[3], pp1[3], pp2[3];
+
+    for(vtkIdType ii = 0; ii < points->GetNumberOfPoints() ; ii++){
+
+        /** input points are in LPS, should be corrected for import offset and mapped onto
+         * the projection geometry origin */
+        pp0[0] = points->GetPoint(ii)[0] + dProjectionLPSOffset [0] - dImportOffsetLPS[0];
+        pp0[1] = points->GetPoint(ii)[1] + dProjectionLPSOffset [1] - dImportOffsetLPS[1];
+        pp0[2] = points->GetPoint(ii)[2] + dProjectionLPSOffset [2] - dImportOffsetLPS[2];
+
+        /** LPS points in projection can now be mapped to the
+         * projection geometry using the inverse of the transform */
+        m_Transform->GetInverse()->TransformPoint(pp0, pp1);
+
+        /** correct x coordinate by the ratio of distance from source
+         *  Basic intercept theorem */
+        if(pp1[2] != 0.){
+            pp1[0] = (pp1[0] * dSCD / abs(pp1[2])) /*+ dPOffset*/;
+        } else {
+            vtkDebugMacro("point at SCD distance from panel. Skipping lateral correction for this point.");
+        }
+        /** Go on the virtual detector */
+        pp1[2] = -dSCD;
+
+        /** Points can be mapped back to LPS coordinates */
+        m_RefTransform->TransformPoint(pp1, pp2);
+
+        /** and back to 'real' LPS */
+        pp2[0] -= dProjectionLPSOffset [0];
+        pp2[1] -= dProjectionLPSOffset [1];
+        pp2[2] -= dProjectionLPSOffset [2];
+
+        PrjPoints->InsertPoint(ii,pp2);
+     }
+
+
+    output->SetPoints( PrjPoints );
+
+    vtkIdType numPoints = points->GetNumberOfPoints();
+
+    VTK_CREATE(vtkCellArray, cells);
+    if (cells->AllocateEstimate(numPoints, 1) == true){
+
+    } else {
+        std::cout << "vtkContourTopogramProjectionFilter Can't allocate cells " << std::endl;
+        PrjPoints->Delete();
+        return 1;
+    }
+
+
+    for (vtkIdType i = 0; i < numPoints; i++)
+    {
+        cells->InsertNextCell(1, &i);
+    }
+
+    output->SetVerts(cells);
+
+    PrjPoints->Delete();
+
+    output->Squeeze();
+
+    return 1;
+
+
+}
+
+void vtkContourTopogramProjectionFilter::SetImportOffsetLPS(const double * dP)
+{
+    memcpy(dImportOffsetLPS,dP,3*sizeof(double));
+}
+
+void vtkContourTopogramProjectionFilter::SetSCD(const double dVal){
+    dSCD = dVal;
+}
+
+void vtkContourTopogramProjectionFilter::SetPanelOffset(const double dVal){
+    dPOffset = dVal;
+}
+
+
+void vtkContourTopogramProjectionFilter::SetLPStoProjectionLPSGeometryOffset(const double * dP)
+{
+    memcpy(dProjectionLPSOffset,dP,3*sizeof(double));
+}
+
+void vtkContourTopogramProjectionFilter::SetTransformMatrix( vtkMatrix4x4 *userMatrix)
+{
+    if(this->m_Transform == nullptr){
+        m_Transform = vtkTransform::New();
+    }
+    m_Transform->SetMatrix(userMatrix);
+    this->Modified();
+}
+
+
+void vtkContourTopogramProjectionFilter::SetReferenceTransformMatrix( vtkMatrix4x4 *userMatrix)
+{
+    if(this->m_RefTransform == nullptr){
+        m_RefTransform = vtkTransform::New();
+    }
+    m_RefTransform ->SetMatrix(userMatrix);
+    this->Modified();
+}
+
+
+void vtkContourTopogramProjectionFilter ::PrintSelf(ostream& os, vtkIndent indent)
+{
+    this->Superclass::PrintSelf(os, indent);
+}
+
+
+int vtkContourTopogramProjectionFilter::FillInputPortInformation(int, vtkInformation* info)
+{
+    info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
+    info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkGraph");
+    info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
+    return 1;
+}

itkReg23DRT/vtkContourTopogramProjectionFilter.h

@@ -0,0 +1,86 @@
+/*=========================================================================
+  Program:   DTR
+  Module:    vtkContourTopogramProjectionFilter.h
+  Description:
+    This filter projects polydata points using a topogram transform.
+    Projection geometry is defined by the input transform matrix.
+    Reference transform is used to trasnform points back into LPS and
+    is expect to be the matching projection trasform without correction
+    offsets and rotations.
+
+  gfattori Feb 2022
+=========================================================================*/
+
+#ifndef VTKCONTOURTOPOGRAMPROJECTIONFILTER_H
+#define VTKCONTOURTOPOGRAMPROJECTIONFILTER_H
+
+#include "vtkFiltersCoreModule.h" // For export macro
+#include "vtkPolyDataAlgorithm.h"
+#include "vtkMatrix4x4.h"
+#include "vtkTransform.h"
+
+//#include "itkMacro.h"
+#include "DRTMetaInformation.h"
+
+class VTKFILTERSCORE_EXPORT vtkContourTopogramProjectionFilter :
+        public vtkPolyDataAlgorithm
+{
+
+public:
+    vtkTypeMacro(vtkContourTopogramProjectionFilter, vtkPolyDataAlgorithm);
+    void PrintSelf(ostream& os, vtkIndent indent) override;
+
+
+    static vtkContourTopogramProjectionFilter * New();
+
+    void SetReferenceTransformMatrix( vtkMatrix4x4 *);
+    void SetTransformMatrix( vtkMatrix4x4 *);
+
+    vtkGetObjectMacro(m_Transform, vtkTransform);
+    vtkGetObjectMacro(m_RefTransform, vtkTransform);
+
+    itkSetEnumMacro(ProjectionOrientation, itk::tProjOrientationType);
+    itkGetEnumMacro(ProjectionOrientation, itk::tProjOrientationType);
+
+    itkSetEnumMacro(PatientOrientation, itk::tPatOrientation);
+    itkGetEnumMacro(PatientOrientation, itk::tPatOrientation);
+
+    void SetLPStoProjectionLPSGeometryOffset (const double * );
+    void SetSCD(const double dVal);
+    void SetImportOffsetLPS(const double *);
+    void SetPanelOffset(const double dVal);
+
+    void cleanUpFilter();
+
+protected:
+     vtkContourTopogramProjectionFilter();
+     ~vtkContourTopogramProjectionFilter() override;
+
+     int RequestData(vtkInformation* request, vtkInformationVector** inputVector,
+            vtkInformationVector* outputVector) override;
+     int FillInputPortInformation(int port, vtkInformation* info) override;
+
+     itk::tProjOrientationType
+     m_ProjectionOrientation;
+
+     itk::tPatOrientation
+     m_PatientOrientation;
+
+private:
+     vtkContourTopogramProjectionFilter (const vtkContourTopogramProjectionFilter &) = delete;
+     void operator=(const vtkContourTopogramProjectionFilter&) = delete;
+
+
+     double
+        dProjectionLPSOffset[3],
+        dImportOffsetLPS[3],
+        dSCD,
+        dPOffset;
+
+     vtkTransform *m_RefTransform;
+     vtkTransform *m_Transform;
+
+};
+
+
+#endif