/****************************************************************************
 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 "metric/onf_metricrastercrown2.h"
#include "ct_shape2ddata/ct_polygon2ddata.h"

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

#include <QDebug>


ONF_MetricRasterCrown2::ONF_MetricRasterCrown2() : CT_AbstractMetric_Raster()
{
    declareAttributes();
}

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

QString ONF_MetricRasterCrown2::getShortDescription() const
{
    return tr("Métriques houppiers (test, H)");
}

QString ONF_MetricRasterCrown2::getDetailledDescription() const
{
    return tr("");
}

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

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

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

void ONF_MetricRasterCrown2::computeMetric()
{
    m_configAndResults.convexArea.value = 0;
    m_configAndResults.maxDistApex.value = 0;
    m_configAndResults.maxDist.value = 0;
    m_configAndResults.minDistApex.value = 0;
    m_configAndResults.minDist.value = 0;
    m_configAndResults.maxExtendHeight.value = 0;
    m_configAndResults.hMoyAboveMaxExtendHeight.value = 0;
    m_configAndResults.volumeAboveMaxExtendHeight.value = 0;

    double inNA = _inRaster->NAAsDouble();
    double inMin = _inRaster->minValueAsDouble();

    QList<Eigen::Vector2d *> allPoints;
    Eigen::Vector3d center;
    Eigen::Vector3d apex;

    apex(0) = -std::numeric_limits<double>::max();
    apex(1) = -std::numeric_limits<double>::max();
    apex(2) = -std::numeric_limits<double>::max();

    CT_Image2D<quint8> mask(NULL, NULL, _inRaster->minX(), _inRaster->minY(), _inRaster->colDim(), _inRaster->linDim(), _inRaster->resolution(), _inRaster->level(), 0, 0);


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

        if (val != inNA && val != inMin)
        {
            mask.setValueAtIndex(index, 1);
            _inRaster->getCellCenterCoordinates(index, center);
            Eigen::Vector2d *point2D = new Eigen::Vector2d(center(0), center(1));
            allPoints.append(point2D);

            if (val > apex(2))
            {
                apex(0) = center(0);
                apex(1) = center(1);
                apex(2) = val;
            }
        }
    }

    // tri par (X,Y) de la liste des points
    CT_Polygon2DData::orderPointsByXY(allPoints);
    CT_Polygon2DData *data = CT_Polygon2DData::createConvexHull(allPoints);

    if (data != NULL)
    {
        m_configAndResults.convexArea.value = data->getArea();

        const QVector<Eigen::Vector2d*>& vertices = data->getVertices();

        double maxDistApex = 0;
        double maxDist = 0;

        double minDistApex = std::numeric_limits<double>::max();
        double minDist = std::numeric_limits<double>::max();


        for (int i = 0 ; i < vertices.size() ; i++)
        {
            Eigen::Vector2d* vert1 = vertices.at(i);

            double distApex = sqrt(pow(apex(0) - (*vert1)(0), 2) + pow(apex(1) - (*vert1)(1), 2));

            if (distApex > maxDistApex) {maxDistApex = distApex;}
            if (distApex < minDistApex) {minDistApex = distApex;}

            for (int j = i + 1 ; j < vertices.size() ; j++)
            {
                Eigen::Vector2d* vert2 = vertices.at(j);
                double dist = sqrt(pow((*vert1)(0) - (*vert2)(0), 2) + pow((*vert1)(1) - (*vert2)(1), 2));

                if (dist > maxDist) {maxDist = dist;}
                if (dist < minDist) {minDist = dist;}
            }
        }


        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 != inMin && val > hThreshold)
            {
                size_t col, lin;
                _inRaster->indexToGrid(index, col, lin);
                if (mask.value(col, lin) <= 0)
                {
                    outHeights.append(val);
                }
            }
        }
        qSort(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)
            {
                sumH += (val - maxExtendHeight);
                ++nH;
            }
        }


        m_configAndResults.maxDistApex.value = maxDistApex;
        m_configAndResults.maxDist.value = maxDist;
        m_configAndResults.minDistApex.value = minDistApex;
        m_configAndResults.minDist.value = minDist;
        m_configAndResults.maxExtendHeight.value = maxExtendHeight;
        if (nH > 0) {m_configAndResults.hMoyAboveMaxExtendHeight.value = sumH / nH;}
        m_configAndResults.volumeAboveMaxExtendHeight.value = sumH * _inRaster->resolution() * _inRaster->resolution();
    }


    setAttributeValueVaB(m_configAndResults.convexArea);
    setAttributeValueVaB(m_configAndResults.maxDistApex);
    setAttributeValueVaB(m_configAndResults.maxDist);
    setAttributeValueVaB(m_configAndResults.minDistApex);
    setAttributeValueVaB(m_configAndResults.minDist);
    setAttributeValueVaB(m_configAndResults.maxExtendHeight);
    setAttributeValueVaB(m_configAndResults.hMoyAboveMaxExtendHeight);
    setAttributeValueVaB(m_configAndResults.volumeAboveMaxExtendHeight);
}

void ONF_MetricRasterCrown2::declareAttributes()
{
    registerAttributeVaB(m_configAndResults.convexArea, CT_AbstractCategory::DATA_NUMBER, tr("convexArea"));
    registerAttributeVaB(m_configAndResults.maxDistApex, CT_AbstractCategory::DATA_NUMBER, tr("maxDistApex"));
    registerAttributeVaB(m_configAndResults.maxDist, CT_AbstractCategory::DATA_NUMBER, tr("maxDist"));
    registerAttributeVaB(m_configAndResults.minDistApex, CT_AbstractCategory::DATA_NUMBER, tr("minDistApex"));
    registerAttributeVaB(m_configAndResults.minDist, CT_AbstractCategory::DATA_NUMBER, tr("minDist"));
    registerAttributeVaB(m_configAndResults.maxExtendHeight, CT_AbstractCategory::DATA_NUMBER, tr("maxExtendHeight"));
    registerAttributeVaB(m_configAndResults.hMoyAboveMaxExtendHeight, CT_AbstractCategory::DATA_NUMBER, tr("hMoyAboveMaxExtendHeight"));
    registerAttributeVaB(m_configAndResults.volumeAboveMaxExtendHeight, CT_AbstractCategory::DATA_NUMBER, tr("volumeAboveMaxExtendHeight"));
}

double ONF_MetricRasterCrown2::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]));
}