#include "onf_stepcomputesloperaster.h"

ONF_StepComputeSlopeRaster::ONF_StepComputeSlopeRaster() : SuperClass()
{
    _computePourc = true;
    _computeDegree = false;

    _nameDeg = tr("Pentes Degrés");
    _namePourc = tr("Pentes Pourcents");

    _maxSlope = 500.0;
}

QString ONF_StepComputeSlopeRaster::description() const
{
    return tr("Pente (paramétrable)");
}

QString ONF_StepComputeSlopeRaster::detailledDescription() const
{
    return tr("Calcul d'un raster de pente (en degrés)"
              "Formula from https://pro.arcgis.com/fr/pro-app/tool-reference/3d-analyst/how-slope-works.htm");
}

QString ONF_StepComputeSlopeRaster::inputDescription() const
{
    return SuperClass::inputDescription() + tr("<br><br>");
}

QString ONF_StepComputeSlopeRaster::outputDescription() const
{
    return SuperClass::outputDescription() + tr("<br><br>");
}

QString ONF_StepComputeSlopeRaster::detailsDescription() const
{
    return tr("");
}

QString ONF_StepComputeSlopeRaster::URL() const
{
    //return tr("STEP URL HERE");
    return SuperClass::URL(); //by default URL of the plugin
}

CT_VirtualAbstractStep* ONF_StepComputeSlopeRaster::createNewInstance() const
{
    return new ONF_StepComputeSlopeRaster();
}

//////////////////// PROTECTED METHODS //////////////////

void ONF_StepComputeSlopeRaster::declareInputModels(CT_StepInModelStructureManager& manager)
{
    manager.addResult(_inResult, tr("MNT"));
    manager.setZeroOrMoreRootGroup(_inResult, _inZeroOrMoreRootGroup);
    manager.addGroup(_inZeroOrMoreRootGroup, _inGroup);
    manager.addItem(_inGroup, _inDEM, tr("MNT"));
}

void ONF_StepComputeSlopeRaster::fillPreInputConfigurationDialog(CT_StepConfigurableDialog* preInputConfigDialog)
{
    preInputConfigDialog->addBool("", tr("Calculer les pentes en pourcents"),"", _computePourc);
    preInputConfigDialog->addBool("", tr("Calculer les pentes en degrés"),"", _computeDegree);
    preInputConfigDialog->addDouble(tr("Pente maximale autorisée"), "%", 0.0, 9999.99, 2, _maxSlope);
}

void ONF_StepComputeSlopeRaster::fillPostInputConfigurationDialog(CT_StepConfigurableDialog* postInputConfigDialog)
{
    if (_computePourc)  {postInputConfigDialog->addString(tr("Nom à donner au raster (pentes en pourcents)"), "", _namePourc, tr("<strong>Attention</strong> : ce paramètre ne peut être pris en compte qu'à l'ajout de l'étape ou au chargement d'un script.<br>Une modification de ce paramètre lors d'une reconfiguration des paramètres se répercutera dans les scripts exportés ultérieurement, <strong>mais pas dans la session en cours</strong>. "));}
    if (_computeDegree) {postInputConfigDialog->addString(tr("Nom à donner au raster (pentes en degrés)"), "", _nameDeg, tr("<strong>Attention</strong> : ce paramètre ne peut être pris en compte qu'à l'ajout de l'étape ou au chargement d'un script.<br>Une modification de ce paramètre lors d'une reconfiguration des paramètres se répercutera dans les scripts exportés ultérieurement, <strong>mais pas dans la session en cours</strong>. "));}
}

void ONF_StepComputeSlopeRaster::declareOutputModels(CT_StepOutModelStructureManager& manager)
{
    manager.addResultCopy(_inResult);
    if (_computePourc) {manager.addItem(_inGroup, _outRasterPourc, _namePourc);}
    if (_computeDegree) {manager.addItem(_inGroup, _outRaster, _nameDeg);}
}

void ONF_StepComputeSlopeRaster::compute()
{
    double maxSlopeRad= std::atan(_maxSlope/100.0);

    for (CT_StandardItemGroup* grp : _inGroup.iterateOutputs(_inResult))
    {
        for (const CT_Image2D<float>* dem : grp->singularItems(_inDEM))
        {
            if (isStopped()) {return;}

            CT_Image2D<float>* slope = nullptr;
            CT_Image2D<float>* slopePourc = nullptr;

            if (_computePourc)
            {
                slopePourc = new CT_Image2D<float>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->level(), dem->NA(), dem->NA());
                grp->addSingularItem(_outRasterPourc, slopePourc);
            }

            if (_computeDegree)
            {
                slope = new CT_Image2D<float>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->level(), dem->NA(), dem->NA());
                grp->addSingularItem(_outRaster, slope);
            }

            for (int xx = 0 ; xx < dem->xdim() ; xx++)
            {
                for (int yy = 0 ; yy < dem->ydim() ; yy++)
                {
                    // formula from https://pro.arcgis.com/fr/pro-app/tool-reference/3d-analyst/how-slope-works.htm
                    //
                    // Original raster
                    //
                    //  a   b   c
                    //  d   e   f
                    //  g   h   i
                    //
                    // slope_radians = ATAN ( sqrt([dz/dx]2 + [dz/dy]2) )
                    // with:
                    // [dz/dx] = ((c + 2f + i) - (a + 2d + g)) / (8 * x_cellsize)
                    // [dz/dy] = ((g + 2h + i) - (a + 2b + c)) / (8 * y_cellsize)

                    float e = dem->value(xx, yy);

                    if (!qFuzzyCompare(e, dem->NA()))
                    {
                        float a = dem->value(xx - 1, yy - 1);
                        float b = dem->value(xx    , yy - 1);
                        float c = dem->value(xx + 1, yy - 1);
                        float d = dem->value(xx - 1, yy    );
                        float f = dem->value(xx + 1, yy    );
                        float g = dem->value(xx - 1, yy + 1);
                        float h = dem->value(xx    , yy + 1);
                        float i = dem->value(xx + 1, yy + 1);

                        if (qFuzzyCompare(a, dem->NA())) {a = e;}
                        if (qFuzzyCompare(b, dem->NA())) {b = e;}
                        if (qFuzzyCompare(c, dem->NA())) {c = e;}
                        if (qFuzzyCompare(d, dem->NA())) {d = e;}
                        if (qFuzzyCompare(f, dem->NA())) {f = e;}
                        if (qFuzzyCompare(g, dem->NA())) {g = e;}
                        if (qFuzzyCompare(h, dem->NA())) {h = e;}
                        if (qFuzzyCompare(i, dem->NA())) {i = e;}

                        float dzdx = ((c + 2.0f*f + i) - (a + 2.0f*d + g)) / (8.0f * float(dem->resolution()));
                        float dzdy = ((g + 2.0f*h + i) - (a + 2.0f*b + c)) / (8.0f * float(dem->resolution()));
                        double slope_radians = std::atan(std::sqrt(double(dzdx)*double(dzdx) + double(dzdy)*double(dzdy)));

                        if (slope_radians > maxSlopeRad) {slope_radians = maxSlopeRad;}

                        double slope_degrees = slope_radians * 180.0 / M_PI;
                        double slope_Pourcent = 100.0 * std::tan(slope_radians);

                        if (_computePourc) {slopePourc->setValue(xx, yy, float(slope_Pourcent));}
                        if (_computeDegree) {slope->setValue(xx, yy, float(slope_degrees));}
                    }
                }
            }

            if (_computePourc) {slopePourc->computeMinMax();}
            if (_computeDegree) {slope->computeMinMax();}
        }
    }
}