/*****************************************************************************
 * $CAMITK_LICENCE_BEGIN$
 *
 * CamiTK - Computer Assisted Medical Intervention ToolKit
 * (c) 2001-2025 Univ. Grenoble Alpes, CNRS, Grenoble INP - UGA, TIMC, 38000 Grenoble, France
 *
 * Visit http://camitk.imag.fr for more information
 *
 * This file is part of CamiTK.
 *
 * CamiTK is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License version 3
 * only, as published by the Free Software Foundation.
 *
 * CamiTK is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License version 3 for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * version 3 along with CamiTK.  If not, see <http://www.gnu.org/licenses/>.
 *
 * $CAMITK_LICENCE_END$
 ****************************************************************************/
// CamiTK includes
#include "DicomComponent.h"

#include <ImageOrientationHelper.h>
#include <FrameOfReference.h>
#include <Transformation.h>
#include <TransformationManager.h>
#include <Log.h>
#include "DicomParser.h"
#include "DicomSeries.h"

// Qt includes
#include <QDir>

// VTK includes
#include <vtkStringArray.h>
#include <vtkPointData.h>
#include <vtkImageChangeInformation.h>

// GDCM includes
#include <vtkGDCMImageReader.h>
#include <gdcmReader.h>
#include <gdcmScanner.h>
#include <gdcmIPPSorter.h>

using namespace camitk;

// --------------- Constructor -------------------
DicomComponent::DicomComponent(DicomSeries* dicomSeries) : ImageComponent("", "Dicom Image") {

    // associated the component's dicom series
    series = dicomSeries;

    std::vector<std::string> stdFileNames = dicomSeries->getStdFileNames();
    // scan files for series description
    gdcm::Scanner scanner;
    gdcm::Tag seriesDescriptionTag = gdcm::Tag(0x0008, 0x103e);
    scanner.AddTag(seriesDescriptionTag);
    scanner.Scan(stdFileNames);
    QString imageName = QString(scanner.GetValue(stdFileNames.at(0).c_str(), seriesDescriptionTag));
    setName(imageName);
    frameOfReference->setName(imageName + " (main)");
    TransformationManager::getFrameOfReferenceOwnership(getDataFrame())->setName(imageName + " (data)");

    // Use Image Position Patient filter (IPPSorter) to correctly sort slices according to their Z spacing
    // Also deduce the actual Z spacing
    std::vector< std::string > files;
    double zSpacing = 0.0;
    if (stdFileNames.size() > 1) {
        gdcm::IPPSorter ippSorter;
        ippSorter.SetComputeZSpacing(true);
        ippSorter.SetZSpacingTolerance(0.001);
        if (!ippSorter.Sort(stdFileNames)) {
            CAMITK_ERROR(tr("IPPSorter sorting failed. Try to adjust Z spacing tolerance."))
            files = stdFileNames;
            zSpacing = DicomParser::getZSpacing(files);
        }
        else {
            files = ippSorter.GetFilenames();
            zSpacing = ippSorter.GetZSpacing();
        }
    }
    else {
        files = stdFileNames;
        zSpacing = DicomParser::getZSpacing(files);
    }


    // convert this list as a vtkStringArray
    vtkSmartPointer<vtkStringArray> fileNamesSorted = vtkSmartPointer<vtkStringArray>::New();
    for (std::string file : files) {
        fileNamesSorted->InsertNextValue(file.c_str());
    }

    // get image orientation information
    // we need to get the rotation matrix from tag "Direct cos angle"
    // ImageOrientationHelper::PossibleImageOrientations orientationFromCos = readDirectCosinesAngle(files);


    // get the image position orientation
    // TODO read the tag (0018, 5100) Patient position
    // https://dicom.nema.org/medical/Dicom/2016e/output/chtml/part03/sect_C.7.3.html#sect_C.7.3.1.1.2

    // create image data corresponding to the component
    imageReader = vtkSmartPointer<vtkGDCMImageReader>::New();
    if (fileNamesSorted->GetSize() == 1) {
        imageReader->SetFileName(fileNamesSorted->GetValue(0).c_str());
    }
    else {
        imageReader->SetFileNames(fileNamesSorted);
    }
    imageReader->Update();
    vtkSmartPointer<vtkImageData> rawData = imageReader->GetOutput();

    /// TODO Store Image position patient as a CamiTK property?

    // Update Z-spacing
    // vtkGDCMImageReader misses this information, see: http://gdcm.sourceforge.net/2.4/html/classvtkGDCMImageReader.html#details
    // Use the value found from IPPSorter or DicomParser::GetZSpacing()
    double* spacing = imageReader->GetDataSpacing();

    // Update Z Spacing using a VTK pipeline
    vtkSmartPointer<vtkImageChangeInformation> imageInfoChanger = vtkSmartPointer<vtkImageChangeInformation>::New();
    imageInfoChanger->SetInputData(rawData); //translatedData);
    double origin[3] = {0.0, 0.0, 0.0}; //We set the image origin at (0,0,0)
    imageInfoChanger->SetOutputOrigin(origin);
    imageInfoChanger->SetOutputSpacing(spacing[0], spacing[1], zSpacing);
    imageInfoChanger->Update();

    // Flip all actors in the Y axis
    // Explanation:
    // DICOM stores the upper left pixel as the first pixel in an image.
    // However, VTK stores the lower left pixel as the first pixel in an image
    // As we based our image frame on RAI coordinate (DICOM LPS located at the bottom left hand corner)
    // with a Radiologist point of view (Axial, Coronal, Sagittal)
    // To ensure that VTK frame is in RAI coordinates, we flip the image in the Y axis
    vtkSmartPointer<vtkImageFlip> flipYFilter = vtkSmartPointer<vtkImageFlip>::New();
    flipYFilter->SetFilteredAxis(1); // flip y axis
    flipYFilter->SetInputConnection(imageInfoChanger->GetOutputPort());
    flipYFilter->Update();

    // Get DICOM patient position and orientation and set it as this component's frame transform
    double* position = imageReader->GetImagePositionPatient();
    double* orientation = imageReader->GetImageOrientationPatient();

    double* xdir = &orientation[0];
    double* ydir = &orientation[3];
    double  zdir[3];
    vtkMath::Cross(xdir, ydir, zdir);
    vtkSmartPointer<vtkMatrix4x4> affineMatrix = vtkSmartPointer<vtkMatrix4x4>::New();

    for (int i = 0; i < 3; i++) {
        affineMatrix->Element[i][0] = xdir[i];
        affineMatrix->Element[i][1] = ydir[i];
        affineMatrix->Element[i][2] = zdir[i];
        affineMatrix->Element[i][3] = position[i];
    }

    // get the processed image (flipped and with correct z-spacing)
    vtkSmartPointer<vtkImageData> image = flipYFilter->GetOutput();
    setImageData(image, false);

    //apply the affine matrix to the image
    updateMainTransformation(affineMatrix);

    // Wait for the LUT update in CamiTK and / or support for color image
    if (getLut()) { // sometimes ImageComponent happens not to have a lut, strange behaviour
        updateLUT();
    }
    else {
        CAMITK_WARNING(tr("Image LUT is null."))
    }

}

// --------------- destructor -------------------
DicomComponent::~DicomComponent() {
    if (series != nullptr) {
        delete series;
    }
}

// --------------- updateLUT -------------------
void DicomComponent::updateLUT() {
    // Update LUT
    // Initialize our lut with vtkGDCMImageReader information found, as our LUT needs repair ...
    double range[2] = {0.0, 0.0};
    imageReader->GetOutput()->GetScalarRange(range);
    getLut()->SetRange(range);  // we need to set up range and table values
    getLut()->SetLevel((range[1] + range[0]) / 2);
    getLut()->SetNumberOfTableValues(abs(range[0]) + abs(range[1]));
}

// --------------- readDirectCosinesAngle -------------------
camitk::ImageOrientationHelper::PossibleImageOrientations DicomComponent::readDirectCosinesAngle(const std::vector< std::string >& fileNames) const {

    // scan files Image Orientation Patient
    gdcm::Scanner scanner;
    gdcm::Tag iopTag = gdcm::Tag(0x0020, 0x0037);
    scanner.AddTag(iopTag);
    // Scan should never failed, since DicomParser::parseDirectory() has already filter files.
    CAMITK_ERROR_IF((!scanner.Scan(fileNames)), tr("Scan failed looking for tag (0x0020, 0x0037) Image Orientation Patient."))

    // Check value tag exists
    gdcm::Scanner::TagToValue const& ttv = scanner.GetMapping(fileNames[0].c_str());
    gdcm::Scanner::TagToValue::const_iterator it = ttv.find(iopTag);
    if (!(it != ttv.end())) {
        CAMITK_WARNING(tr("No tag (0x0020, 0x0037) Image Orientation Patient found on image \"%1\"").arg(QString::fromStdString(fileNames[0])))
        return ImageOrientationHelper::UNKNOWN;
    }

    // Then we know it exists
    std::string value = scanner.GetValue(fileNames[0].c_str(), iopTag);

    // convert the string into the appropriate couple of cosine vectors
    double x[3] = {0.0, 0.0, 0.0};
    double y[3] = {0.0, 0.0, 0.0};
    std::sscanf(value.c_str(), R"(%lf\%lf\%lf\%lf\%lf\%lf)", &x[0], &x[1], &x[2], &y[0], &y[1], &y[2]);

    // get the 90 degrees closest cosines for each vector
    x[0] = roundCosine(x[0]);
    x[1] = roundCosine(x[1]);
    x[2] = roundCosine(x[2]);
    y[0] = roundCosine(y[0]);
    y[1] = roundCosine(y[1]);
    y[2] = roundCosine(y[2]);

    // return the appropriate ImageOrientation corresponding to this cosines
    // we only need to check one component of vector X and Y as they are colinear (other components are 0) to one of the frame vectors

    if ((x[0] == 1.0) && (y[1] == 1.0)) { // identity matrix => Image frame is the scanner frame
        return ImageOrientationHelper::RAI;
    }

    // other cases Image frame is different
    if (x[0] == 1.0) {
        if (y[1] == -1.0) {
            return ImageOrientationHelper::RPS;
        }
        if (y[2] == 1.0) {
            return ImageOrientationHelper::RIP;
        }
        if (y[2] == -1.0) {
            return ImageOrientationHelper::RSA;
        }
    }
    if (x[0] == -1.0) {
        if (y[1] == 1.0) {
            return ImageOrientationHelper::LAS;
        }
        if (y[2] == -1.0) {
            return ImageOrientationHelper::LPI;
        }
        if (y[2] == 1.0) {
            return ImageOrientationHelper::LIA;
        }
    }
    if (x[1] == 1.0) {
        if (y[0] == 1.0) {
            return ImageOrientationHelper::ARS;
        }
        if (y[0] == -1.0) {
            return ImageOrientationHelper::ALI;
        }
        if (y[2] == 1.0) {
            return ImageOrientationHelper::AIR;
        }
        if (y[2] == -1.0) {
            return ImageOrientationHelper::ASL;
        }
    }
    if (x[1] == -1.0) {
        if (y[0] == 1.0) {
            return ImageOrientationHelper::PRI;
        }
        if (y[0] == -1.0) {
            return ImageOrientationHelper::PLS;
        }
        if (y[2] == 1.0) {
            return ImageOrientationHelper::PIL;
        }
        if (y[2] == -1.0) {
            return ImageOrientationHelper::PSR;
        }
    }
    if (x[2] == 1.0) {
        if (y[0] == 1.0) {
            return ImageOrientationHelper::IRA;
        }
        if (y[0] == -1.0) {
            return ImageOrientationHelper::ILP;
        }
        if (y[1] == 1.0) {
            return ImageOrientationHelper::IAL;
        }
        if (y[1] == -1.0) {
            return ImageOrientationHelper::IPR;
        }
    }
    if (x[2] == -1.0) {
        if (y[0] == 1.0) {
            return ImageOrientationHelper::SRP;
        }
        if (y[0] == -1.0) {
            return ImageOrientationHelper::SLA;
        }
        if (y[1] == 1.0) {
            return ImageOrientationHelper::SAR;
        }
        if (y[1] == -1.0) {
            return ImageOrientationHelper::SPL;
        }
    }

    // should never return UNKNOW
    CAMITK_WARNING(tr("No orientation found for this image (direct cosines)."))
    return ImageOrientationHelper::UNKNOWN;

}

// --------------- roundCosine -------------------
double DicomComponent::roundCosine(const double& value) const {
    if (value < -0.5) {
        return -1.0;
    }
    if ((value >= -0.5) && (value <= 0.5)) {
        return 0.0;
    }
    else {
        return 1.0;
    }
}