
//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2016 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#include "MeshProcessTriangleData.hpp"

#include "Progress.hpp"
#include <QFileDialog>

namespace mgx {
  
  //  bool ScaleTriangleValues::run(Mesh *mesh, const QString& scaleMap, float maxValueUsr, float minValueUsr)
  //	{
  //	// Scale with options: manual, auto (percentile 5-95%) or fullRange
  //  // see drawCFMHeatMap
  //    return true;
  //	}
  
  bool StatsTriangleValues::initialize()
  {
    float prc = parm("Percentile").toFloat();
    if(prc < 0)
      setParm("Percentile", QString("0"));
    if(prc > 100)
      setParm("Percentile", QString("100"));
    return true;
  }

  // Compute stats of triangle values on a cell basis (e.g. mean, median, percentile, std etc) and display in heat map
  // For the moment only compute average of triangle values, weighted by triangle area
  bool StatsTriangleValues::run(Mesh *mesh, bool selection, QString stat, float prc)
  {
    if(mesh->meshType() != "MGX3D")
      throw(QString("Mesh type is (%1), it must be converted to (MGX3D)").arg(mesh->meshType()));

    mesh->labelHeat().clear();
    IntFloatAttr& labelHeatMap = mesh->labelHeat();
    TriFloatAttr& triangleValue = mesh->triangleValue();

    // compute average or percentile of triangle values per label
    forall(const cell &c, mesh->cells()) {
      std::vector<Triangle> triList;
      std::vector<double> valueList;
      double sumValues = 0, sumAreas = 0;
      // find all triangles within the cell "c"
      forall(const vertex &v, c->S)
          forall(const vertex &n, c->S.neighbors(v)) {
        vertex m = c->S.nextTo(v, n);
        // check that triangle exist
        if(!c->S.uniqueTri(v, n, m))
          continue;
        // check that triangle selected (if relevant)
        if(selection)
          if(!v->selected or !n->selected or !m->selected)
            continue;
        // check that value stored for this triangle
        TriFloatAttr::iterator it = triangleValue.find(Triangle(v,n,m));
        if(it == triangleValue.end())
          continue;
        float value = it->second;
        float area = triangleArea(v->pos, n->pos, m->pos);

        // for percentile
        triList.push_back(Triangle(v,n,m));
        valueList.push_back(value);

        // for weighted average
        sumValues += value*area;
        sumAreas += area;
      }

      if(valueList.size() == 0)
        continue;

      if(stat == "Percentile"){
        std::sort(valueList.begin(), valueList.end());
        // position of percentile in the list of values
        uint posPrc = 0;
        if (valueList.size() > 0) {
          if (prc > 100) {
            prc = 100;
          }
          uint numPrc = valueList.size()*(prc/100.0);
          posPrc = numPrc - 1;
        }
        float prcValue = valueList[posPrc];
        labelHeatMap[c->label] = prcValue;
        Information::out << "label " << c->label << " percentile " << prcValue << " list size " << valueList.size() << " first value " << valueList[0] << " last value " <<  valueList[valueList.size()-1] << endl;
        Information::out << "label " << c->label << " position_prct " << posPrc << endl;
      }
      else if(stat == "WeightedAverage")
        labelHeatMap[c->label] = sumValues/sumAreas;
      Information::out << "label " << c->label << " sumValues " << sumValues << " sumAreas " << sumAreas << " sumValues/sumAreas " << sumValues/sumAreas << endl;
    }

    // find range of heat map
    float maxValue, minValue;
    bool firstRun = true;
    for(IntFloatAttr::const_iterator h = labelHeatMap.begin(); h != labelHeatMap.end(); h++){
      float value = h->second;
      if(firstRun) {
        maxValue = value;
        minValue = value;
        firstRun = false;
      }
      if(value > maxValue)
        maxValue = value;
      if(value < minValue)
        minValue = value;
    }
    mesh->heatMapBounds() = Point2f(minValue, maxValue);
    mesh->heatMapUnit() = mesh->triangleValueUnit();
    mesh->setShowLabel("Label Heat");
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(StatsTriangleValues);
}
