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

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

#include <QDebug>


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

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

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

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

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.convexArea.value = 0;
    m_configAndResults.maxDistApex.value = 0;
    m_configAndResults.maxDist.value = 0;
    m_configAndResults.maxExtH.value = 0;
    m_configAndResults.hMoyAboveMaxExtH.value = 0;
    m_configAndResults.volAboveMaxExtH.value = 0;
    m_configAndResults.area3dAboveMaxExtH.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;

        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;}

            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;}
            }
        }


        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);

        CT_Image2D<float> upperRaster(NULL, NULL, _inRaster->minX(), _inRaster->minY(), _inRaster->colDim(), _inRaster->linDim(), _inRaster->resolution(), _inRaster->level(), -1, -1);

        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;

                upperRaster.setValueAtIndex(index, deltaH);
            }
        }

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

        double sigma = 0.4 / _inRaster->resolution();
        cv::GaussianBlur(upperRaster.getMat(), gaussianRaster.getMat(), cv::Size2d(0, 0), sigma);


        for (size_t index = 0 ; index < upperRaster.nCells() ; index++)
        {
            double val = upperRaster.valueAtIndexAsDouble(index);
            double valGauss = gaussianRaster.valueAtIndexAsDouble(index);
            size_t xx, yy;

            if (upperRaster.indexToGrid(index, xx, yy) && val != -1)
            {
                // Slope
                size_t index2;
                double a = -1;
                double b = -1;
                double c = -1;
                double d = -1;
                double f = -1;
                double g = -1;
                double h = -1;
                double i = -1;

                if (upperRaster.index(xx - 1, yy - 1, index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {a = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx    , yy - 1, index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {b = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx + 1, yy - 1, index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {c = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx - 1, yy    , index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {d = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx + 1, yy    , index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {f = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx - 1, yy + 1, index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {g = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx    , yy + 1, index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {h = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (upperRaster.index(xx + 1, yy + 1, index2))
                    if (upperRaster.valueAtIndexAsDouble(index2) != -1) {i = gaussianRaster.valueAtIndexAsDouble(index2);}

                if (a == -1) {a = valGauss;}
                if (b == -1) {b = valGauss;}
                if (c == -1) {c = valGauss;}
                if (d == -1) {d = valGauss;}
                if (f == -1) {f = valGauss;}
                if (g == -1) {g = valGauss;}
                if (h == -1) {h = valGauss;}
                if (i == -1) {i = valGauss;}

                double dzdx = ((c + 2.0*f + i) - (a + 2.0*d + g)) / (8.0 * upperRaster.resolution());
                double dzdy = ((g + 2.0*h + i) - (a + 2.0*b + c)) / (8.0 * upperRaster.resolution());

                double slope = std::sqrt(dzdx*dzdx + dzdy*dzdy);

                m_configAndResults.area3dAboveMaxExtH.value += std::sqrt(pow(upperRaster.resolution()*slope, 2) + upperRaster.resolution()*upperRaster.resolution()) * upperRaster.resolution();
            }
        }


        m_configAndResults.maxDistApex.value = maxDistApex;
        m_configAndResults.maxDist.value = maxDist;
        m_configAndResults.maxExtH.value = maxExtendHeight;
        if (nH > 0) {m_configAndResults.hMoyAboveMaxExtH.value = sumH / nH;}
        m_configAndResults.volAboveMaxExtH.value = sumH * _inRaster->resolution() * _inRaster->resolution();
    }


    setAttributeValueVaB(m_configAndResults.convexArea);
    setAttributeValueVaB(m_configAndResults.maxDistApex);
    setAttributeValueVaB(m_configAndResults.maxDist);
    setAttributeValueVaB(m_configAndResults.maxExtH);
    setAttributeValueVaB(m_configAndResults.hMoyAboveMaxExtH);
    setAttributeValueVaB(m_configAndResults.volAboveMaxExtH);
    setAttributeValueVaB(m_configAndResults.area3dAboveMaxExtH);
}

void ONF_MetricRasterCrown4::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.maxExtH, CT_AbstractCategory::DATA_NUMBER, tr("maxExtH"));
    registerAttributeVaB(m_configAndResults.hMoyAboveMaxExtH, CT_AbstractCategory::DATA_NUMBER, tr("hMoyAboveMaxExtH"));
    registerAttributeVaB(m_configAndResults.volAboveMaxExtH, CT_AbstractCategory::DATA_NUMBER, tr("volAboveMaxExtH"));
    registerAttributeVaB(m_configAndResults.area3dAboveMaxExtH, CT_AbstractCategory::DATA_NUMBER, tr("area3dAboveMaxExtH"));
}

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