fix to issue: Disproportionately large alpha in certain geometry cases

Update README.md
2026-01-12 22:35:44 +01:00 · 2025-06-25 00:11:48 +02:00 · 2025-06-25 00:09:20 +02:00 · 2025-06-25 00:01:28 +02:00 · 2025-06-24 23:59:54 +02:00 · 2025-06-24 23:57:26 +02:00
30 changed files with 509 additions and 736 deletions
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -0,0 +1,20 @@
+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
@@ -0,0 +1,29 @@
+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
@@ -0,0 +1,47 @@
+# 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 3.16)
+cmake_minimum_required(VERSION 2.8.12)

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

 add_subdirectory(autoreg)

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

+    
 namespace itk
 {

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

 }

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

-
 TopogramImageMetaInformation
 ::~TopogramImageMetaInformation ()
 {

 }

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

 }

-
 DRTImageMetaInformation::
 DRTImageMetaInformation(){

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

 }

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

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

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

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

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

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

-
 DRTProjectionGeometryImageMetaInformation::
 DRTProjectionGeometryImageMetaInformation(){

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

-
 DRTProjectionGeometryImageMetaInformation
 ::~DRTProjectionGeometryImageMetaInformation    ()
 {

 }

-
-
-
 RTGeometryMetaInformation::
 RTGeometryMetaInformation(){

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

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

 }

-
 R23MetaInformation::
 R23MetaInformation (){

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

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

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


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

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

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

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

-
 class DegreeOfFreedom : public itk::Object {

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

 protected:

-
    tPatOrientation
    m_PatientOrientation;

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

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

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

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

-
-
-
 protected:

    tDegreeOfFreedomEnum
@@ -693,7 +681,6 @@ 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,6 +6,12 @@
 #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])
@@ -24,7 +30,7 @@ const char* computeMethodName(const char (&function)[FL], const char (&prettyFun

 namespace itk {
 /* reference string required for comparison with tag values */
-static const char* ImageOrientationStrings[] = {
+static const char VARIABLE_IS_NOT_USED *ImageOrientationStrings[] = {
    "NotDefined",
    "HFS",
    "FFS",
@@ -40,49 +46,49 @@ static const bool verbose = false;
 /* this is in the end a IEC to HFS...
 * but we keep this for ourselves...
 */
-static double Standard_DRT2LPS[9] = {
+static double VARIABLE_IS_NOT_USED Standard_DRT2LPS[9] = {
    1, 0, 0,
    0, 0, -1,
    0, 1, 0
 };

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

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

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

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

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

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

-static double PAT2IEC[9] = {
+static double VARIABLE_IS_NOT_USED 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,9 +20,7 @@

 #include "itkCovariantVector.h"
 #include "itkMacro.h"
-//#include "itkNormalizedCorrelationImageToImageMetric.hxx"
-#include "itkNormalizedCorrelationImageToImageMetric.h"
-
+#include "itkNormalizedCorrelationImageToImageMetric.hxx"
 #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,7 +164,6 @@ void itkImageProcessor::SetNumberOfWorkingUnits(int iN){
    }
 }

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

 }

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

-
 double itkImageProcessor::GetPanelOffsetPGeo(tProjOrientationType ImgPrj){

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

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

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

 }

-
 void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(
        double dX1,
        double dY1,
@@ -378,7 +371,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)
@@ -399,7 +392,6 @@ 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(){
@@ -424,7 +416,6 @@ int itkImageProcessor::unload3DVolumeAndMeta(){
    m_TransformMetaInfo = NULL;
    m_TransformMetaInfo = TransformMetaInformation::New();

-
    this->ResetROIRegions();

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

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

@@ -450,12 +441,6 @@ 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();

@@ -472,7 +457,6 @@ 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();

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

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

@@ -532,22 +514,13 @@ 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>(
@@ -576,17 +549,9 @@ 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();

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

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

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

-
    /* Calculate projection angle IEC to LPS */
    m_DRTImage2MetaInfo = DRTImageMetaInformation::New();
    m_DRTImage2MetaInfo->SetProjectionAngleLPS(
@@ -679,35 +640,12 @@ 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()
@@ -777,15 +715,12 @@ itkImageProcessor::CalcProjectionAngleLPS(
        tPatOrientation  pOrient,
        double pAngleIEC){

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

    InternalImageType::DirectionType
            currIECtoLPS;

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

-
 void itkImageProcessor::SetProjectionAngle1IEC(double dGantryAngle)
 {
    m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(dGantryAngle);
@@ -952,8 +886,6 @@ 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;
@@ -961,10 +893,8 @@ 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;
@@ -996,9 +926,6 @@ 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++){
@@ -1015,17 +942,10 @@ 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 =
@@ -1039,9 +959,6 @@ 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;

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

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

@@ -1077,10 +993,8 @@ 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;

@@ -1091,9 +1005,7 @@ 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;
@@ -1119,28 +1031,11 @@ 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)
 {

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

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

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


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

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

 }

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

-
 void itkImageProcessor::InitializeProjector()
 {

@@ -1504,7 +1393,6 @@ 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() );
@@ -1514,8 +1402,6 @@ void itkImageProcessor::InitializeProjector()
    filter1->ChangeOriginOn();
    filter1->UpdateOutputInformation();

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

 }

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

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

    m_RTMetaInfo->SetIsocenterIECS(Point);

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

-
    // FIRST

    ImageType3D::PointType NominalIsocenter_wZcorrectionLPS;
@@ -1679,12 +1560,6 @@ 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] +
@@ -1694,7 +1569,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
    CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2];

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

    ImageType3D::PointType CalibratedIsocenterZeroLPS;

@@ -1704,7 +1578,6 @@ 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(
@@ -1877,13 +1750,10 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){

 }

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

-
-
 void itkImageProcessor::GetProjectionImages(){

-    //    std::cout<< "itkImageProcessor::GetProjectionImages" <<std::endl;
    if(m_DRTImage1MetaInfo == NULL ||
            m_DRTImage2MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
@@ -1917,12 +1784,9 @@ 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();

@@ -1948,18 +1812,6 @@ 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.
@@ -2043,8 +1895,6 @@ itkImageProcessor::CalcImportVolumeOffset(
    ImageType3D::PointType IsoSupport_IECPos =
            InputTransform->TransformPoint(rtCouchOffset );

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

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


-
 void itkImageProcessor::Write2DImages(){

    //    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
@@ -2127,59 +1976,18 @@ 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(32768);
-    rescaler1->SetInput(m_PASourceDupli->GetOutput());
-
+    rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
+    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();
 }
@@ -2187,65 +1995,22 @@ 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(32768);
-    rescaler2->SetInput(m_LATSourceDupli ->GetOutput());
-
+    rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
+    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()
 {

@@ -2256,69 +2021,53 @@ 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(32768);
+    rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
    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();
+}

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

-    //    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();
+    if(m_DRTImage1MetaInfo == NULL ||
+            m_DRTGeometryMetaInfo == NULL ||
+            m_TransformMetaInfo == NULL     ){
+        return NULL;
+        //TODO
+    }

-    //    /* 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 ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
+    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
+    imageCalculatorFilter->SetImage(imageDRT1In);
+    imageCalculatorFilter->Compute();

-    //    ImageType3D::PointType LastVoxelPosition =
-    //            origin3D + (imagDirection * Dest);
+    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
+    auto intWindowFilter = IntWindowType::New();
+    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
+    intWindowFilter->SetWindowMaximum(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->SetOutputMinimum(0);
+    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX){
+        intWindowFilter->SetOutputMaximum(SHRT_MAX);
+    }else
+        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());


-    //    double* dBounds = toVTK2D1->GetOutput()->GetBounds();
+    intWindowFilter->SetInput(imageDRT1In);
+    intWindowFilter->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;
+    toVTK2D1->SetInput(intWindowFilter->GetOutput());
+    toVTK2D1->Update();

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

-
-
    return
            m_Projection2VTKTransform;

@@ -2387,56 +2134,53 @@ 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(32768);
+    rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
    rescaler2->SetInput( imageDRT2In );
-
    rescaler2->Update();
+
    toVTK2D2->SetInput(rescaler2->GetOutput());
    toVTK2D2->Update();

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

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

-    //    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 ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
+    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
+    imageCalculatorFilter->SetImage(imageDRT2In);
+    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);
-
-    //    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->SetOutputMinimum(0);
+    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX)
+        intWindowFilter->SetOutputMaximum(SHRT_MAX);
+    else
+        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());


+    intWindowFilter->SetInput(imageDRT2In);
+    intWindowFilter->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;
+    toVTK2D2->SetInput(intWindowFilter->GetOutput());
+    toVTK2D2->Update();

    return
            toVTK2D2->GetOutput();
@@ -2452,14 +2196,11 @@ 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();
@@ -2470,7 +2211,6 @@ void itkImageProcessor::WriteProjectionImages()

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

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

 }

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

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

    m_TransformMetaInfo->SetUserRotations(Rotations);

-
 }

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

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

 }
@@ -2636,21 +2367,12 @@ 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,6 +22,7 @@ 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"
@@ -49,6 +50,8 @@ gfattori 08.11.2021
 #include "itkReg23.h"
 #include "itkReg23MetaInformation.h"

+#define IMG_VIS_MAXIMUM_RANGE 2048
+
 namespace itk
 {

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

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

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


-
    DuplicatorType::Pointer
    m_LATSourceDupli,
    m_PASourceDupli,
--- a/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
@@ -63,8 +63,6 @@ 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,9 +54,6 @@ private:
    bAbortProcessCommand;


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

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


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

    constexpr static unsigned int Dimension = 3;

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

        OptimizerPointer optPow;
        AmoebaOptimizerPointer optAm;
@@ -168,7 +163,6 @@ 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,8 +30,6 @@ itkReg23::itkReg23()
    registration->SetTransform1(Transform1);
    registration->SetTransform2(Transform2);

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

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

@@ -306,8 +270,6 @@ 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);

@@ -323,9 +285,6 @@ 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
@@ -358,7 +317,6 @@ int itkReg23::StartRegistration(std::string extraInfo)
        this->GetCurrentPosition();
    std::cout<<" FinalPars: "<< finalParameters <<std::endl;

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

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


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

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

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

    /* ----  User provided  END */

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

-
 PowellOptimizerMetaInformation
 ::~PowellOptimizerMetaInformation ()
 {

 }

-
 AmoebaOptimizerMetaInformation::
 AmoebaOptimizerMetaInformation(){

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

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

-
 AmoebaOptimizerMetaInformation
 ::~AmoebaOptimizerMetaInformation ()
 {

 }

-
 MIMetricMetaInformation::
 MIMetricMetaInformation(){

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

 }

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

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

 }

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

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

-
 class NCCMetricMetaInformation :
        public itk::Object{

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

    bool m_SubtractMean;

-
-
    /** Default Constructor **/
    NCCMetricMetaInformation  ();
    /** Default Destructor **/
--- a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
+++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
@@ -133,6 +133,7 @@ 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,6 +56,9 @@ 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 =
@@ -113,26 +116,25 @@ template <typename TInputImage, typename TCoordRep>
 typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
 gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
 {
-    float     rayVector[3];
+    // --- Geometry / Siddon working variables use double ---
+    double    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;
-
+    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;

    // Min/max values of the output pixel type AND these values
    // represented as the output type of the interpolator
@@ -140,18 +142,12 @@ 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;
@@ -159,7 +155,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c

    drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);

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

    typename InputImageType::SizeType    sizeCT;
@@ -168,13 +164,12 @@ 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];
@@ -183,22 +178,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
-    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
-            ){
+    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)
+    {
        pixval = static_cast<OutputType>(0.0);
        return pixval;
    }
@@ -216,260 +211,351 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
    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 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 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];
+    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.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
-    //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.
+    // Compute ray from source to "detector" point pDest
    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]);

-    /* 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*/)
+    // just in case (should never really be zero)
+    if (rayLength <= EPSILON)
    {
-        //    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];
+        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];
        alphaXmin = std::min(alphaX1, alphaXN);
        alphaXmax = std::max(alphaX1, alphaXN);
    }
    else
    {
-        alphaXmin = -2;
-        alphaXmax = 2;
+        // 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;
    }

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

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

    /* 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
-  of the ray with the X, Y, and Z-planes after the ray entered the
-  CT volume. */
+       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];

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

-    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(alphaMin, alphaYmin)) // j is special
    {
-        alphaX = 2;
+        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;
    }
    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]);
+            (O[0] + firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+        alphaX  = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
+        alphaUx = ctPixelSpacing[0] / std::fabs(rayVector[0]);
    }

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

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

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

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

-    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. */
+    /* Initialize the current voxel index with the voxel entered at alphaMin. */
    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 */
+    // Current parametric position and previous one
+    double alpha     = alphaMin;
+    double alphaPrev = alphaMin;
+
+    // Siddon traversal
+    while (alpha < alphaMax)
    {
-        /* Store the current ray position */
-        alphaCminPrev = alphaCmin;
+        // Next plane intersection along the ray
+        double alphaNext = std::min(std::min(alphaX, alphaY), alphaZ);
+        if (alphaNext > alphaMax)
+            alphaNext = alphaMax;

-        if ((alphaX <= alphaY) && (alphaX <= alphaZ))
+        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)
        {
-            /* 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. */
+            // advance to the next plane and try again
+            if (nearlyEqual(alphaNext, alphaX))
            {
-                d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
+                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;
+        }
+
+        // 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]))
+        {
+            value = static_cast<float>(inputPtr->GetPixel(cIndex));
+
+            if (value > m_Threshold)
+            {
+                // physical path length = segmentLength (in α) * |rayVector|
+                d12 += segmentLength * rayLength * (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;
@@ -482,7 +568,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,7 +53,6 @@ vtkContourTopogramProjectionFilter::vtkContourTopogramProjectionFilter()
                0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
 }

-
 vtkContourTopogramProjectionFilter::~vtkContourTopogramProjectionFilter()
 {
    if (this->m_Transform)
@@ -111,7 +110,6 @@ 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;
    }
@@ -119,7 +117,6 @@ 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;
    }
@@ -127,7 +124,6 @@ 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;
    }
@@ -198,7 +194,6 @@ 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
@@ -0,0 +1,24 @@
+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
Giovanni Fattori	0e6933461f	fix to issue: Disproportionately large alpha in certain geometry cases	2026-01-12 22:35:44 +01:00
fattori_g	6009befeac	Update README.md	2025-06-25 00:11:48 +02:00
fattori_g	7a6d9abd8b	Update README.md	2025-06-25 00:09:20 +02:00
fattori_g	a003080e33	Update README.md	2025-06-25 00:01:28 +02:00
fattori_g	a16c48a204	Update README.md	2025-06-24 23:59:54 +02:00
fattori_g	dd94a1fb76	Update README.md	2025-06-24 23:57:26 +02:00
Giovanni Fattori	2361354982	A bit of code comments cleanup	2025-05-14 23:20:02 +02:00
fattori_g	06ad957615	Update itkReg23DRT/CMakeLists.txt	2025-05-14 23:01:10 +02:00
Giovanni Fattori	80eedce9cb	clear folder structure	2025-05-14 23:00:14 +02:00
Giovanni Fattori	7e9b1be191	clear folder structure	2025-05-14 22:56:29 +02:00
fattori_g	263da7276f	Update CITATION.cff	2025-05-14 11:02:27 +02:00
fattori_g	6899e563db	Add README.md	2025-05-14 11:01:21 +02:00
fattori_g	8ad613edee	Add third_party_licenses.txt	2025-05-14 10:58:38 +02:00
fattori_g	8703d4b49d	LICENSE.txt	2025-05-14 10:56:52 +02:00
fattori_g	c4f59078f0	CITATION.cff	2025-05-14 10:53:06 +02:00
fattori_g	88ad7118bc	Merge branch 'ScoutRTRelease' into 'master' RTSMarker bug See merge request cpt_bioeng/drt!63	2025-04-17 13:54:52 +02:00
via_r	7b9db5e25d	Merge branch 'ScalarRangeDRT' into 'ScoutRTRelease' Scalar Range of DRT for writing See merge request cpt_bioeng/drt!57	2024-09-23 09:20:46 +02:00
Proton local user	5196241a7b	Solved Issue #149 . Introduced new function for returning projection from imageProcessor with original Scalar Range which is used when saving the DRR. Visualisation is kept fixed (hard-coded upper limits: 0-2048).	2024-08-22 15:54:16 +02:00
Proton local user	affd24b48e	Working on Scalar Range issue. Intorudcing two different functions for writing and rendering Projections	2024-08-22 13:07:00 +02:00
Proton local user	4989d021bd	Strarting playing with Scalar Range of DRT	2024-08-21 09:18:27 +02:00
via_r	66762fac79	Merge branch 'OffRevRenderingAlternative' into 'ScoutRTRelease' DICOM Saving v2.0 See merge request cpt_bioeng/drt!56	2024-08-07 09:31:43 +02:00
via_r	c8b15d2d35	Merge branch 'DRRSavingAndErrorHandling' into 'OffRevRenderingAlternative' Unsigned short checks + DICOM Savnign error handling See merge request cpt_bioeng/drt!55	2024-08-07 09:25:25 +02:00
Proton local user	dd63d59841	Finished work on usigned short checks for DRR images as well as better handling and logging of errors in Data Saving	2024-07-23 14:54:26 +02:00
Proton local user	f15a6f5871	Cleaning unused variables here and there - variables are commented out or warning prevented using a unused_variable macro	2024-05-23 00:10:31 +02:00