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

#include <QDebug>


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

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

QString ONF_MetricRasterCrown::getShortDescription() const
{
    return tr("Rumple, volume et pente (houppier)");
}

QString ONF_MetricRasterCrown::getDetailledDescription() const
{
    return tr("Les valeurs suivantes sont calculées :<br>"
              "- rumple<br>"
              "- area3d<br>"
              "- area2d<br>"
              "- crownThickness<br>"
              "- crownThicknessQ25<br>"
              "- crownThicknessQ50<br>"
              "- crownThicknessQ75<br>"
              "- crownThicknessQ95<br>"
              "- crownEquivRadius<br>"
              "- volume (m3)<br>"
              "- volume_Q25 (m3)<br>"
              "- volume_Q50 (m3)<br>"
              "- volume_Q75 (m3)<br>"
              "- slope_max (degrees)<br>"
              "- slope_min (degrees)<br>"
              "- slope_moy (degrees)<br>"
              "- slope_sd (degrees)<br>"
              "- slope_Q25 (degrees)<br>"
              "- slope_Q50 - median (degrees)<br>"
              "- slope_Q75 (degrees)<br>"
              );
}

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

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

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

void ONF_MetricRasterCrown::computeMetric()
{    
    m_configAndResults.rumple.value = 0;
    m_configAndResults.volume.value = 0;
    m_configAndResults.volume_topQ25.value = 0;
    m_configAndResults.volume_topQ50.value = 0;
    m_configAndResults.volume_topQ75.value = 0;
    m_configAndResults.slope_max.value = 0;
    m_configAndResults.slope_min.value = std::numeric_limits<double>::max();
    m_configAndResults.slope_moy.value = 0;
    m_configAndResults.slope_sd.value = 0;
    m_configAndResults.slope_Q25.value = 0;
    m_configAndResults.slope_Q50.value = 0;
    m_configAndResults.slope_Q75.value = 0;

    m_configAndResults.area3d.value = 0;
    m_configAndResults.area2d.value = 0;
    m_configAndResults.crownThickness.value = 0;
    m_configAndResults.crownThicknessQ25.value = 0;
    m_configAndResults.crownThicknessQ50.value = 0;
    m_configAndResults.crownThicknessQ75.value = 0;
    m_configAndResults.crownThicknessQ95.value = 0;
    m_configAndResults.crownEquivRadius.value = 0;

    // ((crownThickness / crownEquivRadius) + 1)^(2/3) * area3d^(3/2)
    // ((crownThicknessQ95 / crownEquivRadius) + 1)^(2/3) * area3d^(3/2)

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

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

        if (val != inNA && val != inMin)
        {
            zvals.append(val);
        }
    }
    qSort(zvals.begin(), zvals.end());

    double zQ25 = computePercentile(zvals, 0.75); // Quantiles inversés  car depuis le haut
    double zQ50 = computePercentile(zvals, 0.50);
    double zQ75 = computePercentile(zvals, 0.25);
    double zQ95 = computePercentile(zvals, 0.05);

    m_configAndResults.crownThickness.value = zvals.last() - zvals.first();
    m_configAndResults.crownThicknessQ25.value = zvals.last() - zQ25;
    m_configAndResults.crownThicknessQ50.value = zvals.last() - zQ50;
    m_configAndResults.crownThicknessQ75.value = zvals.last() - zQ75;
    m_configAndResults.crownThicknessQ95.value = zvals.last() - zQ95;

    QList<double> slopeValues;
    for (size_t index = 0 ; index < _inRaster->nCells() ; index++)
    {
        double val = _inRaster->valueAtIndexAsDouble(index);
        size_t xx, yy;

        if (_inRaster->indexToGrid(index, xx, yy) && val != inNA)
        {
            // Volume
            m_configAndResults.volume.value += val * _inRaster->resolution() * _inRaster->resolution();

            if (val > zQ25) {m_configAndResults.volume_topQ25.value += (val - zQ25) * _inRaster->resolution() * _inRaster->resolution();}
            if (val > zQ50) {m_configAndResults.volume_topQ50.value += (val - zQ50) * _inRaster->resolution() * _inRaster->resolution();}
            if (val > zQ75) {m_configAndResults.volume_topQ75.value += (val - zQ75) * _inRaster->resolution() * _inRaster->resolution();}

            // Slope
            size_t index2;
            double a = inNA;
            double b = inNA;
            double c = inNA;
            double d = inNA;
            double f = inNA;
            double g = inNA;
            double h = inNA;
            double i = inNA;

            if (_inRaster->index(xx - 1, yy - 1, index2))
                {a = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx    , yy - 1, index2))
                {b = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx + 1, yy - 1, index2))
                {c = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx - 1, yy    , index2))
                {d = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx + 1, yy    , index2))
                {f = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx - 1, yy + 1, index2))
                {g = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx    , yy + 1, index2))
                {h = _inRaster->valueAtIndexAsDouble(index2);}

            if (_inRaster->index(xx + 1, yy + 1, index2))
                {i = _inRaster->valueAtIndexAsDouble(index2);}

            if (a == inNA) {a = val;}
            if (b == inNA) {b = val;}
            if (c == inNA) {c = val;}
            if (d == inNA) {d = val;}
            if (f == inNA) {f = val;}
            if (g == inNA) {g = val;}
            if (h == inNA) {h = val;}
            if (i == inNA) {i = val;}

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

            double slope = std::sqrt(dzdx*dzdx + dzdy*dzdy);
            //double slope_percent = 100.0*slope;
            double slope_radians = std::atan(std::sqrt(dzdx*dzdx + dzdy*dzdy));
            double slope_degrees = slope_radians * 180.0 / M_PI;

            slopeValues.append(slope_degrees);

            n += 1.0;
            if (slope_degrees < m_configAndResults.slope_min.value) {m_configAndResults.slope_min.value = slope_degrees;}
            if (slope_degrees > m_configAndResults.slope_max.value) {m_configAndResults.slope_max.value = slope_degrees;}
            m_configAndResults.slope_moy.value += slope_degrees;
            sumx2 += slope_degrees*slope_degrees;

            // for Rumple
            m_configAndResults.area3d.value += std::sqrt(pow(_inRaster->resolution()*slope, 2) + _inRaster->resolution()*_inRaster->resolution()) * _inRaster->resolution();
        }
    }

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

    m_configAndResults.slope_Q25.value = computePercentile(slopeValues, 0.25);
    m_configAndResults.slope_Q50.value = computePercentile(slopeValues, 0.50);
    m_configAndResults.slope_Q75.value = computePercentile(slopeValues, 0.75);

    // Slope
    if (n > 0)
    {
        m_configAndResults.slope_moy.value /= n;
        m_configAndResults.slope_sd.value = std::sqrt(sumx2/n - m_configAndResults.slope_moy.value*m_configAndResults.slope_moy.value);
        m_configAndResults.area2d.value = n*_inRaster->resolution()*_inRaster->resolution();
        m_configAndResults.rumple.value = m_configAndResults.area3d.value / m_configAndResults.area2d.value;
        m_configAndResults.crownEquivRadius.value = sqrt(m_configAndResults.area2d.value / M_PI);
    }

    setAttributeValueVaB(m_configAndResults.volume);
    setAttributeValueVaB(m_configAndResults.volume_topQ25);
    setAttributeValueVaB(m_configAndResults.volume_topQ50);
    setAttributeValueVaB(m_configAndResults.volume_topQ75);
    setAttributeValueVaB(m_configAndResults.slope_max);
    setAttributeValueVaB(m_configAndResults.slope_min);
    setAttributeValueVaB(m_configAndResults.slope_moy);
    setAttributeValueVaB(m_configAndResults.slope_sd);
    setAttributeValueVaB(m_configAndResults.slope_Q25);
    setAttributeValueVaB(m_configAndResults.slope_Q50);
    setAttributeValueVaB(m_configAndResults.slope_Q75);
    setAttributeValueVaB(m_configAndResults.rumple);
    setAttributeValueVaB(m_configAndResults.area3d);
    setAttributeValueVaB(m_configAndResults.area2d);
    setAttributeValueVaB(m_configAndResults.crownThickness);
    setAttributeValueVaB(m_configAndResults.crownThicknessQ25);
    setAttributeValueVaB(m_configAndResults.crownThicknessQ50);
    setAttributeValueVaB(m_configAndResults.crownThicknessQ75);
    setAttributeValueVaB(m_configAndResults.crownThicknessQ95);
    setAttributeValueVaB(m_configAndResults.crownEquivRadius);
}

void ONF_MetricRasterCrown::declareAttributes()
{
    registerAttributeVaB(m_configAndResults.rumple, CT_AbstractCategory::DATA_NUMBER, tr("rumple"));
    registerAttributeVaB(m_configAndResults.area3d, CT_AbstractCategory::DATA_NUMBER, tr("area3d"));
    registerAttributeVaB(m_configAndResults.area2d, CT_AbstractCategory::DATA_NUMBER, tr("area2d"));
    registerAttributeVaB(m_configAndResults.crownThickness, CT_AbstractCategory::DATA_NUMBER, tr("crownThickness"));
    registerAttributeVaB(m_configAndResults.crownThicknessQ25, CT_AbstractCategory::DATA_NUMBER, tr("crownThicknessQ25"));
    registerAttributeVaB(m_configAndResults.crownThicknessQ50, CT_AbstractCategory::DATA_NUMBER, tr("crownThicknessQ50"));
    registerAttributeVaB(m_configAndResults.crownThicknessQ75, CT_AbstractCategory::DATA_NUMBER, tr("crownThicknessQ75"));
    registerAttributeVaB(m_configAndResults.crownThicknessQ95, CT_AbstractCategory::DATA_NUMBER, tr("crownThicknessQ95"));
    registerAttributeVaB(m_configAndResults.crownEquivRadius, CT_AbstractCategory::DATA_NUMBER, tr("crownEquivRadius"));
    registerAttributeVaB(m_configAndResults.volume, CT_AbstractCategory::DATA_NUMBER, tr("volume"));
    registerAttributeVaB(m_configAndResults.volume_topQ25, CT_AbstractCategory::DATA_NUMBER, tr("volume_topQ25"));
    registerAttributeVaB(m_configAndResults.volume_topQ50, CT_AbstractCategory::DATA_NUMBER, tr("volume_topQ50"));
    registerAttributeVaB(m_configAndResults.volume_topQ75, CT_AbstractCategory::DATA_NUMBER, tr("volume_topQ75"));
    registerAttributeVaB(m_configAndResults.slope_max, CT_AbstractCategory::DATA_NUMBER, tr("slope_max"));
    registerAttributeVaB(m_configAndResults.slope_min, CT_AbstractCategory::DATA_NUMBER, tr("slope_min"));
    registerAttributeVaB(m_configAndResults.slope_moy, CT_AbstractCategory::DATA_NUMBER, tr("slope_moy"));
    registerAttributeVaB(m_configAndResults.slope_sd, CT_AbstractCategory::DATA_NUMBER, tr("slope_sd"));
    registerAttributeVaB(m_configAndResults.slope_Q25, CT_AbstractCategory::DATA_NUMBER, tr("slope_Q25"));
    registerAttributeVaB(m_configAndResults.slope_Q50, CT_AbstractCategory::DATA_NUMBER, tr("slope_Q50"));
    registerAttributeVaB(m_configAndResults.slope_Q75, CT_AbstractCategory::DATA_NUMBER, tr("slope_Q75"));
}

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