/****************************************************************************
 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_stepcomputeprotectionforestmask.h"

#include "ct_log/ct_logmanager.h"

ONF_StepComputeProtectionForestMask::ONF_StepComputeProtectionForestMask() : SuperClass()
{
    _neededFields.append(CT_TextFileConfigurationFields("Class", QRegExp("([cC][lL]|[aA][sS][sS]|[cC][lL]|[aA][sS][sS][eE])"), false));
    _neededFields.append(CT_TextFileConfigurationFields("Snow Height (cm)", QRegExp("[hH][nN]|[sS][hH]|[hH][aA][uU][tT][eE][uU][rR]|[hH][eE][iI][gG][tT][hH]"), false));

    _refFileName = "";
    _refHeader = true;
    _refSeparator = "\t";
    _refDecimal = ".";
    _refLocale = QLocale(QLocale::English, QLocale::UnitedKingdom).name();
    _refSkip = 0;

    _HnZRef = 3.2;
    _ZRef = 2000.0;
    _gradientInfZ = 0.18;
    _gradientSupZ = 0.24;

    _threshold = 0.7;

    _altitudeMin = 1300;
    _altitudeMax = 2300;

    _HsvDefault = 1.5;

}

QString ONF_StepComputeProtectionForestMask::description() const
{
    return tr("Créer masque de forêt de protection");
}

QString ONF_StepComputeProtectionForestMask::detailledDescription() const
{
    return tr("Cette étape permet de créer un masque des forêts de protection, à partir d'un raster type de forêt et d'un fichier de hauteurs de neiges par type.");
}

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

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

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

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

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

void ONF_StepComputeProtectionForestMask::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.addResult(_inResultCHM, tr("MNH"), "", true);
    manager.setZeroOrMoreRootGroup(_inResultCHM, _inZeroOrMoreRootGroupCHM);
    manager.addGroup(_inZeroOrMoreRootGroupCHM, _inGroupCHM);
    manager.addItem(_inGroupCHM, _inCHM, tr("MNH"));

    manager.addResult(_inResultForest, tr("Types de forêts (optionnel)"), "", true);
    manager.setZeroOrMoreRootGroup(_inResultForest, _inZeroOrMoreRootGroupForest);
    manager.addGroup(_inZeroOrMoreRootGroupForest, _inGroupForest);
    manager.addItem(_inGroupForest, _inForest, tr("Types de forêts"));
}

void ONF_StepComputeProtectionForestMask::declareOutputModels(CT_StepOutModelStructureManager& manager)
{    
    manager.addResultCopy(_inResultCHM);
    manager.addItem(_inGroupCHM, _outMask, tr("Forêts protection brut"));
    manager.addItem(_inGroupCHM, _outMask_kernel, tr("Forêts protection kernel"));
}

void ONF_StepComputeProtectionForestMask::fillPostInputConfigurationDialog(CT_StepConfigurableDialog* postInputConfigDialog)
{
    postInputConfigDialog->addDouble(tr("Hauteur de neige de référence"), "cm", 0, 9999, 2, _HnZRef, 100);
    postInputConfigDialog->addDouble(tr("Altitude de référence"), "m", 0, 9999, 2, _ZRef);
    postInputConfigDialog->addDouble(tr("Correction si < altitude de référence"), "cm / 100 m", 0, 9999, 2, _gradientInfZ, 100);
    postInputConfigDialog->addDouble(tr("Correction si > altitude de référence"), "cm / 100 m", 0, 9999, 2, _gradientSupZ, 100);

    postInputConfigDialog->addEmpty();
    postInputConfigDialog->addDouble(tr("Seuil"), "%", 0, 100, 0, _threshold, 100);
    postInputConfigDialog->addEmpty();

    postInputConfigDialog->addDouble(tr("Altitude minimale"), "m", 0, 99999, 2, _altitudeMin);
    postInputConfigDialog->addDouble(tr("Altitude maximale"), "m", 0, 99999, 2, _altitudeMax);

    postInputConfigDialog->addEmpty();

    if (_inForest.hasAtLeastOnePossibilitySelected(_inResultForest))
    {
        postInputConfigDialog->addAsciiFileChoice("Fichier de référence", "Fichier ASCII (*.txt ; *.asc)", true, _neededFields, _refFileName, _refHeader, _refSeparator, _refDecimal, _refLocale, _refSkip, _refColumns);
    } else {
        postInputConfigDialog->addDouble(tr("Hauteur de stabilisation de la végétation (Hsv)"), "cm", 0, 9999, 2, _HsvDefault, 100);
    }
}


void ONF_StepComputeProtectionForestMask::compute()
{
    // get DTM raster
    const CT_Image2D<float>* mnt = nullptr;
    for (const CT_Image2D<float>* constMnt : _inDTM.iterateInputs(_inResultDTM))
    {
        mnt = constMnt;
    }
    if (mnt == nullptr)
    {
        PS_LOG->addMessage(LogInterface::error, LogInterface::step, QString(tr("Pas de MNT fourni")));
        return;
    }


    // get forest types raster
    const CT_AbstractImage2D* forest = nullptr;
    for (const CT_StandardItemGroup* group : _inGroupForest.iterateInputs(_inResultForest))
    {
        forest = group->singularItem(_inForest);
    }
    if (forest == nullptr)
    {
        PS_LOG->addMessage(LogInterface::warning, LogInterface::step, QString(tr("Pas de raster type de forêt.")));
    }

    // Create Snow Height reference by class, from reference file
    QMap<quint8, double> snowHeights;

    if (forest != nullptr)
    {
        QFile fRef(_refFileName);
        if (fRef.exists() && fRef.open(QIODevice::ReadOnly | QIODevice::Text))
        {
            QTextStream stream(&fRef);
            stream.setLocale(_refLocale);

            int colClass        = _refColumns.value("Class", -1);
            int colSnowHeight   = _refColumns.value("Snow Height (cm)", -1);

            if (colClass < 0) {PS_LOG->addMessage(LogInterface::error, LogInterface::step, QString(tr("Champ Class non défini")));}
            if (colSnowHeight < 0) {PS_LOG->addMessage(LogInterface::error, LogInterface::step, QString(tr("Champ Snow Height non défini")));}

            if (colClass >=0 && colSnowHeight >= 0)
            {
                int colMax = colClass;
                if (colSnowHeight > colMax) {colMax = colSnowHeight;}

                size_t cpt = 1;
                for (int i = 0 ; i < _refSkip ; i++) {++cpt; stream.readLine();}
                if (_refHeader) {++cpt; stream.readLine();}

                while (!stream.atEnd())
                {
                    QString line = stream.readLine();
                    cpt++;
                    if (!line.isEmpty())
                    {
                        QStringList values = line.split(_refSeparator);
                        if (values.size() >= colMax)
                        {
                            bool okClass, okSnowHeight;
                            quint8 classe = _refLocale.toUInt(values.at(colClass), &okClass);
                            double snowHeight = _refLocale.toDouble(values.at(colSnowHeight), &okSnowHeight) / 100.0; // convert cm to m

                            if (okClass && okSnowHeight)
                            {
                                snowHeights.insert(classe, snowHeight);
                            } else {
                                PS_LOG->addMessage(LogInterface::warning, LogInterface::step, QString(tr("Valeur non définie, ligne %1").arg(cpt)));
                            }
                        }
                    }
                }
            }
            fRef.close();
        } else {

            PS_LOG->addMessage(LogInterface::warning, LogInterface::step, QString(tr("Fichier de référence invalide")));
            return;
        }
    }

    setProgress(float(10.0));

    // compute protection forest mask
    for (CT_StandardItemGroup* group : _inGroupCHM.iterateOutputs(_inResultCHM))
    {
        for (const CT_Image2D<float>* mnh : group->singularItems(_inCHM))
        {
            if (isStopped()) {return;}

            CT_Image2D<quint8>* mask = new CT_Image2D<quint8>(mnh->minX(), mnh->minY(), mnh->xdim(), mnh->ydim(), mnh->resolution(), mnh->minZ(), 255, 0);
            CT_Image2D<quint8>* mask_kenel = new CT_Image2D<quint8>(mnh->minX(), mnh->minY(), mnh->xdim(), mnh->ydim(), mnh->resolution(), mnh->minZ(), 255, 0);


            // Compute base mask
            size_t ncells = mnh->nCells();
            for (size_t index = 0 ; index < ncells ; index++)
            {
                Eigen::Vector3d center;
                mnh->getCellCenterCoordinates(index, center);

                float mnhVal = mnh->valueAtIndex(index);
                float mntVal = mnt->valueAtCoords(center(0), center(1));


                if (!qFuzzyCompare(mnhVal, mnh->NA()) && mntVal >= _altitudeMin && mntVal <= _altitudeMax)
                {
                    double refSH = _HsvDefault;

                    if (forest != nullptr)
                    {
                        size_t forestIndex;
                        forest->indexAtCoords(center(0), center(1), forestIndex);
                        int forestType = int(forest->valueAtIndexAsDouble(forestIndex));

                        refSH = snowHeights.value(forestType, -1);
                    }

                    if (refSH < 0)
                    {
                        mask->setValueAtIndex(index, 0);
                    } else {
                        double height = mnhVal - refSH;
                        double refHeight = _HnZRef;

                        if (mntVal < _ZRef)
                        {
                            refHeight = _HnZRef - ((_ZRef - mntVal) / 100.0)*_gradientInfZ;
                        } else if (mntVal > _ZRef)
                        {
                            refHeight = _HnZRef + ((mntVal - _ZRef) / 100.0)*_gradientSupZ;
                        }

                        if (height > refHeight)
                        {
                            mask->setValueAtIndex(index, 1);
                        } else {
                            mask->setValueAtIndex(index, 0);
                        }
                    }
                } else {
                    mask->setValueAtIndex(index, 0);
                }

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

            // Compute kernel mask
            for (int xx = 0 ; xx < mask->xdim() ; xx++)
            {
                for (int yy = 0 ; yy < mask->ydim() ; yy++)
                {
                    double p1 = 0;
                    double p2 = 0;

                    // 3x3
                    if (mask->value(xx - 1, yy - 1) == 1) {++p1; ++p2;}
                    if (mask->value(xx    , yy - 1) == 1) {++p1; ++p2;}
                    if (mask->value(xx + 1, yy - 1) == 1) {++p1; ++p2;}

                    if (mask->value(xx - 1, yy    ) == 1) {++p1; ++p2;}
                    if (mask->value(xx    , yy    ) == 1) {++p1; ++p2;}
                    if (mask->value(xx + 1, yy    ) == 1) {++p1; ++p2;}

                    if (mask->value(xx - 1, yy + 1) == 1) {++p1; ++p2;}
                    if (mask->value(xx    , yy + 1) == 1) {++p1; ++p2;}
                    if (mask->value(xx + 1, yy + 1) == 1) {++p1; ++p2;}

                    // adding tests for 5x5
                    if (mask->value(xx - 2, yy - 2) == 1) {++p2;}
                    if (mask->value(xx - 1, yy - 2) == 1) {++p2;}
                    if (mask->value(xx    , yy - 2) == 1) {++p2;}
                    if (mask->value(xx + 1, yy - 2) == 1) {++p2;}
                    if (mask->value(xx + 2, yy - 2) == 1) {++p2;}

                    if (mask->value(xx - 2, yy - 1) == 1) {++p2;}
                    if (mask->value(xx + 2, yy - 1) == 1) {++p2;}

                    if (mask->value(xx - 2, yy    ) == 1) {++p2;}
                    if (mask->value(xx + 2, yy    ) == 1) {++p2;}

                    if (mask->value(xx - 2, yy + 1) == 1) {++p2;}
                    if (mask->value(xx + 2, yy + 1) == 1) {++p2;}

                    if (mask->value(xx - 2, yy + 2) == 1) {++p2;}
                    if (mask->value(xx - 1, yy + 2) == 1) {++p2;}
                    if (mask->value(xx    , yy + 2) == 1) {++p2;}
                    if (mask->value(xx + 1, yy + 2) == 1) {++p2;}
                    if (mask->value(xx + 2, yy + 2) == 1) {++p2;}

                    double meanP = (p1 / 9.0 + p2 / 25.0) / 2.0;

                    if (meanP > _threshold)
                    {
                        mask_kenel->setValue(xx, yy, 1);
                    }
                }
            }


            // ajout des rasters Masques
            group->addSingularItem(_outMask, mask);
            mask->computeMinMax();

            group->addSingularItem(_outMask_kernel, mask_kenel);
            mask->computeMinMax();
        }
        setProgress(100.0f);
    }
}