Merge branch 'cleanUpUIHandler' into 'ScoutRT_Qt6_handler'

working on Handler in qt6 See merge request cpt_bioeng/drt!44
itk metric wrong include
2024-02-22 16:40:26 +01:00 · 2024-02-06 17:06:33 +01:00 · 2024-02-06 16:56:27 +01:00
30 changed files with 733 additions and 506 deletions
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -1,20 +0,0 @@
-cff-version: 1.2.0
-message: "If you use this software, please cite it as below."
-title: "reg23Topograms"
-version: "1.0"
-doi: "10.5281/zenodo.15399663"
-date-released: 2025-05-13
-authors:
-  - family-names: Fattori
-    given-names: Giovanni
-    orcid: "https://orcid.org/0000-0002-7639-7759"
-  - family-names: Via
-    given-names: Riccardo
-    orcid: "https://orcid.org/0000-0003-3310-4802"
-
-repository-code: "https://gitea.psi.ch/fattori_g/reg23Topograms"
-license: "BSD-3-Clause"
-keywords:
-  - 2D/3D registration
-  - Digitally Reconstructed Topograms (DRT)
-  - CT topogram imaging
--- a/LICENSE.txt
+++ b/LICENSE.txt
@@ -1,29 +0,0 @@
-BSD 3-Clause License
-
-Copyright (c) 2025, Giovanni Fattori
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-1. Redistributions of source code must retain the above copyright notice, this
-   list of conditions and the following disclaimer.
-
-2. Redistributions in binary form must reproduce the above copyright notice,
-   this list of conditions and the following disclaimer in the documentation
-   and/or other materials provided with the distribution.
-
-3. Neither the name of the copyright holder nor the names of its
-   contributors may be used to endorse or promote products derived from
-   this software without specific prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
-ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
-WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--- a/README.md
+++ b/README.md
@@ -1,47 +0,0 @@
-# 2D/3D Registration with Digitally Reconstructed Topograms
-
-This project provides a framework for 2D/3D registration using Digitally Reconstructed Topogram (DRT) images, which are generated from 3D CT volumes using a modified version of the Siddon-Jacobs ray casting algorithm.
-
-In addition to DRT generation, the framework extends [ITKTwoProjectionRegistration](https://github.com/InsightSoftwareConsortium/ITKTwoProjectionRegistration) in several key areas, including:
-
- Compliance with IEC 61217  
- Revised and consistent reference coordinate systems  
- This project provides a framework for 2D/3D registration using Digitally Reconstructed Topograms (DRTs), which are generated from 3D CT volumes using a modified version of the Siddon-Jacobs RayCast algorithm. See the [`siddonTril` branch](https://gitea.psi.ch/fattori_g/reg23Topograms/src/branch/siddonTril) for implementation details.
-
-## Requirements
-
-Work in progress — requirements and example code to be added.
-
-## Citation
-
-If you use this code or method in your research, please cite the following paper:
-
-> **Giovanni Fattori, Riccardo Via, et al.**  
-> *Topogram and 3DCT geometry calibration for image-guided proton therapy with in-room CT-on-rails*  
-> *Physics and Imaging in Radiation Oncology*, 2025.  
-> [https://doi.org/10.1016/j.phro.2025.100799](https://doi.org/10.1016/j.phro.2025.100799)
-
-<details>
-<summary>BibTeX</summary>
-
-```bibtex
-@article{Fattori2025,
-  author  = {Giovanni Fattori and Riccardo Via and others},
-  title   = {Topogram and 3DCT geometry calibration for image-guided proton therapy with in-room CT-on-rails},
-  journal = {Physics and Imaging in Radiation Oncology},
-  year    = {2025},
-  doi     = {10.1016/j.phro.2025.100799}
-}
-```
-
-</details>
-
-## Requirements
-
-Work in progress — requirements and example code to be added.
-
-## Contributors
-
- Giovanni Fattori  
- Riccardo Via  
- Reto Ansorge
--- a/reg23Topograms/itkDTRrecon/CMakeLists.txt
+++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt
@@ -1,6 +1,6 @@
-cmake_minimum_required(VERSION 2.8.12)
+cmake_minimum_required(VERSION 3.16)

-SET(LIB_NAME itkReg23DRT)
+SET(LIB_NAME itkDTRrecon)

 add_subdirectory(autoreg)

--- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
+++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
@@ -11,7 +11,6 @@ static double FFS2IEC[9] = {
    0, 0, -1,
    0, -1, 0};

-    
 namespace itk
 {

@@ -32,6 +31,7 @@ TopogramImageMetaInformation(){

 }

+
 void
 TopogramImageMetaInformation
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
@@ -39,12 +39,14 @@ TopogramImageMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 TopogramImageMetaInformation
 ::~TopogramImageMetaInformation ()
 {

 }

+
 void TopogramImageMetaInformation ::
 SetPatientOrientation(tPatOrientation m_orient)
 {
@@ -69,6 +71,7 @@ SetPatientOrientation(tPatOrientation m_orient)

 }

+
 DRTImageMetaInformation::
 DRTImageMetaInformation(){

@@ -109,6 +112,7 @@ DRTImageMetaInformation

 }

+
 DRTImageMetaInformation::PointType
 DRTImageMetaInformation::GetOrigin()
 {
@@ -128,12 +132,67 @@ DRTImageMetaInformation::GetOrigin()
            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;
@@ -143,6 +202,7 @@ void DRTImageMetaInformation::SetSizeWithBounds(
    return;
 }

+
 DRTImageMetaInformation::PointType
 DRTImageMetaInformation::GetLPStoProjectionGeoLPSOffset()
 {
@@ -203,6 +263,10 @@ CTVolumeImageMetaInformation::GetProjectionOriginLPS(
        PointType ProjectionCenter)
 {

+    //    PointType Origin;
+    //    Origin = this->m_OriginLPS;
+    //    Origin = Origin - this->m_ImportOffset;
+
    DirectionType IECtoLPS_Directions;
    IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose();

@@ -296,6 +360,7 @@ CTVolumeImageMetaInformation(){
    this->m_ImportOffset.Fill(0.);
 }

+
 void
 CTVolumeImageMetaInformation
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
@@ -303,6 +368,7 @@ CTVolumeImageMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 CTVolumeImageMetaInformation
 ::~CTVolumeImageMetaInformation   ()
 {
@@ -366,6 +432,7 @@ CTVolumeImageMetaInformation::GetOriginWOffset(){
            this->m_OriginLPS - this->m_ImportOffset;
 }

+
 DRTProjectionGeometryImageMetaInformation::
 DRTProjectionGeometryImageMetaInformation(){

@@ -420,12 +487,16 @@ DRTProjectionGeometryImageMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 DRTProjectionGeometryImageMetaInformation
 ::~DRTProjectionGeometryImageMetaInformation    ()
 {

 }

+
+
+
 RTGeometryMetaInformation::
 RTGeometryMetaInformation(){

@@ -443,6 +514,7 @@ RTGeometryMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 RTGeometryMetaInformation
 ::~RTGeometryMetaInformation()
 {
@@ -472,6 +544,7 @@ InternalTransformMetaInformation

 }

+
 R23MetaInformation::
 R23MetaInformation (){

@@ -480,6 +553,7 @@ R23MetaInformation (){
    m_MetricType = tMetricTypeEnum::NCC;
 }

+
 void
 R23MetaInformation
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
@@ -487,6 +561,7 @@ R23MetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 R23MetaInformation
 ::~R23MetaInformation  ()
 {
@@ -557,4 +632,11 @@ TransformMetaInformation
 }


+
+
+
+
+
+
+
 }
--- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h
+++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h
@@ -24,6 +24,7 @@ typedef enum eImageOrientationType{
    FFS = 2
 } tPatOrientation;

+
 typedef enum eHandlingRotImpTransform {
    ALWAYS_USE = 0,
    NEVER_USE =1,
@@ -41,6 +42,7 @@ typedef enum eDegreeOfFreedomType {
    OTHER =7
 } tDegreeOfFreedomEnum;

+
 typedef enum eOptimizerType{
    POWELL = 0,
    AMOEBA = 1,
@@ -52,6 +54,7 @@ typedef enum eMetricType{
    MI = 1
 } tMetricTypeEnum;

+
 inline int GetNumberOfDOF(eDegreeOfFreedomType dof)
 {
    switch (dof) {
@@ -73,6 +76,7 @@ inline int GetNumberOfDOF(eDegreeOfFreedomType dof)
    return 0;
 };

+
 class DegreeOfFreedom : public itk::Object {

 public:
@@ -130,6 +134,7 @@ public:

 protected:

+
    tPatOrientation
    m_PatientOrientation;

@@ -144,6 +149,7 @@ protected:
    DirectionType
    m_LPS2IECDirections;

+
    /** Default Constructor **/
    TopogramImageMetaInformation  ();
    /** Default Destructor **/
@@ -187,6 +193,9 @@ public:
    itkSetMacro(Spacing,SpacingType);
    itkGetMacro(Spacing,SpacingType);

+    //itkSetMacro(Origin,PointType);
+    //itkGetMacro(Origin,PointType);
+
    itkSetMacro(OriginLPS,PointType);
    itkGetMacro(OriginLPS,PointType);

@@ -636,6 +645,9 @@ public:
     itkSetEnumMacro(MetricType, tMetricTypeEnum);
     itkGetEnumMacro(MetricType, tMetricTypeEnum);

+
+
+
 protected:

    tDegreeOfFreedomEnum
@@ -681,6 +693,7 @@ public:
    /** object information streaming **/
    void PrintSelf(std::ostream& os, itk::Indent indent) const;

+
    itkSetMacro(HiddenTranslations,PointType);
    itkGetMacro(HiddenTranslations,PointType);

--- a/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt
+++ b/reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt
--- a/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h
@@ -6,12 +6,6 @@
 #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])
@@ -30,7 +24,7 @@ const char* computeMethodName(const char (&function)[FL], const char (&prettyFun

 namespace itk {
 /* reference string required for comparison with tag values */
-static const char VARIABLE_IS_NOT_USED *ImageOrientationStrings[] = {
+static const char* ImageOrientationStrings[] = {
    "NotDefined",
    "HFS",
    "FFS",
@@ -46,49 +40,49 @@ 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] = {
+static double Standard_DRT2LPS[9] = {
    1, 0, 0,
    0, 0, -1,
    0, 1, 0
 };

-static double VARIABLE_IS_NOT_USED PAElementsIEC[9] = {
+static double PAElementsIEC[9] = {
    1, 0, 0,
    0, -1, 0,
    0, 0, -1
 };

-static double VARIABLE_IS_NOT_USED LATElementsIEC[9] = {
+static double LATElementsIEC[9] = {
    0, 0, -1,
    0, -1, 0,
    -1, 0, 0
 };

-static double VARIABLE_IS_NOT_USED HFS2IEC[9] = {
+static double HFS2IEC[9] = {
    1, 0, 0,
    0, 0, 1,
    0, -1, 0
 };

-static double VARIABLE_IS_NOT_USED FFS2IEC[9] = {
+static double FFS2IEC[9] = {
    -1, 0, 0,
    0, 0, -1,
    0, -1, 0
 };

-static double VARIABLE_IS_NOT_USED HFP2IEC[9] = {
+static double HFP2IEC[9] = {
    -1, 0, 0,
    0, 0, 1,
    0, 1, 0
 };

-static double VARIABLE_IS_NOT_USED FFP2IEC[9] = {
+static double FFP2IEC[9] = {
    1, 0, 0,
    0, 0, -1,
    0, 1, 0
 };

-static double VARIABLE_IS_NOT_USED PAT2IEC[9] = {
+static double PAT2IEC[9] = {
    1, 0, 0,
    0, 0, 1,
    0, -1, 0
--- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h
--- a/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx
@@ -143,7 +143,7 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
    //std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl;

    auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion();
-    //unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
+    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.
@@ -180,7 +180,7 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
    //std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl;

    movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion();
-    //numberOfPixels = movingImageRegion.GetNumberOfPixels();
+    numberOfPixels = movingImageRegion.GetNumberOfPixels();

    //numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
    //metric2->SetNumberOfSpatialSamples(numberOfSamples);
--- a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h
@@ -20,7 +20,9 @@

 #include "itkCovariantVector.h"
 #include "itkMacro.h"
-#include "itkNormalizedCorrelationImageToImageMetric.hxx"
+//#include "itkNormalizedCorrelationImageToImageMetric.hxx"
+#include "itkNormalizedCorrelationImageToImageMetric.h"
+
 #include "itkPoint.h"
 #include "itkgTwoImageToOneImageMetric.h"

--- a/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx
@@ -262,7 +262,7 @@ NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetVal
    }

    // Calculate the measure value between fixed image 1 and the moving image
-    //auto oldprecision = std::cout.precision();
+    auto oldprecision = std::cout.precision();


 //    std::cout<<"Region "  <<this->GetFixedImageRegion1() <<std::endl;
--- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h
--- a/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx
--- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h
--- a/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx
+++ b/reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx
--- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
@@ -119,7 +119,7 @@ itkImageProcessor::itkImageProcessor()
    Point3D[1]=0.;
    Point3D[2]=175.;
    m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D);
-
+    //    m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D);
    ImageType3D::SizeType
            ImageSize;
    ImageSize[0]=512;
@@ -164,6 +164,7 @@ void itkImageProcessor::SetNumberOfWorkingUnits(int iN){
    }
 }

+
 const CTVolumeImageMetaInformation::Pointer
 itkImageProcessor::GetCTMetaInfo(
        ){
@@ -187,6 +188,7 @@ itkImageProcessor::GetDRTGeoMetaInfo(){

 }

+
 void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
 {
    Superclass::PrintSelf(os, indent);
@@ -219,6 +221,7 @@ double itkImageProcessor::GetSCD2(){
            m_DRTImage2MetaInfo ->GetSCD();
 }

+
 double itkImageProcessor::GetPanelOffsetPGeo(tProjOrientationType ImgPrj){

    if(m_CTMetaInfo == NULL ||
@@ -266,6 +269,8 @@ double itkImageProcessor::GetPanelOffset2(){
            m_DRTGeometryMetaInfo ->GetPanel2Offset();
 }

+
+
 tDegreeOfFreedomEnum
 itkImageProcessor::GetDegreeOfFreedom()
 {
@@ -283,6 +288,7 @@ void itkImageProcessor::SetUseRotationsForImportOffset(bool bVal){
    return;
 }

+
 void itkImageProcessor::SetCustom_ImportTransform(double dTx,
                                                  double dTy,
                                                  double dTz,
@@ -337,6 +343,7 @@ void itkImageProcessor::SetCustom_ProjectionCenterOffsetFixedAxes_IEC(

 }

+
 void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(
        double dX1,
        double dY1,
@@ -371,7 +378,7 @@ void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double
    Spacing [1] = nY2;
    Spacing [2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT2Spacing(Spacing);
-
+    //TODO UPDATE TO FOLLOW
 }

 void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
@@ -392,6 +399,7 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
    Size [2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT2Size(Size);

+    //TODO UPDATE TO FOLLOW
 }

 void itkImageProcessor::SetCustom_UpdateMetaInfo(){
@@ -416,6 +424,7 @@ int itkImageProcessor::unload3DVolumeAndMeta(){
    m_TransformMetaInfo = NULL;
    m_TransformMetaInfo = TransformMetaInformation::New();

+
    this->ResetROIRegions();

    return 1;
@@ -427,8 +436,8 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
    tPatOrientation m_PatOrientation
            = tPatOrientation::NotDefined;

-    /* Since are not sure what we get here,
-     * we run the IPP sorter. */
+    /* Since are not sure DB indexer to sort images,
+     * we run the IPP sorter here. */
    using IppSorterType = gdcm::IPPSorter;
    IppSorterType ZSorter;

@@ -441,6 +450,12 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
        return EXIT_FAILURE;
    }

+    /* Some debugging output */
+    //std::cout << "Sorting succeeded:" << std::endl;
+    //ZSorter.Print( std::cout );
+    //std::cout << "Found z-spacing:" << std::endl;
+    //std::cout << ZSorter.GetZSpacing() << std::endl;
+
    const std::vector<std::string> & v_sortedFnames =
            ZSorter.GetFilenames();

@@ -457,6 +472,7 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
        imageSeriesReader3D->SetFileNames(v_sortedFnames);
        imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits);
        imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt
+        //imageSeriesReader3D->SetSpacingWarningRelThreshold();

        imageSeriesReader3D->Update();

@@ -489,6 +505,8 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
            return -1;
        }

+        ////////////////////////*****************//////////////////////////
+
        CastFilterType3D::Pointer caster3D =
                CastFilterType3D::New();

@@ -514,13 +532,22 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
            const typename ImageType3D::SizeType& inputSize =
                    inputRegion.GetSize();

+            // Inputs Info
+    //        std::cout << "inputSpacing " << std::endl;
+    //        std::cout << inputSpacing << std::endl;
+    //        std::cout << "inputRegion " << std::endl;
+    //        std::cout << inputRegion << std::endl;
+    //        std::cout << "inputSize " << std::endl;
+    //        std::cout << inputSize  << std::endl;
+
+
            ImageType3D::SpacingType outputSpacing;
            outputSpacing[0] = inputSpacing[0];
            outputSpacing[1] = inputSpacing[1];
            outputSpacing[2] = m_SamplingLNG;

            ImageType3D::SizeType   outputSize;
-
+    //        typedef typename ImageType3D::SizeType::SizeValueType SizeValueType;
            outputSize[0] = static_cast<SizeValueType>(
                        inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5);
            outputSize[1] = static_cast<SizeValueType>(
@@ -549,9 +576,17 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){

        }

+        ////////////////////////*****************//////////////////////////
+
        if (m_VolumeSourceDupli == NULL)
            std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl;

+        //        if (m_VolumeSourceDupli) {
+        //            std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl;
+        //            m_VolumeSourceDupli = NULL;
+        //            m_VolumeSourceDupli = DuplicatorType::New();
+        //        }
+
        m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput());
        m_VolumeSourceDupli->Update();

@@ -598,6 +633,7 @@ int itkImageProcessor::fill3DVolumeMeta(
        m_RTMetaInfo = NULL;
    }

+
    /* copy useful meta information into the CT container */
    m_CTMetaInfo = CTVolumeImageMetaInformation::New();
    m_CTMetaInfo->SetPatientOrientation(m_PatOrientation);
@@ -612,6 +648,8 @@ int itkImageProcessor::fill3DVolumeMeta(
    PointOffset.Fill(0.);
    m_CTMetaInfo->SetImportOffset(PointOffset);

+
+
    /* initialise DRT meta */
    m_DRTImage1MetaInfo = DRTImageMetaInformation::New();
    m_DRTImage1MetaInfo->SetProjectionAngleLPS(
@@ -626,6 +664,7 @@ int itkImageProcessor::fill3DVolumeMeta(
                m_DRTGeometryMetaInfo->GetSCD1Offset()
                );

+
    /* Calculate projection angle IEC to LPS */
    m_DRTImage2MetaInfo = DRTImageMetaInformation::New();
    m_DRTImage2MetaInfo->SetProjectionAngleLPS(
@@ -640,12 +679,35 @@ int itkImageProcessor::fill3DVolumeMeta(
                m_DRTGeometryMetaInfo->GetSCD2Offset()
                );

+    //    std::cout<< "///////////////// 3D VOLUME BEG ///////////////" <<std::endl;
+
+    //    std::cout<< "OriginLPS"<< m_CTMetaInfo->GetOriginLPS() <<std::endl;
+    //    std::cout<< "COV"<< m_CTMetaInfo->GetCOV() <<std::endl;
    std::cout<< "ImportOffset"<< m_CTMetaInfo->GetImportOffset() <<std::endl;
+    //    std::cout<< "LPS2IEC"<< m_CTMetaInfo->GetLPS2IECDirections() <<std::endl;
+    //    std::cout<< "Size"<< m_CTMetaInfo->GetSize() <<std::endl;
+    //    std::cout<< "Spacing"<< m_CTMetaInfo->GetSpacing() <<std::endl;
+
+    //    std::cout<< "pAngleLPS1"<< m_DRTImage1MetaInfo->GetProjectionAngleLPS() <<std::endl;
+    //    std::cout<< "SCD1"<< m_DRTImage1MetaInfo->GetSCD() <<std::endl;
+
+    //    std::cout<< "pAngleLPS2"<< m_DRTImage2MetaInfo->GetProjectionAngleLPS() <<std::endl;
+    //    std::cout<< "SCD2"<< m_DRTImage2MetaInfo->GetSCD() <<std::endl;
+
+    //    std::cout<< "///////////////// 3D VOLUME END ///////////////" <<std::endl;

    this->UpdateProjectionGeometryMeta();

+    // To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image
+    // to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
+    // geometry is fully defined. The origin of the CT image becomes irrelavent.
+
+    // this is done internally in the dSiddonFilter.
+
+
    return EXIT_SUCCESS;

+
 }

 const std::vector <double>  itkImageProcessor::GetRTImportOffset()
@@ -715,12 +777,15 @@ itkImageProcessor::CalcProjectionAngleLPS(
        tPatOrientation  pOrient,
        double pAngleIEC){

+    //    std::cout<< "CalcProjectèionAngleLPS :: pAngleIEC " << pAngleIEC <<std::endl;
+
    double dProj = pAngleIEC;

    InternalImageType::DirectionType
            currIECtoLPS;

-    /* NOTE that we transpose on the fly... */
+
+    /* NOTE WE TRANSPOSE ON THE FLY... */
    switch (pOrient) {
    case tPatOrientation :: HFS:
        for(int iIdx = 0 ; iIdx < 3; iIdx++){
@@ -779,6 +844,7 @@ void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2)
    m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2);
 }

+
 void itkImageProcessor::SetProjectionAngle1IEC(double dGantryAngle)
 {
    m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(dGantryAngle);
@@ -886,6 +952,8 @@ int itkImageProcessor::load2D(const char * pcFName){
    }
    tokens.clear();

+    //    std::cout<<"Intensity window center and width: "<< dIntensityWindow[0] <<" "<< dIntensityWindow[1] <<std::endl;
+
    /* Image orientation */
    eImageOrientationType
            currImgOrient;
@@ -893,8 +961,10 @@ int itkImageProcessor::load2D(const char * pcFName){
    *std::remove(sTmpString, sTmpString + strlen(sTmpString), ' ') = 0;

    if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
+        //        std::cout<<"Image orientation: HFS"<<std::endl;
        currImgOrient = eImageOrientationType::HFS;
    } else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
+        //        std::cout<<"Image orientation: FFS"<<std::endl;
        currImgOrient = eImageOrientationType::FFS;
    } else {
        std::cerr<< "Image Orientation: unrecognised"<< sTmpString <<std::endl;
@@ -926,6 +996,9 @@ int itkImageProcessor::load2D(const char * pcFName){
    ImageDirectionType3D LocalizerImagDirectionDCM =
            imageReader2D->GetOutput()->GetDirection();

+//    std::cout<<"LocalizerImagDirectionDCM " <<LocalizerImagDirectionDCM <<std::endl;
+
+
    InternalImageType::DirectionType toIECDirection;
    for(int iIdx = 0 ; iIdx < 3; iIdx++){
        for(int jIdx = 0 ; jIdx < 3; jIdx++){
@@ -942,10 +1015,17 @@ int itkImageProcessor::load2D(const char * pcFName){
        }
    }

+//    std::cout<<"toIECDirection " <<toIECDirection <<std::endl;
+
    /* calculate image orientation with respect to fixed IEC */
    InternalImageType::DirectionType LocalizerImagDirectionIEC =
            toIECDirection * LocalizerImagDirectionDCM;

+//    std::cout<<"toIECDirection " <<toIECDirection <<std::endl;
+//    std::cout<<"LocalizerImagDirectionDCM " <<LocalizerImagDirectionDCM <<std::endl;
+//    std::cout<<"LocalizerImagDirectionIEC " <<LocalizerImagDirectionIEC <<std::endl;
+
+
    /* we calculate dot products between the Z components */
    double dSumPA = 0;
    dSumPA =
@@ -959,6 +1039,9 @@ int itkImageProcessor::load2D(const char * pcFName){
            (LocalizerImagDirectionIEC[1][2] * LATElementsIEC[5]) +
            (LocalizerImagDirectionIEC[2][2] * LATElementsIEC[8]);

+//    std::cout<<"dSumPA " <<dSumPA <<std::endl;
+//    std::cout<<"dSumLAT " <<dSumLAT <<std::endl;
+
    tProjOrientationType
            currProjOrientation = NA;

@@ -982,6 +1065,7 @@ int itkImageProcessor::load2D(const char * pcFName){
    InternalImageType::Pointer MyLocalizerImage =
            caster2DInput->GetOutput();

+    //double* m_ImageIntensity;
    DuplicatorType::Pointer m_Duplicator;
    TopogramImageMetaInformation::Pointer m_TImageMeta;

@@ -993,8 +1077,10 @@ int itkImageProcessor::load2D(const char * pcFName){
            m_PASourceDupli = DuplicatorType::New();
        }

+        //m_ImageIntensity = d_image1IntensityWindow;
        m_Duplicator = m_PASourceDupli;
        m_TImageMeta = m_TImage1MetaInfo;
+        //        m_ImLoaded = &image2D1Loaded;

        break;

@@ -1005,7 +1091,9 @@ int itkImageProcessor::load2D(const char * pcFName){
            m_LATSourceDupli = DuplicatorType::New();
        }

+        //m_ImageIntensity = d_image2IntensityWindow;
        m_Duplicator = m_LATSourceDupli;
+        //        m_ImLoaded = &image2D2Loaded;
        m_TImageMeta = m_TImage2MetaInfo;

        break;
@@ -1031,11 +1119,28 @@ int itkImageProcessor::load2D(const char * pcFName){
    m_Duplicator->SetInputImage(caster2DInput->GetOutput() );
    m_Duplicator->Update();

+//    InternalImageType::IndexType pIndex;
+//    pIndex[0]=200;
+//    pIndex[1]=300;
+//    InternalImageType::PixelType pPixel =
+//            m_Duplicator->GetOutput()->GetPixel(pIndex);
+    //std::cout<<"processro PIXEL: "<<pPixel <<std::endl;
+    //    std::cout<< " ////////////// 2D Topo BEG " <<std::endl;
+
+    //    std::cout<<"GetDirection()"<<m_Duplicator->GetOutput()->GetDirection()<<std::endl;
+    //    std::cout<<"GetOrigin()"<<m_Duplicator->GetOutput()->GetOrigin()<<std::endl;
+    //    std::cout<<"GetSpacing()"<<m_Duplicator->GetOutput()->GetSpacing()<<std::endl;
+    //    std::cout<<"GetSize()"<<m_Duplicator->GetOutput()->GetLargestPossibleRegion().GetSize()<<std::endl;
+
+    //    std::cout<< " ////////////// 2D Topo END " <<std::endl;
+
    return
            (int) currProjOrientation;

 }

+
+
 double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
 {

@@ -1100,6 +1205,10 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
    return 0.;
 }

+
+
+
+
 Optimizer::ParametersType
 itkImageProcessor::GetFinalR23Parameters(){

@@ -1189,7 +1298,6 @@ void itkImageProcessor::SetOptimizer(std::string optimizer)


 }
-
 void itkImageProcessor::SetMetric(std::string metric)
 {

@@ -1217,6 +1325,8 @@ void itkImageProcessor::SetHandleRotationImportOffset(eHandlingRotImpTransform h

 }

+
+
 void itkImageProcessor::SetPowellOptimParameters(
        double dStepT,
        double dValTol,
@@ -1254,6 +1364,7 @@ void itkImageProcessor::SetNCCMetricParameters(double dMaxT,bool bSm){
    m_NCCMetricMetaInfo->SetSubtractMean(bSm);
 }

+
 void itkImageProcessor::InitializeProjector()
 {

@@ -1393,6 +1504,7 @@ void itkImageProcessor::InitializeProjector()
    filter2 =
            ChangeInformationFilterType::New();

+
    /* First of all we need the center of the Projection images in its reference frame */
    resampleFilter1->Update();
    filter1->SetInput( resampleFilter1->GetOutput() );
@@ -1402,6 +1514,8 @@ void itkImageProcessor::InitializeProjector()
    filter1->ChangeOriginOn();
    filter1->UpdateOutputInformation();

+    //    std::cout<< "itkImageProcessor::InitializeProjector() " <<std::endl;
+
    /* Again.. */
    resampleFilter2->Update();
    filter2->SetInput( resampleFilter2 ->GetOutput() );
@@ -1419,6 +1533,8 @@ void itkImageProcessor::InitializeProjector()

 }

+
+
 itkImageProcessor::InternalImageType::DirectionType
 itkImageProcessor::CalcDRTImageDirections(
        double angle,
@@ -1455,6 +1571,7 @@ int itkImageProcessor::unloadRTPlanAndMeta(){
    return 0;
 }

+
 void itkImageProcessor::loadRTPlanInfo(
        double dIsoX, double dIsoY, double dIsoZ,
        double dLAT, double dVRT ,double dLNG){
@@ -1490,6 +1607,7 @@ void itkImageProcessor::loadRTPlanInfo(

    m_RTMetaInfo->SetIsocenterIECS(Point);

+
    ImageType3D::PointType
            pZero (3);
    pZero.Fill(0.);
@@ -1538,6 +1656,7 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
        //TODO.
    }

+
    // FIRST

    ImageType3D::PointType NominalIsocenter_wZcorrectionLPS;
@@ -1560,6 +1679,12 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
    IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint(
                m_DRTGeometryMetaInfo->GetProjectionCenterOffset1());

+    //    std::cout<< "///////////////// PGEOM META BEG ///////////////" <<std::endl;
+//        std::cout<<"NominalIsocenter IEC"<< m_DRTGeometryMetaInfo->GetProjectionCenter() <<std::endl;
+//        std::cout<<"CALIBRATION NominalIsocenter_wZcorrectionLPS "<<NominalIsocenter_wZcorrectionLPS<<std::endl;
+//        std::cout<<"CALIBRATION NominalIsocenterZero_wZcorrectionLPS  "<<NominalIsocenterZero_wZcorrectionLPS<<std::endl;
+//        std::cout<<"CALIBRATION IsocenterOffsetLPS  "<<IsocenterOffsetLPS<<std::endl;
+
    ImageType3D::PointType CalibratedIsocenterLPS;

    CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] +
@@ -1569,6 +1694,7 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
    CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2];

+//        std::cout<<"CALIBRATION CalibratedIsocenterLPS  "<<CalibratedIsocenterLPS<<std::endl;

    ImageType3D::PointType CalibratedIsocenterZeroLPS;

@@ -1578,6 +1704,7 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
            IsocenterOffsetLPS[1] /*+ m_CTMetaInfo->GetSpacing()[1]/2*/;
    CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2] /*+ m_CTMetaInfo->GetSpacing()[2]/2*/;
+//        std::cout<<"CALIBRATION CalibratedIsocenterZeroLPS  "<<CalibratedIsocenterZeroLPS<<std::endl;

    m_DRTImage1MetaInfo->SetProjectionAngleLPS(
                this->CalcProjectionAngleLPS(
@@ -1750,10 +1877,13 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){

 }

+
 /*
    Nobody should go through all this...
    To be very careful editing ...
+
 */
+
 itkImageProcessor::ImageType3D::PointType
 itkImageProcessor::CalcDRTImageOrigin(
        ImageType3D::PointType m_DRTOrigin,
@@ -1773,8 +1903,11 @@ itkImageProcessor::CalcDRTImageOrigin(
            NewOrigin;
 }

+
+
 void itkImageProcessor::GetProjectionImages(){

+    //    std::cout<< "itkImageProcessor::GetProjectionImages" <<std::endl;
    if(m_DRTImage1MetaInfo == NULL ||
            m_DRTImage2MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
@@ -1784,9 +1917,12 @@ void itkImageProcessor::GetProjectionImages(){
        return;
    }

+  //  const double dtr = (atan(1.0) * 4.0) / 180.0;
+
    ImageType3D::PointType ZeroPoint;
    ZeroPoint.Fill(0.);

+
    transform1->SetComputeZYX(true);
    transform1->SetIdentity();

@@ -1812,6 +1948,18 @@ void itkImageProcessor::GetProjectionImages(){
    transform1 ->SetCenter(
                m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());

+//    std::cout<< "CurrTransform Rotations: "<<
+//                CurrTransform->GetAngleX() <<" "<<
+//                CurrTransform->GetAngleY() <<" "<<
+//                CurrTransform->GetAngleZ() <<" "<< std::endl;
+
+//    std::cout<< "CurrTransform Translations: "<<
+//                CurrTransform->GetTranslation()<<std::endl;
+//    std::cout<< "CurrTransform center: "<<
+//                CurrTransform->GetCenter()<<std::endl;
+
+
+    //transform1 ->Print(std::cout);
    // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
    // have been set before registration. Here we only need to replace the initial
    // transform with the final transform.
@@ -1895,6 +2043,8 @@ itkImageProcessor::CalcImportVolumeOffset(
    ImageType3D::PointType IsoSupport_IECPos =
            InputTransform->TransformPoint(rtCouchOffset );

+
+
    InternalImageType::DirectionType VolumeIECtoLPS;
    VolumeIECtoLPS =  VolumeLPStoIEC.GetTranspose();

@@ -1915,6 +2065,7 @@ itkImageProcessor::CalcImportVolumeOffset(
 }


+
 void itkImageProcessor::Write2DImages(){

    //    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
@@ -1976,18 +2127,59 @@ void itkImageProcessor::Write2DImages(){

 vtkImageData* itkImageProcessor::GetLocalizer1VTK()
 {
-
+    // Rescale the intensity of the projection images to 0-32768 for output.
    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
    rescaler1->SetOutputMinimum(0);
-    rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
-    rescaler1->SetInput( m_PASourceDupli->GetOutput() );
+    rescaler1->SetOutputMaximum(32768);
+    rescaler1->SetInput(m_PASourceDupli->GetOutput());
+
    rescaler1->Update();

    toVTKLocalizer1->SetInput(rescaler1->GetOutput());
    toVTKLocalizer1->Update();

+//    if(false) {
+////        using ImageRegionType3D = ImageType3D::RegionType;
+////        using SizeType3D = ImageRegionType3D::SizeType;
+////        using ImageDirectionType3D = ImageType3D::DirectionType;
+
+////        ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion();
+////        SizeType3D        size3D = region3D.GetSize();
+////        //ImageType3D::PointType       origin3D = rescaler1->GetOutput()->GetOrigin();
+////        const itk::Vector<double, 3> resolution3D = rescaler1->GetOutput()->GetSpacing();
+////        ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection();
+
+////        /* calculate image size in 3D Space */
+////        using VectorType = itk::Vector<double, 3>;
+////        VectorType Dest;
+////        Dest[0]=(size3D[0]-1) * resolution3D[0];
+////        Dest[1]=(size3D[1]-1) * resolution3D[1];
+////        Dest[2]=(size3D[2]-1) * resolution3D[2];
+
+
+//        // ImageType3D::PointType LastVoxelPosition =
+//        //         origin3D + (imagDirection * Dest);
+
+//        //        std::cout<<" -------- Localizer 1 ITK --------"<<std::endl;
+//        //        std::cout<<"size3D " <<size3D <<std::endl;
+//        //        std::cout<<"resolution3D " <<resolution3D <<std::endl;
+//        //        std::cout<<"imagDirection " <<imagDirection <<std::endl;
+//        //        std::cout<<"First Voxel position " <<origin3D <<std::endl;
+//        //        std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
+
+
+
+//        // double* dBounds = toVTKLocalizer1->GetOutput()->GetBounds();
+
+//        //        std::cout<< "-------- Localizer 1 VTK --------" <<std::endl;
+//        //        std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
+//        //                <<dBounds[2]<<" "<<dBounds[3]<<" "
+//        //               <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
+//        //        std::cout<< "-------- -------- --------" <<std::endl;
+//    }
+
    return
            toVTKLocalizer1->GetOutput();
 }
@@ -1995,22 +2187,65 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK()
 vtkImageData* itkImageProcessor::GetLocalizer2VTK()
 {

+
+    // Rescale the intensity of the projection images to 0-32768 for output.
    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
    rescaler2->SetOutputMinimum(0);
-    rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
-    rescaler2->SetInput( m_LATSourceDupli->GetOutput() );
+    rescaler2->SetOutputMaximum(32768);
+    rescaler2->SetInput(m_LATSourceDupli ->GetOutput());
+
    rescaler2->Update();

    toVTKLocalizer2->SetInput(rescaler2->GetOutput());
    toVTKLocalizer2->Update();

+
+
+//    if(true) {
+//        using ImageRegionType3D = ImageType3D::RegionType;
+//        using SizeType3D = ImageRegionType3D::SizeType;
+//        using ImageDirectionType3D = ImageType3D::DirectionType;
+
+//        ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion();
+//        SizeType3D        size3D = region3D.GetSize();
+//        ImageType3D::PointType       origin3D = rescaler2->GetOutput()->GetOrigin();
+//        const itk::Vector<double, 3> resolution3D = rescaler2->GetOutput()->GetSpacing();
+//        ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection();
+
+//        /* calculate image size in 3D Space */
+//        using VectorType = itk::Vector<double, 3>;
+//        VectorType Dest;
+//        Dest[0]=(size3D[0]-1) * resolution3D[0];
+//        Dest[1]=(size3D[1]-1) * resolution3D[1];
+//        Dest[2]=(size3D[2]-1) * resolution3D[2];
+
+//        ImageType3D::PointType LastVoxelPosition =
+//                origin3D + (imagDirection * Dest);
+
+//        //        std::cout<<" -------- Localizer 2 ITK --------"<<std::endl;
+//        //        std::cout<<"size3D " <<size3D <<std::endl;
+//        //        std::cout<<"resolution3D " <<resolution3D <<std::endl;
+//        //        std::cout<<"imagDirection " <<imagDirection <<std::endl;
+//        //        std::cout<<"First Voxel position " <<origin3D <<std::endl;
+//        //        std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
+
+//        //        double* dBounds = toVTKLocalizer2->GetOutput()->GetBounds();
+//        //        std::cout<< "-------- Localizer 2 VTK --------" <<std::endl;
+//        //        std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
+//        //                             <<dBounds[2]<<" "<<dBounds[3]<<" "
+//        //                             <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
+//        //        std::cout<< "-------- -------- --------" <<std::endl;
+
+//    }
+
+
+
    return
            toVTKLocalizer2->GetOutput();
 }

-
 vtkImageData* itkImageProcessor::GetProjection1VTK()
 {

@@ -2021,53 +2256,69 @@ vtkImageData* itkImageProcessor::GetProjection1VTK()
        //TODO
    }

+    // Rescale the intensity of the projection images to 0-32768 for output.
    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
    rescaler1->SetOutputMinimum(0);
-    rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
+    rescaler1->SetOutputMaximum(32768);
    rescaler1->SetInput( imageDRT1In );
    rescaler1->Update();

+    //    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
+    //    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
+    //    imageCalculatorFilter->SetImage(imageDRT1In);
+    //    imageCalculatorFilter->Compute();
+
+    //    std::cout<< "itkImageProcessor::imageDRT1In() " <<
+    //                imageCalculatorFilter <<std::endl;
+
+    //    using ImageCalculatorFilterType2 = itk::MinimumMaximumImageCalculator<OutputImageType>;
+    //    auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New();
+    //    imageCalculatorFilter2->SetImage(rescaler1->GetOutput());
+    //    imageCalculatorFilter2->Compute();
+    //    std::cout<< "itkImageProcessor::imageDRT2In() " <<
+    //                imageCalculatorFilter2 <<std::endl;
+
    toVTK2D1->SetInput(rescaler1->GetOutput());
    toVTK2D1->Update();

-    return
-            toVTK2D1->GetOutput();
-}

-vtkImageData* itkImageProcessor::GetProjection1VTKToWrite()
-{
+    //    using ImageRegionType3D = ImageType3D::RegionType;
+    //    using SizeType3D = ImageRegionType3D::SizeType;
+    //    using ImageDirectionType3D = ImageType3D::DirectionType;

-    if(m_DRTImage1MetaInfo == NULL ||
-            m_DRTGeometryMetaInfo == NULL ||
-            m_TransformMetaInfo == NULL     ){
-        return NULL;
-        //TODO
-    }
+    //    ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion();
+    //    SizeType3D        size3D = region3D.GetSize();
+    //    ImageType3D::PointType       origin3D = rescaler1->GetOutput()->GetOrigin();
+    //    const itk::Vector<double, 3> resolution3D = rescaler1->GetOutput()->GetSpacing();
+    //    ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection();

-    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
-    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
-    imageCalculatorFilter->SetImage(imageDRT1In);
-    imageCalculatorFilter->Compute();
+    //    /* calculate image size in 3D Space */
+    //    using VectorType = itk::Vector<double, 3>;
+    //    VectorType Dest;
+    //    Dest[0]=(size3D[0]-1) * resolution3D[0];
+    //    Dest[1]=(size3D[1]-1) * resolution3D[1];
+    //    Dest[2]=(size3D[2]-1) * resolution3D[2];

-    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
-    auto intWindowFilter = IntWindowType::New();
-    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
-    intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());
+    //    ImageType3D::PointType LastVoxelPosition =
+    //            origin3D + (imagDirection * Dest);

-    intWindowFilter->SetOutputMinimum(0);
-    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX){
-        intWindowFilter->SetOutputMaximum(SHRT_MAX);
-    }else
-        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());
+    //    std::cout<<" -------- Projection 1 ITK --------"<<std::endl;
+    //    std::cout<<"size3D " <<size3D <<std::endl;
+    //    std::cout<<"resolution3D " <<resolution3D <<std::endl;
+    //    std::cout<<"imagDirection " <<imagDirection <<std::endl;
+    //    std::cout<<"First Voxel position " <<origin3D <<std::endl;
+    //    std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;


-    intWindowFilter->SetInput(imageDRT1In);
-    intWindowFilter->Update();
+    //    double* dBounds = toVTK2D1->GetOutput()->GetBounds();

-    toVTK2D1->SetInput(intWindowFilter->GetOutput());
-    toVTK2D1->Update();
+    //    std::cout<< "-------- Proj 1 --------" <<std::endl;
+    //    std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
+    //                         <<dBounds[2]<<" "<<dBounds[3]<<" "
+    //                         <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
+    //    std::cout<< "-------- -------- --------" <<std::endl;

    return
            toVTK2D1->GetOutput();
@@ -2119,6 +2370,8 @@ vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform()
    m_Projection2VTKTransform->SetElement(2,3,
                                          interpolator2->GetInverseTransform()->GetTranslation()[2]);

+
+
    return
            m_Projection2VTKTransform;

@@ -2134,53 +2387,56 @@ vtkImageData* itkImageProcessor::GetProjection2VTK()
        //TODO
    }

-
+    // Rescale the intensity of the projection images to 0-32768 for output.
    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
    rescaler2->SetOutputMinimum(0);
-    rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
+    rescaler2->SetOutputMaximum(32768);
    rescaler2->SetInput( imageDRT2In );
-    rescaler2->Update();

+    rescaler2->Update();
    toVTK2D2->SetInput(rescaler2->GetOutput());
    toVTK2D2->Update();

-    return
-            toVTK2D2->GetOutput();
-}
+    //    rescaler2->Print(std::cout);

-vtkImageData* itkImageProcessor::GetProjection2VTKToWrite()
-{
-    if(m_DRTImage2MetaInfo == NULL ||
-            m_DRTGeometryMetaInfo == NULL ||
-            m_TransformMetaInfo == NULL     ){
-        return NULL;
-        //TODO
-    }
+    //    using ImageRegionType3D = ImageType3D::RegionType;
+    //    using SizeType3D = ImageRegionType3D::SizeType;
+    //    using ImageDirectionType3D = ImageType3D::DirectionType;

-    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
-    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
-    imageCalculatorFilter->SetImage(imageDRT2In);
-    imageCalculatorFilter->Compute();
+    //    ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion();
+    //    SizeType3D        size3D = region3D.GetSize();
+    //    ImageType3D::PointType       origin3D = rescaler2->GetOutput()->GetOrigin();
+    //    const itk::Vector<double, 3> resolution3D = rescaler2->GetOutput()->GetSpacing();
+    //    ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection();

-    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
-    auto intWindowFilter = IntWindowType::New();
-    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
-    intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());
+    /* calculate image size in 3D Space */
+    //    using VectorType = itk::Vector<double, 3>;
+    //    VectorType Dest;
+    //    Dest[0]=(size3D[0]-1) * resolution3D[0];
+    //    Dest[1]=(size3D[1]-1) * resolution3D[1];
+    //    Dest[2]=(size3D[2]-1) * resolution3D[2];

-    intWindowFilter->SetOutputMinimum(0);
-    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX)
-        intWindowFilter->SetOutputMaximum(SHRT_MAX);
-    else
-        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());
+    //    ImageType3D::PointType LastVoxelPosition =
+    //            origin3D + (imagDirection * Dest);
+
+    //    std::cout<<" -------- Projection 2 ITK --------"<<std::endl;
+    //    std::cout<<"size3D " <<size3D <<std::endl;
+    //    std::cout<<"resolution3D " <<resolution3D <<std::endl;
+    //    std::cout<<"imagDirection " <<imagDirection <<std::endl;
+    //    std::cout<<"First Voxel position " <<origin3D <<std::endl;
+    //    std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;


-    intWindowFilter->SetInput(imageDRT2In);
-    intWindowFilter->Update();

-    toVTK2D2->SetInput(intWindowFilter->GetOutput());
-    toVTK2D2->Update();
+
+    //double* dBounds = toVTK2D2->GetOutput()->GetBounds();
+    //        std::cout<< "-------- Proj 2 --------" <<std::endl;
+    //        std::cout<<"Bounds: "<<dBounds[0]<<" "<<dBounds[1]<<" "
+    //                             <<dBounds[2]<<" "<<dBounds[3]<<" "
+    //                             <<dBounds[4]<<" "<<dBounds[5]<<std::endl;
+    //        std::cout<< "-------- -------- --------" <<std::endl;

    return
            toVTK2D2->GetOutput();
@@ -2196,11 +2452,14 @@ void itkImageProcessor::WriteProjectionImages()
    rescaler1->SetOutputMinimum(0);
    rescaler1->SetOutputMaximum(255);
    rescaler1->SetInput( imageDRT1In );
+    //                flipFilter1->GetOutput());

    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
    rescaler2->SetOutputMinimum(0);
    rescaler2->SetOutputMaximum(255);
    rescaler2->SetInput( imageDRT2In );
+    //                flipFilter2->GetOutput());
+

    using WriterType = itk::ImageFileWriter<OutputImageType>;
    WriterType::Pointer writer1 = WriterType::New();
@@ -2211,6 +2470,7 @@ void itkImageProcessor::WriteProjectionImages()

    try
    {
+        //        std::cout << "Writing image 1 " << std::endl;
        writer1->Update();
    }
    catch (itk::ExceptionObject & err)
@@ -2224,6 +2484,7 @@ void itkImageProcessor::WriteProjectionImages()

    try
    {
+        //        std::cout << "Writing image 2" << std::endl;
        writer2->Update();
    }
    catch (itk::ExceptionObject & err)
@@ -2234,6 +2495,10 @@ void itkImageProcessor::WriteProjectionImages()

 }

+
+
+
+
 void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ)
 {

@@ -2258,6 +2523,7 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ)

    m_TransformMetaInfo->SetUserRotations(Rotations);

+
 }

 Optimizer::ParametersType
@@ -2345,6 +2611,9 @@ void itkImageProcessor::SetRegionFixed1(
    roiAutoReg1.SetIndex(index1);
    roiAutoReg1.SetSize(size1);

+    //std::cout << "itkImageProcessor " << std::endl;
+    //std::cout << roiAutoReg1 << std::endl;
+
    this->m_R23->SetroiAutoReg1(roiAutoReg1);

 }
@@ -2367,12 +2636,21 @@ void itkImageProcessor::SetRegionFixed2(
    roiAutoReg2.SetIndex(index2);
    roiAutoReg2.SetSize(size2);

+    //std::cout << "itkImageProcessor " << std::endl;
+    //std::cout << roiAutoReg2 << std::endl;
+
    this->m_R23->SetroiAutoReg2(roiAutoReg2);

 }

 void itkImageProcessor::ResetROIRegions(){

+//    auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion();
+//    auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
+
+//    this->m_R23->SetroiAutoReg1(region1temp);
+//    this->m_R23->SetroiAutoReg2(region2temp);
+
    this->m_R23->setROISizeToZero();

 }
--- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h
@@ -22,7 +22,6 @@ gfattori 08.11.2021
 #include "itkResampleImageFilter.h"
 #include "itkCastImageFilter.h"
 #include "itkRescaleIntensityImageFilter.h"
-#include "itkIntensityWindowingImageFilter.h"
 #include "itkChangeInformationImageFilter.h"
 #include "itkGDCMImageIO.h"
 #include "itkMetaDataObject.h"
@@ -50,8 +49,6 @@ gfattori 08.11.2021
 #include "itkReg23.h"
 #include "itkReg23MetaInformation.h"

-#define IMG_VIS_MAXIMUM_RANGE 2048
-
 namespace itk
 {

@@ -245,10 +242,6 @@ public:
    vtkImageData* GetProjection1VTK();
    vtkImageData* GetProjection2VTK();

-    /** Conveniency methods to get VTK images for rendering */
-    vtkImageData* GetProjection1VTKToWrite();
-    vtkImageData* GetProjection2VTKToWrite();
-
    vtkImageData* GetLocalizer1VTK();
    vtkImageData* GetLocalizer2VTK();

@@ -342,6 +335,7 @@ private:
    imageDRT2In;


+
    DuplicatorType::Pointer
    m_LATSourceDupli,
    m_PASourceDupli,
--- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
@@ -63,6 +63,8 @@ MapTransformToNewOrigin(
    return m_OutputTransform;
 }

+
+
 TransformType::Pointer
 CalculateInternalTransformV3(
        ImageType3D::PointType m_Translation,       //IEC
--- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h
--- a/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h
+++ b/reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h
@@ -54,6 +54,9 @@ private:
    bAbortProcessCommand;


+
+
+
 signals:
 void sendRegistrationProgress(int,double,double,double,double);

@@ -64,6 +67,7 @@ namespace itk


 class CommandIterationUpdate : public itk::Command {
+    // TODO: Move to own files.

    constexpr static unsigned int Dimension = 3;

@@ -123,6 +127,7 @@ public:
            objIterUpdate->setAbortFlag(false);
            throw itk::ProcessAborted();
        }
+//        std::cout << "Progress: " << this->m_Process->GetAbortGenerateData() << std::endl;

        OptimizerPointer optPow;
        AmoebaOptimizerPointer optAm;
@@ -163,6 +168,7 @@ public:
 };

 class ExhaustiveCommandIterationUpdate : public itk::Command {
+    // TODO: Move to own files.

 public:
    using Self = ExhaustiveCommandIterationUpdate;
--- a/reg23Topograms/itkDTRrecon/itkReg23.cpp
+++ b/reg23Topograms/itkDTRrecon/itkReg23.cpp
@@ -30,6 +30,8 @@ itkReg23::itkReg23()
    registration->SetTransform1(Transform1);
    registration->SetTransform2(Transform2);

+
+
    // to be provided by the user
    m_r23Meta = nullptr;
    m_Volume = nullptr;
@@ -42,11 +44,14 @@ itkReg23::itkReg23()
    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()
 {

@@ -119,6 +124,21 @@ void itkReg23::InitializeRegistration()
        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:
@@ -177,6 +197,14 @@ void itkReg23::InitializeRegistration()
    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:
@@ -260,9 +288,17 @@ void itkReg23::InitializeRegistration()

    registration->SetIsocIECTransform(IsocTransform);

+    //    if (verbose) {
+    //        registration->DebugOn();
+    //        registration->Print(std::cout);
+    //    }
+
+
    std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
 }

+
+
 int itkReg23::StartRegistration(std::string extraInfo)
 {

@@ -270,6 +306,8 @@ int itkReg23::StartRegistration(std::string extraInfo)
    // 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);

@@ -285,6 +323,9 @@ int itkReg23::StartRegistration(std::string extraInfo)
        fs << extraInfo;
        fs << "Value\tX\tY\tZ " << std::endl;

+        //        if (r23Meta->GetOptimizerType() == tOptimizerTypeEnum::EXHAUSTIVE) {
+        //            ExhaustiveOptimizerObserver->set_stream(fs);
+        //        }
    }

    // Start the registration
@@ -317,6 +358,7 @@ int itkReg23::StartRegistration(std::string extraInfo)
        this->GetCurrentPosition();
    std::cout<<" FinalPars: "<< finalParameters <<std::endl;

+
    std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;

    return 0;
@@ -327,9 +369,11 @@ 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()) {
@@ -390,4 +434,6 @@ double itkReg23::GetCurrentMetricValue()
 }


+
+
 }
--- a/reg23Topograms/itkDTRrecon/itkReg23.h
+++ b/reg23Topograms/itkDTRrecon/itkReg23.h
@@ -161,7 +161,6 @@ private:
    ExhaustiveOptimizer;

    /* ---- User provided */
-    
    R23MetaInformation::Pointer
    m_r23Meta;

@@ -202,6 +201,7 @@ private:

    /* ----  User provided  END */

+
    TransformType::Pointer
    IsocTransform,
    Transform1,
--- a/reg23Topograms/itkDTRrecon/itkReg23MetaInformation.cpp
+++ b/reg23Topograms/itkDTRrecon/itkReg23MetaInformation.cpp
@@ -21,12 +21,14 @@ PowellOptimizerMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 PowellOptimizerMetaInformation
 ::~PowellOptimizerMetaInformation ()
 {

 }

+
 AmoebaOptimizerMetaInformation::
 AmoebaOptimizerMetaInformation(){

@@ -37,6 +39,7 @@ AmoebaOptimizerMetaInformation(){
    this->m_MaxIterations = 100;
 }

+
 void
 AmoebaOptimizerMetaInformation
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
@@ -44,12 +47,14 @@ AmoebaOptimizerMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 AmoebaOptimizerMetaInformation
 ::~AmoebaOptimizerMetaInformation ()
 {

 }

+
 MIMetricMetaInformation::
 MIMetricMetaInformation(){

@@ -58,6 +63,7 @@ MIMetricMetaInformation(){

 }

+
 void
 MIMetricMetaInformation
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
@@ -65,6 +71,7 @@ MIMetricMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 MIMetricMetaInformation
 ::~MIMetricMetaInformation ()
 {
@@ -79,6 +86,7 @@ NCCMetricMetaInformation(){

 }

+
 void
 NCCMetricMetaInformation
 ::PrintSelf(std::ostream& os, itk::Indent indent) const
@@ -86,6 +94,7 @@ NCCMetricMetaInformation
    Superclass::PrintSelf(os, indent);
 }

+
 NCCMetricMetaInformation
 ::~NCCMetricMetaInformation ()
 {
--- a/reg23Topograms/itkDTRrecon/itkReg23MetaInformation.h
+++ b/reg23Topograms/itkDTRrecon/itkReg23MetaInformation.h
@@ -173,6 +173,7 @@ private:
    void operator=(const Self&);
 };

+
 class NCCMetricMetaInformation :
        public itk::Object{

@@ -204,6 +205,8 @@ protected:

    bool m_SubtractMean;

+
+
    /** Default Constructor **/
    NCCMetricMetaInformation  ();
    /** Default Destructor **/
--- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
+++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
@@ -133,7 +133,6 @@ public:

  /** Internal floating point comparison accuracy **/
  static const double EPSILON;
-  static const double TOL;

  /** \brief
   * Interpolate the image at a point position.
--- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
+++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
@@ -56,9 +56,6 @@ CIT 6, 89-94 (1998).
 namespace itk
 {

-template <typename TInputImage, typename TCoordRep>
-const double gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::TOL =
-    1e-9;

 template <typename TInputImage, typename TCoordRep>
 const double gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EPSILON =
@@ -116,25 +113,26 @@ template <typename TInputImage, typename TCoordRep>
 typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
 gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
 {
-    // --- Geometry / Siddon working variables use double ---
-    double    rayVector[3];
+    float     rayVector[3];
    IndexType cIndex;

    PointType  drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
    OutputType pixval;

-    double firstIntersection[3];
-    double alphaX1, alphaXN, alphaXmin, alphaXmax;
-    double alphaY1, alphaYN, alphaYmin, alphaYmax;
-    double alphaZ1, alphaZN, alphaZmin, alphaZmax;
-    double alphaMin, alphaMax;
-    double alphaX, alphaY, alphaZ;
-    double alphaUx, alphaUy, alphaUz;
-    double alphaIntersectionUp[3], alphaIntersectionDown[3];
-    double d12, value;
-    double firstIntersectionIndex[3];
-    int    firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
-    int    iU, jU, kU;
+
+    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
@@ -142,12 +140,18 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
    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 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;
@@ -155,7 +159,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c

    drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);

-    // Get the input pointers
+    // Get ths input pointers
    InputImageConstPointer inputPtr = this->GetInputImage();

    typename InputImageType::SizeType    sizeCT;
@@ -164,12 +168,13 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
    typename InputImageType::PointType   ctOrigin;

    ctPixelSpacing = inputPtr->GetSpacing();
-    ctOrigin       = inputPtr->GetOrigin();
-    regionCT       = inputPtr->GetLargestPossibleRegion();
-    sizeCT         = regionCT.GetSize();
+    ctOrigin = inputPtr->GetOrigin();
+    regionCT = inputPtr->GetLargestPossibleRegion();
+    sizeCT = regionCT.GetSize();
+

    PointType SlidingSourcePoint = m_SourcePoint;
-    SlidingSourcePoint[0]        = 0.;
+    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];
@@ -178,22 +183,22 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
    PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);

    PointType O(3);
-    O[0] = -ctPixelSpacing[0] / 2.;
-    O[1] = -ctPixelSpacing[1] / 2.;
-    O[2] = -ctPixelSpacing[2] / 2.;
+    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
-    double xSizeCT = O[0] + (sizeCT[0]) * ctPixelSpacing[0];
-    double ySizeCT = O[1] + (sizeCT[1]) * ctPixelSpacing[1];
-    double zSizeCT = O[2] + (sizeCT[2]) * ctPixelSpacing[2];
-    double xDrrPix = drrPixelWorld[0];
-    double yDrrPix = drrPixelWorld[1];
-    double zDrrPix = drrPixelWorld[2];
-
-    if (xDrrPix <= O[0] || xDrrPix >= xSizeCT ||
-        yDrrPix <= O[1] || yDrrPix >= ySizeCT ||
-        zDrrPix <= O[2] || zDrrPix >= zSizeCT)
-    {
+    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;
    }
@@ -211,351 +216,260 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
    p2[1] = drrPixelWorld[1];
    p2[2] = drrPixelWorld[2];

-    double P12D = std::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 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 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];

-    double pDest[3];
-    pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2.0 * p12v[0];
-    pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2.0 * p12v[1];
-    pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2.0 * p12v[2];

    // The following is the Siddon-Jacob fast ray-tracing algorithm
-    // Compute ray from source to "detector" point pDest
+    //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];

-    const double rayLength = std::sqrt(rayVector[0]*rayVector[0] +
-                                       rayVector[1]*rayVector[1] +
-                                       rayVector[2]*rayVector[2]);

-    // just in case (should never really be zero)
-    if (rayLength <= EPSILON)
+    /* 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*/)
    {
-        pixval = static_cast<OutputType>(0.0);
-        return pixval;
-    }
-    const double INF = std::numeric_limits<double>::infinity();
-
-    // Convenience bounds for the CT volume in world coords
-    double xMin = O[0];
-    double xMax = O[0] + sizeCT[0] * ctPixelSpacing[0];
-    double yMin = O[1];
-    double yMax = O[1] + sizeCT[1] * ctPixelSpacing[1];
-    double zMin = O[2];
-    double zMax = O[2] + sizeCT[2] * ctPixelSpacing[2];
-
-    /* Calculate the parametric values of the first and last
-       intersection points of the ray with the X, Y, and Z-planes
-       that define the CT volume. */
-
-    // ---- X slab ----
-    if (std::fabs(rayVector[0]) > EPSILON)
-    {
-        alphaX1 = (xMin - SourceWorld[0]) / rayVector[0];
-        alphaXN = (xMax - SourceWorld[0]) / rayVector[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
    {
-        // Ray is parallel to X planes: must lie within [xMin, xMax] or no intersection
-        if (SourceWorld[0] < xMin || SourceWorld[0] > xMax)
-        {
-            pixval = static_cast<OutputType>(0.0);
-            return pixval;
-        }
-        alphaXmin = -INF;
-        alphaXmax =  INF;
+        alphaXmin = -2;
+        alphaXmax = 2;
    }

-    // ---- Y slab ----
-    if (std::fabs(rayVector[1]) > EPSILON)
+    if (fabs(rayVector[1]) > EPSILON/* != 0*/)
    {
-        alphaY1 = (yMin - SourceWorld[1]) / rayVector[1];
-        alphaYN = (yMax - SourceWorld[1]) / rayVector[1];
+        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
    {
-        if (SourceWorld[1] < yMin || SourceWorld[1] > yMax)
-        {
-            pixval = static_cast<OutputType>(0.0);
-            return pixval;
-        }
-        alphaYmin = -INF;
-        alphaYmax =  INF;
+        alphaYmin = -2;
+        alphaYmax = 2;
    }

-    // ---- Z slab ----
-    if (std::fabs(rayVector[2]) > EPSILON)
+    if (fabs(rayVector[2]) > EPSILON/* != 0*/)
    {
-        alphaZ1 = (zMin - SourceWorld[2]) / rayVector[2];
-        alphaZN = (zMax - SourceWorld[2]) / rayVector[2];
+        //    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
    {
-        if (SourceWorld[2] < zMin || SourceWorld[2] > zMax)
-        {
-            pixval = static_cast<OutputType>(0.0);
-            return pixval;
-        }
-        alphaZmin = -INF;
-        alphaZmax =  INF;
+        alphaZmin = -2;
+        alphaZmax = 2;
    }

    /* Get the very first and the last alpha values when the ray
-       intersects with the CT volume. */
+  intersects with the CT volume. */
    alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin);
    alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax);

-    // No intersection at all
-    if (alphaMax <= alphaMin)
-    {
-        pixval = static_cast<OutputType>(0.0);
-        return pixval;
-    }
-
-    // Restrict to the [0,1] segment (from SourceWorld to pDest)
-    alphaMax = std::min(alphaMax, 1.0);
-    alphaMin = std::max(alphaMin, 0.0);
-
-    if (alphaMax <= alphaMin)
-    {
-        pixval = static_cast<OutputType>(0.0);
-        return pixval;
-    }
-
    /* Calculate the parametric values of the first intersection point
-       after the ray entered the CT volume. */
-    firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0] - O[0];
-    firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - O[1];
-    firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2] - O[2];
+  of the ray with the X, Y, and Z-planes after the ray entered the
+  CT volume. */

-    /* Transform world coordinates to continuous indices of 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];

-    // Helper for equality with tolerance (for on-grid cases)
-    auto nearlyEqual = [](double a, double b)
+    // now "correct" the indices (there are special cases):
+    if (alphaMin == alphaYmin) // j is special
    {
-        return std::fabs(a - b) < TOL;
-    };
-
-    // "Correct" the indices for special entry cases
-    if (nearlyEqual(alphaMin, alphaXmin)) // i is special
-    {
-        if (nearlyEqual(alphaXmin, alphaX1))
-            firstIntersectionIndex[0] = 0.0;            // continuous index
+        if (alphaYmin == alphaY1)
+            firstIntersectionIndex [1] = 0;
        else
-            firstIntersectionIndex[0] = sizeCT[0] - 1.; // continuous index
+            // 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;
    }

-    if (nearlyEqual(alphaMin, alphaYmin)) // j is special
+    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*/)
    {
-        if (nearlyEqual(alphaYmin, alphaY1))
-            firstIntersectionIndex[1] = 0.0;
-        else
-            firstIntersectionIndex[1] = sizeCT[1] - 1.;
-    }
-
-    if (nearlyEqual(alphaMin, alphaZmin)) // k is special
-    {
-        if (nearlyEqual(alphaZmin, alphaZ1))
-            firstIntersectionIndex[2] = 0.0;
-        else
-            firstIntersectionIndex[2] = sizeCT[2] - 1.;
-    }
-
-    firstIntersectionIndexUp[0]   = static_cast<int>(std::ceil(firstIntersectionIndex[0]));
-    firstIntersectionIndexUp[1]   = static_cast<int>(std::ceil(firstIntersectionIndex[1]));
-    firstIntersectionIndexUp[2]   = static_cast<int>(std::ceil(firstIntersectionIndex[2]));
-    firstIntersectionIndexDown[0] = static_cast<int>(std::floor(firstIntersectionIndex[0]));
-    firstIntersectionIndexDown[1] = static_cast<int>(std::floor(firstIntersectionIndex[1]));
-    firstIntersectionIndexDown[2] = static_cast<int>(std::floor(firstIntersectionIndex[2]));
-
-    // ---- Compute first plane hits and increments ----
-
-    // X
-    if (std::fabs(rayVector[0]) <= EPSILON)
-    {
-        alphaX  = INF;
-        alphaUx = INF;
+        alphaX = 2;
    }
    else
    {
        alphaIntersectionUp[0] =
-            (O[0] + firstIntersectionIndexUp[0]   * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[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]);
-        alphaUx = ctPixelSpacing[0] / std::fabs(rayVector[0]);
+                (O[0] + firstIntersectionIndexDown[0] * ctPixelSpacing[0]  - SourceWorld[0]) / rayVector[0];
+        alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
    }

-    // Y
-    if (std::fabs(rayVector[1]) <= EPSILON)
+    if (fabs(rayVector[1] ) <= EPSILON/* == 0*/)
    {
-        alphaY  = INF;
-        alphaUy = INF;
+        alphaY = 2;
    }
    else
    {
        alphaIntersectionUp[1] =
-            (O[1] + firstIntersectionIndexUp[1]   * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[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]);
-        alphaUy = ctPixelSpacing[1] / std::fabs(rayVector[1]);
+                (O[1] + firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+        alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
    }

-    // Z
-    if (std::fabs(rayVector[2]) <= EPSILON)
+    if (fabs(rayVector[2]) <= EPSILON/* == 0*/)
    {
-        alphaZ  = INF;
-        alphaUz = INF;
+        alphaZ = 2;
    }
    else
    {
        alphaIntersectionUp[2] =
-            (O[2] + firstIntersectionIndexUp[2]   * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[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]);
-        alphaUz = ctPixelSpacing[2] / std::fabs(rayVector[2]);
+                (O[2] + firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+        alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
    }

-    // Make sure the first plane alphas are not before entry
-
-    // Only advance if this axis is not parallel (alphaUx finite)
-    if (!std::isinf(alphaUx))
+    /* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */
+    if (fabs(rayVector[0]) > EPSILON /*!= 0*/)
    {
-        while (alphaX <= alphaMin + TOL)
-            alphaX += alphaUx;
+        alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
    }
-    if (!std::isinf(alphaUy))
+    else
    {
-        while (alphaY <= alphaMin + TOL)
-            alphaY += alphaUy;
+        alphaUx = 999;
    }
-    if (!std::isinf(alphaUz))
+    if (fabs(rayVector[1]) > EPSILON /*!= 0*/)
    {
-        while (alphaZ <= alphaMin + TOL)
-            alphaZ += alphaUz;
+        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. */
-    // Use the direction of the ray
-    iU = (rayVector[0] > 0.0) ? 1 : -1;
-    jU = (rayVector[1] > 0.0) ? 1 : -1;
-    kU = (rayVector[2] > 0.0) ? 1 : -1;

-    // Initialize line integral
-    d12 = 0.0;
+    iU = (SourceWorld[0] < drrPixelWorld[0]) ? 1 : -1;
+    jU = (SourceWorld[1] < drrPixelWorld[1]) ? 1 : -1;
+    kU = (SourceWorld[2] < drrPixelWorld[2]) ? 1 : -1;

-    /* Initialize the current voxel index with the voxel entered at alphaMin. */
+    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];

-    // Current parametric position and previous one
-    double alpha     = alphaMin;
-    double alphaPrev = alphaMin;
-
-    // Siddon traversal
-    while (alpha < alphaMax)
+    while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
    {
-        // Next plane intersection along the ray
-        double alphaNext = std::min(std::min(alphaX, alphaY), alphaZ);
-        if (alphaNext > alphaMax)
-            alphaNext = alphaMax;
+        /* Store the current ray position */
+        alphaCminPrev = alphaCmin;

-        double segmentLength = alphaNext - alphaPrev;
-        // don’t break immediately on a single zero-length segment. This would cause black pixels!
-        // if (segmentLength <= 0.0)
-        //  break; // FP guard, especially for grid-aligned rays
-        if (segmentLength <= 0.0)
+        if ((alphaX <= alphaY) && (alphaX <= alphaZ))
        {
-            // advance to the next plane and try again
-            if (nearlyEqual(alphaNext, alphaX))
-            {
-                cIndex[0] += iU;
-                alphaX    += alphaUx;
-            }
-            else if (nearlyEqual(alphaNext, alphaY))
-            {
-                cIndex[1] += jU;
-                alphaY    += alphaUy;
-            }
-            else
-            {
-                cIndex[2] += kU;
-                alphaZ    += alphaUz;
-            }
-
-            alphaPrev = alphaNext;
-            alpha     = alphaNext;
-
-            if (alpha >= alphaMax)
-                break;
-
-            continue;
+            /* 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;
        }

-        // Accumulate contribution if inside the volume
-        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 ((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)
+            if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
            {
-                // physical path length = segmentLength (in α) * |rayVector|
-                d12 += segmentLength * rayLength * (value - m_Threshold);
+                d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
            }
        }
-
-        if (alphaNext >= alphaMax)
-            break;
-
-        alphaPrev = alphaNext;
-        alpha     = alphaNext;
-
-        // Advance along the axis whose plane we actually hit
-        if (nearlyEqual(alphaNext, alphaX))
-        {
-            cIndex[0] += iU;
-            alphaX    += alphaUx;
-        }
-        else if (nearlyEqual(alphaNext, alphaY))
-        {
-            cIndex[1] += jU;
-            alphaY    += alphaUy;
-        }
-        else // must be z
-        {
-            cIndex[2] += kU;
-            alphaZ    += alphaUz;
-        }
    }

-    // Clamp output value
    if (d12 < minOutputValue)
    {
        pixval = minOutputValue;
@@ -568,7 +482,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
    {
        pixval = static_cast<OutputType>(d12);
    }
-    return pixval;
+    return (pixval);
 }


--- a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx
+++ b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx
@@ -53,6 +53,7 @@ vtkContourTopogramProjectionFilter::vtkContourTopogramProjectionFilter()
                0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
 }

+
 vtkContourTopogramProjectionFilter::~vtkContourTopogramProjectionFilter()
 {
    if (this->m_Transform)
@@ -110,6 +111,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
    // 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;
    }
@@ -117,6 +119,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
    // 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;
    }
@@ -124,6 +127,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
    // 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;
    }
@@ -194,6 +198,7 @@ int vtkContourTopogramProjectionFilter::RequestData(

    return 1;

+
 }

 void vtkContourTopogramProjectionFilter::SetImportOffsetLPS(const double * dP)
--- a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h
+++ b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h
--- a/third_party_licenses.txt
+++ b/third_party_licenses.txt
@@ -1,24 +0,0 @@
-This project uses the following third-party libraries:
-
---
-
-GDCM (Grassroots DICOM)
-Website: https://github.com/malaterre/GDCM
-License: BSD-style license
-
---
-
-VTK (The Visualization Toolkit)
-Website: https://vtk.org/
-License: BSD 3-Clause License
-
---
-
-ITK (Insight Segmentation and Registration Toolkit)
-Website: https://itk.org/
-License: Apache License 2.0
-
---
-
-These libraries are used as external dependencies and are not redistributed with this source code.
-Please refer to their respective websites for full license texts.
Author	SHA1	Message	Date
via_r	db542cea9a	Merge branch 'cleanUpUIHandler' into 'ScoutRT_Qt6_handler' working on Handler in qt6 See merge request cpt_bioeng/drt!44	2024-02-22 16:40:26 +01:00
fattori_g	74736603c6	itk metric wrong include	2024-02-06 17:06:33 +01:00
fattori_g	ed75b0920d	most of CTK should be clear	2024-02-06 16:56:27 +01:00