/*

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 "itkDRTHelpers.h"
#include "itkVersion.h"
#include <gdcmReader.h>
#include <itkMatrix.h>
#include "itkImageDuplicator.h"
#include "itkPermuteAxesImageFilter.h"
#include "itkMatrix.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkOrientImageFilter.h"
#include <gdcmIPPSorter.h>

#include "itkCommand.h"

#include "itkTimeProbesCollectorBase.h"

#include <string>



namespace itk
{

itkImageProcessor::itkImageProcessor()
{
    iNWUnits = 1;


    transform1 = TransformType::New();
    transform1->SetIdentity();
    transform2 = TransformType::New();
    transform2->SetIdentity();

    m_InternalTransf1 = InternalTransformMetaInformation::New();
    m_InternalTransf2 = InternalTransformMetaInformation::New();


    interpolator1 = InterpolatorType::New();
    interpolator2 = InterpolatorType::New();

    toVTK2D1 = ITKtoVTKFilterType::New();
    toVTK2D2 = ITKtoVTKFilterType::New();
    toVTKLocalizer1 = ITKtoVTKFilterType::New();
    toVTKLocalizer2 = ITKtoVTKFilterType::New();


    m_LATSourceDupli = DuplicatorType::New();
    m_PASourceDupli = DuplicatorType::New();
    m_VolumeSourceDupli = DuplicatorType::New();

    Std_DRT2LPS.SetIdentity();
    for(int iIdx = 0 ; iIdx < 3; iIdx++){
        for(int jIdx = 0 ; jIdx < 3; jIdx++){
            Std_DRT2LPS.GetVnlMatrix()[iIdx][jIdx] = Standard_DRT2LPS [jIdx+iIdx*3];
        }
    }

    m_Projection1VTKTransform = vtkMatrix4x4::New();
    m_Projection1VTKTransform->Identity();
    m_Projection2VTKTransform = vtkMatrix4x4::New();
    m_Projection2VTKTransform->Identity();

    /* Set to NULL the metainfo */
    m_CTMetaInfo = NULL;

    m_TImage1MetaInfo = NULL;
    m_TImage2MetaInfo = NULL;

    m_DRTImage1MetaInfo = NULL;
    m_DRTImage2MetaInfo = NULL;

    m_RTMetaInfo = NULL;

    m_TransformMetaInfo = NULL;
    m_TransformMetaInfo = TransformMetaInformation::New();

    m_r23MetaInfo = NULL;
    m_r23MetaInfo = R23MetaInformation::New();

    m_PowellOptimizerMetaInfo = NULL;
    m_PowellOptimizerMetaInfo = PowellOptimizerMetaInformation::New();

    m_AmoebaOptimizerMetaInfo = NULL;
    m_AmoebaOptimizerMetaInfo = AmoebaOptimizerMetaInformation::New();

    m_MIMetricMetaInfo = NULL;
    m_MIMetricMetaInfo = MIMetricMetaInformation::New();

    m_NCCMetricMetaInfo = NULL;
    m_NCCMetricMetaInfo = NCCMetricMetaInformation::New();

    /* Initialise the projection geoemtry with defaults */
    m_DRTGeometryMetaInfo = NULL;
    m_DRTGeometryMetaInfo
            = DRTProjectionGeometryImageMetaInformation::New();
    m_DRTGeometryMetaInfo->SetSCD1(570.);
    m_DRTGeometryMetaInfo->SetSCD2(570.);
    m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(180.);
    m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(90.);
    m_DRTGeometryMetaInfo->SetIntensityThreshold(-300.);
    ImageType3D::PointType
            Point3D(3);
    Point3D[0]=0.;
    Point3D[1]=0.;
    Point3D[2]=175.;
    m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D);

    ImageType3D::SizeType
            ImageSize;
    ImageSize[0]=512;
    ImageSize[1]=512;
    ImageSize[2]=1;
    m_DRTGeometryMetaInfo->SetDRT1Size(ImageSize);
    m_DRTGeometryMetaInfo->SetDRT2Size(ImageSize);
    ImageType3D::SpacingType
            ImageRes(3);
    ImageRes[0] = 1.;
    ImageRes[1] = 1.;
    ImageRes[2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT1Spacing(ImageRes);
    m_DRTGeometryMetaInfo->SetDRT2Spacing(ImageRes);
    Point3D[0]=0.;
    Point3D[1]=0.;
    Point3D[2]=0.;
    m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Point3D);
    m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Point3D);
    m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D);
    m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D);

    m_DRTGeometryMetaInfo->SetDegreeOfFreedom(eDegreeOfFreedomType::THREE_DOF);
    m_DRTGeometryMetaInfo->SetHandleRotationImpOffset(eHandlingRotImpTransform::ALWAYS_USE);

    optimizerObserver = CommandIterationUpdate::New();
    ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();

    m_R23 = itkReg23::New();
    m_R23->SetOptimizerObserver(optimizerObserver);

}

itkImageProcessor::~itkImageProcessor()
{

}

void itkImageProcessor::SetNumberOfWorkingUnits(int iN){
    if(iN>0 && iN<50){
        iNWUnits = iN;
    }
}

const CTVolumeImageMetaInformation::Pointer
itkImageProcessor::GetCTMetaInfo(
        ){

    return m_CTMetaInfo;

}

const RTGeometryMetaInformation::Pointer
itkImageProcessor::GetRTMetaInfo(
        ){

    return m_RTMetaInfo;

}

const DRTProjectionGeometryImageMetaInformation::Pointer
itkImageProcessor::GetDRTGeoMetaInfo(){

    return m_DRTGeometryMetaInfo;

}

void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
{
    Superclass::PrintSelf(os, indent);
}

void itkImageProcessor::SetIntensityThreshold(double dT)
{
    m_DRTGeometryMetaInfo->SetIntensityThreshold(dT);
}

void itkImageProcessor::SetSCDOffset(double dOff1, double dOff2)
{
    m_DRTGeometryMetaInfo->SetSCD1Offset(dOff1);
    m_DRTGeometryMetaInfo->SetSCD2Offset(dOff2);

}

void itkImageProcessor::SetSCD(double dDist1, double dDist2)
{
    m_DRTGeometryMetaInfo->SetSCD1(dDist1);
    m_DRTGeometryMetaInfo->SetSCD2(dDist2);
}

double itkImageProcessor::GetSCD1(){
    return
            m_DRTImage1MetaInfo->GetSCD();
}
double itkImageProcessor::GetSCD2(){
    return
            m_DRTImage2MetaInfo ->GetSCD();
}

double itkImageProcessor::GetPanelOffsetPGeo(tProjOrientationType ImgPrj){

    if(m_CTMetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL)
        return 0.;


    switch (ImgPrj) {
    case tProjOrientationType::PA:
        return
            this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
                                 -m_DRTGeometryMetaInfo->GetPanel1Offset());
        break;

    case tProjOrientationType::LAT:
        return
            this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
                                 -m_DRTGeometryMetaInfo->GetPanel2Offset());
        break;

    case tProjOrientationType::NA:
        return 0.;
        break;

    default:
        return 0.;
        break;
    }

}

void itkImageProcessor::SetPanelOffsets(double dOff1, double dOff2)
{
    m_DRTGeometryMetaInfo->SetPanel1Offset(dOff1);
    m_DRTGeometryMetaInfo->SetPanel2Offset(dOff2);
}

double itkImageProcessor::GetPanelOffset1(){
    return
            m_DRTGeometryMetaInfo->GetPanel1Offset();
}

double itkImageProcessor::GetPanelOffset2(){
    return
            m_DRTGeometryMetaInfo ->GetPanel2Offset();
}

tDegreeOfFreedomEnum
itkImageProcessor::GetDegreeOfFreedom()
{
    return
            m_r23MetaInfo->GetDegreeOfFreedom();
}

void itkImageProcessor::SetUseRotationsForImportOffset(bool bVal){

    if(m_DRTGeometryMetaInfo == NULL){
        return;
    }
    m_DRTGeometryMetaInfo->SetUseRotationsForHiddenTransform(bVal);

    return;
}

void itkImageProcessor::SetCustom_ImportTransform(double dTx,
                                                  double dTy,
                                                  double dTz,
                                                  double dRx,
                                                  double dRy,
                                                  double dRz)
{

    if(m_DRTGeometryMetaInfo == NULL){
        //todo
    }

    ImageType3D::PointType
            Punto;
    Punto[0] = dTx;
    Punto[1] = dTy;
    Punto[2] = dTz;

    m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Punto);

    Punto[0] = dRx;
    Punto[1] = dRy;
    Punto[2] = dRz;

    m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Punto);

}

void itkImageProcessor::SetCustom_ProjectionCenterOffsetFixedAxes_IEC(
        double dX1,
        double dY1,
        double dZ1,
        double dX2,
        double dY2,
        double dZ2)
{

    typedef itk::Point<double, 3> PointType;
    PointType m_point;

    m_point [0] = dX1;
    m_point [1] = dY1;
    m_point [2] = dZ1;

    m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(m_point);

    m_point [0] = dX2;
    m_point [1] = dY2;
    m_point [2] = dZ2;

    m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(m_point);

}

void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(
        double dX1,
        double dY1,
        double dZ1)
{

    typedef itk::Point<double, 3> PointType;
    PointType m_point;

    m_point [0] = dX1;
    m_point [1] = dY1;
    m_point [2] = dZ1;

    m_DRTGeometryMetaInfo->SetProjectionCenter(m_point);

}

void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double nY2)
{

    if(m_DRTGeometryMetaInfo == NULL) {
        // todo
    }

    ImageType3D::SpacingType Spacing;
    Spacing [0] = nX1;
    Spacing [1] = nY1;
    Spacing [2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT1Spacing(Spacing);

    Spacing [0] = nX2;
    Spacing [1] = nY2;
    Spacing [2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT2Spacing(Spacing);

}

void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
{

    if(m_DRTGeometryMetaInfo == NULL) {
        // todo
    }

    ImageType3D::SizeType Size;
    Size [0] = nX1;
    Size [1] = nY1;
    Size [2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT1Size(Size);

    Size [0] = nX2;
    Size [1] = nY2;
    Size [2] = 1.;
    m_DRTGeometryMetaInfo->SetDRT2Size(Size);

}

void itkImageProcessor::SetCustom_UpdateMetaInfo(){

    this->UpdateProjectionGeometryMeta();
}

int itkImageProcessor::unload3DVolumeAndMeta(){

    std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;

    if (m_VolumeSourceDupli) {
        m_VolumeSourceDupli = NULL;
        m_VolumeSourceDupli = DuplicatorType::New();
    }

    m_CTMetaInfo = NULL;

    m_DRTImage1MetaInfo = NULL;
    m_DRTImage2MetaInfo = NULL;

    m_TransformMetaInfo = NULL;
    m_TransformMetaInfo = TransformMetaInformation::New();

    this->ResetROIRegions();

    return 1;
}


int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){

    tPatOrientation m_PatOrientation
            = tPatOrientation::NotDefined;

    /* Since are not sure what we get here,
     * we run the IPP sorter. */
    using IppSorterType = gdcm::IPPSorter;
    IppSorterType ZSorter;

    ZSorter.SetComputeZSpacing( true );
    ZSorter.SetZSpacingTolerance( 1e-3 );
    bool b = ZSorter.Sort(v_fnames);
    if( !b )
    {
        std::cerr << "itkImageProcessor::load3DSerieFromFiles :: Failed to sort files based on IPP." << std::endl;
        return EXIT_FAILURE;
    }

    const std::vector<std::string> & v_sortedFnames =
            ZSorter.GetFilenames();

    try
    {

        using ImageIOType = itk::GDCMImageIO;
        ImageIOType::Pointer dicomIO = ImageIOType::New();

        ImageSeriesReaderType3D::Pointer
                imageSeriesReader3D = ImageSeriesReaderType3D::New();

        imageSeriesReader3D->SetImageIO(dicomIO);
        imageSeriesReader3D->SetFileNames(v_sortedFnames);
        imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits);
        imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt

        imageSeriesReader3D->Update();

        if(v_fnames.size()>0){

            gdcm::Reader R;
            R.SetFileName(v_sortedFnames.at(0).c_str());
            if(!R.Read())
            {
                cerr<<"ERROR: cannot read file: "<<v_sortedFnames.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])){
                m_PatOrientation = eImageOrientationType::HFS;
            } else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
                m_PatOrientation = eImageOrientationType::FFS;
            } else {
                m_PatOrientation =  eImageOrientationType::NotDefined;
            }

            delete []  (sTmpString);
        } else {
            return -1;
        }

        CastFilterType3D::Pointer caster3D =
                CastFilterType3D::New();

        if(m_ResamplepCT == false){

            caster3D->SetInput(imageSeriesReader3D->GetOutput());
            caster3D->Update();
        }else{
            std::cout << "Resampling pCT @ spacing in Z " << m_SamplingLNG << std::endl;

            LinearInterpolatorType::Pointer
                    linearInterpolator = LinearInterpolatorType::New();

            // Set identity transform
            TransformType::Pointer
                    transform = TransformType::New();
            transform->SetIdentity();

            const typename ImageType3D::SpacingType& inputSpacing =
                    imageSeriesReader3D->GetOutput()->GetSpacing();
            const typename ImageType3D::RegionType& inputRegion =
                    imageSeriesReader3D->GetOutput()->GetLargestPossibleRegion();
            const typename ImageType3D::SizeType& inputSize =
                    inputRegion.GetSize();

            ImageType3D::SpacingType outputSpacing;
            outputSpacing[0] = inputSpacing[0];
            outputSpacing[1] = inputSpacing[1];
            outputSpacing[2] = m_SamplingLNG;

            ImageType3D::SizeType   outputSize;

            outputSize[0] = static_cast<SizeValueType>(
                        inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5);
            outputSize[1] = static_cast<SizeValueType>(
                        inputSize[1] * inputSpacing[1] / outputSpacing[1] + .5);
            outputSize[2] = static_cast<SizeValueType>(
                        inputSize[2] * inputSpacing[2] / outputSpacing[2] + .5);

            ResampleInputFilterType::Pointer
                    resampler = ResampleInputFilterType::New();

            resampler->SetInput( imageSeriesReader3D->GetOutput() );
            resampler->SetTransform( transform );

            resampler->SetInterpolator( linearInterpolator );
            resampler->SetDefaultPixelValue(-1024);

            resampler->SetOutputOrigin( imageSeriesReader3D->GetOutput()->GetOrigin() );
            resampler->SetOutputSpacing( outputSpacing );
            resampler->SetOutputDirection( imageSeriesReader3D->GetOutput()->GetDirection() );
            resampler->SetSize( outputSize );
            resampler->Update();


            caster3D->SetInput(resampler->GetOutput());
            caster3D->Update();

        }

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

        m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput());
        m_VolumeSourceDupli->Update();

        caster3D = NULL;


    }
    catch (itk::ExceptionObject & ex)
    {
        //        std::cout << ex << std::endl;
        return EXIT_FAILURE;
    }

    InternalImageType::Pointer m_InputImage =
            m_VolumeSourceDupli->GetOutput();
    return
            this->fill3DVolumeMeta(m_InputImage, m_PatOrientation);

}

/** Fill Meta after 3D volume load */
int itkImageProcessor::fill3DVolumeMeta(
        InternalImageType::Pointer m_InputImage,
        tPatOrientation m_PatOrientation){


    if(m_CTMetaInfo != NULL){
        std::cout << "NEVER HERE m_CTMetaInfo" << std::endl;
        m_CTMetaInfo = NULL;
    }

    if(m_DRTImage1MetaInfo != NULL){
        std::cout << "NEVER HERE m_DRTImage1MetaInfo" << std::endl;
        m_DRTImage1MetaInfo = NULL;
    }

    if(m_DRTImage2MetaInfo != NULL){
        std::cout << "NEVER HERE m_DRTImage2MetaInfo" << std::endl;
        m_DRTImage2MetaInfo = NULL;
    }

    if(m_RTMetaInfo != NULL){
        std::cout << "NEVER HERE m_RTMetaInfo" << std::endl;
        m_RTMetaInfo = NULL;
    }

    /* copy useful meta information into the CT container */
    m_CTMetaInfo = CTVolumeImageMetaInformation::New();
    m_CTMetaInfo->SetPatientOrientation(m_PatOrientation);
    m_CTMetaInfo->SetSpacing(m_InputImage->GetSpacing());
    m_CTMetaInfo->SetSize(
                m_InputImage->GetBufferedRegion().GetSize() );
    std::cout<<"ORIGIN: "<<m_InputImage->GetOrigin()<<std::endl;
    m_CTMetaInfo->SetOriginLPS(m_InputImage->GetOrigin());
    m_CTMetaInfo->SetImageDirections(m_InputImage->GetDirection());
    ImageType3D::PointType
            PointOffset;
    PointOffset.Fill(0.);
    m_CTMetaInfo->SetImportOffset(PointOffset);

    /* initialise DRT meta */
    m_DRTImage1MetaInfo = DRTImageMetaInformation::New();
    m_DRTImage1MetaInfo->SetProjectionAngleLPS(
                this->CalcProjectionAngleLPS(
                    m_CTMetaInfo->GetPatientOrientation(),
                    m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() +
                    m_DRTGeometryMetaInfo->GetProjectionAngle1OffsetIEC())
                );

    m_DRTImage1MetaInfo->SetSCD(
                m_DRTGeometryMetaInfo->GetSCD1() +
                m_DRTGeometryMetaInfo->GetSCD1Offset()
                );

    /* Calculate projection angle IEC to LPS */
    m_DRTImage2MetaInfo = DRTImageMetaInformation::New();
    m_DRTImage2MetaInfo->SetProjectionAngleLPS(
                this->CalcProjectionAngleLPS(
                    m_CTMetaInfo->GetPatientOrientation(),
                    m_DRTGeometryMetaInfo->GetProjectionAngle2IEC() +
                    m_DRTGeometryMetaInfo->GetProjectionAngle2OffsetIEC())
                );

    m_DRTImage2MetaInfo->SetSCD(
                m_DRTGeometryMetaInfo->GetSCD2()+
                m_DRTGeometryMetaInfo->GetSCD2Offset()
                );

    std::cout<< "ImportOffset"<< m_CTMetaInfo->GetImportOffset() <<std::endl;

    this->UpdateProjectionGeometryMeta();

    return EXIT_SUCCESS;

}

const std::vector <double>  itkImageProcessor::GetRTImportOffset()
{
    std::vector <double>  vOffset;
    vOffset.clear();

    vOffset.push_back(
                m_CTMetaInfo->GetImportOffset()[0]
            );

    vOffset.push_back(
                m_CTMetaInfo->GetImportOffset()[1]
            );

    vOffset.push_back(
                m_CTMetaInfo->GetImportOffset()[2]
            );
    return
            vOffset;
}

const std::vector <double> itkImageProcessor::GetLPStoProjectionGeoLPSOffset()
{
    std::vector <double>  vOffset;
    vOffset.clear();
    vOffset.push_back(
                m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[0]
            );
    vOffset.push_back(
                m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[1]
            );
    vOffset.push_back(
                m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[2]
            );
    return
            vOffset;
}

double
itkImageProcessor::CalcPanelOffsetPGeo(
        tPatOrientation pOrient,
        double dPOffset){

    double dVal = dPOffset;//0.;

    switch (pOrient) {
        case tPatOrientation :: HFS:

        break;

    case tPatOrientation :: FFS:
        dVal *= -1;
        break;

    default:
        break;

    }

    return
            dVal;
}

double
itkImageProcessor::CalcProjectionAngleLPS(
        tPatOrientation  pOrient,
        double pAngleIEC){

    double dProj = pAngleIEC;

    InternalImageType::DirectionType
            currIECtoLPS;

    /* NOTE that we transpose on the fly... */
    switch (pOrient) {
    case tPatOrientation :: HFS:
        for(int iIdx = 0 ; iIdx < 3; iIdx++){
            for(int jIdx = 0 ; jIdx < 3; jIdx++){
                currIECtoLPS.GetVnlMatrix()[jIdx][iIdx] = HFS2IEC [jIdx+iIdx*3];
            }
        }
        break;
    case tPatOrientation :: FFS:
        for(int iIdx = 0 ; iIdx < 3; iIdx++){
            for(int jIdx = 0 ; jIdx < 3; jIdx++){
                currIECtoLPS.GetVnlMatrix()[jIdx][iIdx] = FFS2IEC [jIdx+iIdx*3];
            }
        }
        break;
    default:
        break;

    }

    InternalImageType::DirectionType
            DRTDirsIEC;

    DRTDirsIEC.SetIdentity();
    DRTDirsIEC. GetVnlMatrix()[0][0] =  cos (pAngleIEC * itk::Math::pi/180);
    DRTDirsIEC. GetVnlMatrix()[0][2] =  sin (pAngleIEC * itk::Math::pi/180);
    DRTDirsIEC. GetVnlMatrix()[2][0] = -sin (pAngleIEC * itk::Math::pi/180);
    DRTDirsIEC. GetVnlMatrix()[2][2] =  cos (pAngleIEC * itk::Math::pi/180);


    InternalImageType::DirectionType
            DRTtoCurrPatientOrient;

    DRTtoCurrPatientOrient =
            currIECtoLPS * DRTDirsIEC;

    ImageType3D::PointType SourcePositionIEC;
    SourcePositionIEC [0] = 0;
    SourcePositionIEC [1] = 0;
    SourcePositionIEC [2] = 1;

    ImageType3D::PointType SourcePositionPatOrient =
            DRTtoCurrPatientOrient * SourcePositionIEC;

    dProj =
            (atan2(SourcePositionPatOrient [1], SourcePositionPatOrient [0]) - atan2(-1,0)) * 180 / itk::Math::pi;

    return
            dProj;

}

void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2)
{
    m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1);
    m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2);
}

void itkImageProcessor::SetProjectionAngle1IEC(double dGantryAngle)
{
    m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(dGantryAngle);
}

double itkImageProcessor::GetProjectionAngle1LPS(){

    return
            m_DRTImage1MetaInfo->GetProjectionAngleLPS();
}

void itkImageProcessor::SetProjectionAngle2IEC(double dGantryAngle)
{
    m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(dGantryAngle);
}

double itkImageProcessor::GetProjectionAngle2LPS(){
    return
            m_DRTImage2MetaInfo->GetProjectionAngleLPS();
}

int itkImageProcessor::unload2DAndMeta(int iImgType){


    switch (iImgType) {
    case eProjectionOrientationType::PA:

        std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << " PA "<< std::endl;
        m_TImage1MetaInfo = NULL;

        if (m_PASourceDupli) {
            m_PASourceDupli = NULL;
            m_PASourceDupli = DuplicatorType::New();
        }
        break;

    case eProjectionOrientationType::LAT:

        std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << " LAT "<< std::endl;
        m_TImage2MetaInfo = NULL;

        if (m_LATSourceDupli) {
            m_LATSourceDupli = NULL;
            m_LATSourceDupli = DuplicatorType::New();
        }

        break;
    }

    return 1;

}

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();

    /* Image orientation */
    eImageOrientationType
            currImgOrient;
    strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0018,0x5100), ds));
    *std::remove(sTmpString, sTmpString + strlen(sTmpString), ' ') = 0;

    if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
        currImgOrient = eImageOrientationType::HFS;
    } else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
        currImgOrient = eImageOrientationType::FFS;
    } else {
        std::cerr<< "Image Orientation: unrecognised"<< sTmpString <<std::endl;
    }

    free (sTmpString);

    /* Actual file read */
    ImageReaderType2D::Pointer imageReader2D;
    imageReader2D = ImageReaderType2D::New();
    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;
    }

    /* identifing a projection to be lateral or posterior requires the interpretation of direction
     * and patient orientation */

    using ImageDirectionType3D = ImageType3D::DirectionType;

    ImageDirectionType3D LocalizerImagDirectionDCM =
            imageReader2D->GetOutput()->GetDirection();

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

    /* calculate image orientation with respect to fixed IEC */
    InternalImageType::DirectionType LocalizerImagDirectionIEC =
            toIECDirection * LocalizerImagDirectionDCM;

    /* we calculate dot products between the Z components */
    double dSumPA = 0;
    dSumPA =
            (LocalizerImagDirectionIEC[0][2] * PAElementsIEC[2]) +
            (LocalizerImagDirectionIEC[1][2] * PAElementsIEC[5]) +
            (LocalizerImagDirectionIEC[2][2] * PAElementsIEC[8]);

    double dSumLAT = 0;
    dSumLAT=
            (LocalizerImagDirectionIEC[0][2] * LATElementsIEC[2]) +
            (LocalizerImagDirectionIEC[1][2] * LATElementsIEC[5]) +
            (LocalizerImagDirectionIEC[2][2] * LATElementsIEC[8]);

    tProjOrientationType
            currProjOrientation = NA;

    /* dot shoudl be 1 or -1 if the normal is parallel */
    if(dSumPA == 1 || dSumPA == -1){
        currProjOrientation = PA;
        //        std::cout<<"It's a PA image"<<std::endl;
    } else if (dSumLAT == 1 || dSumLAT == -1){
        currProjOrientation = LAT;
        //        std::cout<<"It's a LAT image"<<std::endl;
    } else {
        return -1;
    }

    /* cast type to internal processing type */
    CastFilterType3D::Pointer caster2DInput;
    caster2DInput = CastFilterType3D::New();
    caster2DInput->SetInput(imageReader2D->GetOutput());
    caster2DInput->Update();

    InternalImageType::Pointer MyLocalizerImage =
            caster2DInput->GetOutput();

    DuplicatorType::Pointer m_Duplicator;
    TopogramImageMetaInformation::Pointer m_TImageMeta;

    switch (currProjOrientation) {
    case eProjectionOrientationType::PA:

        if (m_PASourceDupli) {
            m_PASourceDupli = NULL;
            m_PASourceDupli = DuplicatorType::New();
        }

        m_Duplicator = m_PASourceDupli;
        m_TImageMeta = m_TImage1MetaInfo;

        break;

    case eProjectionOrientationType::LAT:

        if (m_LATSourceDupli) {
            m_LATSourceDupli = NULL;
            m_LATSourceDupli = DuplicatorType::New();
        }

        m_Duplicator = m_LATSourceDupli;
        m_TImageMeta = m_TImage2MetaInfo;

        break;

    default:
        return -1;
        break;
    }

    if(m_TImageMeta != NULL){
        //TODO
    }


    m_TImageMeta = TopogramImageMetaInformation::New();
    m_TImageMeta->SetWLLevel(dIntensityWindow[0]);
    m_TImageMeta->SetWLWindow(dIntensityWindow[1]);
    m_TImageMeta->SetPatientOrientation(
                currImgOrient);
    m_TImageMeta->SetProjectionOrientation(
                currProjOrientation);

    m_Duplicator->SetInputImage(caster2DInput->GetOutput() );
    m_Duplicator->Update();

    return
            (int) currProjOrientation;

}

double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
{


    TopogramImageMetaInformation::Pointer m_TImageMeta;

    switch (iImg) {

    case 0:
        m_TImageMeta = m_TImage1MetaInfo;
        break;

    case 1:
        m_TImageMeta = m_TImage2MetaInfo;

        break;

    default:
        return 0;
        break;

    }

    if(m_TImageMeta == NULL){
        return 0;
    } else {
        return
                m_TImageMeta->GetWLLevel();
    }

    return 0.;
}

double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
{
    TopogramImageMetaInformation::Pointer m_TImageMeta;

    switch (iImg) {

    case 0:
        m_TImageMeta = m_TImage1MetaInfo;
        break;

    case 1:
        m_TImageMeta = m_TImage2MetaInfo;

        break;

    default:
        return 0;
        break;

    }

    if(m_TImageMeta == NULL){
        return 0;
    } else {
        return
                m_TImageMeta->GetWLWindow();
    }

    return 0.;
}

Optimizer::ParametersType
itkImageProcessor::GetFinalR23Parameters(){


    if (!m_TransformMetaInfo) {
        itkExceptionMacro(<< "m_TransformMetaInfo not present");
    }

    Optimizer::ParametersType pPars(6);
    pPars.fill(0.);

    itk::Optimizer::ParametersType currentPosition =
            m_R23->GetCurrentPosition();

    switch (m_r23MetaInfo->GetDegreeOfFreedom()) {
    case THREE_DOF:
        pPars[3] = currentPosition[0]
                - m_TransformMetaInfo->GetHiddenTranslations()[0];
        pPars[4] = currentPosition[1]
                - m_TransformMetaInfo->GetHiddenTranslations()[1];
        pPars[5] = currentPosition[2]
                - m_TransformMetaInfo->GetHiddenTranslations()[2];
        break;

    case SIX_DOF:
        pPars[3] = currentPosition[0]
                - m_TransformMetaInfo->GetHiddenTranslations()[0];
        pPars[4] = currentPosition[1]
                - m_TransformMetaInfo->GetHiddenTranslations()[1];
        pPars[5] = currentPosition[2]
                - m_TransformMetaInfo->GetHiddenTranslations()[2];

        pPars[0] = currentPosition[3]
                - m_TransformMetaInfo->GetHiddenRotations()[0];
        pPars[1] = currentPosition[4]
                - m_TransformMetaInfo->GetHiddenRotations()[1];
        pPars[2] = currentPosition[5]
                - m_TransformMetaInfo->GetHiddenRotations()[2];

        break;

    default:
        pPars.fill(0.);
        break;
    }


    return
            pPars;

}

void itkImageProcessor::SetOptimizer(tOptimizerTypeEnum eOpti){
    if(eOpti != tOptimizerTypeEnum::AMOEBA &&
       eOpti != tOptimizerTypeEnum::POWELL){
        itkExceptionMacro(<< "Unkown optimzer : " << eOpti);
        return;
    }
    std::cout<<"Setting Optimizer: "<<eOpti<<std::endl;
    m_r23MetaInfo->SetOptimizerType(eOpti);

}

void itkImageProcessor::SetMetric(tMetricTypeEnum eMetric){
    if(eMetric != tMetricTypeEnum::NCC &&
       eMetric != tMetricTypeEnum::MI){
        itkExceptionMacro(<< "Unkown metric string : " << eMetric);
        return;
    }
    std::cout<<"Setting Metric: "<<eMetric<<std::endl;
    m_r23MetaInfo->SetMetricType(eMetric);

}

void itkImageProcessor::SetOptimizer(std::string optimizer)
{

    if (optimizer.compare("Powell") == 0) {
        m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::POWELL);
    } else if (optimizer.compare("Amoeba") == 0) {
        m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::AMOEBA);
    } else if (optimizer.compare("Exhaustive") == 0) {
        m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::EXHAUSTIVE);
    } else {
        itkExceptionMacro(<< "Unkown optimzer string : " << optimizer);
    }


}

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

    if (metric.compare("NCC") == 0) {
        m_r23MetaInfo->SetMetricType(tMetricTypeEnum::NCC);
    } else if (metric.compare("MI") == 0) {
        m_r23MetaInfo->SetMetricType(tMetricTypeEnum::MI);
    } else {
        itkExceptionMacro(<< "Unkown metric string : " << metric);
    }


}

void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof){

    m_DRTGeometryMetaInfo->SetDegreeOfFreedom(dof);
    m_r23MetaInfo->SetDegreeOfFreedom(dof);

}

void itkImageProcessor::SetHandleRotationImportOffset(eHandlingRotImpTransform hrio){

    m_DRTGeometryMetaInfo->SetHandleRotationImpOffset(hrio);

}

void itkImageProcessor::SetPowellOptimParameters(
        double dStepT,
        double dValTol,
        double dStepL,
        int iMaxLinI,
        int iMaxIter){

    m_PowellOptimizerMetaInfo->SetStepTolerance(dStepT);
    m_PowellOptimizerMetaInfo->SetValueTolerance(dValTol);
    m_PowellOptimizerMetaInfo->SetStepLength(dStepL);
    m_PowellOptimizerMetaInfo->SetMaximumLineInteration(iMaxLinI);
    m_PowellOptimizerMetaInfo->SetMaxIterations(iMaxIter);

}

void itkImageProcessor::SetAmoebaOptimParameters(
        double dParConvTol,
        double dFuntConvTol,
        double dSimplex,
        int iMaxIter){
    m_AmoebaOptimizerMetaInfo->SetParametersConvergenceTolerance(dParConvTol);
    m_AmoebaOptimizerMetaInfo->SetFunctionConvergenceTolerance(dFuntConvTol);
    m_AmoebaOptimizerMetaInfo->SetSimplexDelta(dSimplex);
    m_AmoebaOptimizerMetaInfo->SetMaxIterations(iMaxIter);

}

void itkImageProcessor::SetMIMetricParameters(double dMaxT,int iNhb){
    m_MIMetricMetaInfo->SetMaxTranslation(dMaxT);
    m_MIMetricMetaInfo->SetNumberOfHistogramBins(iNhb);
}

void itkImageProcessor::SetNCCMetricParameters(double dMaxT,bool bSm){
    m_NCCMetricMetaInfo->SetMaxTranslation(dMaxT);
    m_NCCMetricMetaInfo->SetSubtractMean(bSm);
}

void itkImageProcessor::InitializeProjector()
{

    if(m_DRTImage1MetaInfo == NULL ||
            m_DRTImage2MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
            m_TransformMetaInfo == NULL     ){

        //TODO
        return;
    }

    const double dtr = (atan(1.0) * 4.0) / 180.0;


    ImageType3D::PointType ZeroPoint;
    ZeroPoint.Fill(0.);

    TransformType::Pointer CurrTransform;

    CurrTransform=    CalculateInternalTransformV3(
                m_TransformMetaInfo->GetT(),
                m_TransformMetaInfo->GetR(),
                m_InternalTransf1->GetpCalProjCenter(),
                m_InternalTransf1->GetpRTIsocenter(),
                m_InternalTransf1->GetIECtoLPSDirs());


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

    transform1->SetTranslation(
                CurrTransform->GetTranslation());
    transform1->SetRotation(
                CurrTransform->GetAngleX(),
                CurrTransform->GetAngleY(),
                CurrTransform->GetAngleZ()
                );

    transform1->SetCenter(
                m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );

    // 2D Image 1
    interpolator1->SetProjectionAngle(
                dtr * m_DRTImage1MetaInfo->GetProjectionAngleLPS() );
    interpolator1->SetFocalPointToIsocenterDistance(
                m_DRTImage1MetaInfo->GetSCD());
    interpolator1->SetThreshold(
                m_DRTGeometryMetaInfo->GetIntensityThreshold() );
    interpolator1->SetPanelOffset(
                m_DRTImage1MetaInfo->GetPanelOffsetPGeo());

    interpolator1->SetTransform(transform1);
    interpolator1->Initialize();

    CurrTransform=    CalculateInternalTransformV3(
                m_TransformMetaInfo->GetT(),
                m_TransformMetaInfo->GetR(),
                m_InternalTransf2->GetpCalProjCenter(),
                m_InternalTransf2->GetpRTIsocenter(),
                m_InternalTransf2->GetIECtoLPSDirs());


    transform2->SetComputeZYX(true);
    transform2->SetIdentity();

    transform2->SetTranslation(
                CurrTransform->GetTranslation());
    transform2->SetRotation(
                CurrTransform->GetAngleX(),
                CurrTransform->GetAngleY(),
                CurrTransform->GetAngleZ()
                );

    transform2->SetCenter(
                m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() );

    // 2D Image 2
    interpolator2->SetProjectionAngle(
                dtr * m_DRTImage2MetaInfo->GetProjectionAngleLPS() );
    interpolator2->SetFocalPointToIsocenterDistance(
                m_DRTImage2MetaInfo->GetSCD() );
    interpolator2->SetThreshold(
                m_DRTGeometryMetaInfo->GetIntensityThreshold() );
    interpolator2->SetPanelOffset(
                m_DRTImage2MetaInfo->GetPanelOffsetPGeo());

    interpolator2->SetTransform(transform2);
    interpolator2->Initialize();


    resampleFilter1 = ResampleFilterType::New();
    resampleFilter1->SetInput(
                m_VolumeSourceDupli->GetOutput());

    ImageType3D::RegionType lagerReg =
        m_VolumeSourceDupli->GetOutput()->GetLargestPossibleRegion();

    resampleFilter1->SetDefaultPixelValue(0);
    resampleFilter1->SetNumberOfWorkUnits(iNWUnits);

    // 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(transform1);
    interpolator1->Initialize();
    resampleFilter1->SetInterpolator(interpolator1);


    resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize());
    resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing());
    resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin());


    resampleFilter2 = ResampleFilterType::New();
    resampleFilter2->SetInput(
                m_VolumeSourceDupli->GetOutput());

    resampleFilter2->SetDefaultPixelValue(0);
    resampleFilter2->SetNumberOfWorkUnits(iNWUnits);


    // 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(transform2);
    interpolator2->Initialize();
    resampleFilter2->SetInterpolator(interpolator2);

    resampleFilter2->SetSize(m_DRTImage2MetaInfo->GetSize());
    resampleFilter2->SetOutputSpacing(m_DRTImage2MetaInfo->GetSpacing());
    resampleFilter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOrigin());

    filter1 =
            ChangeInformationFilterType::New();

    filter2 =
            ChangeInformationFilterType::New();

    /* First of all we need the center of the Projection images in its reference frame */
    resampleFilter1->Update();
    filter1->SetInput( resampleFilter1->GetOutput() );
    filter1->SetOutputDirection(m_DRTImage1MetaInfo->GetImageDirectionsLPS() );//IECtoLPS_Directions);
    filter1->ChangeDirectionOn();
    filter1->SetOutputOrigin( m_DRTImage1MetaInfo->GetOriginLPS() );
    filter1->ChangeOriginOn();
    filter1->UpdateOutputInformation();

    /* Again.. */
    resampleFilter2->Update();
    filter2->SetInput( resampleFilter2 ->GetOutput() );
    filter2->SetOutputDirection(m_DRTImage2MetaInfo->GetImageDirectionsLPS() );
    filter2->ChangeDirectionOn();
    filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() );
    filter2->ChangeOriginOn();
    filter2->UpdateOutputInformation();

    filter1->Update();
    filter2->Update();

    imageDRT1In = filter1->GetOutput();
    imageDRT2In = filter2->GetOutput();

}

itkImageProcessor::InternalImageType::DirectionType
itkImageProcessor::CalcDRTImageDirections(
        double angle,
        InternalImageType::DirectionType DRT2LPS){


    InternalImageType::DirectionType
            DRTDirs;

    DRTDirs.SetIdentity();
    DRTDirs. GetVnlMatrix()[0][0] =  cos (angle*itk::Math::pi/180);
    DRTDirs. GetVnlMatrix()[0][1] = -sin (angle*itk::Math::pi/180);
    DRTDirs. GetVnlMatrix()[1][0] =  sin (angle*itk::Math::pi/180);
    DRTDirs. GetVnlMatrix()[1][1] =  cos (angle*itk::Math::pi/180);

    DRTDirs *=  DRT2LPS;

    return
            DRTDirs;

}

int itkImageProcessor::unloadRTPlanAndMeta(){

    std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
    m_RTMetaInfo = NULL;

    ImageType3D::PointType pZero;
    pZero.Fill(0.);

    m_CTMetaInfo->SetImportOffset(pZero);
    m_TransformMetaInfo->SetHiddenRotations(pZero);

    return 0;
}

void itkImageProcessor::loadRTPlanInfo(
        double dIsoX, double dIsoY, double dIsoZ,
        double dLAT, double dVRT ,double dLNG){

    if(m_RTMetaInfo != NULL){
        std::cout << " NEVER HERE loadRTPlanInfo m_RTMetaInfo" << std::endl;
        //TODO
    }

    if(m_CTMetaInfo != NULL){
        std::cout << " ALWAYS HERE loadRTPlanInfo m_CTMetaInfo" << std::endl;
        //TODO
    }

    if(m_DRTGeometryMetaInfo != NULL){
        std::cout << " ALWAYS HERE loadRTPlanInfo m_DRTGeometryMetaInfo" << std::endl;
        //TODO
    }

    m_RTMetaInfo = RTGeometryMetaInformation::New();

    ImageType3D::PointType
            Point;
    Point[0] = dIsoX;
    Point[1] = dIsoY;
    Point[2] = dIsoZ;

    m_RTMetaInfo->SetIsocenterLPS(Point);

    Point[0] = dLAT;
    Point[1] = dLNG;
    Point[2] = dVRT;

    m_RTMetaInfo->SetIsocenterIECS(Point);

    ImageType3D::PointType
            pZero (3);
    pZero.Fill(0.);

    m_CTMetaInfo->SetImportOffset(
                CalcImportVolumeOffset(
                    m_RTMetaInfo->GetIsocenterIECS(),
                    m_CTMetaInfo->GetLPS2IECDirections(),
                    m_RTMetaInfo->GetIsocenterLPS(),
                    m_DRTGeometryMetaInfo->GetIECS2IECScannerT(),
                    (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ?
                     m_DRTGeometryMetaInfo->GetIECS2IECScannerR()  : pZero)
                    )
                );

    m_TransformMetaInfo->SetHiddenRotations(
                (m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ?
                 m_DRTGeometryMetaInfo->GetIECS2IECScannerR()  : pZero)
                );


    this->UpdateProjectionGeometryMeta();
    this->InitializeProjector();

}

void itkImageProcessor::UpdateProjectionGeometryMeta(){

    if(m_CTMetaInfo == NULL){
        return;
        //TODO.
    }

    if(m_DRTGeometryMetaInfo == NULL){
        return;
        //TODO.
    }

    if(m_DRTImage1MetaInfo == NULL){
        return;
        //TODO.
    }

    if(m_DRTImage2MetaInfo == NULL){
        return;
        //TODO.
    }

    // FIRST

    ImageType3D::PointType NominalIsocenter_wZcorrectionLPS;

    NominalIsocenter_wZcorrectionLPS =
            m_CTMetaInfo->GetProjectionOriginLPS(
                m_DRTGeometryMetaInfo->GetProjectionCenter()
                );

    ImageType3D::PointType NominalIsocenterZero_wZcorrectionLPS;

    NominalIsocenterZero_wZcorrectionLPS =
            m_CTMetaInfo->GetProjectionOriginLPSZero(
                m_DRTGeometryMetaInfo->GetProjectionCenter()
                );


    ImageType3D::PointType IsocenterOffsetLPS;

    IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint(
                m_DRTGeometryMetaInfo->GetProjectionCenterOffset1());

    ImageType3D::PointType CalibratedIsocenterLPS;

    CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] +
            IsocenterOffsetLPS[0];
    CalibratedIsocenterLPS[1] = NominalIsocenter_wZcorrectionLPS[1] +
            IsocenterOffsetLPS[1];
    CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2];


    ImageType3D::PointType CalibratedIsocenterZeroLPS;

    CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] +
            IsocenterOffsetLPS[0] /*+ m_CTMetaInfo->GetSpacing()[0]/2*/;
    CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] +
            IsocenterOffsetLPS[1] /*+ m_CTMetaInfo->GetSpacing()[1]/2*/;
    CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2] /*+ m_CTMetaInfo->GetSpacing()[2]/2*/;

    m_DRTImage1MetaInfo->SetProjectionAngleLPS(
                this->CalcProjectionAngleLPS(
                    m_CTMetaInfo->GetPatientOrientation(),
                    m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() +
                    m_DRTGeometryMetaInfo->GetProjectionAngle1OffsetIEC())
                );

    m_DRTImage1MetaInfo->SetSCD(
                m_DRTGeometryMetaInfo->GetSCD1() +
                m_DRTGeometryMetaInfo->GetSCD1Offset()
                );

    m_DRTImage1MetaInfo->SetProjectionOriginLPS(
                CalibratedIsocenterLPS);

    m_DRTImage1MetaInfo->SetProjectionOriginLPSZero(
                CalibratedIsocenterZeroLPS);

    m_DRTImage1MetaInfo->SetSizeWithBounds(
                NULL,//vBounds.data(),
                m_DRTGeometryMetaInfo->GetDRT1Size(),
                m_DRTGeometryMetaInfo->GetDRT1Spacing()
                );

    m_DRTImage1MetaInfo->SetPanelOffsetPGeo(
                this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
                                         -m_DRTGeometryMetaInfo->GetPanel1Offset())
                );

    // This HAS TO be calculated from the nominal isocenter.
    m_DRTImage1MetaInfo->SetOriginLPS(
                CalcDRTImageOrigin(
                    m_DRTImage1MetaInfo->GetOrigin(),
                    m_DRTImage1MetaInfo->GetCOV(),
                    NominalIsocenter_wZcorrectionLPS,
                    this->CalcProjectionAngleLPS(
                        m_CTMetaInfo->GetPatientOrientation(),
                        m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ),
                    Std_DRT2LPS
                    )

                );

    m_DRTImage1MetaInfo->SetImageDirectionsLPS(
                CalcDRTImageDirections(
                    this->CalcProjectionAngleLPS(
                        m_CTMetaInfo->GetPatientOrientation(),
                        m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()),
                    Std_DRT2LPS)
                );

    if(m_RTMetaInfo == NULL)
    {
        m_InternalTransf1->SetpCalProjCenter(
                    m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
        m_InternalTransf1->SetpRTIsocenter(
                    m_DRTImage1MetaInfo->GetProjectionOriginLPS());
        InternalImageType::DirectionType IECtoLPS_Directions;
        IECtoLPS_Directions =
                m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
        m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions);

    } else {
        m_InternalTransf1->SetpCalProjCenter(
                    m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
        m_InternalTransf1->SetpRTIsocenter(
                    m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS());
        InternalImageType::DirectionType IECtoLPS_Directions;
        IECtoLPS_Directions =
                m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
        m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions);

    }

    // SECOND


    IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint(
                m_DRTGeometryMetaInfo->GetProjectionCenterOffset2());

    CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] +
            IsocenterOffsetLPS[0];
    CalibratedIsocenterLPS[1] = NominalIsocenter_wZcorrectionLPS[1] +
            IsocenterOffsetLPS[1];
    CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2];


    CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] +
            IsocenterOffsetLPS[0];
    CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] +
            IsocenterOffsetLPS[1];
    CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] +
            IsocenterOffsetLPS[2];


    m_DRTImage2MetaInfo->SetProjectionAngleLPS(
                this->CalcProjectionAngleLPS(
                    m_CTMetaInfo->GetPatientOrientation(),
                    m_DRTGeometryMetaInfo->GetProjectionAngle2IEC() +
                    m_DRTGeometryMetaInfo->GetProjectionAngle2OffsetIEC())
                );

    m_DRTImage2MetaInfo->SetSCD(
                m_DRTGeometryMetaInfo->GetSCD2()+
                m_DRTGeometryMetaInfo->GetSCD2Offset()
                );


    m_DRTImage2MetaInfo->SetProjectionOriginLPS(
                CalibratedIsocenterLPS);

    m_DRTImage2MetaInfo->SetProjectionOriginLPSZero(
                CalibratedIsocenterZeroLPS);


    m_DRTImage2MetaInfo->SetSizeWithBounds(
                NULL,//vBounds.data(),
                m_DRTGeometryMetaInfo->GetDRT2Size(),
                m_DRTGeometryMetaInfo->GetDRT2Spacing()
                );

    m_DRTImage2MetaInfo->SetPanelOffsetPGeo(
                this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
                                         -m_DRTGeometryMetaInfo->GetPanel2Offset()));

    m_DRTImage2MetaInfo->SetOriginLPS(
                CalcDRTImageOrigin(
                    m_DRTImage2MetaInfo->GetOrigin(),
                    m_DRTImage2MetaInfo->GetCOV(),
                    NominalIsocenter_wZcorrectionLPS,
                    this->CalcProjectionAngleLPS(
                        m_CTMetaInfo->GetPatientOrientation(),
                        m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ,
                    Std_DRT2LPS
                    ));


    m_DRTImage2MetaInfo->SetImageDirectionsLPS(
                CalcDRTImageDirections(
                    this->CalcProjectionAngleLPS(
                        m_CTMetaInfo->GetPatientOrientation(),
                        m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ,
                    Std_DRT2LPS)
                );

    if(m_RTMetaInfo == NULL)
    {
        m_InternalTransf2->SetpCalProjCenter(
                    m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
        m_InternalTransf2->SetpRTIsocenter(
                    m_DRTImage2MetaInfo->GetProjectionOriginLPS());
        InternalImageType::DirectionType IECtoLPS_Directions;
        IECtoLPS_Directions =
                m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
        m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions);

    } else {
        m_InternalTransf2->SetpCalProjCenter(
                    m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
        m_InternalTransf2->SetpRTIsocenter(
                    m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS());
        InternalImageType::DirectionType IECtoLPS_Directions;
        IECtoLPS_Directions =
                m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
        m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions);

    }

}

/*
    Nobody should go through all this...
    To be very careful editing ...
*/
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcDRTImageOrigin(
        ImageType3D::PointType m_DRTOrigin,
        ImageType3D::PointType m_DRTCOV,
        ImageType3D::PointType m_ProjectionCenterLPS,
        double dAngle,
        InternalImageType::DirectionType stdDRT2LPS
        ){

    itkImageProcessor::InternalImageType::DirectionType DRT2LPS
            = this->CalcDRTImageDirections(dAngle, stdDRT2LPS);

    ImageType3D::PointType NewOrigin =
            m_ProjectionCenterLPS + DRT2LPS * (m_DRTOrigin - m_DRTCOV);

    return
            NewOrigin;
}

void itkImageProcessor::GetProjectionImages(){

    if(m_DRTImage1MetaInfo == NULL ||
            m_DRTImage2MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
            m_TransformMetaInfo == NULL     ){
        return;
        //TODO
        return;
    }

    ImageType3D::PointType ZeroPoint;
    ZeroPoint.Fill(0.);

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

    TransformType::Pointer CurrTransform;
    CurrTransform=    CalculateInternalTransformV3(
                m_TransformMetaInfo->GetT(),
                m_TransformMetaInfo->GetR(),
                m_InternalTransf1->GetpCalProjCenter(),
                m_InternalTransf1->GetpRTIsocenter(),
                m_InternalTransf1->GetIECtoLPSDirs());


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

    transform1->SetTranslation(
                CurrTransform->GetTranslation());
    transform1->SetRotation(
                CurrTransform->GetAngleX(),
                CurrTransform->GetAngleY(),
                CurrTransform->GetAngleZ()
                );
    transform1 ->SetCenter(
                m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());

    // 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(transform1);
    interpolator1->Initialize();


    transform2->SetComputeZYX(true);
    transform2->SetIdentity();
    CurrTransform = CalculateInternalTransformV3(
                m_TransformMetaInfo->GetT(),
                m_TransformMetaInfo->GetR(),
                m_InternalTransf2->GetpCalProjCenter(),
                m_InternalTransf2->GetpRTIsocenter(),
                m_InternalTransf2->GetIECtoLPSDirs());


    transform2->SetComputeZYX(true);
    transform2->SetIdentity();

    transform2->SetTranslation(
                CurrTransform->GetTranslation());
    transform2->SetRotation(
                CurrTransform->GetAngleX(),
                CurrTransform->GetAngleY(),
                CurrTransform->GetAngleZ()
                );



    transform2 ->SetCenter(
                m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());


    //    // 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(transform2); //finalTransform);
    interpolator2->Initialize();

    filter1->Update();
    filter2->Update();

    imageDRT1In = filter1->GetOutput();
    imageDRT2In = filter2->GetOutput();

}

itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcImportVolumeOffset(
        ImageType3D::PointType rtCouchOffset,
        InternalImageType::DirectionType VolumeLPStoIEC,
        ImageType3D::PointType rtIsocenterLPS,
        ImageType3D::PointType IEC2DCMMapT,
        ImageType3D::PointType IEC2DCMMapR){



    TransformType::Pointer InputTransform = TransformType::New();
    InputTransform->SetComputeZYX(true);
    InputTransform->SetIdentity();

    TransformType::OutputVectorType translation;
    translation[0] = IEC2DCMMapT[0];
    translation[1] = IEC2DCMMapT[1];
    translation[2] = IEC2DCMMapT[2];

    InputTransform->SetTranslation(translation);

    const double dtr = (atan(1.0) * 4.0) / 180.0;
    InputTransform->SetRotation(
                dtr * IEC2DCMMapR[0],
            dtr * IEC2DCMMapR[1],
            dtr * IEC2DCMMapR[2]);

    ImageType3D::PointType m_TransformOrigin;
    m_TransformOrigin.Fill(0.);
    InputTransform->SetCenter(
                m_TransformOrigin );

    ImageType3D::PointType IsoSupport_IECPos =
            InputTransform->TransformPoint(rtCouchOffset );

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


    ImageType3D::PointType IsoSupport_LPSPos =
            VolumeIECtoLPS * IsoSupport_IECPos;

    ImageType3D::PointType Offset =
            rtIsocenterLPS - IsoSupport_LPSPos;

    std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ <<
                 " " << Offset << std::endl;

    return
            Offset;


}


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(m_PASourceDupli->GetOutput());

    //        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(m_LATSourceDupli->GetOutput());

    //        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()
{

    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

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

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

    return
            toVTKLocalizer1->GetOutput();
}

vtkImageData* itkImageProcessor::GetLocalizer2VTK()
{

    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

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

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

    return
            toVTKLocalizer2->GetOutput();
}


vtkImageData* itkImageProcessor::GetProjection1VTK()
{

    if(m_DRTImage1MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
            m_TransformMetaInfo == NULL     ){
        return NULL;
        //TODO
    }

    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

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

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

    return
            toVTK2D1->GetOutput();
}

vtkImageData* itkImageProcessor::GetProjection1VTKToWrite()
{

    if(m_DRTImage1MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
            m_TransformMetaInfo == NULL     ){
        return NULL;
        //TODO
    }

    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
    imageCalculatorFilter->SetImage(imageDRT1In);
    imageCalculatorFilter->Compute();

    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
    auto intWindowFilter = IntWindowType::New();
    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
    intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());

    intWindowFilter->SetOutputMinimum(0);
    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX){
        intWindowFilter->SetOutputMaximum(SHRT_MAX);
    }else
        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());


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

    toVTK2D1->SetInput(intWindowFilter->GetOutput());
    toVTK2D1->Update();

    return
            toVTK2D1->GetOutput();
}

vtkMatrix4x4 * itkImageProcessor::GetProjection1VTKTransform()
{

    m_Projection1VTKTransform->Identity();

    for(int iIdx = 0; iIdx<3 ; iIdx++){
        for(int jIdx = 0; jIdx<3 ; jIdx++){
            m_Projection1VTKTransform->SetElement(iIdx,jIdx,
                                                  interpolator1->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]);
        }
    }

    m_Projection1VTKTransform->SetElement(0,3,
                                          interpolator1->GetInverseTransform()->GetTranslation()[0]);

    m_Projection1VTKTransform->SetElement(1,3,
                                          interpolator1->GetInverseTransform()->GetTranslation()[1]);

    m_Projection1VTKTransform->SetElement(2,3,
                                          interpolator1->GetInverseTransform()->GetTranslation()[2]);

    return
            m_Projection1VTKTransform;
}

vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform()
{

    m_Projection2VTKTransform->Identity();

    for(int iIdx = 0; iIdx<3 ; iIdx++){
        for(int jIdx = 0; jIdx<3 ; jIdx++){
            m_Projection2VTKTransform->SetElement(iIdx,jIdx,
                                                  interpolator2->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]);
        }
    }

    m_Projection2VTKTransform->SetElement(0,3,
                                          interpolator2->GetInverseTransform()->GetTranslation()[0]);

    m_Projection2VTKTransform->SetElement(1,3,
                                          interpolator2->GetInverseTransform()->GetTranslation()[1]);

    m_Projection2VTKTransform->SetElement(2,3,
                                          interpolator2->GetInverseTransform()->GetTranslation()[2]);

    return
            m_Projection2VTKTransform;

}


vtkImageData* itkImageProcessor::GetProjection2VTK()
{
    if(m_DRTImage2MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
            m_TransformMetaInfo == NULL     ){
        return NULL;
        //TODO
    }


    using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

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

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

    return
            toVTK2D2->GetOutput();
}

vtkImageData* itkImageProcessor::GetProjection2VTKToWrite()
{
    if(m_DRTImage2MetaInfo == NULL ||
            m_DRTGeometryMetaInfo == NULL ||
            m_TransformMetaInfo == NULL     ){
        return NULL;
        //TODO
    }

    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
    imageCalculatorFilter->SetImage(imageDRT2In);
    imageCalculatorFilter->Compute();

    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
    auto intWindowFilter = IntWindowType::New();
    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
    intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());

    intWindowFilter->SetOutputMinimum(0);
    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX)
        intWindowFilter->SetOutputMaximum(SHRT_MAX);
    else
        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());


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

    toVTK2D2->SetInput(intWindowFilter->GetOutput());
    toVTK2D2->Update();

    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( imageDRT1In );

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

    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
    {
        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
    {
        writer2->Update();
    }
    catch (itk::ExceptionObject & err)
    {
        std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
        std::cerr << err << std::endl;
    }

}

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

    if(m_TransformMetaInfo == NULL){
        //todo
    }

    ImageType3D::PointType Translations;
    Translations[0]= dX;
    Translations[1]= dY;
    Translations[2]= dZ;

    m_TransformMetaInfo->SetUserTranslations(Translations);
}

void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ)
{
    ImageType3D::PointType Rotations;
    Rotations[0]= dX;
    Rotations[1]= dY;
    Rotations[2]= dZ;

    m_TransformMetaInfo->SetUserRotations(Rotations);

}

Optimizer::ParametersType
itkImageProcessor::GetUserIsocentricTransform(){

    Optimizer::ParametersType
            m_Pars(6);
    m_Pars.fill(0.);

    if(!m_TransformMetaInfo){
        return m_Pars;
    }

    m_Pars[0] = m_TransformMetaInfo->GetUserRotations()[0];
    m_Pars[1] = m_TransformMetaInfo->GetUserRotations()[1];
    m_Pars[2] = m_TransformMetaInfo->GetUserRotations()[2];

    m_Pars[3] = m_TransformMetaInfo->GetUserTranslations()[0];
    m_Pars[4] = m_TransformMetaInfo->GetUserTranslations()[1];
    m_Pars[5] = m_TransformMetaInfo->GetUserTranslations()[2];

    return
            m_Pars;

}



TransformType::Pointer
itkImageProcessor::GetCompleteIsocentricTransform(){

    if(m_TransformMetaInfo == nullptr ||
       m_CTMetaInfo == nullptr){
        return nullptr;
    }

    ImageType3D::PointType mISORTIEC =
    m_CTMetaInfo->GetLPS2IECDirections() * m_RTMetaInfo->GetIsocenterLPS();

    mISORTIEC[0] *= -1;
    mISORTIEC[1] *= -1;
    mISORTIEC[2] *= -1;

    TransformType::Pointer mTransf=
        MapTransformToNewOrigin(
                mISORTIEC,
                m_TransformMetaInfo->GetT(),
                m_TransformMetaInfo->GetR()
                );

    itk::TransformMetaInformation::PointType
            pTM =  m_CTMetaInfo->GetLPS2IECDirections() * m_CTMetaInfo->GetImportOffset();

    TransformType::OutputVectorType pTransl;

    pTransl[0] = mTransf->GetTranslation()[0] - pTM[0];
    pTransl[1] = mTransf->GetTranslation()[1] - pTM[1];
    pTransl[2] = mTransf->GetTranslation()[2] - pTM[2];

    mTransf->SetTranslation(pTransl);
    itk::Point<double,3>pZero;
    pZero.Fill(0.);
    mTransf->SetCenter(pZero);

    return mTransf;
}



void itkImageProcessor::SetRegionFixed1(
        int iIdx0, int iIdx1, int iSz0, int iSz1){

    auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion();

    auto index1 = region1temp.GetIndex();
    auto size1 = region1temp.GetSize();

    index1[0] = iIdx0;
    index1[1] = iIdx1;

    size1[0] = iSz0;
    size1[1] = iSz1;
    size1[2] = 1;

    roiAutoReg1.SetIndex(index1);
    roiAutoReg1.SetSize(size1);

    this->m_R23->SetroiAutoReg1(roiAutoReg1);

}

void itkImageProcessor::SetRegionFixed2(
        int iIdx0, int iIdx1, int iSz0, int iSz1){

    auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion();

    auto index2 = region2temp.GetIndex();
    auto size2 = region2temp.GetSize();

    index2[0] = iIdx0;
    index2[1] = iIdx1;

    size2[0] = iSz0;
    size2[1] = iSz1;
    size2[2] = 1;

    roiAutoReg2.SetIndex(index2);
    roiAutoReg2.SetSize(size2);

    this->m_R23->SetroiAutoReg2(roiAutoReg2);

}

void itkImageProcessor::ResetROIRegions(){

    this->m_R23->setROISizeToZero();

}

void itkImageProcessor::InizializeRegistration(){

    m_R23->SetPA(m_PASourceDupli->GetOutput());
    m_R23->SetLAT(m_LATSourceDupli->GetOutput());
    m_R23->SetVolume(m_VolumeSourceDupli->GetOutput());

    m_R23->Setr23Meta(m_r23MetaInfo);
    m_R23->SetPowellMeta(m_PowellOptimizerMetaInfo);
    m_R23->SetAmoebaMeta(m_AmoebaOptimizerMetaInfo);
    m_R23->SetMIMeta(m_MIMetricMetaInfo);
    m_R23->SetNCCMeta(m_NCCMetricMetaInfo);

    m_R23->SetInternalTransf1(m_InternalTransf1);
    m_R23->SetInternalTransf2(m_InternalTransf2);
    m_R23->Setfilter1(filter1);
    m_R23->Setfilter2(filter2);
    m_R23->Setinterpolator1(interpolator1);
    m_R23->Setinterpolator2(interpolator2);
    m_R23->SetTransformMetaInfo(m_TransformMetaInfo);

    m_R23->InitializeRegistration();

}



} // end namespace itk