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

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

#include <QDebug>


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

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

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

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

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

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

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

void ONF_MetricRasterCrown3::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.maxExtH_50.value = 0;
    m_configAndResults.hMoyMaxExtH_50.value = 0;
    m_configAndResults.volAboveMaxExtH_50.value = 0;

    m_configAndResults.maxExtH_75.value = 0;
    m_configAndResults.hMoyMaxExtH_75.value = 0;
    m_configAndResults.volAboveMaxExtH_75.value = 0;

    m_configAndResults.maxExtH_90.value = 0;
    m_configAndResults.hMoyMaxExtH_90.value = 0;
    m_configAndResults.volAboveMaxExtH_90.value = 0;

    m_configAndResults.convexArea_50.value = 0;
    m_configAndResults.maxDist_50.value = 0;

    m_configAndResults.convexArea_75.value = 0;
    m_configAndResults.maxDist_75.value = 0;

    m_configAndResults.convexArea_90.value = 0;
    m_configAndResults.maxDist_90.value = 0;

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

    QList<Eigen::Vector2d *> allPoints;
    Eigen::Vector3d cellCenter;
    Eigen::Vector3d apex;
    double angle = 0;
    double val = 0;
    double dist = 0;
    double asinx = 0;
    double acosy = 0;
    double dx = 0;
    double dy = 0;

    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++)
    {
        val = _inRaster->valueAtIndexAsDouble(index);

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

            if (val > apex(2))
            {
                apex(0) = cellCenter(0);
                apex(1) = cellCenter(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);
                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> outHeightsNNE;
        QList<double> outHeightsNEE;
        QList<double> outHeightsESE;
        QList<double> outHeightsSSE;
        QList<double> outHeightsSSO;
        QList<double> outHeightsOSO;
        QList<double> outHeightsNOO;
        QList<double> outHeightsNNO;

        double M_1PIs4 = 1.0*M_PI / 4.0;
        double M_2PIs4 = 2.0*M_PI / 4.0;
        double M_3PIs4 = 3.0*M_PI / 4.0;
        double M_4PIs4 = 4.0*M_PI / 4.0;
        double M_5PIs4 = 5.0*M_PI / 4.0;
        double M_6PIs4 = 6.0*M_PI / 4.0;
        double M_7PIs4 = 7.0*M_PI / 4.0;
        double M_8PIs4 = 8.0*M_PI / 4.0;

        for (size_t index = 0 ; index < _inRaster->nCells() ; index++)
        {
            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)
                {
                    _inRaster->getCellCenterCoordinates(index, cellCenter);
                    dx =  cellCenter(0) - apex(0);
                    dy =  cellCenter(1) - apex(1);
                    dist = sqrt(pow(dx, 2) + pow(dy, 2));

                    if (dist > 0)
                    {
                        asinx = std::asin(dx/dist);
                        acosy = std::acos(dy/dist);
                        if (asinx >= 0)
                        {
                            angle = acosy;
                        } else {
                            angle = 2*M_PI-acosy;
                        }
                    } else {
                        angle = 0;
                    }

                    if (angle >= 0       && angle < M_1PIs4) {outHeightsNNE.append(val);}
                    if (angle >= M_1PIs4 && angle < M_2PIs4) {outHeightsNEE.append(val);}
                    if (angle >= M_2PIs4 && angle < M_3PIs4) {outHeightsESE.append(val);}
                    if (angle >= M_3PIs4 && angle < M_4PIs4) {outHeightsSSE.append(val);}
                    if (angle >= M_4PIs4 && angle < M_5PIs4) {outHeightsSSO.append(val);}
                    if (angle >= M_5PIs4 && angle < M_6PIs4) {outHeightsOSO.append(val);}
                    if (angle >= M_6PIs4 && angle < M_7PIs4) {outHeightsNOO.append(val);}
                    if (angle >= M_7PIs4 && angle < M_8PIs4) {outHeightsNNO.append(val);}
                }
            }
        }

        qSort(outHeightsNNE.begin(), outHeightsNNE.end());
        qSort(outHeightsNEE.begin(), outHeightsNEE.end());
        qSort(outHeightsESE.begin(), outHeightsESE.end());
        qSort(outHeightsSSE.begin(), outHeightsSSE.end());
        qSort(outHeightsSSO.begin(), outHeightsSSO.end());
        qSort(outHeightsOSO.begin(), outHeightsOSO.end());
        qSort(outHeightsNOO.begin(), outHeightsNOO.end());
        qSort(outHeightsNNO.begin(), outHeightsNNO.end());

        QList<double> extendsHeights;

        if (outHeightsNNE.size() > 0) {extendsHeights.append(computePercentile(outHeightsNNE, percentile));}
        if (outHeightsNEE.size() > 0) {extendsHeights.append(computePercentile(outHeightsNEE, percentile));}
        if (outHeightsESE.size() > 0) {extendsHeights.append(computePercentile(outHeightsESE, percentile));}
        if (outHeightsSSE.size() > 0) {extendsHeights.append(computePercentile(outHeightsSSE, percentile));}
        if (outHeightsSSO.size() > 0) {extendsHeights.append(computePercentile(outHeightsSSO, percentile));}
        if (outHeightsOSO.size() > 0) {extendsHeights.append(computePercentile(outHeightsOSO, percentile));}
        if (outHeightsNOO.size() > 0) {extendsHeights.append(computePercentile(outHeightsNOO, percentile));}
        if (outHeightsNNO.size() > 0) {extendsHeights.append(computePercentile(outHeightsNNO, percentile));}

        qSort(extendsHeights.begin(), extendsHeights.end());

        double maxExtendHeight_50 = computePercentile(extendsHeights, 0.50);
        double maxExtendHeight_75 = computePercentile(extendsHeights, 0.75);
        double maxExtendHeight_90 = computePercentile(extendsHeights, 0.90);


        double sumH_50 = 0;
        double nH_50 = 0;
        double sumH_75 = 0;
        double nH_75 = 0;
        double sumH_90 = 0;
        double nH_90 = 0;

        QList<Eigen::Vector2d *> allPoints_50;
        QList<Eigen::Vector2d *> allPoints_75;
        QList<Eigen::Vector2d *> allPoints_90;

        for (size_t index = 0 ; index < _inRaster->nCells() ; index++)
        {
            val = _inRaster->valueAtIndexAsDouble(index);
            _inRaster->getCellCenterCoordinates(index, cellCenter);
            Eigen::Vector2d *point2D = new Eigen::Vector2d(cellCenter(0), cellCenter(1));


            if (val != inNA && val >= maxExtendHeight_50)
            {
                sumH_50 += (val - maxExtendHeight_50);
                ++nH_50;
                allPoints_50.append(point2D);
            }

            if (val != inNA && val >= maxExtendHeight_75)
            {
                sumH_75 += (val - maxExtendHeight_75);
                ++nH_75;
                allPoints_75.append(point2D);
            }


            if (val != inNA && val >= maxExtendHeight_90)
            {
                sumH_90 += (val - maxExtendHeight_90);
                ++nH_90;
                allPoints_90.append(point2D);
            }
        }


        m_configAndResults.maxDistApex.value = maxDistApex;
        m_configAndResults.maxDist.value = maxDist;
        m_configAndResults.minDistApex.value = minDistApex;
        m_configAndResults.minDist.value = minDist;

        m_configAndResults.maxExtH_50.value = maxExtendHeight_50;
        if (nH_50 > 0) {m_configAndResults.hMoyMaxExtH_50.value = sumH_50 / nH_50;}
        m_configAndResults.volAboveMaxExtH_50.value = sumH_50 * _inRaster->resolution() * _inRaster->resolution();

        m_configAndResults.maxExtH_75.value = maxExtendHeight_75;
        if (nH_75 > 0) {m_configAndResults.hMoyMaxExtH_75.value = sumH_75 / nH_75;}
        m_configAndResults.volAboveMaxExtH_75.value = sumH_75 * _inRaster->resolution() * _inRaster->resolution();


        m_configAndResults.maxExtH_90.value = maxExtendHeight_90;
        if (nH_90 > 0) {m_configAndResults.hMoyMaxExtH_90.value = sumH_90 / nH_90;}
        m_configAndResults.volAboveMaxExtH_90.value = sumH_90 * _inRaster->resolution() * _inRaster->resolution();




        CT_Polygon2DData::orderPointsByXY(allPoints_50);
        CT_Polygon2DData::orderPointsByXY(allPoints_75);
        CT_Polygon2DData::orderPointsByXY(allPoints_90);
        CT_Polygon2DData *data_50 = CT_Polygon2DData::createConvexHull(allPoints_50);
        CT_Polygon2DData *data_75 = CT_Polygon2DData::createConvexHull(allPoints_75);
        CT_Polygon2DData *data_90 = CT_Polygon2DData::createConvexHull(allPoints_90);

        if (data_50 != NULL)
        {
            m_configAndResults.convexArea_50.value = data_50->getArea();

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

            double maxDist = 0;

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

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

                    if (dist > maxDist) {maxDist = dist;}
                }
            }
            m_configAndResults.maxDist_50.value = maxDist;
        }


        if (data_75 != NULL)
        {
            m_configAndResults.convexArea_75.value = data_75->getArea();

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

            double maxDist = 0;

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

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

                    if (dist > maxDist) {maxDist = dist;}
                }
            }
            m_configAndResults.maxDist_75.value = maxDist;
        }

        if (data_90 != NULL)
        {
            m_configAndResults.convexArea_90.value = data_90->getArea();

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

            double maxDist = 0;

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

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

                    if (dist > maxDist) {maxDist = dist;}
                }
            }
            m_configAndResults.maxDist_90.value = maxDist;
        }
    }


    setAttributeValueVaB(m_configAndResults.convexArea);
    setAttributeValueVaB(m_configAndResults.maxDistApex);
    setAttributeValueVaB(m_configAndResults.maxDist);
    setAttributeValueVaB(m_configAndResults.minDistApex);
    setAttributeValueVaB(m_configAndResults.minDist);

    setAttributeValueVaB(m_configAndResults.maxExtH_50);
    setAttributeValueVaB(m_configAndResults.hMoyMaxExtH_50);
    setAttributeValueVaB(m_configAndResults.volAboveMaxExtH_50);

    setAttributeValueVaB(m_configAndResults.maxExtH_75);
    setAttributeValueVaB(m_configAndResults.hMoyMaxExtH_75);
    setAttributeValueVaB(m_configAndResults.volAboveMaxExtH_75);

    setAttributeValueVaB(m_configAndResults.maxExtH_90);
    setAttributeValueVaB(m_configAndResults.hMoyMaxExtH_90);
    setAttributeValueVaB(m_configAndResults.volAboveMaxExtH_90);

    setAttributeValueVaB(m_configAndResults.convexArea_50);
    setAttributeValueVaB(m_configAndResults.maxDist_50);

    setAttributeValueVaB(m_configAndResults.convexArea_75);
    setAttributeValueVaB(m_configAndResults.maxDist_75);

    setAttributeValueVaB(m_configAndResults.convexArea_90);
    setAttributeValueVaB(m_configAndResults.maxDist_90);

}

void ONF_MetricRasterCrown3::declareAttributes()
{
    registerAttributeVaB(m_configAndResults.convexArea, CT_AbstractCategory::DATA_NUMBER, tr("3convexArea"));
    registerAttributeVaB(m_configAndResults.maxDistApex, CT_AbstractCategory::DATA_NUMBER, tr("3maxDistApex"));
    registerAttributeVaB(m_configAndResults.maxDist, CT_AbstractCategory::DATA_NUMBER, tr("3maxDist"));
    registerAttributeVaB(m_configAndResults.minDistApex, CT_AbstractCategory::DATA_NUMBER, tr("3minDistApex"));
    registerAttributeVaB(m_configAndResults.minDist, CT_AbstractCategory::DATA_NUMBER, tr("3minDist"));

    registerAttributeVaB(m_configAndResults.maxExtH_50, CT_AbstractCategory::DATA_NUMBER, tr("maxExtH_50"));
    registerAttributeVaB(m_configAndResults.hMoyMaxExtH_50, CT_AbstractCategory::DATA_NUMBER, tr("hMoyMaxExtH_50"));
    registerAttributeVaB(m_configAndResults.volAboveMaxExtH_50, CT_AbstractCategory::DATA_NUMBER, tr("volAboveMaxExtH_50"));

    registerAttributeVaB(m_configAndResults.maxExtH_75, CT_AbstractCategory::DATA_NUMBER, tr("maxExtH_75"));
    registerAttributeVaB(m_configAndResults.hMoyMaxExtH_75, CT_AbstractCategory::DATA_NUMBER, tr("hMoyMaxExtH_75"));
    registerAttributeVaB(m_configAndResults.volAboveMaxExtH_75, CT_AbstractCategory::DATA_NUMBER, tr("volAboveMaxExtH_75"));

    registerAttributeVaB(m_configAndResults.maxExtH_90, CT_AbstractCategory::DATA_NUMBER, tr("maxExtH_90"));
    registerAttributeVaB(m_configAndResults.hMoyMaxExtH_90, CT_AbstractCategory::DATA_NUMBER, tr("hMoyMaxExtH_90"));
    registerAttributeVaB(m_configAndResults.volAboveMaxExtH_90, CT_AbstractCategory::DATA_NUMBER, tr("volAboveMaxExtH_90"));

    registerAttributeVaB(m_configAndResults.convexArea_50, CT_AbstractCategory::DATA_NUMBER, tr("convexArea_50"));
    registerAttributeVaB(m_configAndResults.maxDist_50, CT_AbstractCategory::DATA_NUMBER, tr("maxDist_50"));

    registerAttributeVaB(m_configAndResults.convexArea_75, CT_AbstractCategory::DATA_NUMBER, tr("convexArea_75"));
    registerAttributeVaB(m_configAndResults.maxDist_75, CT_AbstractCategory::DATA_NUMBER, tr("maxDist_75"));

    registerAttributeVaB(m_configAndResults.convexArea_90, CT_AbstractCategory::DATA_NUMBER, tr("convexArea_90"));
    registerAttributeVaB(m_configAndResults.maxDist_90, CT_AbstractCategory::DATA_NUMBER, tr("maxDist_90"));

}

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