/**************************************************************************** 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 . *****************************************************************************/ #include "onf_metricrastercrown.h" #include "ct_shape2ddata/ct_polygon2ddata.h" #include ONF_MetricRasterCrown::ONF_MetricRasterCrown(QString pluginName) : CT_AbstractMetric_Raster(pluginName) { declareAttributes(); } ONF_MetricRasterCrown::ONF_MetricRasterCrown(const ONF_MetricRasterCrown &other) : CT_AbstractMetric_Raster(other) { declareAttributes(); m_configAndResults = other.m_configAndResults; } QString ONF_MetricRasterCrown::getShortDisplayableName() const { return tr("Métriques de houppiers (Alti/Ht)"); } QString ONF_MetricRasterCrown::getShortDescription() const { return tr("Ces métriques décrivent la taille et la forme globale d'un houppier.
" "Elles sont conçues pour être calculées sur un modèle numérique de surface ou de hauteur l'échelle de l'arbre.
" "Elles peuvent fonctionner indifféremment en altitude ou en hauteur.
" "Attention cependant, en zone de forte pente la soustraction de l'altitude du sol pour passer à un modèle numérique de hauteur, conduit à une déformation du houppier. En cas de fort relief, il est donc conseillé de calculer ces métriques avec des nuages de points par houppier codés en altitude."); } QString ONF_MetricRasterCrown::getDetailledDescription() const { return tr("Les valeurs suivantes peuvent être calculées :
" "" "
" "Pour les métriques de pentes, un raster de pente est calculé à partir du modèles numérique de surface / de hauteur du houppier. Les statistiques suivantes sont calculées sur les valeurs de pentes des pixels du houppier.
" "" ); } 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.slope_max.value = 0; m_configAndResults.slope_min.value = std::numeric_limits::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.crownEquivRadius.value = 0; m_configAndResults.convexArea.value = 0; m_configAndResults.maxDistApex.value = 0; m_configAndResults.maxDist.value = 0; m_configAndResults.eqR_MxDst.value = 0; m_configAndResults.eqR_CrTh.value = 0; double n = 0.0; double sumx2 = 0.0; double inNA = _inRaster->NAAsDouble(); QList allPoints; Eigen::Vector3d center = Eigen::Vector3d(0.0, 0.0, 0.0); Eigen::Vector3d apex; apex(0) = -std::numeric_limits::max(); apex(1) = -std::numeric_limits::max(); apex(2) = -std::numeric_limits::max(); QList zvals; for (size_t index = 0 ; index < _inRaster->nCells() ; index++) { double val = _inRaster->valueAtIndexAsDouble(index); if (val != inNA) { zvals.append(val); _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; } } } std::sort(zvals.begin(), zvals.end()); m_configAndResults.crownThickness.value = zvals.last() - zvals.first(); QList slopeValues; for (size_t index = 0 ; index < _inRaster->nCells() ; index++) { double val = _inRaster->valueAtIndexAsDouble(index); int xx, yy; if (_inRaster->indexToGrid(index, xx, yy) && val != inNA) { // 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(); } } std::sort(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); } // ConvexArea // tri par (X,Y) de la liste des points CT_Polygon2DData::orderPointsByXY(allPoints); CT_Polygon2DData *data = CT_Polygon2DData::createConvexHull(allPoints); if (data != nullptr) { m_configAndResults.convexArea.value = data->getArea(); const auto& vertices = data->getVertices(); double maxDistApex = 0; double maxDist = 0; for (int i = 0 ; i < vertices.size() ; i++) { const auto& 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++) { const auto& vert2 = vertices.at(j); double dist = sqrt(pow((vert1)(0) - (vert2)(0), 2) + pow((vert1)(1) - (vert2)(1), 2)); if (dist > maxDist) {maxDist = dist;} } } m_configAndResults.maxDistApex.value = maxDistApex; m_configAndResults.maxDist.value = maxDist; delete data; } if (m_configAndResults.maxDist.value > 0) { m_configAndResults.eqR_MxDst.value = m_configAndResults.crownEquivRadius.value / m_configAndResults.maxDist.value; } if (m_configAndResults.crownThickness.value > 0) { m_configAndResults.eqR_CrTh.value = m_configAndResults.crownEquivRadius.value / m_configAndResults.crownThickness.value; } setAttributeValueVaB(m_configAndResults.rumple); setAttributeValueVaB(m_configAndResults.area3d); setAttributeValueVaB(m_configAndResults.area2d); setAttributeValueVaB(m_configAndResults.crownThickness); setAttributeValueVaB(m_configAndResults.crownEquivRadius); setAttributeValueVaB(m_configAndResults.convexArea); setAttributeValueVaB(m_configAndResults.maxDistApex); setAttributeValueVaB(m_configAndResults.maxDist); setAttributeValueVaB(m_configAndResults.eqR_MxDst); setAttributeValueVaB(m_configAndResults.eqR_CrTh); 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); qDeleteAll(allPoints); } void ONF_MetricRasterCrown::declareAttributes() { registerAttributeVaB(m_configAndResults.rumple, CT_AbstractCategory::DATA_NUMBER, "Rumple"); registerAttributeVaB(m_configAndResults.area3d, CT_AbstractCategory::DATA_NUMBER, "Area3d"); registerAttributeVaB(m_configAndResults.area2d, CT_AbstractCategory::DATA_NUMBER, "Area2d"); registerAttributeVaB(m_configAndResults.crownThickness, CT_AbstractCategory::DATA_NUMBER, "CrThick"); registerAttributeVaB(m_configAndResults.crownEquivRadius, CT_AbstractCategory::DATA_NUMBER, "CrEqRad"); registerAttributeVaB(m_configAndResults.convexArea, CT_AbstractCategory::DATA_NUMBER, "CvxArea"); registerAttributeVaB(m_configAndResults.maxDistApex, CT_AbstractCategory::DATA_NUMBER, "MaxDiApx"); registerAttributeVaB(m_configAndResults.maxDist, CT_AbstractCategory::DATA_NUMBER, "MaxDist"); registerAttributeVaB(m_configAndResults.eqR_MxDst, CT_AbstractCategory::DATA_NUMBER, "eqR_MxDst"); registerAttributeVaB(m_configAndResults.eqR_CrTh, CT_AbstractCategory::DATA_NUMBER, "eqR_CrTh"); registerAttributeVaB(m_configAndResults.slope_max, CT_AbstractCategory::DATA_NUMBER, "Slope_max"); registerAttributeVaB(m_configAndResults.slope_min, CT_AbstractCategory::DATA_NUMBER, "Slope_min"); registerAttributeVaB(m_configAndResults.slope_moy, CT_AbstractCategory::DATA_NUMBER, "Slope_moy"); registerAttributeVaB(m_configAndResults.slope_sd, CT_AbstractCategory::DATA_NUMBER, "Slope_sd"); registerAttributeVaB(m_configAndResults.slope_Q25, CT_AbstractCategory::DATA_NUMBER, "Slope_Q25"); registerAttributeVaB(m_configAndResults.slope_Q50, CT_AbstractCategory::DATA_NUMBER, "Slope_Q50"); registerAttributeVaB(m_configAndResults.slope_Q75, CT_AbstractCategory::DATA_NUMBER, "Slope_Q75"); } double ONF_MetricRasterCrown::computePercentile(const QList &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])); }