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

#include "ct_itemdrawable/ct_image2d.h"
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"

#include <QDebug>


ONF_MetricRasterCrown4::ONF_MetricRasterCrown4(QString pluginName) : CT_AbstractMetric_Raster(pluginName)
{
    declareAttributes();
}

ONF_MetricRasterCrown4::ONF_MetricRasterCrown4(const ONF_MetricRasterCrown4 &other) : CT_AbstractMetric_Raster(other)
{
    declareAttributes();
    m_configAndResults = other.m_configAndResults;
}

QString ONF_MetricRasterCrown4::getShortDisplayableName() const
{
    return tr("Métriques de houppiers (Ht)");
}

QString ONF_MetricRasterCrown4::getShortDescription() const
{
    return tr("Ces métriques décrivent la taille et la forme globale d'un houppier.<br>"
              "Elles sont conçues pour être calculées sur un modèle numérique de hauteur l'échelle de l'arbre.<br>"
              "Elles nécessitent que le raster soit en hauteur. Un calcul sur des altitudes sera faux.<br>");
}

QString ONF_MetricRasterCrown4::getDetailledDescription() const
{
    return tr("Les valeurs suivantes peuvent être calculées. A noter que seuls les pixels situés au dessus de 2 m sont pris en compte.<br>"
              "<ul>"
              "<li><strong>Volume</strong> : volume du houppier (en m3). Volume situé entre les pixels composant le houppier et le sol (somme des volumes des colonnes).<br></li>"
              "<li><strong>MaxExtH</strong> : hauteur de plus grande largeur de houppier (en m). Une analyse détermine les pixels de bordure (à l'aide d'une érosion du masque du houppier d'un pixel). MaxExtH est égal au percentile 5% des hauteurs des pixels de bordure.<br></li>"
              "<li><strong>MaxH</strong> : hauteur (en m) du pixel le plus haut, l'apex<br></li>"
              "<li><strong>VAbMEH</strong> : volume de la partie supérieure du houppier (en m3). Calculé comme le volume entre le plan horizontal situé à MaxExtH et la surface supérieure du houppier.<br></li>"
              "<li><strong>MEH_H</strong> : rapport entre la hauteur de plus grande largeur de houppier (MaxExtH) et la hauteur de l'apex (MaxH)<br></li>"
              "</ul>"
              );
}

ONF_MetricRasterCrown4::Config ONF_MetricRasterCrown4::metricConfiguration() const
{
    return m_configAndResults;
}

void ONF_MetricRasterCrown4::setMetricConfiguration(const ONF_MetricRasterCrown4::Config &conf)
{
    m_configAndResults = conf;
}

CT_AbstractConfigurableElement *ONF_MetricRasterCrown4::copy() const
{
    return new ONF_MetricRasterCrown4(*this);
}

void ONF_MetricRasterCrown4::computeMetric()
{
    m_configAndResults.volume.value = 0;
    m_configAndResults.maxExtH.value = 0;
    m_configAndResults.maxH.value = 0;
    m_configAndResults.volAboveMaxExtH.value = 0;
    m_configAndResults.maxExtH_H.value = 0;

    double inNA = _inRaster->NAAsDouble();

    // create binary mask
    CT_Image2D<quint8> mask(_inRaster->minX(), _inRaster->minY(), _inRaster->xdim(), _inRaster->ydim(), _inRaster->resolution(), _inRaster->level(), 0, 0);

    for (size_t index = 0 ; index < _inRaster->nCells() ; index++)
    {
        double val = _inRaster->valueAtIndexAsDouble(index);


        if (val != inNA)
        {
            m_configAndResults.volume.value += val * _inRaster->resolution() * _inRaster->resolution();
            mask.setValueAtIndex(index, 1);

            if (val > m_configAndResults.maxH.value) {m_configAndResults.maxH.value = val;}
        }
    }

    double percentile = 0.05;
    double hThreshold = 2.0;

    cv::Mat_<quint8> maskMat = mask.getMat();
    cv::Mat element = cv::getStructuringElement(cv::MORPH_CROSS, cv::Size(3, 3));
    cv::erode(maskMat, maskMat, element, cv::Point(-1,-1), 2);

    QList<double> outHeights;
    for (size_t index = 0 ; index < _inRaster->nCells() ; index++)
    {
        double val = _inRaster->valueAtIndexAsDouble(index);

        if (val != inNA && val > hThreshold)
        {
            int col = 0, lin = 0;
            _inRaster->indexToGrid(index, col, lin);
            if (mask.value(col, lin) <= 0)
            {
                outHeights.append(val);
            }
        }
    }
    std::sort(outHeights.begin(), outHeights.end());
    double maxExtendHeight = computePercentile(outHeights, percentile);

    double sumH = 0;
    double nH = 0;
    for (size_t index = 0 ; index < _inRaster->nCells() ; index++)
    {
        double val = _inRaster->valueAtIndexAsDouble(index);

        if (val != inNA && val >= maxExtendHeight)
        {
            double deltaH = (val - maxExtendHeight);
            sumH += deltaH;
            ++nH;
        }
    }


    m_configAndResults.maxExtH.value = maxExtendHeight;
    m_configAndResults.volAboveMaxExtH.value = sumH * _inRaster->resolution() * _inRaster->resolution();

    if (m_configAndResults.maxH.value > 0) {
        m_configAndResults.maxExtH_H.value = m_configAndResults.maxExtH.value / m_configAndResults.maxH.value;
    }


    setAttributeValueVaB(m_configAndResults.volume);
    setAttributeValueVaB(m_configAndResults.maxExtH);
    setAttributeValueVaB(m_configAndResults.maxH);
    setAttributeValueVaB(m_configAndResults.volAboveMaxExtH);
    setAttributeValueVaB(m_configAndResults.maxExtH_H);
}

void ONF_MetricRasterCrown4::declareAttributes()
{
    registerAttributeVaB(m_configAndResults.volume, CT_AbstractCategory::DATA_NUMBER, "Volume");
    registerAttributeVaB(m_configAndResults.maxExtH, CT_AbstractCategory::DATA_NUMBER, "MaxExtH");
    registerAttributeVaB(m_configAndResults.maxH, CT_AbstractCategory::DATA_NUMBER, "MaxH");
    registerAttributeVaB(m_configAndResults.volAboveMaxExtH, CT_AbstractCategory::DATA_NUMBER, "VAbMEH");
    registerAttributeVaB(m_configAndResults.maxExtH_H, CT_AbstractCategory::DATA_NUMBER, "MEH_H");
}

double ONF_MetricRasterCrown4::computePercentile(const QList<double> &array, const double &p)
{
    int arraySize = array.size();

    // Second Variant, show wikipedia "Percentile"
    double v = ((double)(arraySize-1)) * p;
    int ip1 = (int)v;
    double f = (v-ip1); // (arraySize-1)*p = ip1+f

    int ip2 = ip1 + 1;

    if(ip2 == arraySize)
        return array[ip1];

    if(f == 0)
        return array[ip1];

    return array[ip1] + (f * (array[ip2] - array[ip1]));
}