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

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"

#include <QtConcurrent>
#include <math.h>
#include <stdlib.h>

ONF_StepComputeGapsJMM::ONF_StepComputeGapsJMM() : SuperClass()
{
    _gapMaxHeight = 10.0;
    _minGapSurface = 25.0;
    _maxGapSurface = 99999.0;
    _closingHeightBin = 1.0;
    _ratio = 2;
    _gapReconstruct = false;
}

QString ONF_StepComputeGapsJMM::description() const
{
    return tr("Créer une carte des trouées");
}

QString ONF_StepComputeGapsJMM::detailledDescription() const
{
    return tr("Cette étape créée une carte des trouées selon la méthode proposée par Jean-Matthieu Monnet dans le package lidaRtRee.");
}

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

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

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

QString ONF_StepComputeGapsJMM::detailsDescription() const
{
    return tr("Idéalement, fournir en entrée de cette étape un Modèle Numérique de Hauteur ayant bénéficié d'un remplissage de trous, ainsi que d'un léger lissage. ");
}

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

void ONF_StepComputeGapsJMM::declareInputModels(CT_StepInModelStructureManager& manager)
{
    manager.addResult(_inResult, tr("Scene(s)"));
    manager.setZeroOrMoreRootGroup(_inResult, _inZeroOrMoreRootGroup);
    manager.addGroup(_inZeroOrMoreRootGroup, _inGroup);
    manager.addItem(_inGroup, _inDEM, tr("MNH (hauteurs)"));
}



void ONF_StepComputeGapsJMM::fillPostInputConfigurationDialog(CT_StepConfigurableDialog* postInputConfigDialog)
{
    postInputConfigDialog->addDouble(tr("Seuil de hauteur de trouée"), "m", 0, 99999, 2, _gapMaxHeight, 1, tr("Hauteur maximale de canopée pour être considéré dans une trouée."));
    postInputConfigDialog->addDouble(tr("Surface minimale de trouée"), "m²", 0, 99999, 2, _minGapSurface);
    postInputConfigDialog->addDouble(tr("Surface maximale de trouée"), "m²", 0, 99999, 2, _maxGapSurface);
    postInputConfigDialog->addDouble(tr("incréments de hauteur (fermetures mrophologiques)"), "m", 0, 99999, 4, _closingHeightBin);
    postInputConfigDialog->addDouble(tr("Ratio hauteur / distance"), "", 0, 99999, 2, _ratio, 1, tr("rapport maximal entre la hauteur de la canopée environnante et la distance de la trouée (un pixel appartient à la trouée uniquement si, pour tout pixel de végétation autour de lui, la distance au pixel de végétation est supérieure à la hauteur du pixel/ratio). Si le rapport est fixé à 0, ce critère n'est pas pris en compte."));
    postInputConfigDialog->addBool(tr("Reconstruire la trouée"), "m", "", _gapReconstruct, tr(("Par défaut, les zones qui ne remplissent pas le critère de ratio sont supprimées des trouées. Si cette option est activée, si certains pixels d'un espace remplissent le critère de distance, les pixels connectés qui remplissent le critère de hauteur sont également intégrés à l'espace.")));
}

void ONF_StepComputeGapsJMM::declareOutputModels(CT_StepOutModelStructureManager& manager)
{
    manager.addResultCopy(_inResult);
    manager.addItem(_inGroup, _outGapIDs, tr("Trouées"));
}

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

            CT_Image2D<qint32>* gapIDs = new CT_Image2D<int>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->minZ(), -9999, -1);
            CT_Image2D<quint8>* tmpMask = new CT_Image2D<quint8>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->minZ(), 0, 0);
            CT_Image2D<quint8>* tmpMaskClosing = new CT_Image2D<quint8>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->minZ(), 0, 0);
            CT_Image2D<quint8>* tmpMaskUnion = new CT_Image2D<quint8>(dem->minX(), dem->minY(), dem->xdim(), dem->ydim(), dem->resolution(), dem->minZ(), 0, 0);

            cv::Mat_<quint8> tmpMaskMat = tmpMask->getMat();
            cv::Mat_<quint8> tmpMaskClosingMat = tmpMaskClosing->getMat();

            size_t ncells = dem->nCells();

            // _ratio == 0 case: gap mask is simply h <= _gapMaxHeight
            if (_gapReconstruct || qFuzzyCompare(_ratio, 0))
            {
                for (size_t index = 0 ; index < ncells ; index++)
                {
                    float hDEM = dem->valueAtIndex(index);

                        if (qFuzzyCompare(hDEM, dem->NA()) || (hDEM > float(_gapMaxHeight)))
                        {
                            tmpMask->setValueAtIndex(index, 0);
                        } else {
                            tmpMask->setValueAtIndex(index, 1);
                        }
                }

                cv::Mat labels = cv::Mat::zeros(tmpMaskMat.rows, tmpMaskMat.cols, CV_32S); // Nécessaire car compare ne prend pas la version template Mat_<Tp> en output !!!
                cv::connectedComponents(tmpMaskMat, labels);
                gapIDs->getMat() = labels;

                gapIDs->computeMinMax();
            }

            if (_ratio > 0)
            {
                float hmaxMNH = std::fmin(dem->dataMax(), 60.0);

                for (float hmin = _gapMaxHeight ; hmin <= hmaxMNH ; hmin += _closingHeightBin)
                {
                    for (size_t index = 0 ; index < ncells ; index++)
                    {
                        float hDEM = dem->valueAtIndex(index);

                        if (!qFuzzyCompare(hDEM, dem->NA()) && hDEM > hmin)
                        {
                            tmpMask->setValueAtIndex(index, 1);
                        } else {
                            tmpMask->setValueAtIndex(index, 0);
                        }
                    }

                    int size = int(std::floor(hmin / _ratio / dem->resolution() / 2.0) * 2.0 + 1.0);

                    // Replaced by home made structuring Element because don't do the same shape than R::lidaRtRee
                    // cv::Mat element = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(size, size));

                    cv::Mat element;
                    createDisk(size, element);
                    cv::morphologyEx(tmpMaskMat, tmpMaskClosingMat, cv::MORPH_CLOSE, element);

                    for (size_t index = 0 ; index < ncells ; index++)
                    {
                        quint8 val = tmpMaskClosing->valueAtIndex(index);
                        float hDEM = dem->valueAtIndex(index);

                        if (!qFuzzyCompare(hDEM, dem->NA()) && val > 0)
                        {
                            tmpMaskUnion->setValueAtIndex(index, 1);
                        }
                    }
                }

                if (_gapReconstruct)
                {
                    // detect closed gaps
                    QVector<bool> closed(gapIDs->dataMax() + 1);
                    closed.fill(true);

                    for (size_t index = 0 ; index < ncells ; index++)
                    {
                        quint8 val = tmpMaskUnion->valueAtIndex(index);
                        qint32 label = gapIDs->valueAtIndex(index);

                        if (label > 0 && val == 0)
                        {
                            closed[label] = false;
                        }
                    }

                    // remove closed gaps
                    for (size_t index = 0 ; index < ncells ; index++)
                    {
                        qint32 label = gapIDs->valueAtIndex(index);

                        if (label > 0 && closed[label])
                        {
                            gapIDs->setValueAtIndex(index, 0);
                        }
                    }

                } else {
                    // invert canopy mask to get gap mask
                    for (size_t index = 0 ; index < ncells ; index++)
                    {
                        quint8 val = tmpMaskUnion->valueAtIndex(index);
                        float hDEM = dem->valueAtIndex(index);

                        if (qFuzzyCompare(hDEM, dem->NA()) || val == 1)
                        {
                            tmpMask->setValueAtIndex(index, 0);
                        } else {
                            tmpMask->setValueAtIndex(index, 1);
                        }
                    }

                    cv::Mat labels = cv::Mat::zeros(tmpMaskMat.rows, tmpMaskMat.cols, CV_32S); // Nécessaire car compare ne prend pas la version template Mat_<Tp> en output !!!
                    cv::connectedComponents(tmpMaskMat, labels);
                    gapIDs->getMat() = labels;
                }
            }

            gapIDs->computeMinMax();

            // compute gaps areas
            QVector<float> areas(gapIDs->dataMax() + 1);
            areas.fill(0);

            double res2 = dem->resolution()*dem->resolution();
            for (size_t index = 0 ; index < ncells ; index++)
            {
                qint32 label = gapIDs->valueAtIndex(index);

                if (label > 0)
                {
                    areas[label] += res2;
                }
            }

            // Select gaps to keep (with area in good area range)
            QVector<bool> keepAreas(gapIDs->dataMax() + 1);
            keepAreas.fill(false);
            for (int i = 1 ; i < keepAreas.size() ; i++)
            {
                if (areas[i] >= _minGapSurface && areas[i] <= _maxGapSurface)
                {
                    keepAreas[i] = true;
                }
            }

            for (size_t index = 0 ; index < ncells ; index++)
            {
                qint32 label = gapIDs->valueAtIndex(index);

                if (label > 0 && keepAreas[label])
                {
                    tmpMask->setValueAtIndex(index, 1);
                } else {
                    tmpMask->setValueAtIndex(index, 0);
                }
            }

            cv::Mat labels = cv::Mat::zeros(tmpMaskMat.rows, tmpMaskMat.cols, CV_32S); // Nécessaire car compare ne prend pas la version template Mat_<Tp> en output !!!
            cv::connectedComponents(tmpMaskMat, labels);
            gapIDs->getMat() = labels;


            delete tmpMask;
            delete tmpMaskClosing;
            delete tmpMaskUnion;

            // ajout du raster MNS
            group->addSingularItem(_outGapIDs, gapIDs);
            gapIDs->computeMinMax();
        }

    }

    setProgress(100);
}

void ONF_StepComputeGapsJMM::createDisk(int ksize, cv::Mat &element)
{
    element.create(ksize, ksize, CV_8U);

    int radius = ksize / 2;
    int radius2 = radius * radius;

    for (int yy = 0 ; yy < element.rows ; yy++)
    {
        for (int xx = 0 ; xx < element.cols ; xx++)
        {
            int dist = pow(xx - radius, 2) + pow(yy - radius, 2);

            if (dist <= radius2)
            {
                element.at<quint8>(yy,xx) = 1;
            } else {
                element.at<quint8>(yy,xx) = 0;
            }
        }
    }
}