#include "onf_stepcomputehillshaderaster.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_StepComputeHillShadeRaster::ONF_StepComputeHillShadeRaster(CT_StepInitializeData &dataInit) : CT_AbstractStep(dataInit)
{
    _azimut = 315.0;
    _altitude = 45.0;

}

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

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

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

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

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

// Creation and affiliation of IN models
void ONF_StepComputeHillShadeRaster::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_StepComputeHillShadeRaster::createOutResultModelListProtected()
{
    CT_OutResultModelGroupToCopyPossibilities *resCpy_res = createNewOutResultModelToCopy(DEFin_res);
    if (resCpy_res != NULL)
    {
        resCpy_res->addItemModel(DEFin_grp, _hillshadeRasterModelName, new CT_Image2D<float>(), tr("HillShade"));
    }
}

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

    configDialog->addDouble(tr("Azimut"), "°", 0, 360, 2, _azimut);
    configDialog->addDouble(tr("Altitude"), "°", 0, 90, 2, _altitude);
}

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


   double zenith_rad = (90.0 - _altitude) * M_PI / 180.0;

   double azimuth_math = 360.0 - _azimut + 90.0;
   if (azimuth_math >= 360.0) {azimuth_math -= 360.0;}
   double azimuth_rad = azimuth_math * M_PI / 180.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<uchar>* hillshade = new CT_Image2D<uchar>(_hillshadeRasterModelName.completeName(), res, dem->minX(), dem->minY(), dem->colDim(), dem->linDim(), dem->resolution(), dem->level(), 0, 0);
            grp->addItemDrawable(hillshade);

            for (size_t xx = 0 ; xx < hillshade->colDim() ; xx++)
            {
                for (size_t yy = 0 ; yy < hillshade->linDim() ; yy++)
                {
                    // formula from https://pro.arcgis.com/fr/pro-app/tool-reference/3d-analyst/how-hillshade-works.htm

                    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 * hillshade->resolution());
                        double dzdy = ((g + 2.0*h + i) - (a + 2.0*b + c)) / (8.0 * hillshade->resolution());
                        double slope_radians = std::atan(std::sqrt(dzdx*dzdx + dzdy*dzdy));

                        double aspect_rad = 0;
                        if (dzdx != 0)
                        {
                            aspect_rad = std::atan2(dzdy, - dzdx);
                              if (aspect_rad < 0)
                              {
                                aspect_rad = 2.0 * M_PI + aspect_rad;
                              }
                        } else {
                            if (dzdy > 0)
                            {
                              aspect_rad = M_PI / 2.0;
                            }
                            else if (dzdy < 0)
                            {
                              aspect_rad = 2.0 * M_PI - M_PI / 2.0;
                            } else
                            {
                              aspect_rad = 0;
                            }
                        }

                        double hillshadeVal = 255.0 * ((cos(zenith_rad) * cos(slope_radians)) + (sin(zenith_rad) * sin(slope_radians) * cos(azimuth_rad - aspect_rad)));

                        hillshade->setValue(xx, yy, (uchar)hillshadeVal);
                    }
                }
            }

            hillshade->computeMinMax();
        }
    }
}