commit c196ae3dbbbd9fd3509bb1cc883ea6acca056da6
Author: fattori_g <giovanni.fattori@psi.ch>
Date:   Fri Dec 17 10:08:47 2021 +0100

    Initial commit.
    Digitally Reconstructed Topogram library and simple test application.

diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt
new file mode 100644
index 0000000..bf4e5e2
--- /dev/null
+++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt
@@ -0,0 +1,40 @@
+cmake_minimum_required(VERSION 2.8.8)
+
+SET(LIB_NAME itkDTRrecon)
+
+find_package(ITK REQUIRED)
+include(${ITK_USE_FILE})
+INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
+
+find_package(VTK REQUIRED)
+INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS})
+
+FIND_PACKAGE(GDCM REQUIRED)
+include(${GDCM_USE_FILE})
+INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS})
+
+SET(SRCS
+    itkImageProcessor.cpp  
+)
+
+SET(HDR
+    itkImageProcessor.h  itkgSiddonJacobsRayCastInterpolateImageFunction.h  itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
+)
+
+ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR})
+
+SET(LINK_LIBS
+    ${LIB_NAME}
+    ${VTK_LIBRARIES}
+    ${ITK_LIBRARIES}
+    ${GDCM_LIBRARIES}
+)
+
+TARGET_LINK_LIBRARIES(
+    ${LINK_LIBS}
+)
+
+
+target_include_directories (${LIB_NAME}
+    PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+)
\ No newline at end of file
diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
new file mode 100644
index 0000000..ee03f7c
--- /dev/null
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
@@ -0,0 +1,1933 @@
+
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+
+#include "itkImageProcessor.h"
+
+#include "itkCreateObjectFunction.h"
+#include "itkVersion.h"
+#include <gdcmReader.h>
+#include <itkMatrix.h>
+#include "itkImageDuplicator.h"
+#include "itkPermuteAxesImageFilter.h"
+#include "itkMatrix.h"
+#include "itkGDCMSeriesFileNames.h"
+
+#include <string>
+
+namespace itk
+{
+
+
+static double coronalElements[9] = {
+    1, 0, 0,
+    0, 0, 1,
+    0,-1, 0 };
+
+static double sagittalElements[9] = {
+    0, 0,-1,
+    1, 0, 0,
+    0,-1, 0};
+
+static double HFS2IEC[9] = {
+    1, 0, 0,
+    0, 0, 1,
+    0, -1, 0};
+
+static double FFS2IEC[9] = {
+    -1, 0, 0,
+    0, 0, -1,
+    0, -1, 0};
+
+
+//FUNCTION DECLARATION NECESSARY FOR COPY/PASTE FROM vtkGDCMImageReader
+// const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds);
+
+const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds)
+{
+    static std::string buffer;
+    buffer = "";  // cleanup previous call
+    if( ds.FindDataElement( t ) )
+    {
+        const gdcm::DataElement& de = ds.GetDataElement( t );
+        const gdcm::ByteValue *bv = de.GetByteValue();
+        if( bv ) // Can be Type 2
+        {
+            buffer = std::string( bv->GetPointer(), bv->GetLength() );
+            // Will be padded with at least one \0
+        }
+    }
+    // Since return is a const char* the very first \0 will be considered
+    return buffer.c_str();
+};
+
+
+itkImageProcessor::itkImageProcessor()
+{
+
+    std::cout<<"itkImageProcessor::itkImageProcessor contructor"<<std::endl;
+    metric = MetricType::New();
+    transform = TransformType::New();
+    optimizer = OptimizerType::New();
+    interpolator1 = InterpolatorType::New();
+    interpolator2 = InterpolatorType::New();
+    registration = RegistrationType::New();
+
+    metric->ComputeGradientOff();
+    metric->SetSubtractMean(true);
+
+    registration->SetMetric(metric);
+    registration->SetOptimizer(optimizer);
+    registration->SetTransform(transform);
+    registration->SetInterpolator1(interpolator1);
+    registration->SetInterpolator2(interpolator2);
+
+    imageReader2D = ImageReaderType2D::New();
+    imageReader3D = ImageReaderType3D::New();
+    imageSeriesReader3D = ImageSeriesReaderType3D::New();
+
+    scd = 1000.; // Source to isocenter distance
+    threshold = 0.;
+    projAngle1 = -999.;
+    projAngle2 = -999.;
+
+    isocenter[0]=isocenter[1]=isocenter[2] = 0.;
+    image1res[0]=image2res[1]= 2.;
+    image1center[0]=image1center[1]= 0.;
+    image2center[0]=image2center[1]= 0.;
+    image1Size[0] = image1Size[1] = 512;
+    image2Size[0] = image2Size[1] = 512;
+    image3DCOV[0] = image3DCOV[1] = image3DCOV[2] = 0.;
+
+    image1IntensityWindow[0] = image1IntensityWindow[1] = 0.;
+    image2IntensityWindow[0] = image2IntensityWindow[1] = 0.;
+
+    TZero[0]=TZero[1]=TZero[2]=0.;
+    RZero[0]=RZero[1]=RZero[2]=0.;
+    customized_iso = false;
+    customized_2DCX = false;
+    customized_2DRES = false;
+    customized_2DSIZE = false;
+
+    image2D1Loaded= false;
+    image2D2Loaded= false;
+    image3DLoaded= false;
+
+    flipFilter1 = FlipFilterType::New();
+    flipFilter2 = FlipFilterType::New();
+    flipFilterLocalizer1 = FlipFilterType::New();
+    flipFilterLocalizer2 = FlipFilterType::New();
+
+    using FlipAxesArrayType = FlipFilterType::FlipAxesArrayType;
+    FlipAxesArrayType flipArray;
+    flipArray[0] = true;
+    flipArray[1] = true;
+    flipArray[2] = false;
+    flipFilter1->SetFlipAxes(flipArray);
+    //flipFilterLocalizer1->SetFlipAxes(flipArray);
+    flipArray[0] = false;
+    flipArray[1] = true;
+    flipArray[2] = false;
+    flipFilter2->SetFlipAxes(flipArray);
+    //flipFilterLocalizer2->SetFlipAxes(flipArray);
+
+    caster3D = CastFilterType3D::New();
+    caster2D1 = CastFilterType3D::New();
+    caster2D2 = CastFilterType3D::New();
+    caster2DInput = CastFilterType3D::New();
+
+    toVTK2D1 = ITKtoVTKFilterType::New();
+    toVTK2D2 = ITKtoVTKFilterType::New();
+    toVTKLocalizer1 = ITKtoVTKFilterType::New();
+    toVTKLocalizer2 = ITKtoVTKFilterType::New();
+
+    image2D1Rescaler = Input2DRescaleFilterType::New();
+    image2D2Rescaler = Input2DRescaleFilterType::New();
+
+
+    image3DIECFilt   = ChangeInformationFilterInputType::New();
+
+}
+
+itkImageProcessor::~itkImageProcessor()
+{
+
+}
+
+void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
+{
+    Superclass::PrintSelf(os, indent);
+}
+
+void itkImageProcessor::SetIntensityThreshold(double dT)
+{
+    threshold = dT;
+}
+
+void itkImageProcessor::SetSCD(double dDist)
+{
+    scd = dDist;
+}
+
+void itkImageProcessor::SetCustom_Isocenter(double dX,double dY,double dZ)
+{
+    isocenter[0] = dX;
+    isocenter[1] = dY;
+    isocenter[2] = dZ;
+    customized_iso = true;
+}
+
+void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double nY2)
+{
+    image1res[0] = nX1;
+    image1res[1] = nY1;
+    image2res[0] = nX2;
+    image2res[1] = nY2;
+    customized_2DRES = true;
+}
+
+void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
+{
+    image1Size[0] = nX1;
+    image1Size[1] = nY1;
+    image2Size[0] = nX2;
+    image2Size[1] = nY2;
+    customized_2DSIZE = true;
+}
+
+void itkImageProcessor::SetCustom_2Dcenter(double dX1,double dY1, double dX2,double dY2)
+{
+    image1center[0] = dX1;
+    image1center[1] = dY1;
+    image2center[0] = dX2;
+    image2center[1] = dY2;
+    customized_2DCX = true;
+}
+
+int itkImageProcessor::load3DSerie(const char * pcDirName)
+{
+
+    std::cout<<"itkImageProcessor::load3DSerie"<<std::endl;
+
+    using NamesGeneratorType = itk::GDCMSeriesFileNames;
+    NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
+
+    nameGenerator->SetUseSeriesDetails(true);
+    nameGenerator->AddSeriesRestriction("0008|0021");
+    nameGenerator->SetGlobalWarningDisplay(false);
+    nameGenerator->SetDirectory(pcDirName);
+
+    try
+    {
+        using SeriesIdContainer = std::vector<std::string>;
+        const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();
+        auto                      seriesItr = seriesUID.begin();
+        auto                      seriesEnd = seriesUID.end();
+
+        if (seriesItr != seriesEnd)
+        {
+            std::cout << "The directory: ";
+            std::cout << pcDirName << std::endl;
+            std::cout << "Contains the following DICOM Series: ";
+            std::cout << std::endl;
+        }
+        else
+        {
+            std::cout << "No DICOMs in: " << pcDirName << std::endl;
+            return EXIT_SUCCESS;
+        }
+
+        while (seriesItr != seriesEnd)
+        {
+            std::cout << seriesItr->c_str() << std::endl;
+            ++seriesItr;
+        }
+
+        /* we expect only one serie at this point... */
+        if(seriesUID.size() != 1){
+            std::cout<<"More than one serie in the selected folder. returning..."<<std::endl;
+            return EXIT_FAILURE;
+        }
+
+
+
+        seriesItr = seriesUID.begin();
+        while (seriesItr != seriesUID.end())
+        {
+          std::string seriesIdentifier;
+          seriesIdentifier = seriesItr->c_str();
+          seriesItr++;
+
+          std::cout << "\nReading: ";
+          std::cout << seriesIdentifier << std::endl;
+          using FileNamesContainer = std::vector<std::string>;
+          FileNamesContainer fileNames = nameGenerator->GetFileNames(seriesIdentifier);
+
+          using ImageIOType = itk::GDCMImageIO;
+          ImageIOType::Pointer dicomIO = ImageIOType::New();
+          imageSeriesReader3D->SetImageIO(dicomIO);
+          imageSeriesReader3D->SetFileNames(fileNames);
+          imageSeriesReader3D->SetNumberOfWorkUnits(8);
+          imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt
+
+          imageSeriesReader3D->Update();
+          image3DIECFilt->SetInput(imageSeriesReader3D->GetOutput());
+          //image3DIn = imageSeriesReader3D->GetOutput();
+
+          eImageOrientationType
+                currImgOrient;
+          /* check patient orientation */
+          if(fileNames.size()>0){
+
+              gdcm::Reader R;
+              R.SetFileName(fileNames.at(0).c_str());
+              if(!R.Read())
+              {
+                  cerr<<"ERROR: cannot read file: "<<fileNames.at(0).c_str()<<endl;
+                  return -1;
+              }
+              const gdcm::File &file = R.GetFile();
+              const gdcm::DataSet &ds = file.GetDataSet();
+
+              char* sTmpString = new char [255];
+              strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0018,0x5100), ds));
+              *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;
+              }
+
+              free (sTmpString);
+          }
+
+          /* transform everything to IEC */
+
+          itk::Matrix<float, 3, 3> toIECMatrix;
+          toIECMatrix.Fill(0.0);
+          InternalImageType::DirectionType toIECDirection;
+
+          for(int iIdx = 0 ; iIdx < 3; iIdx++){
+              for(int jIdx = 0 ; jIdx < 3; jIdx++){
+
+                  switch (currImgOrient) {
+                    case eImageOrientationType::HFS:
+                      toIECDirection.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3];
+                      break;
+                  case eImageOrientationType::FFS:
+                      toIECDirection.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3];
+                      break;
+                  default:
+                      break;
+                  }
+              }
+          }
+
+
+//          std::cout<< "Orignal Direction: "
+//                   << imageSeriesReader3D->GetOutput()->GetDirection()<<std::endl;
+//          std::cout<< "Correctrion Matrix: "
+//                   << toIECDirection<<std::endl;
+//          std::cout<< "Final Direction: "
+//                   << imageSeriesReader3D->GetOutput()->GetDirection() * toIECDirection<<std::endl;
+//          image3DIECFilt->SetOutputDirection(
+//                imageSeriesReader3D->GetOutput()->GetDirection() * toIECDirection
+//                );
+//          image3DIECFilt->ChangeDirectionOn();
+
+          image3DIECFilt->Update();
+          image3DIn = image3DIECFilt->GetOutput();
+
+          /* before changing origin for projection, let's save the center of volume */
+          using ImageRegionType3D = ImageType3D::RegionType;
+          using SizeType3D = ImageRegionType3D::SizeType;
+          using ImageDirectionType3D = ImageType3D::DirectionType;
+
+          ImageRegionType3D region3D = image3DIn->GetBufferedRegion();
+          SizeType3D        size3D = region3D.GetSize();
+          ImageType3D::PointType       origin3D = image3DIn->GetOrigin();
+          const itk::Vector<double, 3> resolution3D = image3DIn->GetSpacing();
+          ImageDirectionType3D imagDirection = image3DIn->GetDirection();
+
+          /* calculate image size in 3D Space by fiding the last voxel position */
+          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);
+
+          image3DCOV[0] = (origin3D[0]+LastVoxelPosition[0])/2;
+          image3DCOV[1] = (origin3D[1]+LastVoxelPosition[1])/2;
+          image3DCOV[2] = (origin3D[2]+LastVoxelPosition[2])/2;
+
+          if(true){
+              std::cout<<" -------- 3D IMAGE --------"<<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;
+              std::cout<<"Image 3D COV: "<<image3DCOV[0]<<" "
+                                      <<image3DCOV[1]<<" "
+                                      <<image3DCOV[2]<<" "<<std::endl;
+              std::cout<<" -------- -------- --------"<<std::endl;
+
+          }
+
+
+          // To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image
+          // to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
+          // geometry is fully defined. The origin of the CT image becomes irrelavent.
+          ImageType3D::PointType image3DOrigin;
+          image3DOrigin[0] = 0.0;
+          image3DOrigin[1] = 0.0;
+          image3DOrigin[2] = 0.0;
+          image3DIn->SetOrigin(image3DOrigin);
+
+
+          caster3D->SetInput(image3DIn);
+          caster3D->Update();
+
+
+        }
+    }
+    catch (itk::ExceptionObject & ex)
+    {
+        std::cout << ex << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    return EXIT_SUCCESS;
+
+}
+
+
+void itkImageProcessor::load3DVolume(const char *pcFName)
+{
+
+
+    std::cout<<"itkImageProcessor::load3DVolume"<<std::endl;
+
+    imageReader3D->SetFileName(pcFName);
+    imageReader3D->Update();
+
+    image3DIn = imageReader3D->GetOutput();
+
+
+    /* before changing origin for projection, let's save the center of volume */
+    using ImageRegionType3D = ImageType3D::RegionType;
+    using SizeType3D = ImageRegionType3D::SizeType;
+    using ImageDirectionType3D = ImageType3D::DirectionType;
+
+    ImageRegionType3D region3D = image3DIn->GetBufferedRegion();
+    SizeType3D        size3D = region3D.GetSize();
+    ImageType3D::PointType       origin3D = image3DIn->GetOrigin();
+    const itk::Vector<double, 3> resolution3D = image3DIn->GetSpacing();
+    ImageDirectionType3D imagDirection = image3DIn->GetDirection();
+
+    /* calculate image size in 3D Space by fiding the last voxel position */
+    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);
+
+    image3DCOV[0] = (origin3D[0]+LastVoxelPosition[0])/2;
+    image3DCOV[1] = (origin3D[1]+LastVoxelPosition[1])/2;
+    image3DCOV[2] = (origin3D[2]+LastVoxelPosition[2])/2;
+
+    if(true){
+        std::cout<<" -------- 3D IMAGE --------"<<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;
+        std::cout<<"Image 3D COV: "<<image3DCOV[0]<<" "
+                                <<image3DCOV[1]<<" "
+                                <<image3DCOV[2]<<" "<<std::endl;
+        std::cout<<" -------- -------- --------"<<std::endl;
+
+    }
+
+
+    // To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image
+    // to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
+    // geometry is fully defined. The origin of the CT image becomes irrelavent.
+    ImageType3D::PointType image3DOrigin;
+    image3DOrigin[0] = 0.0;
+    image3DOrigin[1] = 0.0;
+    image3DOrigin[2] = 0.0;
+    image3DIn->SetOrigin(image3DOrigin);
+
+
+    caster3D->SetInput(image3DIn);
+    caster3D->Update();
+
+
+}
+
+void itkImageProcessor::SetProjectionAngle1(double dGantryAngle) 
+{
+    projAngle1 = dGantryAngle;
+}
+
+void itkImageProcessor::SetProjectionAngle2(double dGantryAngle)
+{
+    projAngle2 = dGantryAngle;
+}
+
+int itkImageProcessor::load2D(const char * pcFName){
+
+    /* Check if we can open the file */
+    gdcm::Reader R;
+    R.SetFileName(pcFName);
+    if(!R.Read())
+    { cerr<<"ERROR: cannot read file: "<<pcFName<<endl;
+        return -1;
+    }
+
+    /* Check if it's a localizer image */
+    const gdcm::File &file = R.GetFile();
+    const gdcm::DataSet &ds = file.GetDataSet();
+
+    char* sTmpString = new char [255];
+    strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0008,0x0008), ds));
+    std::string sLocalizerString = "LOCALIZER";
+    if (std::string(sTmpString).find(sLocalizerString) != std::string::npos) {
+        std::cout << "Loading Localizer: "<<pcFName << '\n';
+    } else {
+        cerr << "ERROR: the image "<<pcFName <<" is not a LOCALIZER" << '\n';
+        free (sTmpString);
+        return -1;
+    }
+
+    /* read some useful tag values */
+    double
+            dIntensityWindow[2];
+    std::vector<std::string> tokens;
+
+
+    /* Intensity window center */
+    strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0028,0x1050), ds));
+    for (auto i = strtok(&sTmpString[0], "\\"); i != NULL; i = strtok(NULL, "\\"))
+        tokens.push_back(i);
+
+    if(tokens.size()>0){
+        dIntensityWindow[0] = std::atof(tokens.at(0).c_str() );
+    } else {
+        dIntensityWindow[0] = 0.;
+    }
+    tokens.clear();
+
+    /* Intensity window width */
+    strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0028,0x1051), ds));
+    for (auto i = strtok(&sTmpString[0], "\\"); i != NULL; i = strtok(NULL, "\\"))
+        tokens.push_back(i);
+
+    if(tokens.size()>0){
+        dIntensityWindow[1] = std::atof(tokens.at(0).c_str() );
+    } else {
+        dIntensityWindow[1] = 0.;
+    }
+    tokens.clear();
+
+    std::cout<<"Intensity window center and width: "<< dIntensityWindow[0] <<" "<< dIntensityWindow[1] <<std::endl;
+
+    /* Image orientation */
+    strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0018,0x5100), ds));
+    *std::remove(sTmpString, sTmpString + strlen(sTmpString), ' ') = 0;
+
+    if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
+        std::cout<<"Image orientation: HFS"<<std::endl;
+    } else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
+        std::cout<<"Image orientation: FFS"<<std::endl;
+    } else {
+        std::cerr<< "Image Orientation: unrecognised"<< sTmpString <<std::endl;
+    }
+
+    free (sTmpString);
+
+    /* Actual file read */
+    ImageIOType::Pointer dicomIO = ImageIOType::New();
+    imageReader2D->SetFileName(pcFName);
+    imageReader2D->SetImageIO(dicomIO);
+    try
+    {
+        imageReader2D->Update();
+    }
+    catch (const itk::ExceptionObject & ex)
+    {
+        std::cout << ex << std::endl;
+        return -1;
+    }
+
+    using ImageRegionType3D = ImageType3D::RegionType;
+    using SizeType3D = ImageRegionType3D::SizeType;
+    using ImageDirectionType3D = ImageType3D::DirectionType;
+
+
+    ImageRegionType3D region3D = imageReader2D->GetOutput()->GetBufferedRegion();
+    SizeType3D        size3D = region3D.GetSize();
+    ImageType3D::PointType       origin3D = imageReader2D->GetOutput()->GetOrigin();
+    const itk::Vector<double, 3> resolution3D = imageReader2D->GetOutput()->GetSpacing();
+    ImageDirectionType3D imagDirection = imageReader2D->GetOutput()->GetDirection();
+
+    /* Identify Coronal or Sagittal image plane */
+    double dSumCoronal = 0;
+    for(int iIdx = 0 ; iIdx < 3; iIdx++){
+        for(int jIdx = 0 ; jIdx < 3; jIdx++){
+            dSumCoronal += fabs(imagDirection[iIdx][jIdx] - coronalElements[jIdx+iIdx*3]);
+        }
+    }
+    double dSumSagittal = 0;
+    for(int iIdx = 0 ; iIdx < 3; iIdx++){
+        for(int jIdx = 0 ; jIdx < 3; jIdx++){
+            dSumSagittal += fabs(imagDirection[iIdx][jIdx] - sagittalElements[jIdx+iIdx*3]);
+        }
+    }
+
+    /* 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);
+
+    if(false){
+        std::cout<<" -------- LOAD OF 2D IMAGE --------"<<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;
+    }
+
+    ImageType2D::Pointer image_tmp1 = imageReader2D->GetOutput();
+
+    //    for(int iIdx =0 ; iIdx <100; iIdx++){
+    //        image_tmp1.GetPointer()->GetBufferPointer()[iIdx] = 2000;
+    //    }
+
+    /* cast type to internal processing type */
+    caster2DInput->SetInput(image_tmp1);
+    caster2DInput->Update();
+
+
+    using DuplicatorType = itk::ImageDuplicator<InternalImageType>;
+    DuplicatorType::Pointer duplicator = DuplicatorType::New();
+    duplicator->SetInputImage(caster2DInput->GetOutput());
+    duplicator->Update();
+
+    ImageType3D::PointType image3DOrigin;
+    image3DOrigin[0]=image3DCOV[0];
+    image3DOrigin[1]=image3DCOV[1];
+    image3DOrigin[2]=image3DCOV[2];
+
+    if(dSumCoronal < 1e-10){
+        //std::cout<< "Coronal" <<std::endl;
+
+        image2D1In = duplicator->GetOutput();
+
+        const InternalImageType::DirectionType image2D1InDirection = image2D1In->GetDirection();
+        //        std::cout<<"image2D1InDirection "<<image2D1InDirection<<std::endl;
+
+        /* center topogram at the COV of 3D image
+           this works only if user does not set custom isocenter
+           to be fixed.
+
+           This is the distance between the topogram origin and the 3DVolume center (which
+           are given in DICOM FOR) mapped into the topogram frame of reference
+
+        */
+        ImageType3D::PointType image2Dcenter_3DCOV;
+        InternalImageType::DirectionType image2D1InDirectionInv;
+        image2D1InDirectionInv.GetVnlMatrix() = image2D1InDirection.GetTranspose();
+
+        image2Dcenter_3DCOV = -
+                (
+                    image2D1InDirectionInv * (image3DOrigin - image2D1In->GetOrigin())
+                    );
+
+        /*
+            We have to set the z coordinate to the focal distance, it would be 0 otherwise...
+        */
+        image2Dcenter_3DCOV[2] = -scd;
+
+        ChangeInformationFilterType::Pointer filter = ChangeInformationFilterType::New();
+        filter->SetInput(image2D1In);
+
+        const InternalImageType::DirectionType newDirection = image2D1InDirection * image2D1InDirectionInv;//direction * rotation.GetMatrix();
+
+        filter->SetOutputDirection(newDirection);
+        filter->ChangeDirectionOn();
+
+        filter->SetOutputOrigin(image2Dcenter_3DCOV);
+        filter->ChangeOriginOn();
+
+        try
+        {
+            filter->UpdateOutputInformation();
+        }
+        catch (itk::ExceptionObject & error)
+        {
+            std::cerr << "Error: " << error << std::endl;
+            return EXIT_FAILURE;
+        }
+
+
+        memcpy(image1IntensityWindow,dIntensityWindow, 2*sizeof(double));
+
+
+        flipFilterLocalizer1->SetInput(filter->GetOutput());
+        flipFilterLocalizer1->Update();
+
+        image2D1Rescaler->SetInput(flipFilterLocalizer1->GetOutput());
+        image2D1Rescaler->SetOutputMinimum(0);
+        image2D1Rescaler->SetOutputMaximum(255);
+
+
+        image2D1Loaded = true;
+        return 1;
+
+    } else if (dSumSagittal < 1e-10){
+
+        //std::cout<< "Sagittal" <<std::endl;
+        image2D2In = duplicator->GetOutput();
+
+        const InternalImageType::DirectionType image2D2InDirection = image2D2In->GetDirection();
+
+        /* center topogram at the COV of 3D image
+           this works only if user does not set custom isocenter
+           to be fixed.
+
+           This is the distance between the topogram origin and the 3DVolume center (which
+           are given in DICOM FOR) mapped into the topogram frame of reference
+
+        */
+        ImageType3D::PointType image2Dcenter_3DCOV;
+        InternalImageType::DirectionType image2D2InDirectionInv;
+        image2D2InDirectionInv.GetVnlMatrix() = image2D2InDirection.GetTranspose();
+
+        image2Dcenter_3DCOV = -
+                (
+                    image2D2InDirectionInv * (image3DOrigin - image2D2In->GetOrigin())
+                    );
+
+        /*
+            We have to set the z coordinate to the focal distance, it would be 0 otherwise...
+        */
+
+        image2Dcenter_3DCOV[2] = -scd;
+
+
+        using FilterType = itk::ChangeInformationImageFilter<InternalImageType>;
+        FilterType::Pointer filter = FilterType::New();
+        filter->SetInput(image2D2In);
+
+        const InternalImageType::DirectionType newDirection = image2D2InDirection * image2D2InDirectionInv;//image2D2InDirectionInv;
+        filter->SetOutputDirection(newDirection);
+        filter->ChangeDirectionOn();
+
+        filter->SetOutputOrigin(image2Dcenter_3DCOV);
+        filter->ChangeOriginOn();
+
+        try
+        {
+            filter->UpdateOutputInformation();
+        }
+        catch (itk::ExceptionObject & error)
+        {
+            std::cerr << "Error: " << error << std::endl;
+            return EXIT_FAILURE;
+        }
+
+        memcpy(image2IntensityWindow,dIntensityWindow, 2*sizeof(double));
+
+
+        flipFilterLocalizer2->SetInput(filter->GetOutput());
+        flipFilterLocalizer2->Update();
+
+        image2D2Rescaler->SetInput(flipFilterLocalizer2->GetOutput());
+        image2D2Rescaler->SetOutputMinimum(0);
+        image2D2Rescaler->SetOutputMaximum(255);
+
+
+        image2D2Loaded = true;
+        return 2;
+    } else {
+        std::cout<< "Not Accepted" <<std::endl;
+        return -1;
+    };
+
+}
+
+
+double* itkImageProcessor::GetImageIntensityWindow(int iImg)
+{
+
+    switch (iImg) {
+
+    case 0:
+        return image1IntensityWindow;
+        break;
+
+    case 1:
+        return image2IntensityWindow;
+        break;
+
+    default:
+        return nullptr;
+        break;
+
+    }
+
+}
+
+
+void itkImageProcessor::InitializeProjector()
+{
+    caster3D->Update();
+    transform->SetComputeZYX(true);
+    TransformType::OutputVectorType translation;
+    translation[0] = TZero[0];
+    translation[1] = TZero[1];
+    translation[2] = TZero[2];
+
+    transform->SetTranslation(translation);
+    const double dtr = (atan(1.0) * 4.0) / 180.0;
+    transform->SetRotation(dtr * RZero[0], dtr * RZero[1], dtr * RZero[2]);
+
+    // The centre of rotation is set by default to the centre of the 3D
+    // volume but can be offset from this position using a command
+    // line specified translation [cx,cy,cz]
+
+    ImageType3D::PointType       origin3D = caster3D->GetOutput()->GetOrigin();
+    const itk::Vector<double, 3> resolution3D = caster3D->GetOutput()->GetSpacing();
+
+    using ImageRegionType3D = ImageType3D::RegionType;
+    using SizeType3D = ImageRegionType3D::SizeType;
+
+    ImageRegionType3D region3D = caster3D->GetOutput()->GetBufferedRegion();
+    SizeType3D        size3D = region3D.GetSize();
+
+    if (customized_iso)
+    {
+        // Isocenter location given by the user.
+        isocenter3D[0] = origin3D[0] + resolution3D[0] * isocenter[0];
+        isocenter3D[1] = origin3D[1] + resolution3D[1] * isocenter[1];
+        isocenter3D[2] = origin3D[2] + resolution3D[2] * isocenter[2];
+    }
+    else
+    {
+        // Set the center of the image as the isocenter.
+        isocenter3D[0] = origin3D[0] + resolution3D[0] * static_cast<double>(size3D[0]) / 2.0;
+        isocenter3D[1] = origin3D[1] + resolution3D[1] * static_cast<double>(size3D[1]) / 2.0;
+        isocenter3D[2] = origin3D[2] + resolution3D[2] * static_cast<double>(size3D[2]) / 2.0;
+    }
+
+    std::cout<<" -------- PROJECTOR GEOMETRY --------"<<std::endl;
+    std::cout<<"Projector Isocenter "<<isocenter3D[0]<<" "<<isocenter3D[1]<<" "<<isocenter3D[2] <<std::endl;
+    std::cout<<"Projector scd "<<scd <<std::endl;
+    std::cout<<" -------- -------- --------"<<std::endl;
+
+    transform->SetCenter(isocenter3D);
+    // 2D Image 1
+    interpolator1->SetProjectionAngle(dtr * projAngle1);
+    interpolator1->SetFocalPointToIsocenterDistance(scd);
+    interpolator1->SetThreshold(threshold);
+    interpolator1->SetTransform(transform);
+
+    interpolator1->Initialize();
+
+    // 2D Image 2
+    interpolator2->SetProjectionAngle(dtr * projAngle2);
+    interpolator2->SetFocalPointToIsocenterDistance(scd);
+    interpolator2->SetThreshold(threshold);
+    interpolator2->SetTransform(transform);
+
+    interpolator2->Initialize();
+
+}
+
+
+void itkImageProcessor::GetProjectionImages(){
+
+    //    std::cout<< "itkImageProcessor::GetProjectionImages" <<std::endl;
+
+    const double dtr = (atan(1.0) * 4.0) / 180.0;
+    TransformType::Pointer finalTransform = TransformType::New();
+
+    finalTransform->SetComputeZYX(true);
+    TransformType::OutputVectorType translation;
+    translation[0] = TZero[0];
+    translation[1] = TZero[1];
+    translation[2] = TZero[2];
+
+    finalTransform->SetTranslation(translation);
+    finalTransform->SetRotation(dtr * RZero[0], dtr * RZero[1], dtr * RZero[2]);
+
+    // The transform is determined by the parameters and the rotation center.
+    finalTransform->SetCenter(isocenter3D);
+
+    using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
+
+    // The ResampleImageFilter is the driving force for the projection image generation.
+    ResampleFilterType::Pointer resampleFilter1 = ResampleFilterType::New();
+    resampleFilter1->SetInput(caster3D->GetOutput()); // Link the 3D volume.
+    resampleFilter1->SetDefaultPixelValue(0);
+
+    // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
+    // have been set before registration. Here we only need to replace the initial
+    // transform with the final transform.
+    interpolator1->SetTransform(finalTransform);
+    interpolator1->Initialize();
+    resampleFilter1->SetInterpolator(interpolator1);
+
+    // The output 2D projection image has the same image size, origin, and the pixel spacing as
+    // those of the input 2D image.
+    InternalImageType::SizeType size;
+    double spacing[Dimension];
+    double o2Dx, o2Dy, origin[Dimension];
+
+    if(false){
+        //    if(image2D1Loaded){
+        resampleFilter1->SetSize(image2D1In->GetLargestPossibleRegion().GetSize());
+        resampleFilter1->SetOutputOrigin(image2D1In->GetOrigin());
+        resampleFilter1->SetOutputSpacing(image2D1In->GetSpacing());
+
+        //        std::cout<<"Size "<<caster2D1->GetOutput()->GetLargestPossibleRegion().GetSize()<<std::endl;
+        //        std::cout<<"Origin "<<caster2D1->GetOutput()->GetOrigin()<<std::endl;
+        //        std::cout<<"Spacing "<<caster2D1->GetOutput()->GetSpacing()<<std::endl;
+    } else {
+
+        size[0] = image1Size[0];       // number of pixels along X of the 2D DRR image
+        size[1] = image1Size[1];       // number of pixels along X of the 2D DRR image
+        size[2] = 1;        // only one slice
+        spacing[0] = image1res[0]; // pixel spacing along X of the 2D DRR image [mm]
+        spacing[1] = image1res[1]; // pixel spacing along Y of the 2D DRR image [mm]
+        spacing[2] = 1.0;   // slice thickness of the 2D DRR image [mm]
+
+        o2Dx = ((double)image1Size[0] - 1.) / 2.;
+        o2Dy = ((double)image1Size[1] - 1.) / 2.;
+        // Compute the origin (in mm) of the 2D Image
+        origin[0] = -image1res[0] * o2Dx;
+        origin[1] = -image1res[1] * o2Dy;
+        origin[2] = -scd;
+
+        resampleFilter1->SetSize(size);
+        resampleFilter1->SetOutputSpacing(spacing);
+        resampleFilter1->SetOutputOrigin(origin);
+
+    }
+
+
+    // Do the same thing for the output image 2.
+    ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New();
+    resampleFilter2->SetInput(caster3D->GetOutput());
+    resampleFilter2->SetDefaultPixelValue(0);
+
+    // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
+    // have been set before registration. Here we only need to replace the initial
+    // transform with the final transform.
+    interpolator2->SetTransform(finalTransform);
+    interpolator2->Initialize();
+    resampleFilter2->SetInterpolator(interpolator2);
+
+    if(false){
+        //    if(image2D2Loaded){
+        resampleFilter2->SetSize(image2D2In->GetLargestPossibleRegion().GetSize());
+        resampleFilter2->SetOutputOrigin(image2D2In->GetOrigin());
+        resampleFilter2->SetOutputSpacing(image2D2In->GetSpacing());
+
+    } else {
+
+        size[0] = image2Size[0];       // number of pixels along X of the 2D DRR image
+        size[1] = image2Size[1];       // number of pixels along X of the 2D DRR image
+        size[2] = 1;        // only one slice
+        spacing[0] = image2res[0]; // pixel spacing along X of the 2D DRR image [mm]
+        spacing[1] = image2res[1]; // pixel spacing along Y of the 2D DRR image [mm]
+        spacing[2] = 1.0;   // slice thickness of the 2D DRR image [mm]
+
+        o2Dx = ((double)image2Size[0] - 1.) / 2.;
+        o2Dy = ((double)image2Size[1] - 1.) / 2.;
+        // Compute the origin (in mm) of the 2D Image
+        origin[0] = -image2res[0] * o2Dx;
+        origin[1] = -image2res[1] * o2Dy;
+        origin[2] = -scd;
+
+        resampleFilter2->SetSize(size);
+        resampleFilter2->SetOutputSpacing(spacing);
+        resampleFilter2->SetOutputOrigin(origin);
+
+    }
+
+
+    // As explained before, the computed projection is upsided-down.
+    // Here we use a FilpImageFilter to flip the images in y-direction.
+    flipFilter1->SetInput(resampleFilter1->GetOutput());
+    flipFilter2->SetInput(resampleFilter2->GetOutput());
+
+    flipFilter1->Update();
+    flipFilter2->Update();
+
+    //    std::cout<< "itkImageProcessor::GetProjectionImages DONE" <<std::endl;
+}
+
+
+void itkImageProcessor::Write2DImages(){
+
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+    using WriterType = itk::ImageFileWriter<OutputImageType>;
+
+
+    if(image2D1Loaded)
+    {
+
+        // Rescale the intensity of the projection images to 0-255 for output.
+        RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+        rescaler1->SetOutputMinimum(0);
+        rescaler1->SetOutputMaximum(255);
+        rescaler1->SetInput(image2D1In);
+
+        WriterType::Pointer writer1 = WriterType::New();
+
+        writer1->SetFileName("/mnt/gdrive/Scratch/gfattori/2D1.tif");
+        writer1->SetInput(rescaler1->GetOutput());
+
+        try
+        {
+            std::cout << "Writing image 2D1" << std::endl;
+            writer1->Update();
+        }
+        catch (itk::ExceptionObject & err)
+        {
+            std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+            std::cerr << err << std::endl;
+        }
+
+    }
+
+    if(image2D2Loaded){
+        // Rescale the intensity of the projection images to 0-255 for output.
+        RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
+        rescaler2->SetOutputMinimum(0);
+        rescaler2->SetOutputMaximum(255);
+        rescaler2->SetInput(image2D2In);
+
+        WriterType::Pointer writer2 = WriterType::New();
+
+        writer2->SetFileName("/mnt/gdrive/Scratch/gfattori/2D2.tif");
+        writer2->SetInput(rescaler2->GetOutput());
+
+        try
+        {
+            std::cout << "Writing image 2D2" << std::endl;
+            writer2->Update();
+        }
+        catch (itk::ExceptionObject & err)
+        {
+            std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+            std::cerr << err << std::endl;
+        }
+    }
+}
+
+//vtkImageData* itkImageProcessor::GetLocalizer1VTK()
+//{
+////    // Rescale the intensity of the projection images to 0-255 for output.
+////    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+////    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+////    rescaler1->SetOutputMinimum(0);
+////    rescaler1->SetOutputMaximum(255);
+////    rescaler1->SetInput(flipFilterLocalizer1->GetOutput());
+
+////    rescaler1->Update();
+//    using ImageRegionType3D = ImageType3D::RegionType;
+//    using SizeType3D = ImageRegionType3D::SizeType;
+//    using ImageDirectionType3D = ImageType3D::DirectionType;
+
+//    ImageRegionType3D region3D = image2D1In->GetBufferedRegion();
+//    SizeType3D        size3D = region3D.GetSize();
+//    ImageType3D::PointType       origin3D = image2D1In->GetOrigin();
+//    const itk::Vector<double, 3> resolution3D = image2D1In->GetSpacing();
+//    ImageDirectionType3D imagDirection = image2D1In->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;
+
+//    toVTKLocalizer1->SetInput("image2D1In",image2D1In);
+//    toVTKLocalizer1->Update();
+
+//    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();
+//}
+
+//vtkImageData* itkImageProcessor::GetLocalizer2VTK()
+//{
+////    // Rescale the intensity of the projection images to 0-255 for output.
+////    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+////    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
+////    rescaler2->SetOutputMinimum(0);
+////    rescaler2->SetOutputMaximum(255);
+////    rescaler2->SetInput(flipFilterLocalizer2->GetOutput());
+
+////    rescaler2->Update();
+//    using ImageRegionType3D = ImageType3D::RegionType;
+//    using SizeType3D = ImageRegionType3D::SizeType;
+//    using ImageDirectionType3D = ImageType3D::DirectionType;
+
+//    ImageRegionType3D region3D = image2D2In->GetBufferedRegion();
+//    SizeType3D        size3D = region3D.GetSize();
+//    ImageType3D::PointType       origin3D = image2D2In->GetOrigin();
+//    const itk::Vector<double, 3> resolution3D = image2D2In->GetSpacing();
+//    ImageDirectionType3D imagDirection = image2D2In->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;
+
+
+
+//    toVTKLocalizer2->SetInput("image2D2In",image2D2In);
+//    toVTKLocalizer2->Update();
+
+//    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::GetLocalizer1VTK()
+{
+    // Rescale the intensity of the projection images to 0-255 for output.
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+    rescaler1->SetOutputMinimum(0);
+    rescaler1->SetOutputMaximum(255);
+    rescaler1->SetInput(flipFilterLocalizer1->GetOutput());
+
+    rescaler1->Update();
+    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;
+
+    toVTKLocalizer1->SetInput(rescaler1->GetOutput());
+    toVTKLocalizer1->Update();
+
+    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();
+}
+
+vtkImageData* itkImageProcessor::GetLocalizer2VTK()
+{
+    // Rescale the intensity of the projection images to 0-255 for output.
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
+    rescaler2->SetOutputMinimum(0);
+    rescaler2->SetOutputMaximum(255);
+    rescaler2->SetInput(flipFilterLocalizer2->GetOutput());
+
+    rescaler2->Update();
+    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;
+
+
+
+    toVTKLocalizer2->SetInput(rescaler2->GetOutput());
+    toVTKLocalizer2->Update();
+
+    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()
+{
+    // Rescale the intensity of the projection images to 0-255 for output.
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+    rescaler1->SetOutputMinimum(0);
+    rescaler1->SetOutputMaximum(255);
+    rescaler1->SetInput(flipFilter1->GetOutput());
+
+    rescaler1->Update();
+    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<<" -------- 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;
+
+
+
+    toVTK2D1->SetInput(rescaler1->GetOutput());
+    toVTK2D1->Update();
+
+    double* dBounds = toVTK2D1->GetOutput()->GetBounds();
+
+    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();
+}
+
+vtkImageData* itkImageProcessor::GetProjection2VTK()
+{
+    // Rescale the intensity of the projection images to 0-255 for output.
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
+    rescaler2->SetOutputMinimum(0);
+    rescaler2->SetOutputMaximum(255);
+    rescaler2->SetInput(flipFilter2->GetOutput());
+
+    rescaler2->Update();
+    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<<" -------- 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;
+
+
+    toVTK2D2->SetInput(rescaler2->GetOutput());
+    toVTK2D2->Update();
+
+    double* dBounds = toVTK2D2->GetOutput()->GetBounds();
+
+    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
+            toVTK2D2->GetOutput();
+}
+
+void itkImageProcessor::WriteProjectionImages()
+{
+
+    // Rescale the intensity of the projection images to 0-255 for output.
+    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+    rescaler1->SetOutputMinimum(0);
+    rescaler1->SetOutputMaximum(255);
+    rescaler1->SetInput(flipFilter1->GetOutput());
+
+    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
+    rescaler2->SetOutputMinimum(0);
+    rescaler2->SetOutputMaximum(255);
+    rescaler2->SetInput(flipFilter2->GetOutput());
+
+
+    using WriterType = itk::ImageFileWriter<OutputImageType>;
+    WriterType::Pointer writer1 = WriterType::New();
+    WriterType::Pointer writer2 = WriterType::New();
+
+    writer1->SetFileName("/mnt/gdrive/Scratch/gfattori/Projection1.tif");
+    writer1->SetInput(rescaler1->GetOutput());
+
+    try
+    {
+        std::cout << "Writing image 1 " << std::endl;
+        writer1->Update();
+    }
+    catch (itk::ExceptionObject & err)
+    {
+        std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+        std::cerr << err << std::endl;
+    }
+
+    writer2->SetFileName("/mnt/gdrive/Scratch/gfattori/Projection2.tif");
+    writer2->SetInput(rescaler2->GetOutput());
+
+    try
+    {
+        std::cout << "Writing image 2" << std::endl;
+        writer2->Update();
+    }
+    catch (itk::ExceptionObject & err)
+    {
+        std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+        std::cerr << err << std::endl;
+    }
+
+}
+
+
+//int itkImageProcessor::InitializeRegistration(){
+
+
+//    //    image2D1Rescaler->Update();
+//    //    image2D2Rescaler->Update();
+
+//    using FlipAxesArrayType = FlipFilterType::FlipAxesArrayType;
+
+//    FlipAxesArrayType flipArray;
+//    flipArray[0] = true;
+//    flipArray[1] = true;
+//    flipArray[2] = false;
+//    flipFilterLocalizer1->SetFlipAxes(flipArray);
+//    flipFilterLocalizer1->Update();
+//    flipArray[0] = false;
+//    flipArray[1] = true;
+//    flipArray[2] = false;
+//    flipFilterLocalizer2->SetFlipAxes(flipArray);
+//    flipFilterLocalizer2->Update();
+
+//    caster3D->Update();
+//    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
+
+//    RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
+//    rescaler1->SetOutputMinimum(0);
+//    rescaler1->SetOutputMaximum(255);
+//    rescaler1->SetInput(flipFilterLocalizer1->GetOutput());
+
+//    rescaler1->Update();
+
+//    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
+
+//    RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
+//    rescaler2->SetOutputMinimum(0);
+//    rescaler2->SetOutputMaximum(255);
+//    rescaler2->SetInput(flipFilterLocalizer2->GetOutput());
+
+//    rescaler2->Update();
+
+
+
+//    //        metric->SetFixedImage1(rescaler1->GetOutput());
+//    //        metric->SetFixedImage2(rescaler2->GetOutput());
+//    //        metric->SetMovingImage(caster3D->GetOutput());
+//    //        metric->SetInterpolator1(interpolator1);
+//    //        metric->SetInterpolator2(interpolator2);
+
+
+//    //        TransformType::OutputVectorType translation;
+
+//    //        translation[0] = 0.;
+//    //        translation[1] = 0.;
+//    //        translation[2] = 0.;
+
+//    //        transform->SetTranslation(translation);
+
+//    //        // constant for converting degrees to radians
+//    //        const double dtr = (atan(1.0) * 4.0) / 180.0;
+//    //        transform->SetRotation(dtr * 0., dtr * 0., dtr * 0.);
+
+//    //        metric->SetTransform(transform);
+
+//    //        std::cout<<" Curr VALUE: " <<
+//    //        metric->GetValue(transform->GetParameters())
+//    //                <<std::endl;
+//    //        return 0;
+
+
+//    std::cout<<
+//                "2d1 "<<
+//                rescaler1->GetOutput()->GetDirection()
+//             <<std::endl;
+//    std::cout<<
+//                "2d2 "<<
+//                rescaler2->GetOutput()->GetDirection()
+//             <<std::endl;
+//    std::cout<<
+//                "3d "<<
+//                caster3D->GetOutput()->GetDirection()
+//             <<std::endl;
+
+
+//    registration->SetFixedImage1(rescaler1->GetOutput());
+//    registration->SetFixedImage2(rescaler2->GetOutput());
+//    registration->SetMovingImage(caster3D->GetOutput());
+
+//    // Initialise the transform
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+//    // Set the order of the computation. Default ZXY
+//    transform->SetComputeZYX(true);
+
+
+//    // The transform is initialised with the translation [tx,ty,tz] and
+//    // rotation [rx,ry,rz] specified on the command line
+
+//    TransformType::OutputVectorType translation;
+
+//    translation[0] = 0.;
+//    translation[1] = 0.;
+//    translation[2] = 0.;
+
+//    transform->SetTranslation(translation);
+
+//    // constant for converting degrees to radians
+//    const double dtr = (atan(1.0) * 4.0) / 180.0;
+//    transform->SetRotation(dtr * 0., dtr * 0., dtr * 0.);
+
+//    // The centre of rotation is set by default to the centre of the 3D
+//    // volume but can be offset from this position using a command
+//    // line specified translation [cx,cy,cz]
+
+//    ImageType3D::PointType       origin3D = image3DIn->GetOrigin();
+//    const itk::Vector<double, 3> resolution3D = image3DIn->GetSpacing();
+
+//    using ImageRegionType3D = ImageType3D::RegionType;
+//    using SizeType3D = ImageRegionType3D::SizeType;
+
+//    ImageRegionType3D region3D = caster3D->GetOutput()->GetBufferedRegion();
+//    SizeType3D        size3D = region3D.GetSize();
+
+//    TransformType::InputPointType isocenter;
+//    if (customized_iso)
+//    {
+//        //      // Isocenter location given by the user.
+//        //      isocenter[0] = origin3D[0] + resolution3D[0] * cx;
+//        //      isocenter[1] = origin3D[1] + resolution3D[1] * cy;
+//        //      isocenter[2] = origin3D[2] + resolution3D[2] * cz;
+//    }
+//    else
+//    {
+//        // Set the center of the image as the isocenter.
+//        isocenter[0] = origin3D[0] + resolution3D[0] * static_cast<double>(size3D[0]) / 2.0;
+//        isocenter[1] = origin3D[1] + resolution3D[1] * static_cast<double>(size3D[1]) / 2.0;
+//        isocenter[2] = origin3D[2] + resolution3D[2] * static_cast<double>(size3D[2]) / 2.0;
+//    }
+
+//    transform->SetCenter(isocenter);
+
+
+//    if (true)
+//    {
+//        std::cout << "3D image size: " << size3D[0] << ", " << size3D[1] << ", " << size3D[2] << std::endl
+//                  << "   resolution: " << resolution3D[0] << ", " << resolution3D[1] << ", " << resolution3D[2] << std::endl
+//                  << "Transform: " << transform << std::endl;
+//    }
+
+
+//    // Set the origin of the 2D image
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+//    // For correct (perspective) projection of the 3D volume, the 2D
+//    // image needs to be placed at a certain distance (the source-to-
+//    // isocenter distance {scd} ) from the focal point, and the normal
+//    // from the imaging plane to the focal point needs to be specified.
+//    //
+//    // By default, the imaging plane normal is set by default to the
+//    // center of the 2D image but may be modified from this using the
+//    // command line parameters [image1centerX, image1centerY,
+//    // image2centerX, image2centerY].
+
+//    double origin2D1[Dimension];
+//    double origin2D2[Dimension];
+
+//    // Note: Two 2D images may have different image sizes and pixel dimensions, although
+//    // scd are the same.
+
+//    const itk::Vector<double, 3> resolution2D1 = rescaler1->GetOutput()->GetSpacing();
+//    const itk::Vector<double, 3> resolution2D2 = rescaler2->GetOutput()->GetSpacing();
+
+//    using ImageRegionType2D = InternalImageType::RegionType;
+//    using SizeType2D = ImageRegionType2D::SizeType;
+
+//    ImageRegionType2D region2D1 = rescaler1->GetOutput()->GetBufferedRegion();
+//    ImageRegionType2D region2D2 = rescaler2->GetOutput()->GetBufferedRegion();
+//    SizeType2D        size2D1 = region2D1.GetSize();
+//    SizeType2D        size2D2 = region2D2.GetSize();
+
+//    double image1centerX = 0.0;
+//    double image1centerY = 0.0;
+//    double image2centerX = 0.0;
+//    double image2centerY = 0.0;
+
+//    double image1resX = 1.0;
+//    double image1resY = 1.0;
+//    double image2resX = 1.0;
+//    double image2resY = 1.0;
+
+//    if (!customized_2DCX)
+//    { // Central axis positions are not given by the user. Use the image centers
+//        // as the central axis position.
+//        image1centerX = ((double)size2D1[0] - 1.) / 2.;
+//        image1centerY = ((double)size2D1[1] - 1.) / 2.;
+//        image2centerX = ((double)size2D2[0] - 1.) / 2.;
+//        image2centerY = ((double)size2D2[1] - 1.) / 2.;
+//    }
+
+//    // 2D Image 1
+//    origin2D1[0] = -resolution2D1[0] * image1centerX;
+//    origin2D1[1] = -resolution2D1[1] * image1centerY;
+//    origin2D1[2] = -scd;
+
+//    flipFilterLocalizer1->GetOutput()->SetOrigin(origin2D1);
+//    rescaler1->GetOutput()->SetOrigin(origin2D1);
+
+//    // 2D Image 2
+//    origin2D2[0] = -resolution2D2[0] * image2centerX;
+//    origin2D2[1] = -resolution2D2[1] * image2centerY;
+//    origin2D2[2] = -scd;
+
+//    flipFilterLocalizer2->GetOutput()->SetOrigin(origin2D2);
+//    rescaler2->GetOutput()->SetOrigin(origin2D2);
+
+//    registration->SetFixedImageRegion1(rescaler1->GetOutput()->GetBufferedRegion());
+//    registration->SetFixedImageRegion2(rescaler2->GetOutput()->GetBufferedRegion());
+
+//    std::cout<<"BuffRegion ROI 1: "<<rescaler1->GetOutput()->GetBufferedRegion()<<std::endl;
+//    std::cout<<"BuffRegion ROI 2: "<<rescaler2->GetOutput()->GetBufferedRegion()<<std::endl;
+
+//    if (true)
+//    {
+//        std::cout << "2D image 1 size: " << size2D1[0] << ", " << size2D1[1] << ", " << size2D1[2] << std::endl
+//                  << "   resolution: " << resolution2D1[0] << ", " << resolution2D1[1] << ", " << resolution2D1[2]
+//                  << std::endl
+//                  << "   and position: " << origin2D1[0] << ", " << origin2D1[1] << ", " << origin2D1[2] << std::endl
+//                  << "2D image 2 size: " << size2D2[0] << ", " << size2D2[1] << ", " << size2D2[2] << std::endl
+//                  << "   resolution: " << resolution2D2[0] << ", " << resolution2D2[1] << ", " << resolution2D2[2]
+//                  << std::endl
+//                  << "   and position: " << origin2D2[0] << ", " << origin2D2[1] << ", " << origin2D2[2] << std::endl;
+//    }
+
+//    // Initialize the ray cast interpolator
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+//    // The ray cast interpolator is used to project the 3D volume. It
+//    // does this by casting rays from the (transformed) focal point to
+//    // each (transformed) pixel coordinate in the 2D image.
+//    //
+//    // In addition a threshold may be specified to ensure that only
+//    // intensities greater than a given value contribute to the
+//    // projected volume. This can be used, for instance, to remove soft
+//    // tissue from projections of CT data and force the registration
+//    // to find a match which aligns bony structures in the images.
+
+//    // 2D Image 1
+//    interpolator1->SetProjectionAngle(dtr * projAngle1);
+//    interpolator1->SetFocalPointToIsocenterDistance(scd);
+//    interpolator1->SetThreshold(threshold);
+//    interpolator1->SetTransform(transform);
+
+//    interpolator1->Initialize();
+
+//    // 2D Image 2
+//    interpolator2->SetProjectionAngle(dtr * projAngle2);
+//    interpolator2->SetFocalPointToIsocenterDistance(scd);
+//    interpolator2->SetThreshold(threshold);
+//    interpolator2->SetTransform(transform);
+
+//    interpolator2->Initialize();
+
+
+//    // Set up the transform and start position
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+//    // The registration start position is intialised using the
+//    // transformation parameters.
+
+//    registration->SetInitialTransformParameters(transform->GetParameters());
+
+//    // We wish to minimize the negative normalized correlation similarity measure.
+
+//    // optimizer->SetMaximize( true );  // for GradientDifferenceTwoImageToOneImageMetric
+//    optimizer->SetMaximize(false); // for NCC
+
+//    optimizer->SetMaximumIteration(10);
+//    optimizer->SetMaximumLineIteration(4); // for Powell's method
+//    optimizer->SetStepLength(4.0);
+//    optimizer->SetStepTolerance(0.02);
+//    optimizer->SetValueTolerance(0.001);
+
+//    // The optimizer weightings are set such that one degree equates to
+//    // one millimeter.
+
+//    itk::Optimizer::ScalesType weightings(transform->GetNumberOfParameters());
+
+//    weightings[0] = 1. / dtr;
+//    weightings[1] = 1. / dtr;
+//    weightings[2] = 1. / dtr;
+//    weightings[3] = 1.;
+//    weightings[4] = 1.;
+//    weightings[5] = 1.;
+
+//    optimizer->SetScales(weightings);
+
+//    if (true)
+//    {
+//        optimizer->Print(std::cout);
+//    }
+
+
+//    // Create the observers
+//    // ~~~~~~~~~~~~~~~~~~~~
+
+//    CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
+
+//    optimizer->AddObserver(itk::IterationEvent(), observer);
+
+
+//    // Start the registration
+//    // ~~~~~~~~~~~~~~~~~~~~~~
+
+//    // Create a timer to record calculation time.
+//    itk::TimeProbesCollectorBase timer;
+
+//    if (true)
+//    {
+//        std::cout << "Starting the registration now..." << std::endl;
+//    }
+
+//    try
+//    {
+//        timer.Start("Registration");
+//        // Start the registration.
+//        registration->StartRegistration();
+//        timer.Stop("Registration");
+//    }
+//    catch (itk::ExceptionObject & err)
+//    {
+//        std::cout << "ExceptionObject caught !" << std::endl;
+//        std::cout << err << std::endl;
+//        return -1;
+//    }
+
+//    using ParametersType = RegistrationType::ParametersType;
+//    ParametersType finalParameters = registration->GetLastTransformParameters();
+
+//    const double RotationAlongX = finalParameters[0] / dtr; // Convert radian to degree
+//    const double RotationAlongY = finalParameters[1] / dtr;
+//    const double RotationAlongZ = finalParameters[2] / dtr;
+//    const double TranslationAlongX = finalParameters[3];
+//    const double TranslationAlongY = finalParameters[4];
+//    const double TranslationAlongZ = finalParameters[5];
+
+//    const int numberOfIterations = optimizer->GetCurrentIteration();
+
+//    const double bestValue = optimizer->GetValue();
+
+//    std::cout << "Result = " << std::endl;
+//    std::cout << " Rotation Along X = " << RotationAlongX << " deg" << std::endl;
+//    std::cout << " Rotation Along Y = " << RotationAlongY << " deg" << std::endl;
+//    std::cout << " Rotation Along Z = " << RotationAlongZ << " deg" << std::endl;
+//    std::cout << " Translation X = " << TranslationAlongX << " mm" << std::endl;
+//    std::cout << " Translation Y = " << TranslationAlongY << " mm" << std::endl;
+//    std::cout << " Translation Z = " << TranslationAlongZ << " mm" << std::endl;
+//    std::cout << " Number Of Iterations = " << numberOfIterations << std::endl;
+//    std::cout << " Metric value  = " << bestValue << std::endl;
+
+
+//    // Write out the projection images at the registration position
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+//    TransformType::Pointer finalTransform = TransformType::New();
+//    // The transform is determined by the parameters and the rotation center.
+//    finalTransform->SetParameters(finalParameters);
+//    finalTransform->SetCenter(isocenter);
+
+//    using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
+
+//    // The ResampleImageFilter is the driving force for the projection image generation.
+//    ResampleFilterType::Pointer resampleFilter1 = ResampleFilterType::New();
+
+//    resampleFilter1->SetInput(caster3D->GetOutput()); // Link the 3D volume.
+//    resampleFilter1->SetDefaultPixelValue(0);
+
+//    // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
+//    // have been set before registration. Here we only need to replace the initial
+//    // transform with the final transform.
+//    interpolator1->SetTransform(finalTransform);
+//    interpolator1->Initialize();
+//    resampleFilter1->SetInterpolator(interpolator1);
+
+//    // The output 2D projection image has the same image size, origin, and the pixel spacing as
+//    // those of the input 2D image.
+//    resampleFilter1->SetSize(rescaler1->GetOutput()->GetLargestPossibleRegion().GetSize());
+//    resampleFilter1->SetOutputOrigin(rescaler1->GetOutput()->GetOrigin());
+//    resampleFilter1->SetOutputSpacing(rescaler1->GetOutput()->GetSpacing());
+
+//    // Do the same thing for the output image 2.
+//    ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New();
+//    resampleFilter2->SetInput(caster3D->GetOutput());
+//    resampleFilter2->SetDefaultPixelValue(0);
+
+//    // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
+//    // have been set before registration. Here we only need to replace the initial
+//    // transform with the final transform.
+//    interpolator2->SetTransform(finalTransform);
+//    interpolator2->Initialize();
+//    resampleFilter2->SetInterpolator(interpolator2);
+
+//    resampleFilter2->SetSize(rescaler2->GetOutput()->GetLargestPossibleRegion().GetSize());
+//    resampleFilter2->SetOutputOrigin(rescaler2->GetOutput()->GetOrigin());
+//    resampleFilter2->SetOutputSpacing(rescaler2->GetOutput()->GetSpacing());
+
+//    /////////////////////////////---DEGUG--START----////////////////////////////////////
+//    if (true)
+//    {
+//        InternalImageType::PointType outputorigin2D1 = rescaler1->GetOutput()->GetOrigin();
+//        std::cout << "Output image 1 origin: " << outputorigin2D1[0] << ", " << outputorigin2D1[1] << ", "
+//                << outputorigin2D1[2] << std::endl;
+//        InternalImageType::PointType outputorigin2D2 = rescaler2->GetOutput()->GetOrigin();
+//        std::cout << "Output image 2 origin: " << outputorigin2D2[0] << ", " << outputorigin2D2[1] << ", "
+//                << outputorigin2D2[2] << std::endl;
+//    }
+//    /////////////////////////////---DEGUG--END----//////////////////////////////////////
+
+
+//    // As explained before, the computed projection is upsided-down.
+//    // Here we use a FilpImageFilter to flip the images in y-direction.
+//    flipFilter1->SetInput(resampleFilter1->GetOutput());
+//    flipFilter2->SetInput(resampleFilter2->GetOutput());
+
+//    // Rescale the intensity of the projection images to 0-255 for output.
+//    //using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
+
+//    RescaleFilterType::Pointer rescaler1b = RescaleFilterType::New();
+//    rescaler1b->SetOutputMinimum(0);
+//    rescaler1b->SetOutputMaximum(255);
+//    rescaler1b->SetInput(flipFilter1->GetOutput());
+
+//    RescaleFilterType::Pointer rescaler2b = RescaleFilterType::New();
+//    rescaler2b->SetOutputMinimum(0);
+//    rescaler2b->SetOutputMaximum(255);
+//    rescaler2b->SetInput(flipFilter2->GetOutput());
+
+
+//    using WriterType = itk::ImageFileWriter<OutputImageType>;
+//    WriterType::Pointer writer1 = WriterType::New();
+//    WriterType::Pointer writer2 = WriterType::New();
+
+//    writer1->SetFileName("Out1.tif");
+//    writer1->SetInput(0,rescaler1b->GetOutput());
+
+//    try
+//    {
+//        //std::cout << "Writing image: " << fileOutput1 << std::endl;
+//        writer1->Update();
+//    }
+//    catch (itk::ExceptionObject & err)
+//    {
+//        std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+//        std::cerr << err << std::endl;
+//    }
+
+//    writer2->SetFileName("Out2.tif");
+//    writer2->SetInput(0,rescaler2b->GetOutput());
+
+//    try
+//    {
+//        writer2->Update();
+//    }
+//    catch (itk::ExceptionObject & err)
+//    {
+//        std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
+//        std::cerr << err << std::endl;
+//    }
+//    timer.Report();
+
+//    return 0;
+
+
+
+
+
+//}
+
+//void itkImageProcessor::RunRegistration()
+//{
+
+//}
+
+
+
+
+void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ)
+{
+    TZero[0]= dX;
+    TZero[1]= dY;
+    TZero[2]= dZ;
+}
+
+void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ)
+{
+    RZero[0]=dX;
+    RZero[1]=dY;
+    RZero[2]=dZ;
+}
+
+
+} // end namespace itk
+
diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h
new file mode 100644
index 0000000..99c88d8
--- /dev/null
+++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h
@@ -0,0 +1,273 @@
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+#ifndef ITKIMAGEPROCESSOR_H
+#define ITKIMAGEPROCESSOR_H
+
+
+#include "itkCommand.h"
+
+#include "itkPowellOptimizer.h"
+#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
+#include "itkEuler3DTransform.h"
+#include "itkTwoProjectionImageRegistrationMethod.h"
+#include "itkResampleImageFilter.h"
+#include "itkCastImageFilter.h"
+#include "itkRescaleIntensityImageFilter.h"
+#include "itkFlipImageFilter.h"
+#include "itkChangeInformationImageFilter.h"
+#include "itkTimeProbesCollectorBase.h"
+#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
+#include "itkGDCMImageIO.h"
+#include "itkMetaDataObject.h"
+
+#include "gdcmGlobal.h"
+
+#include "itkImageToVTKImageFilter.h"
+
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkImageSeriesReader.h"
+
+
+namespace itk
+{
+
+
+class CommandIterationUpdate : public itk::Command
+{
+public:
+  using Self = CommandIterationUpdate;
+  using Superclass = itk::Command;
+  using Pointer = itk::SmartPointer<Self>;
+  itkNewMacro(Self);
+
+protected:
+  CommandIterationUpdate() = default;
+
+public:
+  using OptimizerType = itk::PowellOptimizer;
+  using OptimizerPointer = const OptimizerType *;
+
+  void
+  Execute(itk::Object * caller, const itk::EventObject & event) override
+  {
+    Execute((const itk::Object *)caller, event);
+  }
+
+  void
+  Execute(const itk::Object * object, const itk::EventObject & event) override
+  {
+    auto optimizer = dynamic_cast<OptimizerPointer>(object);
+    if (typeid(event) != typeid(itk::IterationEvent))
+    {
+      return;
+    }
+    std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
+    std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
+    std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
+  }
+};
+
+
+/* reference string required for comparison with tag values */
+static const char *ImageOrientationStrings[] = {
+    "HFS",
+    "FFS"
+};
+
+class ITK_EXPORT itkImageProcessor : public Object
+{
+
+    constexpr static unsigned int Dimension = 3;
+
+public:
+  /** Standard "Self" typedef. */
+  typedef itkImageProcessor  Self;
+  /** Standard "Superclass" typedef. */
+  typedef Object Superclass;
+  /** Smart pointer typedef support */
+  typedef SmartPointer<Self>  Pointer;
+  typedef SmartPointer<const Self>  ConstPointer;
+  /** Method of creation through the object factory. */
+  itkNewMacro(Self);
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(itkImageProcessor, Object);
+
+  /** Input data load methods*/
+  void load3DVolume(const char *);
+  int load3DSerie(const char* );
+  int load2D(const char *);
+
+  /** Projection angles - Gantry angle */
+  void SetProjectionAngle1(double);
+  void SetProjectionAngle2(double);
+
+  /** Source to axis distance - SAD */
+  void SetSCD(double);
+
+  /** Intensity threshold used for image projection,
+   *  any voxel below threshold is ignored */
+  void SetIntensityThreshold(double);
+
+  /** Custom settings of the projection images */
+  void SetCustom_Isocenter(double,double,double);
+  void SetCustom_2Dres(double,double,double,double);
+  void SetCustom_2Dsize(int,int,int,int);
+  void SetCustom_2Dcenter(double dX1, double dY1, double dX2, double dY2);
+
+  /** Set transform parameters for 3D Volume */
+  void SetInitialTranslations(double, double, double);
+  void SetInitialRotations(double, double, double);
+
+  /** Get the Localizer intensity window and center for rendering */
+  double* GetImageIntensityWindow(int );
+
+  /** 2D-3D registration routines */
+//  int InitializeRegistration();
+  void InitializeProjector();
+//  void RunRegistration();
+
+  /** Compute Digital Localizers */
+  void GetProjectionImages();
+
+  /** Conveniency methods to get VTK images for rendering */
+  vtkImageData* GetProjection1VTK();
+  vtkImageData* GetProjection2VTK();
+
+  vtkImageData* GetLocalizer1VTK();
+  vtkImageData* GetLocalizer2VTK();
+
+  /** Debug writers */
+  void WriteProjectionImages();
+  void Write2DImages();
+
+protected:
+  itkImageProcessor();
+  virtual ~itkImageProcessor();
+  void PrintSelf(std::ostream& os, Indent indent) const;
+
+
+private:
+  itkImageProcessor(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+
+  void Finalize();
+
+/** Various pixel types */
+  using InternalPixelType = float;
+  using PixelType2D = short;
+  using PixelType3D = short;
+  using OutputPixelType = unsigned char;
+
+  /** Image types */
+  using ImageType2D = itk::Image<PixelType3D, Dimension>;
+  using ImageType3D = itk::Image<PixelType3D, Dimension>;
+  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
+
+  /** The following lines define each of the components used in the
+      registration: The transform, optimizer, metric, interpolator and
+      the registration method itself. */
+  using InternalImageType = itk::Image<InternalPixelType, Dimension>;
+  using TransformType = itk::Euler3DTransform<double>;
+  using OptimizerType = itk::PowellOptimizer;
+  using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
+  using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
+  using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
+  using ParametersType = RegistrationType::ParametersType;
+
+  /** Image reader types */
+  using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
+  using ImageReaderType3D = itk::ImageFileReader<ImageType3D>;
+  using ImageSeriesReaderType3D = itk::ImageSeriesReader<ImageType3D>;
+  using ImageIOType = itk::GDCMImageIO;
+
+  /** widely used filters: flip, cast, rescale, ... */
+  using FlipFilterType = itk::FlipImageFilter<InternalImageType>;
+  using Input2DRescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
+  using CastFilterType3D = itk::CastImageFilter<ImageType3D, InternalImageType>;
+  using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
+  using ChangeInformationFilterInputType = itk::ChangeInformationImageFilter<ImageType3D>;
+
+  /** ITK to VTK filter */
+  using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
+
+  /** Enum type for Orientation */
+  enum eImageOrientationType { HFS = 0, FFS = 1};
+
+
+  MetricType::Pointer       metric;
+  TransformType::Pointer    transform;
+  OptimizerType::Pointer    optimizer;
+  InterpolatorType::Pointer interpolator1;
+  InterpolatorType::Pointer interpolator2;
+  RegistrationType::Pointer registration;
+
+  ImageReaderType2D::Pointer imageReader2D;
+  ImageReaderType3D::Pointer imageReader3D;
+  ImageSeriesReaderType3D::Pointer imageSeriesReader3D;
+
+  CommandIterationUpdate::Pointer observer;
+
+  FlipFilterType::Pointer flipFilter1;
+  FlipFilterType::Pointer flipFilter2;
+  FlipFilterType::Pointer flipFilterLocalizer1;
+  FlipFilterType::Pointer flipFilterLocalizer2;
+
+  CastFilterType3D::Pointer caster3D;
+  CastFilterType3D::Pointer caster2D1;
+  CastFilterType3D::Pointer caster2D2;
+  CastFilterType3D::Pointer caster2DInput;
+
+  TransformType::InputPointType isocenter3D;
+
+  ITKtoVTKFilterType::Pointer toVTK2D1;
+  ITKtoVTKFilterType::Pointer toVTK2D2;
+  ITKtoVTKFilterType::Pointer toVTKLocalizer1;
+  ITKtoVTKFilterType::Pointer toVTKLocalizer2;
+
+  ImageType3D::Pointer image3DIn;
+  InternalImageType::Pointer image2D1In;
+  InternalImageType::Pointer image2D2In;
+
+  Input2DRescaleFilterType::Pointer image2D1Rescaler;
+  Input2DRescaleFilterType::Pointer image2D2Rescaler;
+
+  InternalImageType::Pointer TMPImg1;
+  InternalImageType::Pointer TMPImg2;
+
+
+  ChangeInformationFilterInputType::Pointer image3DIECFilt;
+
+  double image1res[2], image2res[2];
+  double image1center[2], image2center[2];
+  int image1Size[2], image2Size[2];
+  double isocenter[3];
+  double scd;
+  double threshold;
+  double projAngle1, projAngle2;
+
+  double TZero[3], RZero[3];
+  double image3DCOV[3];
+
+  bool customized_iso,customized_2DCX,customized_2DRES, customized_2DSIZE;
+  bool image2D1Loaded, image2D2Loaded, image3DLoaded;
+
+  double image1IntensityWindow[2], image2IntensityWindow[2];
+
+};
+
+
+}
+
+
+#endif
diff --git a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
new file mode 100644
index 0000000..666c0ea
--- /dev/null
+++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
@@ -0,0 +1,248 @@
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+/*=========================================================================
+Calculate DRR from a CT dataset using incremental ray-tracing algorithm
+The algorithm was initially proposed by Robert Siddon and improved by
+Filip Jacobs etc.
+
+-------------------------------------------------------------------------
+References:
+
+R. L. Siddon, "Fast calculation of the exact radiological path for a
+threedimensional CT array," Medical Physics 12, 252-55 (1985).
+
+F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
+"A fast algorithm to calculate the exact radiological path through a pixel
+or voxel space," Journal of Computing and Information Technology ?
+CIT 6, 89-94 (1998).
+
+=========================================================================*/
+#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_h
+#define itkgSiddonJacobsRayCastInterpolateImageFunction_h
+
+#include "itkInterpolateImageFunction.h"
+#include "itkTransform.h"
+#include "itkVector.h"
+#include "itkEuler3DTransform.h"
+
+namespace itk
+{
+
+/** \class gSiddonJacobsRayCastInterpolateImageFunction
+ * \brief Projective interpolation of an image at specified positions.
+ *
+ * SiddonJacobsRayCastInterpolateImageFunction casts rays through a 3-dimensional
+ * image
+ * \warning This interpolator works for 3-dimensional images only.
+ *
+ * \ingroup ImageFunctions
+ * \ingroup TwoProjectionRegistration
+ */
+template <typename TInputImage, typename TCoordRep = float>
+class gSiddonJacobsRayCastInterpolateImageFunction : public InterpolateImageFunction<TInputImage, TCoordRep>
+{
+public:
+  ITK_DISALLOW_COPY_AND_ASSIGN(gSiddonJacobsRayCastInterpolateImageFunction);
+
+  /** Standard class type alias. */
+  using Self = gSiddonJacobsRayCastInterpolateImageFunction;
+  using Superclass = InterpolateImageFunction<TInputImage, TCoordRep>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** Constants for the image dimensions */
+  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
+
+
+  using TransformType = Euler3DTransform<TCoordRep>;
+
+  using TransformPointer = typename TransformType::Pointer;
+  using InputPointType = typename TransformType::InputPointType;
+  using OutputPointType = typename TransformType::OutputPointType;
+  using TransformParametersType = typename TransformType::ParametersType;
+  using TransformJacobianType = typename TransformType::JacobianType;
+
+  using PixelType = typename Superclass::InputPixelType;
+
+  using SizeType = typename TInputImage::SizeType;
+
+  using DirectionType = Vector<TCoordRep, 3>;
+
+  /**  Type of the Interpolator Base class */
+  using InterpolatorType = InterpolateImageFunction<TInputImage, TCoordRep>;
+
+  using InterpolatorPointer = typename InterpolatorType::Pointer;
+
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(gSiddonJacobsRayCastInterpolateImageFunction, InterpolateImageFunction);
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** OutputType type alias support. */
+  using OutputType = typename Superclass::OutputType;
+
+  /** InputImageType type alias support. */
+  using InputImageType = typename Superclass::InputImageType;
+
+  /** InputImageConstPointer type alias support. */
+  using InputImageConstPointer = typename Superclass::InputImageConstPointer;
+
+  /** RealType type alias support. */
+  using RealType = typename Superclass::RealType;
+
+  /** Dimension underlying input image. */
+  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
+
+  /** Point type alias support. */
+  using PointType = typename Superclass::PointType;
+
+  /** Index type alias support. */
+  using IndexType = typename Superclass::IndexType;
+
+  /** ContinuousIndex type alias support. */
+  using ContinuousIndexType = typename Superclass::ContinuousIndexType;
+
+
+  /** \brief
+   * Interpolate the image at a point position.
+   *
+   * Returns the interpolated image intensity at a
+   * specified point position. No bounds checking is done.
+   * The point is assume to lie within the image buffer.
+   *
+   * ImageFunction::IsInsideBuffer() can be used to check bounds before
+   * calling the method.
+   */
+  OutputType
+  Evaluate(const PointType & point) const override;
+
+  /** Interpolate the image at a continuous index position
+   *
+   * Returns the interpolated image intensity at a
+   * specified index position. No bounds checking is done.
+   * The point is assume to lie within the image buffer.
+   *
+   * Subclasses must override this method.
+   *
+   * ImageFunction::IsInsideBuffer() can be used to check bounds before
+   * calling the method.
+   */
+  OutputType
+  EvaluateAtContinuousIndex(const ContinuousIndexType & index) const override;
+
+  virtual void
+  Initialize();
+
+  /** Connect the Transform. */
+  itkSetObjectMacro(Transform, TransformType);
+  /** Get a pointer to the Transform.  */
+  itkGetConstObjectMacro(Transform, TransformType);
+
+  /** Set and get the focal point to isocenter distance in mm */
+  itkSetMacro(FocalPointToIsocenterDistance, double);
+  itkGetMacro(FocalPointToIsocenterDistance, double);
+
+  /** Set and get the Lianc grantry rotation angle in radians */
+  itkSetMacro(ProjectionAngle, double);
+  itkGetMacro(ProjectionAngle, double);
+
+  /** Set and get the Threshold */
+  itkSetMacro(Threshold, double);
+  itkGetMacro(Threshold, double);
+
+  /** Check if a point is inside the image buffer.
+   * \warning For efficiency, no validity checking of
+   * the input image pointer is done. */
+  inline bool
+  IsInsideBuffer(const PointType &) const override
+  {
+    return true;
+  }
+  bool
+  IsInsideBuffer(const ContinuousIndexType &) const override
+  {
+    return true;
+  }
+  bool
+  IsInsideBuffer(const IndexType &) const override
+  {
+    return true;
+  }
+
+#if !defined(ITKV4_COMPATIBILITY)
+  SizeType
+  GetRadius() const override
+  {
+    const InputImageType * input = this->GetInputImage();
+    if (!input)
+    {
+      itkExceptionMacro("Input image required!");
+    }
+    return input->GetLargestPossibleRegion().GetSize();
+  }
+#endif
+
+protected:
+  gSiddonJacobsRayCastInterpolateImageFunction();
+
+  ~gSiddonJacobsRayCastInterpolateImageFunction() override = default;
+
+  void
+  PrintSelf(std::ostream & os, Indent indent) const override;
+
+  /// Transformation used to calculate the new focal point position
+  TransformPointer m_Transform; // Displacement of the volume
+  // Overall inverse transform used to calculate the ray position in the input space
+  TransformPointer m_InverseTransform;
+
+  // The threshold above which voxels along the ray path are integrated
+  double m_Threshold;
+  double m_FocalPointToIsocenterDistance; // Focal point to isocenter distance
+  double m_ProjectionAngle;               // Linac gantry rotation angle in radians
+
+private:
+  void
+                   ComputeInverseTransform() const;
+  TransformPointer m_GantryRotTransform; // Gantry rotation transform
+  TransformPointer m_CamShiftTransform;  // Camera shift transform camRotTransform
+  TransformPointer m_CamRotTransform;    // Camera rotation transform
+  TransformPointer m_ComposedTransform;  // Composed transform
+  PointType        m_SourcePoint;        // Coordinate of the source in the standard Z projection geometry
+  PointType        m_SourceWorld;        // Coordinate of the source in the world coordinate system
+};
+
+} // namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkgSiddonJacobsRayCastInterpolateImageFunction.hxx"
+#endif
+
+#endif
diff --git a/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
new file mode 100644
index 0000000..3c65087
--- /dev/null
+++ b/reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
@@ -0,0 +1,488 @@
+
+/*
+
+Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
+version of incremental ray-tracing algorithm
+
+gfattori 08.11.2021
+
+*/
+
+
+
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+/*=========================================================================
+The algorithm was initially proposed by Robert Siddon and improved by
+Filip Jacobs etc.
+
+-------------------------------------------------------------------------
+References:
+
+R. L. Siddon, "Fast calculation of the exact radiological path for a
+threedimensional CT array," Medical Physics 12, 252-55 (1985).
+
+F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
+"A fast algorithm to calculate the exact radiological path through a pixel
+or voxel space," Journal of Computing and Information Technology ?
+CIT 6, 89-94 (1998).
+
+=========================================================================*/
+
+#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
+#define itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
+
+#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
+
+#include "itkMath.h"
+#include <cstdlib>
+#include <unistd.h>
+
+namespace itk
+{
+
+template <typename TInputImage, typename TCoordRep>
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::gSiddonJacobsRayCastInterpolateImageFunction()
+{
+  m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm.
+  m_ProjectionAngle = 0.;                  // Angle in radians betweeen projection central axis and reference axis
+  m_Threshold = 0.;                        // Intensity threshold, below which is ignored.
+
+  m_SourcePoint[0] = 0.;
+  m_SourcePoint[1] = 0.;
+  m_SourcePoint[2] = 0.;
+
+  m_InverseTransform = TransformType::New();
+  m_InverseTransform->SetComputeZYX(true);
+
+  m_ComposedTransform = TransformType::New();
+  m_ComposedTransform->SetComputeZYX(true);
+
+  m_GantryRotTransform = TransformType::New();
+  m_GantryRotTransform->SetComputeZYX(true);
+  m_GantryRotTransform->SetIdentity();
+
+  m_CamShiftTransform = TransformType::New();
+  m_CamShiftTransform->SetComputeZYX(true);
+  m_CamShiftTransform->SetIdentity();
+
+  m_CamRotTransform = TransformType::New();
+  m_CamRotTransform->SetComputeZYX(true);
+  m_CamRotTransform->SetIdentity();
+  // constant for converting degrees into radians
+  const float dtr = (atan(1.0) * 4.0) / 180.0;
+  m_CamRotTransform->SetRotation(dtr * (-90.0), 0.0, 0.0);
+
+  m_Threshold = 0;
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+void
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
+{
+  this->Superclass::PrintSelf(os, indent);
+
+  os << indent << "Threshold: " << m_Threshold << std::endl;
+  os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
+{
+  float     rayVector[3];
+  IndexType cIndex;
+
+  PointType  drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
+  OutputType pixval;
+
+
+  float firstIntersection[3];
+  float alphaX1, alphaXN, alphaXmin, alphaXmax;
+  float alphaY1, alphaYN, alphaYmin, alphaYmax;
+  float alphaZ1, alphaZN, alphaZmin, alphaZmax;
+  float alphaMin, alphaMax;
+  float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
+  float alphaUx, alphaUy, alphaUz;
+  float alphaIntersectionUp[3], alphaIntersectionDown[3];
+  float d12, value;
+  float firstIntersectionIndex[3];
+  int   firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
+  int   iU, jU, kU;
+
+
+  // Min/max values of the output pixel type AND these values
+  // represented as the output type of the interpolator
+  const OutputType minOutputValue = itk::NumericTraits<OutputType>::NonpositiveMin();
+  const OutputType maxOutputValue = itk::NumericTraits<OutputType>::max();
+
+  // If the volume was shifted, recalculate the overall inverse transform
+  unsigned long int interpMTime = this->GetMTime();
+  unsigned long int vTransformMTime = m_Transform->GetMTime();
+
+  if (interpMTime < vTransformMTime)
+  {
+    this->ComputeInverseTransform();
+    // The m_SourceWorld should be computed here to avoid the repeatedly calculation
+    // for each projection ray. However, we are in a const function, which prohibits
+    // the modification of class member variables. So the world coordiate of the source
+    // point is calculated for each ray as below. Performance improvement may be made
+    // by using a static variable?
+    // m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
+  }
+
+  //PointType        m_SourcePointSliding = m_SourcePoint;
+  //m_SourcePointSliding[2] = point[2];
+  PointType SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
+
+  // Get ths input pointers
+  InputImageConstPointer inputPtr = this->GetInputImage();
+
+  typename InputImageType::SizeType    sizeCT;
+  typename InputImageType::RegionType  regionCT;
+  typename InputImageType::SpacingType ctPixelSpacing;
+  typename InputImageType::PointType   ctOrigin;
+
+  ctPixelSpacing = inputPtr->GetSpacing();
+  ctOrigin = inputPtr->GetOrigin();
+  regionCT = inputPtr->GetLargestPossibleRegion();
+  sizeCT = regionCT.GetSize();
+
+  drrPixelWorld = m_InverseTransform->TransformPoint(point);
+
+  //slide source along z axis...
+  SourceWorld[2] = drrPixelWorld[2];
+
+  // calculate the detector position for this pixel center by moving 
+  // 2*m_FocalPointToIsocenterDistance from the source in the pixel
+  // directions
+  double p1[2], p2[2];
+
+  p1[0] = SourceWorld[0];
+  p1[1] = SourceWorld[1];
+
+  p2[0] = drrPixelWorld[0];
+  p2[1] = drrPixelWorld[1];
+
+  double P12D = sqrt( (p2[0]-p1[0]) * (p2[0]-p1[0]) + (p2[1]-p1[1]) * (p2[1]-p1[1]) );
+  double p12V[2];
+  p12V[0] = p2[0] - p1[0];
+  p12V[1] = p2[1] - p1[1];
+
+  double p12v [2];
+  p12v[0] = p12V[0]/P12D;
+  p12v[1] = p12V[1]/P12D;
+
+  double pDest [2];
+  pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
+  pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
+
+
+  // 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
+  rayVector[0] = pDest[0] - SourceWorld[0];
+  rayVector[1] = pDest[1] - SourceWorld[1];
+  rayVector[2] = drrPixelWorld[2] - SourceWorld[2];
+
+
+  /* Calculate the parametric  values of the first  and  the  last
+  intersection points of  the  ray  with the X,  Y, and Z-planes  that
+  define  the  CT volume. */
+  if (rayVector[0] != 0)
+  {
+    alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
+    alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+    alphaXmin = std::min(alphaX1, alphaXN);
+    alphaXmax = std::max(alphaX1, alphaXN);
+  }
+  else
+  {
+    alphaXmin = -2;
+    alphaXmax = 2;
+  }
+
+  if (rayVector[1] != 0)
+  {
+    alphaY1 = (0.0 - SourceWorld[1]) / rayVector[1];
+    alphaYN = (sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+    alphaYmin = std::min(alphaY1, alphaYN);
+    alphaYmax = std::max(alphaY1, alphaYN);
+  }
+  else
+  {
+    alphaYmin = -2;
+    alphaYmax = 2;
+  }
+
+  if (rayVector[2] != 0)
+  {
+    alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
+    alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+    alphaZmin = std::min(alphaZ1, alphaZN);
+    alphaZmax = std::max(alphaZ1, alphaZN);
+  }
+  else
+  {
+    alphaZmin = -2;
+    alphaZmax = 2;
+  }
+
+  /* Get the very first and the last alpha values when the ray
+  intersects with the CT volume. */
+  alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin);
+  alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax);
+
+  /* Calculate the parametric values of the first intersection point
+  of the ray with the X, Y, and Z-planes after the ray entered the
+  CT volume. */
+
+  firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0];
+  firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1];
+  firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[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];
+
+  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 (rayVector[0] == 0)
+  {
+    alphaX = 2;
+  }
+  else
+  {
+    alphaIntersectionUp[0] = (firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+    alphaIntersectionDown[0] = (firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
+    alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
+  }
+
+  if (rayVector[1] == 0)
+  {
+    alphaY = 2;
+  }
+  else
+  {
+    alphaIntersectionUp[1] = (firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+    alphaIntersectionDown[1] = (firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
+    alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
+  }
+
+  if (rayVector[2] == 0)
+  {
+    alphaZ = 2;
+  }
+  else
+  {
+    alphaIntersectionUp[2] = (firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+    alphaIntersectionDown[2] = (firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
+    alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
+  }
+
+  /* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */
+  if (rayVector[0] != 0)
+  {
+    alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
+  }
+  else
+  {
+    alphaUx = 999;
+  }
+  if (rayVector[1] != 0)
+  {
+    alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]);
+  }
+  else
+  {
+    alphaUy = 999;
+  }
+  if (rayVector[2] != 0)
+  {
+    alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]);
+  }
+  else
+  {
+    alphaUz = 999;
+  }
+  /* Calculate voxel index incremental values along the ray path. */
+  if (SourceWorld[0] < drrPixelWorld[0])
+  {
+    iU = 1;
+  }
+  else
+  {
+    iU = -1;
+  }
+  if (SourceWorld[1] < drrPixelWorld[1])
+  {
+    jU = 1;
+  }
+  else
+  {
+    jU = -1;
+  }
+
+  if (SourceWorld[2] < drrPixelWorld[2])
+  {
+    kU = 1;
+  }
+  else
+  {
+    kU = -1;
+  }
+
+  d12 = 0.0; /* Initialize the sum of the voxel intensities along the ray path to zero. */
+
+
+  /* Initialize the current ray position. */
+  alphaCmin = std::min(std::min(alphaX, alphaY), alphaZ);
+
+  /* Initialize the current voxel index. */
+  cIndex[0] = firstIntersectionIndexDown[0];
+  cIndex[1] = firstIntersectionIndexDown[1];
+  cIndex[2] = firstIntersectionIndexDown[2];
+
+  while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
+  {
+    /* Store the current ray position */
+    alphaCminPrev = alphaCmin;
+
+    if ((alphaX <= alphaY) && (alphaX <= alphaZ))
+    {
+      /* Current ray front intercepts with x-plane. Update alphaX. */
+      alphaCmin = alphaX;
+      cIndex[0] = cIndex[0] + iU;
+      alphaX = alphaX + alphaUx;
+    }
+    else if ((alphaY <= alphaX) && (alphaY <= alphaZ))
+    {
+      /* Current ray front intercepts with y-plane. Update alphaY. */
+      alphaCmin = alphaY;
+      cIndex[1] = cIndex[1] + jU;
+      alphaY = alphaY + alphaUy;
+    }
+    else
+    {
+      /* Current ray front intercepts with z-plane. Update alphaZ. */
+      alphaCmin = alphaZ;
+      cIndex[2] = cIndex[2] + kU;
+      alphaZ = alphaZ + alphaUz;
+    }
+
+    if ((cIndex[0] >= 0) && (cIndex[0] < static_cast<IndexValueType>(sizeCT[0])) && (cIndex[1] >= 0) &&
+        (cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
+        (cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
+    {
+      /* If it is a valid index, get the voxel intensity. */
+      value = static_cast<float>(inputPtr->GetPixel(cIndex));
+      if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
+      {
+        d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
+      }
+    }
+  }
+
+  if (d12 < minOutputValue)
+  {
+    pixval = minOutputValue;
+  }
+  else if (d12 > maxOutputValue)
+  {
+    pixval = maxOutputValue;
+  }
+  else
+  {
+    pixval = static_cast<OutputType>(d12);
+  }
+  return (pixval);
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateAtContinuousIndex(
+  const ContinuousIndexType & index) const
+{
+  OutputPointType point;
+  this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
+
+  return this->Evaluate(point);
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+void
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInverseTransform() const
+{
+  m_ComposedTransform->SetIdentity();
+  m_ComposedTransform->Compose(m_Transform, 0);
+
+  typename TransformType::InputPointType isocenter;
+  isocenter = m_Transform->GetCenter();
+  // An Euler 3D transform is used to rotate the volume to simulate the roation of the linac gantry.
+  // The rotation is about z-axis. After the transform, a AP projection geometry (projecting
+  // towards positive y direction) is established.
+  m_GantryRotTransform->SetRotation(0.0, 0.0, -m_ProjectionAngle);
+  m_GantryRotTransform->SetCenter(isocenter);
+  m_ComposedTransform->Compose(m_GantryRotTransform, 0);
+
+  // An Euler 3D transfrom is used to shift the source to the origin.
+  typename TransformType::OutputVectorType focalpointtranslation;
+  focalpointtranslation[0] = -isocenter[0];
+  focalpointtranslation[1] = m_FocalPointToIsocenterDistance - isocenter[1];
+  focalpointtranslation[2] = -isocenter[2];
+  m_CamShiftTransform->SetTranslation(focalpointtranslation);
+  m_ComposedTransform->Compose(m_CamShiftTransform, 0);
+
+  // A Euler 3D transform is used to establish the standard negative z-axis projection geometry. (By
+  // default, the camera is situated at the origin, points down the negative z-axis, and has an up-
+  // vector of (0, 1, 0).)
+
+  m_ComposedTransform->Compose(m_CamRotTransform, 0);
+
+  // The overall inverse transform is computed. The inverse transform will be used by the interpolation
+  // procedure.
+  m_ComposedTransform->GetInverse(m_InverseTransform);
+  this->Modified();
+}
+
+
+template <typename TInputImage, typename TCoordRep>
+void
+gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Initialize()
+{
+  this->ComputeInverseTransform();
+  m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
+}
+
+} // namespace itk
+
+#endif