#include "onf_stepcomputesloperaster.h"

#include "ct_itemdrawable/ct_image2d.h"
#include "ct_itemdrawable/tools/iterator/ct_groupiterator.h"
#include "ct_result/ct_resultgroup.h"
#include "ct_result/model/inModel/ct_inresultmodelgrouptocopy.h"
#include "ct_result/model/outModel/tools/ct_outresultmodelgrouptocopypossibilities.h"
#include "ct_view/ct_stepconfigurabledialog.h"

// Alias for indexing models
#define DEFin_res "res"
#define DEFin_grp "grp"
#define DEFin_DEM "dem"


// Constructor : initialization of parameters
ONF_StepComputeSlopeRaster::ONF_StepComputeSlopeRaster(CT_StepInitializeData &dataInit) : CT_AbstractStep(dataInit)
{
}

// Step description (tooltip of contextual menu)
QString ONF_StepComputeSlopeRaster::getStepDescription() const
{
    return tr("Compute slope raster");
}

// Step detailled description
QString ONF_StepComputeSlopeRaster::getStepDetailledDescription() const
{
    return tr("No detailled description for this step");
}

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

// Step copy method
CT_VirtualAbstractStep* ONF_StepComputeSlopeRaster::createNewInstance(CT_StepInitializeData &dataInit)
{
    return new ONF_StepComputeSlopeRaster(dataInit);
}

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

// Creation and affiliation of IN models
void ONF_StepComputeSlopeRaster::createInResultModelListProtected()
{
    CT_InResultModelGroupToCopy *resIn_res = createNewInResultModelForCopy(DEFin_res, tr("DEM"));
    resIn_res->setZeroOrMoreRootGroup();
    resIn_res->addGroupModel("", DEFin_grp, CT_AbstractItemGroup::staticGetType(), tr("Groupe"));
    resIn_res->addItemModel(DEFin_grp, DEFin_DEM, CT_Image2D<float>::staticGetType(), tr("DEM"));

}

// Creation and affiliation of OUT models
void ONF_StepComputeSlopeRaster::createOutResultModelListProtected()
{
    CT_OutResultModelGroupToCopyPossibilities *resCpy_res = createNewOutResultModelToCopy(DEFin_res);
    if (resCpy_res != NULL)
    {
        resCpy_res->addItemModel(DEFin_grp, _slopeRasterModelName, new CT_Image2D<float>(), tr("Slope"));
    }
}

// Semi-automatic creation of step parameters DialogBox
void ONF_StepComputeSlopeRaster::createPostConfigurationDialog()
{
//    CT_StepConfigurableDialog *configDialog = newStandardPostConfigurationDialog();

}

void ONF_StepComputeSlopeRaster::compute()
{
    QList<CT_ResultGroup*> outResultList = getOutResultList();
    CT_ResultGroup* res = outResultList.at(0);

    // COPIED results browsing
    CT_ResultGroupIterator itCpy_grp(res, this, DEFin_grp);
    while (itCpy_grp.hasNext() && !isStopped())
    {
        CT_StandardItemGroup* grp = (CT_StandardItemGroup*) itCpy_grp.next();
        
        CT_Image2D<float>* dem = (CT_Image2D<float>*)grp->firstItemByINModelName(this, DEFin_DEM);
        if (dem != NULL)
        {
            CT_Image2D<float>* slope = new CT_Image2D<float>(_slopeRasterModelName.completeName(), res, dem->minX(), dem->minY(), dem->colDim(), dem->linDim(), dem->resolution(), dem->level(), dem->NA(), dem->NA());
            grp->addItemDrawable(slope);

            for (size_t xx = 0 ; xx < slope->colDim() ; xx++)
            {
                for (size_t yy = 0 ; yy < slope->linDim() ; 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 (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 (a == dem->NA()) {a = e;}
                        if (b == dem->NA()) {b = e;}
                        if (c == dem->NA()) {c = e;}
                        if (d == dem->NA()) {d = e;}
                        if (f == dem->NA()) {f = e;}
                        if (g == dem->NA()) {g = e;}
                        if (h == dem->NA()) {h = e;}
                        if (i == dem->NA()) {i = e;}

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

                        slope->setValue(xx, yy, slope_degrees);
                    }
                }
            }

            slope->computeMinMax();
        }
    }
}