#include "onf_stepcomputehillshaderaster.h"

ONF_StepComputeHillShadeRaster::ONF_StepComputeHillShadeRaster() : SuperClass()
{
    _azimut = 315.0;
    _altitude = 45.0;
}

QString ONF_StepComputeHillShadeRaster::description() const
{
    return tr("Créer un raster ombré");
}

QString ONF_StepComputeHillShadeRaster::detailledDescription() const
{
    return tr("Formula from https://pro.arcgis.com/fr/pro-app/tool-reference/3d-analyst/how-hillshade-works.htm");
}

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

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

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

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

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

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

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

void ONF_StepComputeHillShadeRaster::declareOutputModels(CT_StepOutModelStructureManager& manager)
{
    manager.addResultCopy(_inResult);
    manager.addItem(_inGroup, _outHillShade, tr("Ombrage"));
}

void ONF_StepComputeHillShadeRaster::fillPostInputConfigurationDialog(CT_StepConfigurableDialog* postInputConfigDialog)
{
    postInputConfigDialog->addDouble(tr("Azimut"), "°", 0, 360, 2, _azimut);
    postInputConfigDialog->addDouble(tr("Altitude"), "°", 0, 90, 2, _altitude);
}

void ONF_StepComputeHillShadeRaster::compute()
{
    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;

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

            CT_Image2D<uchar>* hillshade = new CT_Image2D<uchar>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->level(), 0, 0);
            grp->addSingularItem(_outHillShade, hillshade);

            for (int xx = 0 ; xx < hillshade->xdim() ; xx++)
            {
                for (int yy = 0 ; yy < hillshade->ydim() ; 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 (!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(hillshade->resolution()));
                        float dzdy = ((g + 2.0f*h + i) - (a + 2.0f*b + c)) / (8.0f * float(hillshade->resolution()));
                        double slope_radians = std::atan(std::sqrt(double(dzdx)*double(dzdx) + double(dzdy)*double(dzdy)));

                        double aspect_rad = 0;
                        if (!qFuzzyCompare(dzdx, 0))
                        {
                            aspect_rad = std::atan2(double(dzdy), double(- 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, static_cast<uchar>(hillshadeVal));
                    }
                }
            }

            hillshade->computeMinMax();
        }
    }
}