/****************************************************************************
 Copyright (C) 2010-2012 the Office National des Forêts (ONF), France
                         All rights reserved.

 Contact : alexandre.piboule@onf.fr

 Developers : Alexandre PIBOULE (ONF)

 This file is part of PluginONF library.

 PluginONF is free library: you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 PluginONF 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 General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with PluginONF.  If not, see <http://www.gnu.org/licenses/lgpl.html>.
*****************************************************************************/

#include "onf_stepcomputepotentialstartarea.h"

#include "ct_log/ct_logmanager.h"

ONF_StepComputePotentialStartArea::ONF_StepComputePotentialStartArea() : SuperClass()
{
    _altitudeMin = 1300;
    _slopeMin = 30;
    _slopeMax = 55;
    _distanceMax = 600;
}

QString ONF_StepComputePotentialStartArea::description() const
{
    return tr("Délimiter les aires de Départ Potentielles");
}

QString ONF_StepComputePotentialStartArea::detailledDescription() const
{
    return tr("Cette étape permet de délimiter les Zones de Départ Potentielles.");
}

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

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

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

CT_VirtualAbstractStep* ONF_StepComputePotentialStartArea::createNewInstance() const
{
    // cree une copie de cette etape
    return new ONF_StepComputePotentialStartArea();
}

/////////////////////// PROTECTED ///////////////////////

void ONF_StepComputePotentialStartArea::declareInputModels(CT_StepInModelStructureManager& manager)
{
    manager.addResult(_inResultDTM, tr("MNT"), "", true);
    manager.setZeroOrMoreRootGroup(_inResultDTM, _inZeroOrMoreRootGroupDTM);
    manager.addGroup(_inZeroOrMoreRootGroupDTM, _inGroupDTM);
    manager.addItem(_inGroupDTM, _inDTM, tr("MNT"));
    manager.addItem(_inGroupDTM, _inGeomorphon, tr("Géomorphons"));

    manager.addResult(_inResultProtectionForests, tr("Forêt de protection (optionnel)"), "", true);
    manager.setZeroOrMoreRootGroup(_inResultProtectionForests, _inZeroOrMoreRootGroupProtectionForests);
    manager.addGroup(_inZeroOrMoreRootGroupProtectionForests, _inGroupProtectionForests);
    manager.addItem(_inGroupProtectionForests, _inProtectionForests, tr("Forêt de protection"));
}

void ONF_StepComputePotentialStartArea::declareOutputModels(CT_StepOutModelStructureManager& manager)
{
    manager.addResultCopy(_inResultDTM);
    manager.addItem(_inGroupDTM, _outMaskCloseToRidge, tr("Proche Crêtes"));
    manager.addItem(_inGroupDTM, _outMaskAreaOfInterest, tr("Zones d'intérêt brutes"));
    manager.addItem(_inGroupDTM, _outMaskAreaOfInterestForest, tr("Zones d'intérêt"));
}

void ONF_StepComputePotentialStartArea::fillPostInputConfigurationDialog(CT_StepConfigurableDialog* postInputConfigDialog)
{
    postInputConfigDialog->addDouble(tr("Altitude minimale"), "m", 0, 99999, 2, _altitudeMin);
    postInputConfigDialog->addDouble(tr("Pente minimale"), "°", 0, 90, 2, _slopeMin);
    postInputConfigDialog->addDouble(tr("Pente maximale"), "°", 0, 90, 2, _slopeMax);
    postInputConfigDialog->addDouble(tr("Distance maximale aux crètes"), "m", 0, 9999, 2, _distanceMax);
}


void ONF_StepComputePotentialStartArea::compute()
{
    // get Protection forest raster
    const CT_Image2D<quint8>* pforest = nullptr;

    for (const CT_StandardItemGroup* group : _inGroupProtectionForests.iterateInputs(_inResultProtectionForests))
    {
        pforest = group->singularItem(_inProtectionForests);
    }
    if (pforest == nullptr)
    {
        PS_LOG->addMessage(LogInterface::error, LogInterface::step, QString(tr("Pas de raster forêt de protecion fourni")));
    }

    setProgress(float(10.0));

    // compute area of interest mask
    for (CT_StandardItemGroup* group : _inGroupDTM.iterateOutputs(_inResultDTM))
    {
        const CT_Image2D<float>* mnt = group->singularItem(_inDTM);
        const CT_Image2D<qint8>* geomorphon = group->singularItem(_inGeomorphon);

        if (mnt == nullptr) {return;}
        if (geomorphon == nullptr) {return;}

        if (isStopped()) {return;}

        CT_Image2D<quint8>* maskAOI = new CT_Image2D<quint8>(mnt->minX(), mnt->minY(), mnt->xdim(), mnt->ydim(), mnt->resolution(), mnt->minZ(), 255, 0);
        CT_Image2D<quint8>* maskAOIForest = new CT_Image2D<quint8>(mnt->minX(), mnt->minY(), mnt->xdim(), mnt->ydim(), mnt->resolution(), mnt->minZ(), 255, 0);
        CT_Image2D<quint8>* maskCloseToRidge = new CT_Image2D<quint8>(mnt->minX(), mnt->minY(), mnt->xdim(), mnt->ydim(), mnt->resolution(), mnt->minZ(), 255, 0);

        computeCloseRoRidgeMask(geomorphon, maskCloseToRidge);

        size_t ncells = maskAOI->nCells();
        for (size_t index = 0 ; index < ncells ; index++)
        {
            float mntVal = mnt->valueAtIndex(index);

            int xx, yy;
            mnt->indexToGrid(index, xx, yy);

            Eigen::Vector3d center;
            mnt->getCellCenterCoordinates(index, center);

            if (mntVal >= _altitudeMin)
            {
                double slopeVal = computeSlope(xx, yy, mnt);

                if (slopeVal >= _slopeMin && slopeVal <= _slopeMax)
                {
                    if (maskCloseToRidge->value(xx, yy) == 1)
                    {
                        maskAOI->setValueAtIndex(index, 1);

                        int valForest = 0;
                        if (pforest != nullptr)
                        {
                            valForest = pforest->valueAtCoords(center(0), center(1));
                        }

                        if (valForest != 1)
                        {
                            maskAOIForest->setValueAtIndex(index, 1);
                        }
                    }
                }
            }

            setProgress(float(10.0 + 60.0*index/ncells));
        }

        // ajout du raster Masque
        group->addSingularItem(_outMaskCloseToRidge, maskCloseToRidge);
        maskCloseToRidge->computeMinMax();

        group->addSingularItem(_outMaskAreaOfInterest, maskAOI);
        maskAOI->computeMinMax();

        group->addSingularItem(_outMaskAreaOfInterestForest, maskAOIForest);
        maskAOIForest->computeMinMax();

        setProgress(100.0f);
    }
}


double ONF_StepComputePotentialStartArea::computeSlope(int xx, int yy, const CT_Image2D<float>* mnt)
{
    // formula from https://pro.arcgis.com/fr/pro-app/tool-reference/3d-analyst/how-slope-works.htm

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

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

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

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

        return slope_radians * 180.0 / M_PI;
    }

    return 0;
}


void ONF_StepComputePotentialStartArea::computeCloseRoRidgeMask(const CT_Image2D<qint8>* geomorphon, CT_Image2D<quint8>* mask)
{
    int distPixel = std::ceil(_distanceMax / geomorphon->resolution());
    double distanceMax2 = _distanceMax * _distanceMax;

    for (int xx = 0 ; xx < geomorphon->xdim() ; xx++)
    {
        for (int yy = 0 ; yy < geomorphon->ydim() ; yy++)
        {
            if (geomorphon->value(xx, yy) == 3)
            {
                int xnmin = xx - distPixel; if (xnmin < 0) {xnmin = 0;}
                int ynmin = yy - distPixel; if (ynmin < 0) {ynmin = 0;}

                int xnmax = xx + distPixel; if (xnmax >= geomorphon->xdim()) {xnmax = geomorphon->xdim() - 1;}
                int ynmax = yy + distPixel; if (ynmax >= geomorphon->ydim()) {ynmax = geomorphon->ydim() - 1;}

                Eigen::Vector3d center, centerN;
                geomorphon->getCellCenterCoordinates(xx, yy, center);

                for (int xn = xnmin ; xn <= xnmax ; xn++)
                {
                    for (int yn = ynmin ; yn <= ynmax ; yn++)
                    {
                        if (mask->value(xn, yn) < 1)
                        {
                            geomorphon->getCellCenterCoordinates(xn, yn, centerN);

                            double distance2 = pow(center(0) - centerN(0), 2) + pow(center(1) - centerN(1), 2);

                            if (distance2 < distanceMax2)
                            {
                                mask->setValue(xn, yn, 1);
                            }
                        }
                    }
                }
            }
        }
    }
}