//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 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 "MeshProcessHeatMap.hpp"
#include <MeshProcessCellMesh.hpp>
#include <Triangulate.hpp>
#include <ThirdParty/triangle.h>
#include <GraphUtils.hpp>
#include <limits>
#include "Progress.hpp"
#include <QFileDialog>
#include <QMessageBox>

#include "MeshProcessMeasures.hpp"
#include "MeshProcessLineage.hpp"
#include "MeshProcessCellAxis.hpp"
#include "MeshProcessPDG.hpp"

#include <cmath>

namespace mgx
{

  void scaleHeatmap(Mesh* mesh, QString type, double parameter)
  {
    IntFloatAttr& labelHeat = mesh->labelHeat();
    if(labelHeat.empty() or type == "linear")
      return;


    std::map<float, std::set<int> > valueLabelMap;

    if(type == "ordinal"){
      //sort the data
      for(IntFloatAttr::iterator it = labelHeat.begin(); it != labelHeat.end(); ++it) {
        valueLabelMap[it->second].insert(it->first);
      }
      int counter = 1;
      for(std::map<float, std::set<int> >::iterator it = valueLabelMap.begin(); it != valueLabelMap.end(); ++it) {
        forall(int cellLabel, it->second){
          labelHeat[cellLabel] = counter;
        }
        counter++;
      }

    } else {


      for(IntFloatAttr::iterator it = labelHeat.begin(); it != labelHeat.end(); ++it) {
        float& v = it->second;
        if(type == "log"){
          v = log(v)/log(10);
        } else if(type == "square"){
          v = v*v;
        } else if(type == "sqrt"){
          v = sqrt(v);
        } else if(type == "inverse"){
          v = v==0? 0 : 1./v;
        } else if(type == "revert log"){
          v = pow(10,v);
        } else if(type == "pow"){
          v = pow(v,parameter);
        } else if(type == "add"){
          v += parameter;
        } else if(type == "mult"){
          v *= parameter;
        }
      }

    }

    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->setShowLabel("Label Heat");
    mesh->heatMapScaling() = type;

    mesh->updateTriangles();

  }

  bool ScaleHeatMap::run(Mesh* mesh, QString type, double parameter)
  {
    scaleHeatmap(mesh, type, parameter);
    return true;
  }
  REGISTER_PROCESS(ScaleHeatMap);


  bool NormalizeHeatMap::run(Mesh* mesh, QString type)
  {

    double minValue = HUGE_VAL;
    double maxValue = -HUGE_VAL;

    if(type == "Max Min"){



      forall(auto &p, mesh->labelHeat()){
        if(p.second < minValue) minValue = p.second;
        if(p.second > maxValue) maxValue = p.second;
      }

      forall(auto &p, mesh->labelHeat()){
        p.second = (p.second - minValue) / (maxValue - minValue);
      }


    } else { // Selected Cells

      // find min and max of selected cells
      forall(vertex v, mesh->graph()){
        // TODO
      }


    }

    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->setShowLabel("Label Heat");

    mesh->updateTriangles();


    return true;
  }
  REGISTER_PROCESS(NormalizeHeatMap);

  bool ComputeHeatMap::initialize(QWidget* parent)
  {
    if(!parent) 
      return true;

    dlg = new QDialog(parent);
    Ui_HeatMapDialog ui;
    ui.setupUi(dlg);

    connect(ui.HeatMapType, SIGNAL(currentIndexChanged(const QString &)), this, SLOT(changeMapType(const QString &)));
    connect(ui.SelectSpreadsheetFile, SIGNAL(clicked()), this, SLOT(selectSpreadsheetFile()));

    for(int i = 0; i < ui.HeatMapType->count(); ++i) {
      if(ui.HeatMapType->itemText(i) == parm("Type")) {
        ui.HeatMapType->setCurrentIndex(i);
        break;
      }
    }
    for(int i = 0; i < ui.HeatMapSignal->count(); ++i) {
      if(ui.HeatMapSignal->itemText(i) == parm("Visualize")) {
        ui.HeatMapSignal->setCurrentIndex(i);
        break;
      }
    }
    ui.SpreadsheetFile->setText(parm("FileName"));
    if(parm("ReportFields").isEmpty()) {
      ui.CreateSpreadsheet->setChecked(false);
    } else {
      QStringList spread = parm("ReportFields").split(",");
      ui.HeatMapRepGeom->setChecked(spread.contains("Geometry"));
      ui.HeatMapRepSig->setChecked(spread.contains("Signal"));
      ui.HeatMapRepBord->setChecked(spread.contains("Border"));
    }
    ui.HeatMapRange->setChecked(stringToBool(parm("Man. Range")));
    ui.RangeMin->setValue(parm("Range Min").toFloat());
    ui.RangeMax->setValue(parm("Range Max").toFloat());
    ui.HeatMapAvg->setChecked(stringToBool(parm("Signal Avg")));
    ui.GlobalCoordinates->setChecked(stringToBool(parm("Global Coord")));
    ui.HeatMapPolarity->setChecked(parm("Polarity Type") != "None");
    if(ui.HeatMapPolarity->isChecked()) {
      for(int i = 0; i < ui.HeatMapPolarityType->count(); ++i) {
        if(ui.HeatMapPolarityType->itemText(i) == parm("Polarity Type")) {
          ui.HeatMapPolarityType->setCurrentIndex(i);
          break;
        }
      }
    }

    ui.HeatMapChange->setChecked(stringToBool(parm("Change Map")));
    if(parm("Increasing") == "Increasing")
      ui.HeatMapIncr->setCurrentIndex(0);
    else
      ui.HeatMapIncr->setCurrentIndex(1);
    if(parm("Diff Type") == "Ratio")
      ui.HeatMapDiffType->setCurrentIndex(0);
    else if(parm("Diff Type") == "Difference")
      ui.HeatMapDiffType->setCurrentIndex(1);
    else if(parm("Diff Type") == "Growth")
      ui.HeatMapDiffType->setCurrentIndex(2);

    ui.HeatMapGrowthTime->setValue(parm("Growth Time").toFloat());
    borderSize = parm("Border Size(µm)").toFloat();

    bool go = dlg->exec() == QDialog::Accepted;
    if(go) {
      QStringList spread;
      if(ui.CreateSpreadsheet->isChecked()) {
        if(ui.HeatMapRepGeom->isChecked())
          spread << "Geometry";
        if(ui.HeatMapRepSig->isChecked())
          spread << "Signal";
        if(ui.HeatMapRepBord->isChecked())
          spread << "Border";
      }
       setParm("Type" , ui.HeatMapType->currentText());
       setParm("Visualize" , ui.HeatMapSignal->currentText());
       setParm("FileName" , ui.SpreadsheetFile->text());
       setParm("ReportFields" , spread.join(","));
       setParm( "Man. Range" , (ui.HeatMapRange->isChecked() ? "Yes" : "No"));
       setParm("Range Min" , QString::number(ui.RangeMin->value())) ;
       setParm("Range Max" , QString::number(ui.RangeMax->value()));
       setParm("Signal Avg" , (ui.HeatMapAvg->isChecked() ? "Yes" : "No"));
       setParm("Global Coord" , (ui.GlobalCoordinates->isChecked() ? "Yes" : "No"));
       setParm("Polarity Type" , (ui.HeatMapPolarity->isChecked() ? ui.HeatMapPolarityType->currentText() : QString("None")));
       setParm("Change Map" , (ui.HeatMapChange->isChecked() ? "Yes" : "No"));
       setParm("Increasing" , ui.HeatMapIncr->currentText());
       setParm("Diff Type" , ui.HeatMapDiffType->currentText());
       setParm("Growth Time", QString::number(ui.HeatMapGrowthTime->value()));
       setParm("Border Size(µm)", QString::number(borderSize));
    }
    delete dlg;
    dlg = 0;
    return go;
  }

  bool ComputeHeatMap::run()
  {
    MapType map;
    if(parm("Type") == "Area")
      map = AREA;
    else if(parm("Type") == "Volume")
      map = VOLUME;
    else if(parm("Type") == "Walls")
      map = WALL;
    else {
      setErrorMessage("Error, first string must be one of 'Area', 'Volume' or 'Walls'");
      return false;
    }
    SignalType signal;
    if(parm("Visualize") == "Geometry")
      signal = GEOMETRY;
    else if(parm("Visualize") == "Border signal")
      signal = BORDER;
    else if(parm("Visualize") == "Interior signal")
      signal = INTERIOR;
    else if(parm("Visualize") == "Border/total")
      signal = BORDER_TOTAL;
    else if(parm("Visualize") == "Interior/total")
      signal = INTERIOR_TOTAL;
    else if(parm("Visualize") == "Total signal")
      signal = TOTAL;
    else {
      setErrorMessage("Error, second string must be one of 'Geometry', 'Border signal',"
                      " 'Interior signal', 'Border/total', 'Interior/total' or 'Total signal'");
      return false;
    }
    QString reportFile = parm("FileName");
    if(!reportFile.isEmpty() and !reportFile.endsWith(".csv", Qt::CaseInsensitive))
      reportFile += ".csv";
    QStringList spread = parm("ReportFields").split(",");
    int report = 0;
    if(spread.contains("Geometry"))
      report |= GEOMETRY;
    if(spread.contains("Signal"))
      report |= TOTAL;
    if(spread.contains("Border"))
      report |= BORDER;
    bool manualRange = stringToBool(parm("Man. Range"));
    float rangeMin = parm("Range Min").toFloat();
    float rangeMax = parm("Range Max").toFloat();
    if(manualRange and rangeMin >= rangeMax) {
      setErrorMessage("If range is manual, then rangeMin (2nd value) must be less than rangeMax (3rd value)");
      return false;
    }
    bool signalAverage = stringToBool(parm("Signal Avg"));
    bool globalCoordinates = stringToBool(parm("Global Coord"));
    PolarityType polarity = NONE;
    if(parm("Polarity Type") == "None")
      polarity = NONE;
    else if(parm("Polarity Type") == "Cell Average")
      polarity = CELL_AVERAGE;
    else if(parm("Polarity Type") == "Wall/Min")
      polarity = WALL_MIN;
    else {
      setErrorMessage("Error, 9th string must be one of 'None', 'Cell Average or 'Wall/Min'");
      return false;
    }

    MultiMapType multiMapType;
    float growthTime = parm("Growth Time").toFloat();
    bool hasChange = stringToBool(parm("Change Map"));
    if(!hasChange)
      multiMapType = SINGLE;
    else {
      bool incr = parm("Increasing") == "Increasing";
      bool ratio = parm("Diff Type") == "Ratio";
      bool diff = parm("Diff Type") == "Difference";
      bool growth = parm("Diff Type") == "Growth";
      if(growth and growthTime <= 0) {
        setErrorMessage("Time for growth rate cannot be <= 0)");
        return false;
      }
      if(incr) {
        if(ratio)
          multiMapType = INCREASING_RATIO;
        else if(diff)
          multiMapType = INCREASING_DIFF;
        else
          multiMapType = INCREASING_GROWTH;
      } else {
        if(ratio)
          multiMapType = DECREASING_RATIO;
        else if(diff)
          multiMapType = DECREASING_DIFF;
        else
          multiMapType = DECREASING_GROWTH;
      }
    }
    float borderSize = parm("Border Size(µm)").toFloat();

    // Check the mesh state
    Mesh* mesh1 = currentMesh(), *mesh2 = 0;
    if(!checkState().mesh(MESH_NON_EMPTY, mesh1->id())) {
      setErrorMessage("Current mesh is empty");
      return false;
    }

    if(multiMapType != SINGLE) {
      if(currentMesh() == mesh(0))
        mesh2 = mesh(1);
      else if(currentMesh() == mesh(1))
        mesh2 = mesh(0);
      else
        return false;

      if(!checkState().mesh(MESH_NON_EMPTY, mesh2->id())) {
        setErrorMessage("Second mesh is empty");
        return false;
      }
    }

    bool res = run(mesh1, mesh2, map, signal, reportFile, report, manualRange, rangeMin, rangeMax, signalAverage,
                   globalCoordinates, polarity, multiMapType, growthTime, borderSize);
    return res;
  }

  bool ComputeHeatMap::run(Mesh* mesh1, Mesh* mesh2, MapType _map, SignalType _signal, const QString& reportFile,
                           int _report, bool _manualRange, float _rangeMin, float _rangeMax, bool _signalAverage,
                           bool _globalCoordinates, PolarityType _polarity, MultiMapType _multiMapType,
                           float _growthTime, float _borderSize)
  {
    mapType = _map;
    signal = _signal;
    report = _report;
    manualRange = _manualRange;
    rangeMin = _rangeMin;
    rangeMax = _rangeMax;
    signalAverage = _signalAverage;
    globalCoordinates = _globalCoordinates;
    polarity = _polarity;
    multiMapType = _multiMapType;
    growthTime = _growthTime;
    borderSize = _borderSize;

    bool multiMap = multiMapType != SINGLE;
    if(multiMap and (!mesh2 or mesh2->graph().empty())) {
      setErrorMessage("Error, to compute a heat map based on two maps, you must specify two meshes");
      return false;
    }

    // Check if the report/signal can actually be computed

    if(mapType == VOLUME) {
      int anysignal = (BORDER | INTERIOR | TOTAL | BORDER_TOTAL | INTERIOR_TOTAL);
      if((mesh1 and (mesh1->stack()->empty() or mesh1->stack()->currentStore()->labels()))
         or (mesh2 and (mesh2->stack()->empty() or mesh2->stack()->currentStore()->labels()))) {
        if((signal & anysignal)or (report & anysignal))
          return setErrorMessage("Error, trying to measure 3D signal without a stack being loaded");
      }
    }

    QString sigStr = "";
    if(signalAverage and signal != GEOMETRY)
      sigStr += "Average ";

    // Heat map signal type
    switch(signal) {
    case GEOMETRY:
      switch(mapType) {
      case AREA:
        sigStr += "Area";
        break;
      case VOLUME:
        sigStr += "Volume";
        break;
      case WALL:
        sigStr += "Wall";
      }
      break;
    case BORDER:
      sigStr += "Border Signal";
      break;
    case INTERIOR:
      sigStr += "Interior Signal";
      break;
    case BORDER_TOTAL:
      sigStr += "Border/Total Signal";
      break;
    case INTERIOR_TOTAL:
      sigStr += "Interior/Total Signal";
      break;
    case TOTAL:
      sigStr += "Total Signal";
    }

    ImageHeatMaps map1, map2;
    if(!getLabelMaps(map1, mesh1))
      return false;

    IntFloatAttr& LabelHeat = mesh1->labelHeat();
    IntIntFloatAttr& WallHeat = mesh1->wallHeat();

    // If we are doing a change map, get it from the other mesh and take the ratio as the signal
    // Otherwise just copy the info
    LabelHeat.clear();
    WallHeat.clear();
    if(!multiMap) {
      if(mapType == AREA or mapType == VOLUME)
        LabelHeat = map1.LabelHeatAmt;
      else if(mapType == WALL)
        WallHeat = map1.WallHeatAmt;
    } else {
      switch(multiMapType) {
      case INCREASING_RATIO:
        sigStr += " increase ratio";
        break;
      case INCREASING_DIFF:
        sigStr += " increase difference";
        break;
      case INCREASING_GROWTH:
        sigStr += QString(" increase growth(%1)").arg(growthTime);
        break;
      case DECREASING_RATIO:
        sigStr += " decrease ratio";
        break;
      case DECREASING_DIFF:
        sigStr += " decrease difference";
        break;
      case DECREASING_GROWTH:
        sigStr += QString(" decrease growth(%1)").arg(growthTime);
        break;
      case SINGLE:
        break;
      }

      // Get values from other mesh
      if(!getLabelMaps(map2, mesh2))
        return false;

      // Change map for area and volume
      switch(mapType) {
      case AREA:
      case VOLUME: {
        std::vector<int> LabelDel;
        forall(const IntFloatPair& p, map1.LabelHeatAmt) {
          int i = p.first;
          // Calculate change
          if(map2.LabelHeatAmt.find(i) != map2.LabelHeatAmt.end()) {
            switch(multiMapType) {
            case INCREASING_RATIO:
              LabelHeat[i] = (map1.LabelHeatAmt[i] == 0 ? 0 : map2.LabelHeatAmt[i] / map1.LabelHeatAmt[i]);
              break;
            case DECREASING_RATIO:
              LabelHeat[i] = (map2.LabelHeatAmt[i] == 0 ? 0 : map1.LabelHeatAmt[i] / map2.LabelHeatAmt[i]);
              break;
            case INCREASING_DIFF:
              LabelHeat[i] = map2.LabelHeatAmt[i] - map1.LabelHeatAmt[i];
              break;
            case DECREASING_DIFF:
              LabelHeat[i] = map1.LabelHeatAmt[i] - map2.LabelHeatAmt[i];
              break;
            case INCREASING_GROWTH:
              if(growthTime == 0 or map1.LabelHeatAmt[i] == 0)
                LabelHeat[i] = 0;
              else
                LabelHeat[i] = (map2.LabelHeatAmt[i] / map1.LabelHeatAmt[i] - 1) / growthTime;
              break;
            case DECREASING_GROWTH:
              if(growthTime == 0 or map2.LabelHeatAmt[i] == 0)
                LabelHeat[i] = 0;
              else
                LabelHeat[i] = (map1.LabelHeatAmt[i] / map2.LabelHeatAmt[i] - 1) / growthTime;
              break;
            case SINGLE:
              break;
            }
          } else
            LabelDel.push_back(i);
        }
        // Delete labels not in both meshes
        forall(int i, LabelDel) {
          LabelHeat.erase(i);
          SETSTATUS("MakeHeatMap Warning mesh " << mesh1->userId() << " erasing label " << i);
        }
      } break;
      case WALL: {     // Calculate wall change map
        std::vector<IntIntPair> WallDel;
        forall(const IntIntFloatPair& p, map1.WallHeatAmt) {
          IntIntPair i = p.first;
          // Calculate change
          if(map2.WallHeatAmt.find(i) != map2.WallHeatAmt.end()) {
            switch(multiMapType) {
            case INCREASING_RATIO:
              WallHeat[i] = (map1.WallHeatAmt[i] == 0 ? 0 : map2.WallHeatAmt[i] / map1.WallHeatAmt[i]);
              break;
            case DECREASING_RATIO:
              WallHeat[i] = (map2.WallHeatAmt[i] == 0 ? 0 : map1.WallHeatAmt[i] / map2.WallHeatAmt[i]);
              break;
            case INCREASING_DIFF:
              WallHeat[i] = map2.WallHeatAmt[i] - map1.WallHeatAmt[i];
              break;
            case DECREASING_DIFF:
              WallHeat[i] = map1.WallHeatAmt[i] - map2.WallHeatAmt[i];
            case INCREASING_GROWTH:
              if(growthTime == 0 or map1.WallHeatAmt[i] == 0)
                WallHeat[i] = 0;
              else
                WallHeat[i] = (map2.WallHeatAmt[i] / map1.WallHeatAmt[i] - 1) / growthTime;
              break;
            case DECREASING_GROWTH:
              if(growthTime == 0 or map2.WallHeatAmt[i] == 0)
                WallHeat[i] = 0;
              else
                WallHeat[i] = (map1.WallHeatAmt[i] / map2.WallHeatAmt[i] - 1) / growthTime;
              break;
            case SINGLE:
              break;
            }
          } else
            WallDel.push_back(i);
        }
        // Delete labels not in both meshes
        forall(IntIntPair i, WallDel) {
          WallHeat.erase(i);
          SETSTATUS("MakeHeatMap::Warning:Mesh" << mesh1->userId() << " erasing wall " << i.first << "->"
                    << i.second);
        }
      }
      }
    }

    // Calculate ranges if not using user defined range
    float heatMinAmt = rangeMin;
    float heatMaxAmt = rangeMax;
    if(!manualRange) {
      heatMinAmt = HUGE_VAL, heatMaxAmt = -HUGE_VAL;

      // Find max and min
      switch(mapType) {
      case AREA:
      case VOLUME:
        forall(const IntFloatPair& p, LabelHeat) {
          float amt = LabelHeat[p.first];
          if(amt < heatMinAmt)
            heatMinAmt = amt;
          if(amt > heatMaxAmt)
            heatMaxAmt = amt;
        }
        break;
      case WALL:
        forall(const IntIntFloatPair& p, WallHeat) {
          float amt = WallHeat[p.first];
          if(amt < heatMinAmt)
            heatMinAmt = amt;
          if(amt > heatMaxAmt)
            heatMaxAmt = amt;
        }
      }
    }
    mesh1->heatMapBounds() = Point2f(heatMinAmt, heatMaxAmt);

    // Find units
    QString gunits(""), hunits("");
    switch(mapType) {
    case WALL:
      gunits = UM;
      break;
    case AREA:
      gunits = UM2;
      break;
    case VOLUME:
      gunits = UM3;
      break;
    }
    if(!multiMap and signal == GEOMETRY)
      hunits = gunits;

    mesh1->heatMapUnit() = hunits;

    // Write to file
    if(report != 0) {
      QFile file(reportFile);
      if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) {
        SETSTATUS("makeHeatMap::Error:Cannot open output file:" << reportFile);
      } else {
        QTextStream out(&file);

        // Area heat map
        switch(mapType) {
        case AREA:
        case VOLUME: {
          // Write header
          out << "Label,Value";
          if(mapType == AREA) {
            out << QString(",Area (%1)").arg(UM2);
            if(report & BORDER) {
              out << QString(",Border Area (%1),Interior Area (%1)").arg(UM2);
            }
          } else {
            out << QString(",Volume (%1)").arg(UM3);
            if(report & BORDER) {
              out << QString(",Border Volume (%1),Interior Volume (%1)").arg(UM3);
            }
          }
          if((report & BORDER)or (report & TOTAL)) {
            out << ",Total Signal";
            if(report & BORDER) {
              out << ",Border Signal,Interior Signal";
            }
          }
          out << ",Center_X,Center_Y,Center_Z";
          out << endl;

          // Write data
          IntPoint3fAttr& LabelCenter = mesh1->labelCenter();
          const qglviewer::Frame& frame = mesh1->stack()->getFrame();
          forall(const IntFloatPair& p, LabelHeat) {
            int label = p.first;
            out << label << "," << LabelHeat[label];
            out << "," << map1.LabelGeom[label];
            if(report & BORDER) {
              out << "," << map1.LabelBordGeom[label];
              out << "," << map1.LabelIntGeom[label];
            }
            if((report & SINGLE)or (report & TOTAL)) {
              out << "," << map1.LabelTotalSig[label];
              if(report & BORDER) {
                out << "," << map1.LabelBordSig[label];
                out << "," << map1.LabelIntSig[label];
              }
            }
            // Write cell centers
            if(globalCoordinates) {
              Point3f tr = Point3f(frame.inverseCoordinatesOf(qglviewer::Vec(LabelCenter[label])));
              out << "," << tr.x();
              out << "," << tr.y();
              out << "," << tr.z();
            } else {
              out << "," << LabelCenter[label].x();
              out << "," << LabelCenter[label].y();
              out << "," << LabelCenter[label].z();
            }

            out << endl;
          }
        } break;
        case WALL: {
          IntIntFloatAttr & WallGeom = mesh1->wallGeom();
          // Write header
          out << QString("Label,Neighbor,Value,Wall Length (%1)").arg(UM);
          if(signal == TOTAL or (report & TOTAL))
            out << QString(",Wall Border Area (%1),Wall Signal").arg(UM2);
          out << endl;

          // Write data
          forall(const IntIntFloatPair& p, WallHeat) {
            IntIntPair wall = p.first;
            out << wall.first << "," << wall.second << ",";
            out << WallHeat[wall] << "," << WallGeom[wall];
            if(signal == TOTAL or (report & TOTAL))
              out << "," << map1.WallBordGeom[wall] << "," << map1.WallBordSig[wall];
            out << endl;
          }
        }
        }
        // Footer same for all heat maps
        QString msg = QString("%1 HeatMap results for mesh %2").arg(sigStr).arg(mesh1->userId());
        out << endl << endl << qPrintable(msg) << endl;

        out << "Signal range: " << heatMinAmt << "-" << heatMaxAmt << " " << hunits << endl;
        file.close();
      }
    }

    // Print out stats
    int ncells = map1.LabelGeom.size();
    int acells = map1.LabelTotalSig.size();
    double totageom = 0, totgeom = 0;
    if(mapType == AREA or mapType == VOLUME) {
      forall(const IntFloatPair& p, map1.LabelGeom) {
        totgeom += map1.LabelGeom[p.first];
        if(map1.LabelTotalSig.find(p.first) != map1.LabelTotalSig.end())
          totageom += map1.LabelGeom[p.first];
      }
      SETSTATUS("Mesh " << mesh1->userId() << " " << sigStr << " analysis, Cells: " << acells << "/" << ncells << ", "
                << totageom << "/" << totgeom << " " << gunits);
    }
    Information::out << "Range min: " << heatMinAmt << " max: " << heatMaxAmt << " " << hunits << endl;

    if(mapType == WALL) {
      mesh1->setShowLabel("Wall Heat");
      mesh1->updateLines();
    } else
      mesh1->setShowLabel("Label Heat");

    mesh1->updateTriangles();
    return true;
  }

  void ComputeHeatMap::selectSpreadsheetFile()
  {
    QLineEdit* edit = dlg->findChild<QLineEdit*>("SpreadsheetFile");
    if(edit) {
      QString filename = edit->text();
      if(filename.isEmpty())
        filename = "spreadsheet.csv";
      filename = QFileDialog::getSaveFileName(dlg, "Select file to save report in", filename, "CSV files (*.csv)");
      if(!filename.isEmpty()) {
        if(!filename.endsWith(".csv", Qt::CaseInsensitive))
          filename += ".csv";
        edit->setText(filename);
      }
    } else
      SETSTATUS("SpreadsheetFile not found ...");
  }

  void ComputeHeatMap::changeMapType(const QString& type)
  {
    QCheckBox* pol = dlg->findChild<QCheckBox*>("HeatMapPolarity");
    if(pol)
      pol->setEnabled(type == "Walls");
    else
      SETSTATUS("HeatMapPolarity not found");
    QComboBox* signal = dlg->findChild<QComboBox*>("HeatMapSignal");
    if(signal) {
      QString curtext = signal->currentText();
      signal->clear();
      if(type == "Area") {
        signal->addItem("Geometry");
        signal->addItem("Border signal");
        signal->addItem("Interior signal");
        signal->addItem("Border/total");
        signal->addItem("Interior/total");
        signal->addItem("Total signal");
      } else if(type == "Volume") {
        signal->addItem("Geometry");
        signal->addItem("Border signal");
        signal->addItem("Interior signal");
        signal->addItem("Border/total");
        signal->addItem("Interior/total");
        signal->addItem("Total signal");
      } else {
        signal->addItem("Geometry");
        signal->addItem("Total signal");
      }
      for(int i = 0; i < signal->count(); ++i) {
        if(signal->itemText(i) == curtext) {
          signal->setCurrentIndex(i);
          break;
        }
      }
    } else
      SETSTATUS("HeatMapSignal not found");
    QCheckBox* border = dlg->findChild<QCheckBox*>("HeatMapRepBord");
    if(border)
      border->setEnabled(type == "Area" || type == "Volume");
    else
      SETSTATUS("HeatMapRepBord not found");
  }

  void getBBox(const HVec3F& pts, Point3f* bBox)
  {
    bBox[0] = Point3f(1e10, 1e10, 1e10);
    bBox[1] = -bBox[0];
    forall(const Point3f& pt, pts)
        for(int i = 0; i < 3; i++) {
      if(pt[i] < bBox[0][i])
        bBox[0][i] = pt[i];
      if(pt[i] > bBox[1][i])
        bBox[1][i] = pt[i];
    }
  }

  bool ComputeHeatMap::getLabelMaps(ImageHeatMaps& map, Mesh* mesh)
  {
    const Stack* stack = mesh->stack();
    int label;

    IntFloatAttr& LabelGeom = map.LabelGeom;
    IntFloatAttr& LabelBordGeom = map.LabelBordGeom;
    IntFloatAttr& LabelIntGeom = map.LabelIntGeom;
    IntFloatAttr& LabelTotalSig = map.LabelTotalSig;
    IntFloatAttr& LabelBordSig = map.LabelBordSig;
    IntFloatAttr& LabelIntSig = map.LabelIntSig;
    IntFloatAttr& LabelHeatAmt = map.LabelHeatAmt;
    IntHVec3uMap& LabelTris = map.LabelTris;
    IntVIdSetMap& LabelPts = map.LabelPts;

    IntIntFloatAttr& WallBordGeom = map.WallBordGeom;
    IntIntFloatAttr& WallBordSig = map.WallBordSig;
    IntIntFloatAttr& WallHeatAmt = map.WallHeatAmt;

    IntPoint3fAttr& LabelCenter = mesh->labelCenter();

    LabelGeom.clear();
    LabelBordGeom.clear();
    LabelIntGeom.clear();
    LabelTotalSig.clear();
    LabelBordSig.clear();
    LabelIntSig.clear();
    LabelTris.clear();
    LabelPts.clear();
    LabelHeatAmt.clear();
    LabelCenter.clear();

    WallBordGeom.clear();
    WallBordSig.clear();
    WallHeatAmt.clear();

    // Get the graph
    vvGraph& S = mesh->graph();

    // Get the parents
    IntIntAttr & parents = mesh->parents();

    // Find selected labels
    std::set<int> SelectedLabels;
    std::vector<vertex> av = mesh->activeVertices();
    bool DoSelected = false;
    if(av.size() != S.size()) {
      DoSelected = true;
      forall(const vertex v, av)
          forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;

        if(v->label > 0)
          SelectedLabels.insert(v->label);
      }
      SETSTATUS("Notice: Computing heatmap only for " << SelectedLabels.size() << " selected labels.");
    }

    // Find border vertices for wall maps or area maps that require it
    bool doborder = false;
    QList<SignalType> needBorder = QList<SignalType>() << BORDER << BORDER_TOTAL << INTERIOR << INTERIOR_TOTAL;
    if((mapType == AREA or mapType == VOLUME)and ((report & BORDER) or needBorder.contains(signal)))
      doborder = true;

    // Find all vertices within a limited distance from the walls
    if(doborder or mapType == WALL)
      mesh->updateWallGeometry(borderSize);

    if(mapType == AREA or mapType == VOLUME) {
      // For volumes, delete labels with open meshes
      std::set<int> LabelDel;

      // Loop over triangles
      forall(const vertex& v, S) {
        forall(const vertex& n, S.neighbors(v)) {
          vertex m = S.nextTo(v, n);
          if(!S.uniqueTri(v, n, m))
            continue;
          label = mesh->getLabel(v, n, m, parents);
          if(label <= 0)
            continue;
          if(DoSelected and SelectedLabels.count(label) == 0)
            continue;

          if(mapType == AREA) {
            // Total area, signal, and label center
            float area = triangleArea(v->pos, n->pos, m->pos);
            LabelGeom[label] += area;
            LabelCenter[label] += area * Point3f(v->pos + n->pos + m->pos) / 3.0;
            LabelTotalSig[label] += area * (v->signal + n->signal + m->signal) / 3.0;

            // Calculate how much of triangle is in border area
            if(doborder) {
              int border = 0;
              if(v->minb != 0)
                border++;
              if(n->minb != 0)
                border++;
              if(m->minb != 0)
                border++;
              // Border area and signal
              float bamt = float(border) / 3.0;
              LabelBordSig[label] += area * bamt * (v->signal + n->signal + m->signal) / 3.0;
              LabelBordGeom[label] += area * bamt;

              // Interior area and signal
              float iamt = float(3 - border) / 3.0;
              LabelIntSig[label] += area * iamt * (v->signal + n->signal + m->signal) / 3.0;
              LabelIntGeom[label] += area * iamt;
            }
          } else if(mapType == VOLUME) {
            // For volume just grab the triangle list
            // If any vertices of triangle are not the same as the label, then the label
            // does not define a closed volume, add it to deletion list
            // RSS: Fix this for parents
            if(v->label != label or n->label != label or m->label != label)
              LabelDel.insert(label);

            // Calculate volume enclosed by mesh and save triangles
            float volume = signedTetraVolume(v->pos, n->pos, m->pos);
            LabelGeom[label] += volume;
            LabelCenter[label] += volume * Point3f(v->pos + n->pos + m->pos) / 4.0;

            // Save triangles and points if processing volume signal
            if(signal != GEOMETRY or (report & TOTAL) or (report & BORDER)) {
              LabelTris[label].push_back(Point3u(v.id(), n.id(), m.id()));
              LabelPts[label].insert(v.id());
              LabelPts[label].insert(n.id());
              LabelPts[label].insert(m.id());
            }
          }
        }
      }

      // Do the volume heat map signal calculation
      if(mapType == VOLUME) {
        // Remove labels which define volume calculations on open meshes
        forall(int i, LabelDel) {
          LabelTotalSig.erase(i);
          SETSTATUS("MakeHeatMap::Warning: Mesh " << mesh->userId() << " open mesh for Vol calc, erasing label "
                    << i);
        }

        // Get mesh volume fluorescence
        if(signal != GEOMETRY or (report & TOTAL) or (report & BORDER)) {
          // Get current store
          const Store* store = mesh->stack()->currentStore();
          if(!store) {
            setErrorMessage("Volume signal calculation requires associated store");
            return false;
          }
          const HVecUS& data = store->data();

          // Clear previous geometries and use pixel based volumes, keep the centers
          LabelGeom.clear();
          LabelBordGeom.clear();
          LabelIntGeom.clear();

          // Loop over labels
          uint step = 0;
          progressStart(QString("Finding signal inside mesh %1").arg(mesh->userId()), LabelGeom.size());

          forall(const IntPoint3fPair& p, LabelCenter) {
            int label = p.first;

            // Process triangle and point lists
            IntIntMap vmap;
            HVec3F pts(LabelPts[label].size());
            int idx = 0;
            forall(const intptr_t u, LabelPts[label]) {
              vertex v(u);
              pts[idx] = Point3f(v->pos);
              vmap[u] = idx++;
            }
            // Calculate triangle normals.
            idx = 0;
            HVec3F nrmls(LabelTris[label].size());
            // Convert triangle index list to new vertex list
            forall(Point3u& tri, LabelTris[label]) {
              tri = Point3u(vmap[tri.x()], vmap[tri.y()], vmap[tri.z()]);
              if(tri.x() > pts.size() or tri.y() > pts.size() or tri.z() > pts.size()) {
                setErrorMessage("Error, bad triangle index");
                return false;
              }
              Point3f a = pts[tri.x()], b = pts[tri.y()], c = pts[tri.z()];
              nrmls[idx++] = ((b - a) ^ (c - a)).normalize();
            }

            // Get bounding box
            Point3f bBoxW[2];
            getBBox(pts, bBoxW);
            BoundingBox3i bBoxI(stack->worldToImagei(bBoxW[0]),
                stack->worldToImagei(bBoxW[1]) + Point3i(1, 1, 1));
            Point3i imgSize(stack->size());
            bBoxI &= BoundingBox3i(Point3i(0, 0, 0), imgSize);
            Point3i size = bBoxI.size();

            // call cuda
            if(!progressAdvance(step++))
              userCancel();
            HVecUS totMask(size_t(size.x()) * size.y() * size.z(), 0),
                bordMask(size.x() * size.y() * size.z(), 0);
            if(insideMeshGPU(bBoxI[0], size, stack->step(), stack->origin(), pts, LabelTris[label], nrmls,
                             totMask))
              return false;
            if(!progressAdvance(0))
              userCancel();

            if(doborder) {
              if(nearMeshGPU(bBoxI[0], size, stack->step(), stack->origin(), pts, nrmls, 0, borderSize,
                             bordMask))
                return false;
              if(!progressAdvance(0))
                userCancel();
            }

            // Count interior, border, and total pixels
            uint borderPix = 0, interiorPix = 0, totalPix = 0;

            // sum results
            for(int z = bBoxI[0].z(); z < bBoxI[1].z(); ++z)
              for(int y = bBoxI[0].y(); y < bBoxI[1].y(); ++y) {
                int xIdx = stack->offset(bBoxI[0].x(), y, z);
                for(int x = bBoxI[0].x(); x < bBoxI[1].x(); ++x, ++xIdx) {
                  size_t mIdx = ((size_t(z) - bBoxI[0].z()) * size.y() + y - bBoxI[0].y()) * size.x() + x
                      - bBoxI[0].x();
                  if(totMask[mIdx]) {
                    LabelTotalSig[label] += data[xIdx];
                    totalPix++;
                    if(doborder and bordMask[mIdx]) {
                      borderPix++;
                      LabelBordSig[label] += data[xIdx];
                    } else {
                      interiorPix++;
                      LabelIntSig[label] += data[xIdx];
                    }
                  }
                }
              }

            // Calculate volume and average signal per label
            float pixVolume = stack->step().x() * stack->step().y() * stack->step().z();

            LabelGeom[label] = totalPix * pixVolume;
            LabelBordGeom[label] = borderPix * pixVolume;
            LabelIntGeom[label] = interiorPix * pixVolume;

            LabelTotalSig[label] *= pixVolume;
            LabelBordSig[label] *= pixVolume;
            LabelIntSig[label] *= pixVolume;
          }
        }
      }

      // Fill in heatmap color
      forall(const IntFloatPair& p, LabelGeom) {
        int label = p.first;

        // Calculate label center
        LabelCenter[label] = (LabelGeom[label] == 0 ? Point3f(0, 0, 0) : LabelCenter[label] / LabelGeom[label]);

        float totalSignal, bordSignal = 0, intSignal = 0;

        // Set signal amounts, average over area if required
        totalSignal = LabelTotalSig[label];
        if(doborder) {
          bordSignal = LabelBordSig[label];
          intSignal = LabelIntSig[label];
        }
        if(signalAverage) {
          totalSignal = (LabelGeom[label] == 0 ? 0 : totalSignal / LabelGeom[label]);
          if(doborder) {
            bordSignal = (LabelBordGeom[label] == 0 ? 0 : bordSignal / LabelBordGeom[label]);
            intSignal = (LabelIntGeom[label] == 0 ? 0 : intSignal / LabelIntGeom[label]);
          }
        }

        // Fill in heat depending on user choice
        switch(signal) {
        case GEOMETRY:
          LabelHeatAmt[label] = LabelGeom[label];
          break;
        case BORDER:
          LabelHeatAmt[label] = bordSignal;
          break;
        case BORDER_TOTAL:
          LabelHeatAmt[label] = (totalSignal == 0 ? 0 : bordSignal / totalSignal);
          break;
        case INTERIOR:
          LabelHeatAmt[label] = intSignal;
          break;
        case INTERIOR_TOTAL:
          LabelHeatAmt[label] = (totalSignal == 0 ? 0 : intSignal / totalSignal);
          break;
        case TOTAL:
          LabelHeatAmt[label] = totalSignal;
          break;
        }
      }
    } else if(mapType == WALL) {   // Do Wall heatmap
      IntIntFloatAttr& WallGeom = mesh->wallGeom();

      forall(const IntIntFloatPair& wall, WallGeom)
          LabelGeom[wall.first.first] += wall.second;

      // Find area, and signal on each wall
      if(signal == TOTAL or (report & TOTAL)) {
        forall(const vertex& v, S) {
          forall(const vertex& n, S.neighbors(v)) {
            vertex m = S.nextTo(v, n);
            if(!S.uniqueTri(v, n, m))
              continue;

            // Sum signal and area
            IntIntPair wall;
            int label = mesh->getLabel(v, n, m, parents);
            if(mesh->isBordTriangle(label, v, n, m, wall)) {
              float area = triangleArea(v->pos, n->pos, m->pos);
              WallBordGeom[wall] += area;
              WallBordSig[wall] += (v->signal + n->signal + m->signal) / 3.0 * area;
              if(polarity != NONE) {
                LabelBordGeom[wall.first] += area;
                LabelBordSig[wall.first] += (v->signal + n->signal + m->signal) / 3.0 * area;
              }
            }
          }
        }
      }

      forall(const IntIntFloatPair& wall, WallGeom) {
        // Fill in heat depending on user choice
        switch(signal) {
        case GEOMETRY:
          // Heat based on wall length
          if(polarity == CELL_AVERAGE)
            WallHeatAmt[wall.first]
                = (LabelGeom[wall.first.first] == 0 ? 0 : WallGeom[wall.first] / LabelGeom[wall.first.first]);
          else{
            if(mesh->useParents()){
              IntInt pP = std::make_pair(mesh->parents()[wall.first.first], mesh->parents()[wall.first.second]);
              WallHeatAmt[pP] += WallGeom[wall.first];
            } else
              WallHeatAmt[wall.first] = WallGeom[wall.first];
          }
          break;
        case TOTAL:
          // Heat based on wall signal
          WallHeatAmt[wall.first] = WallBordSig[wall.first];
          // Divide by area if required
          if(signalAverage)
            WallHeatAmt[wall.first]
                = (WallBordGeom[wall.first] == 0 ? 0 : WallHeatAmt[wall.first] / WallBordGeom[wall.first]);
          // Do relative to cell average
          if(polarity == CELL_AVERAGE) {
            if(signalAverage)
              WallHeatAmt[wall.first]
                  = (LabelBordSig[wall.first.first] == 0 or LabelBordGeom[wall.first.first] == 0
                  ? 0
                  : WallHeatAmt[wall.first]
                  / (LabelBordSig[wall.first.first] / LabelBordGeom[wall.first.first]));
            else
              WallHeatAmt[wall.first]
                  = (LabelBordSig[wall.first.first] == 0 ? 0 : WallHeatAmt[wall.first]
                  / LabelBordSig[wall.first.first]);
          }
          break;
        case BORDER:
        case INTERIOR:
        case BORDER_TOTAL:
        case INTERIOR_TOTAL:
          break;
        }
      }

      // Maps for polarity calc relative to minimum wall
      if(polarity == WALL_MIN) {
        IntFloatAttr LabelMin;
        // Find min heat per cell
        forall(const IntIntFloatPair& wall, WallHeatAmt)
            if(LabelMin.find(wall.first.first) == LabelMin.end())
            LabelMin[wall.first.first] = wall.second;
        else if(LabelMin[wall.first.first] > wall.second)
          LabelMin[wall.first.first] = wall.second;
        // Divide heat amounts
        forall(const IntIntFloatPair& wall, WallHeatAmt)
            WallHeatAmt[wall.first]
            = (LabelMin[wall.first.first] == 0 ? 0 : WallHeatAmt[wall.first] / LabelMin[wall.first.first]);
      }
    }
    return true;
  }
  REGISTER_PROCESS(ComputeHeatMap);

  bool DeleteHeatRangeLabels::run(Mesh* mesh, bool rescale, bool deleteCells, float min, float max)
  {
    Point2f& heatMapBounds = mesh->heatMapBounds();
    if(isNan(min))
      min = heatMapBounds[0];
    if(isNan(max))
      max = heatMapBounds[1];
    if(min >= max) {
      setErrorMessage("Error, the minimum must be strictly less than the maximum");
      return false;
    }
    IntFloatAttr& labelHeat = mesh->labelHeat();
    if(labelHeat.empty())
      return true;
    vvGraph& S = mesh->graph();
    std::vector<vertex> T;

    float vmin = 1, vmax = 0;

    IntSet deleteLabel;
    for(IntFloatAttr::iterator it = labelHeat.begin(); it != labelHeat.end(); ++it) {
      float v = it->second;
      if(v >= min and v <= max)
        deleteLabel.insert(it->first);
      else if(rescale) {
        if(v < vmin)
          vmin = v;
        if(v > vmax)
          vmax = v;
      }
    }
    if(deleteLabel.empty())
      return true;

    if(rescale and vmin < vmax)
      heatMapBounds = Point2f(vmin, vmax);

    forall(const vertex& v, S) {
      if(deleteLabel.count(v->label)) {
        if(deleteCells)
          T.push_back(v);
      }
    }
    forall(const vertex& v, T)
        S.erase(v);
    forall(int label, deleteLabel)
        labelHeat.erase(label);

    if(deleteCells) {
      // At last, erase all vertices not part of a triangle
      do {
        T.clear();
        forall(const vertex& v, S) {
          if(S.valence(v) < 2)
            T.push_back(v);
          else if(S.valence(v) == 2) {
            const vertex& n = S.anyIn(v);
            const vertex& m = S.nextTo(v, n);
            if(!S.edge(n, m))
              T.push_back(v);
          }
        }
        forall(const vertex& v, T)
            S.erase(v);
      } while(!T.empty());
      mesh->updateLines();
    }
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(DeleteHeatRangeLabels);

  bool RescaleHeatMap::run(Mesh* m, float minValue, float maxValue)
  {
    if(isNan(minValue))
      minValue = m->heatMapBounds()[0];
    if(isNan(maxValue))
      maxValue = m->heatMapBounds()[1];

    if(maxValue <= minValue)
      maxValue = minValue + 1;

    m->heatMapBounds() = Point2f(minValue, maxValue);

    return true;
  }
  REGISTER_PROCESS(RescaleHeatMap);

  bool SaveHeatMap::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(parent)
      filename = QFileDialog::getSaveFileName(parent, "Choose spreadsheet file to save", QDir::currentPath(),
                                              "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Filename", filename);
    return true;
  }

  bool SaveHeatMap::run(Mesh* mesh, const QString& filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);

    out << "Label,Value" << endl;
    forall(const IntFloatPair& p, mesh->labelHeat())
        out << p.first << "," << p.second << endl;

    SETSTATUS(QString("Wrote heat map file with %1 labels").arg(mesh->labelHeat().size()));
    return true;
  }
  REGISTER_PROCESS(SaveHeatMap);

  bool LoadHeatMap::initialize(QWidget* parent)
  {
    QString fileName = parm("File Name");
    ui = 0;
    // Return if file name filled in
    if(!fileName.isEmpty())
      return selectFile(parm("File Name"));

    if(!parent) 
      return true;

    ui = new Ui_LoadHeatMap;
    dlg = new QDialog(parent);
    ui->setupUi(dlg);
    connect(ui->SelectHeatMapFile, SIGNAL(clicked()), this, SLOT(selectFile()));
    connect(ui->HeatMapColumn, SIGNAL(currentIndexChanged(int)), this, SLOT(selectColumn(int)));
    ui->HeatMapColumn->setSizeAdjustPolicy(QComboBox::AdjustToContents);
    if(dlg->exec() == QDialog::Accepted) {
      setParm("File Name", ui->HeatMapFile->text());
      setParm("Column", QString::number(ui->HeatMapColumn->currentIndex()));
      return true;
    }
    return false;
  }

  bool LoadHeatMap::selectFile(const QString &filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly))
      return setErrorMessage(QString("'%1' cannot be opened for reading").arg(filename));

    QTextStream ss(&file);
    QString line = ss.readLine();
    QStringList fields = line.split(",");
    if(fields.size() < 2)
      return setErrorMessage(("CSV file must have at least 2 columns"));

    QRegExp unit("(.*)\\((.*)\\)\\s*");
    if(ui)
      ui->HeatMapColumn->clear();
    columnUnits.clear();
    QStringList columnNames;
    for(int i = 0; i < fields.size(); ++i) {
      // Trim off whitespace
      QString txt = fields[i].trimmed();
      // Trim off quotes
      if(txt[0] == '"' and txt[txt.size() - 1] == '"')
        txt = txt.mid(1, txt.size() - 2).trimmed();
      // Units can't be at the front
      int index = unit.indexIn(txt);
      if(index == 0) {
        columnNames << unit.cap(1).trimmed();
        columnUnits << unit.cap(2).trimmed();
      } else {
        columnNames << txt;
        columnUnits << "";
      }
      if(ui)
        ui->HeatMapColumn->addItem(txt);
    }

    QString columnName = parm("Column");
    int columnIndex = columnName.toInt();
    // If not a number, search by name
    if(columnIndex == 0) {
      for(int i = 1; i < columnNames.size(); i++)
        if(columnName == columnNames[i]) {
          columnIndex = i;
          break;
      }
    }
    if(ui)
      ui->HeatMapFile->setText(filename);
    if(columnIndex == 0) 
      return setErrorMessage(QString("Cannot find column called '%1'").arg(columnName));

    setParm("Column", QString::number(columnIndex));

    if(ui)
      ui->HeatMapColumn->setCurrentIndex(columnIndex);

    return true;
  }

  void LoadHeatMap::selectFile()
  {
    QString filename = QFileDialog::getOpenFileName(dlg, "Choose spreadsheet to load", ui->HeatMapFile->text(),
                                                    "CSV files (*.csv);;All files (*.*)");
    if(!filename.isEmpty())
      selectFile(filename);
  }

  void LoadHeatMap::selectColumn(int i)
  {
    if(i >= 0 and i < columnUnits.size())
      ui->unitLabel->setText(columnUnits[i]);
  }

  bool LoadHeatMap::run(Mesh* mesh, const QString& filename, int column, float border_size)
  {
    if(column < 0) {
      setErrorMessage(QString("Invalid column '%1'").arg(column));
      return false;
    }

    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream ss(&file);
    QString line = ss.readLine();
    QStringList fields = line.split(",");
    if(fields.size() < column + 1) {
      setErrorMessage(QString("File '%1' should be a CSV file with at least %2 columns").arg(filename).arg(column + 1));
      return false;
    }
    QRegExp re_unit(".*\\((.*)\\)\\s*");
    QString unit;
    QString txt = fields[column].trimmed();
    int nb_fields = fields.size();
    if(re_unit.indexIn(txt) == 0)
      unit = re_unit.cap(1);
    QString col2 = fields[1];
    if(col2[0] == '"' and col2[col2.size() - 1] == '"')
      col2 = col2.mid(1, col2.size() - 2).trimmed();
    bool is_wall = col2.trimmed() == "Neighbor";
    if(is_wall)
      mesh->updateWallGeometry(border_size);
    float minValue = HUGE_VAL, maxValue = -HUGE_VAL;
    int line_nb = 1;
    IntFloatAttr& labelMap = mesh->labelHeat();
    IntIntFloatAttr& wallMap = mesh->wallHeat();
    while(ss.status() == QTextStream::Ok) {
      ++line_nb;
      fields = ss.readLine().split(",");
      if(fields.size() != nb_fields or fields[0].isEmpty())
        break;
      bool ok;
      int neighbor = 0;
      int label = fields[0].toInt(&ok);
      if(!ok) {
        setErrorMessage(QString("Error line %1: the first column is not an integer number but '%2'")
                        .arg(line_nb)
                        .arg(fields[0].trimmed()));
        return false;
      }
      if(is_wall) {
        neighbor = fields[1].toInt(&ok);
        if(!ok) {
          setErrorMessage(QString("Error line %1: the second column is not an integer number but '%2'")
                          .arg(line_nb)
                          .arg(fields[1].trimmed()));
          return false;
        }
      }
      float value = fields[column].toFloat(&ok);
      if(!ok and fields[column]!="") {
        setErrorMessage(QString("Error line %1: the column %2 is not an floating point number but '%3'")
                        .arg(line_nb)
                        .arg(column)
                        .arg(fields[column].trimmed()));
        return false;
      }
      if(is_wall)
        wallMap[std::make_pair(label, neighbor)] = value;
      else
        labelMap[label] = value;
      if(value < minValue)
        minValue = value;
      if(value > maxValue)
        maxValue = value;
    }
    if(minValue >= maxValue)
      return setErrorMessage("The file contains 0 or 1 value only. We cannot create a heat map from that.");
    mesh->heatMapUnit() = unit;
    mesh->heatMapBounds() = Point2f(minValue, maxValue);
    mesh->updateTriangles();
    mesh->setShowLabel("Label Heat");
    SETSTATUS(QString("Loaded heat map with values ranging from %1%3 to %2%3").arg(minValue).arg(maxValue).arg(unit));
    return true;
  }
  REGISTER_PROCESS(LoadHeatMap);

  bool HeatMapGrowth::run(Mesh* m, Mesh* m2, QString name, QString t1, QString t2, bool createCustom)
  {

    QString hmap_name;
    CleanHeatMap *clean;
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
     throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process: %1"));
    CreateAttrMap *save_attr;
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
     throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process: %1"));

    if(createCustom){
      //Save heatmap to attr
      std::cout<<"Clean and save generating heatmap"<<std::endl;

      clean->run();

      save_attr->setParm("Name", name + "CustomDirHeatmap" + t1 + t2);
      hmap_name = save_attr->parm("Name");
      save_attr->run();

      //Compute custom directions
      std::cout<<"Reload generative heatmap and compute custom growth directions"<<std::endl;
      HeatAndChangeMap *GenMap;
      if(!getProcess("Mesh/Heat Map/Heat Map", GenMap))
        throw QString("MeshProcessHeatMap:: Unable to make HeatAndChangeMap process");
      GenMap->parm("Selection") = name + "CustomDirHeatmap" + t1 + t2;
      //GenMap->parm("Max valid") = "No";
      //GenMap->parm("Parent based") = "No";
      GenMap->run();

      CreateCustomDirectionsHeat *custom_dir;
	    if(!getProcess("Mesh/Cell Axis/Custom/Create Heatmap Directions", custom_dir))
        throw QString("MeshProcessHeatMap:: Unable to make CreateCustomDirectionsHeat process");
        custom_dir->run();
    }

    //Save T1 area
    std::cout<<"Create clean and save AreaT1"<<std::endl;
    MeasureArea *Area1;
    if(!getProcess("Mesh/Heat Map/Measures/Geometry/Area", Area1))
      throw(QString("MeshProcessHeatMap:: Unable to make MeasureArea process"));
    //std::cout<<"got parms"<<std::endl;
    Area1->run();
    std::cout<<"got area"<<std::endl;

    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    std::cout<<"cleaned map"<<std::endl;

    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name+"CellArea" + t1);
    hmap_name = save_attr->parm("Name Time Point 1");
    save_attr->run();
    //std::cout<<"saved parms"<<std::endl;

    //Save T2 area
    std::cout<<"Create clean and save AreaT2 on the other mesh"<<std::endl;
    HeatAndChangeMap *Area2;
    if(!getProcess("Mesh/Heat Map/Heat Map", Area2))
      throw(QString("MeshProcessHeatMap:: Unable to make HeatAndChangeMap process"));
    Area2->setParm("Identical Label Behavior", "Sum"); // parent based heat map
    Area2->setParm("Measure", "Geometry/Area"); // parent based heat map
    Area2->setParm("Use Labels Active Mesh", "No");

    // switch to the other mesh
    if(m == mesh(0)){
      Area2->setCurrentMeshId(1);
    } else {
      Area2->setCurrentMeshId(0);
    }
    m2->setUseParents(true);
    m2->updateTriangles();
    Area2->run();

    //std::cout << __LINE__ << std::endl;

    //Export to attr map
    //default + save to other mesh (= the first mesh)
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "CellArea" + t2);
    save_attr->setParm("Key", "Parent");
    save_attr->setParm("Heat Map Type", "Parent Heat");
    save_attr->setParm("Save to", "Active Mesh");
    if(m == mesh(0))
      save_attr->setCurrentMeshId(1);
    else
      save_attr->setCurrentMeshId(0);
    save_attr->run();

    //std::cout << __LINE__ << std::endl;

    save_attr->setParm("Save to", "Other Mesh");
    save_attr->run();
    //std::cout << __LINE__ << std::endl;
    // switch back to the first mesh
    if(m == mesh(0))
      save_attr->setCurrentMeshId(0);
    else
      save_attr->setCurrentMeshId(1);

    //Load using HeatAndChangeMap
    HeatAndChangeMap *Area2M1;
    if(!getProcess("Mesh/Heat Map/Heat Map", Area2M1))
      throw(QString("MeshProcessHeatMap:: Unable to make HeatAndChangeMap process"));
    Area2M1->setParm("Measure", name + "CellArea" + t1);
    //Area2M1->setParm("Parent based","No");
    Area2M1->run();
    //Clean
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    // save_attr = getProcessParms<CreateAttrMap>(this,parms);
    // parms[1] = name + "CellArea" + t2;
    // parms[2] = "Label";
    // parms[3] = "Label Heat";
    // parms[4] = "Active Mesh";
    // save_attr->run(parms);
//    return true;
    //Save T2 growth
    std::cout<<"Create clean and save areal increase"<<std::endl;
    CombineAttrMaps *Growth2;
    if(!getProcess("Mesh/Heat Map/Operators/Heat Map Transform", Growth2))
      throw(QString("MeshProcessHeatMap:: Unable to make CombineAttrMaps process"));

    Growth2->setParm("Attr Map 1", name + "CellArea" + t2);
    Growth2->setParm("Attr Map 2", name + "CellArea" + t1);
    Growth2->setParm("Combination", "Yes");
    Growth2->setParm("Combination Type", "Ratio (A/B)");
    Growth2->setParm("Transformation", "No");
    Growth2->setParm("New Attr Map", name + "GrowthArea" + t1 + t2);
    Growth2->run(); //Should save growth

    /* should be clean as T1 and T2 are clean - so try to skip
    clean = getProcessParms<CleanHeatMap>(this,parms);
    clean->run(parms);
    //Also - should already be saved by combine - so skip this as well
    //Resave using saved name
    save_attr = getProcessParms<CreateAttrMap>(this,parms);
    parms[1] = name + "GrowthT2";
    parms[4] = "Active Mesh";
    save_attr->run(parms);
*/
    //Save T2 prolif
    std::cout<<"Create, clean and save Proliferation"<<std::endl;
    HeatMapProliferation *Prolif2;
    if(!getProcess("Mesh/Lineage Tracking/Heat Map Proliferation", Prolif2))
      throw(QString("MeshProcessHeatMap:: Unable to make HeatMapProliferation process"));

    if(m == mesh(0))
      Prolif2->setCurrentMeshId(1);
    else
      Prolif2->setCurrentMeshId(0);

    Prolif2->run();

    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "Proliferation" + t1 + t2);
    save_attr->setParm("Key", "Parent");
    save_attr->setParm("Heat Map Type", "Parent Heat");
    save_attr->setParm("Save to", "Other Mesh");
    if(m == mesh(0))
      save_attr->setCurrentMeshId(1);
    else
      save_attr->setCurrentMeshId(0);
    save_attr->run();
    if(m == mesh(0))
      save_attr->setCurrentMeshId(0);
    else
      save_attr->setCurrentMeshId(1);

    //Load using HeatAndChangeMap
    HeatAndChangeMap *Prolif2M1;
    if(!getProcess("Mesh/Heat Map/Heat Map", Prolif2M1))
      throw(QString("MeshProcessHeatMap:: Unable to make HeatAndChangeMap process"));

    Prolif2M1->setParm("Measure", name + "Proliferation" + t1 + t2);
    //parms[0] = saved name
    //Prolif2M1->setParm("Parent based", "No");
    Prolif2M1->run();
    //Clean
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "Proliferation" + t1 + t2);
    save_attr->setParm("Key", "Parent");
    save_attr->setParm("Heat Map Type", "Parent Heat");
    save_attr->setParm("Save to", "Other Mesh");
    save_attr->run();

    //Check correspondences and such
    std::cout<<"Compute correspondence and growth directions"<<std::endl;

    CorrespondenceJunctions *checkCor;
    if(!getProcess("Mesh/Cell Axis/PDG/Check Correspondence", checkCor))
      throw(QString("MeshProcessHeatMap:: Unable to make CorrespondenceJunctions process"));
    checkCor->run();

    GrowthDirections *comp_growth;
    if(!getProcess("Mesh/Cell Axis/PDG/Compute Growth Directions", comp_growth))
      throw(QString("MeshProcessHeatMap:: Unable to make GrowthDirections process"));
    comp_growth->run();

    CellAxisAttrMapExport *cellAxisExport;
    if(!getProcess("Mesh/Cell Axis/Cell Axis Export To Attr Map", cellAxisExport))
      throw(QString("MeshProcessHeatMap:: Unable to make CellAxisAttrMapExport process"));
    cellAxisExport->setParm("Name", name + "CellAxis" + t1 + t2);
    cellAxisExport->run();

    //Save max growth
    std::cout<<"Compute clean and save max growth"<<std::endl;

    DisplayPDGs *maxG;
    if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", maxG))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPDGs process"));
    maxG->setParm("Heatmap", "StretchMax");
    maxG->run();
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthMax" + t1 + t2);
    save_attr->setParm("Save to", "Active Mesh");
    save_attr->run();


    //Save min growth
    std::cout<<"Compute clean and save min growth"<<std::endl;
    DisplayPDGs *minG;
    if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", minG))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPDGs process"));
    minG->setParm("Heatmap", "StretchMin");
    minG->run();

    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthMin" + t1 + t2);
    save_attr->setParm("Save to" , "Active Mesh");
    save_attr->run();

    //Save aniso
    std::cout<<"Compute clean and save anisotropy"<<std::endl;
    DisplayPDGs *aniso;
    if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", aniso))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPDGs process"));
    aniso->setParm("Heatmap", "Aniso=StretchMax/StretchMin");
    aniso->run();
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthAnisotropy" + t1 + t2);
    save_attr->setParm("Save to", "Active Mesh");
    save_attr->run();

    // check if there are any custom directions, if not finish
    AttrMap<int, SymmetricTensor>& customDirections = m->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    if(customDirections.empty()){
      std::cout << "No Custom Directions found." << std::endl;
      return true;
    }

    //Save growth X
    std::cout<<"Compute clean and save customX"<<std::endl;
    DisplayPDGs *growthX;
    if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", growthX))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPDGs process"));
    growthX->setParm("Heatmap", "StretchCustomX");
    growthX->run();
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthCustomX" + t1 + t2);
    save_attr->setParm("Save to", "Active Mesh");
    save_attr->run();

    //Save growth Y
    std::cout<<"Compute clean and save customY"<<std::endl;
    DisplayPDGs *growthY;
    if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", growthY))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPDGs process"));
    growthY->setParm("Heatmap", "StretchCustomY");
    growthY->run();
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthCustomY" + t1 + t2);
    save_attr->setParm("Save to", "Active Mesh");
    save_attr->run();


    //Save ratio of X/Y
    std::cout<<"Compute clean and save custom growth ratio (Y/X)"<<std::endl;
    CombineAttrMaps *RatioCustom;
    if(!getProcess("Mesh/Heat Map/Operators/Heat Map Transform", RatioCustom))
      throw(QString("MeshProcessHeatMap:: Unable to make CombineAttrMaps process"));

	  RatioCustom->setParm("Attr Map 1", name + "GrowthCustomY" + t1 + t2);
    RatioCustom->setParm("Attr Map 2", name + "GrowthCustomX" + t1 + t2);
    RatioCustom->setParm("Combination", "Yes");
    RatioCustom->setParm("Combination Type", "Ratio (A/B)");
    RatioCustom->setParm("Transformation", "No");
    RatioCustom->setParm("New Attr Map", name + "GrowthCustomRatio" + t1 + t2);
    RatioCustom->run(); //Should save growth

    return true;

    //Save angle X
    std::cout<<"Compute clean and save angleX"<<std::endl;
    CustomDirectionsAngle *angleX;
    if(!getProcess("Mesh/Cell Axis/Custom/Custom Direction Angle", angleX))
      throw(QString("MeshProcessHeatMap:: Unable to make CustomDirectionsAngle process"));
    angleX->setParm("Direction Custom", "X");
    angleX->run();
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthCustomAngleX" + t1 + t2);
    save_attr->setParm("Save to", "Active Mesh");
    save_attr->run();


    //Save angle Y
    std::cout<<"Compute clean and save angleY"<<std::endl;
    CustomDirectionsAngle *angleY;
    if(!getProcess("Mesh/Cell Axis/Custom/Custom Direction Angle", angleY))
      throw(QString("MeshProcessHeatMap:: Unable to make CustomDirectionsAngle process"));

    angleY->setParm("Direction Custom", "Y");
    angleY->run();
    if(!getProcess("Mesh/Heat Map/Operators/Clean Heat Map", clean))
      throw(QString("MeshProcessHeatMap:: Unable to make CleanHeatMap process"));
    clean->run();
    //Resave using saved name
    if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
    save_attr->setParm("Name", name + "GrowthCustomAngleY" + t1 + t2);
    save_attr->setParm("Save to", "Active Mesh");
    save_attr->run();

    return true;
  }
  REGISTER_PROCESS(HeatMapGrowth);


  bool CleanHeatMap::run(Mesh* m, Mesh* m2, /*QString selectedMap,*/ bool useSelection, QString minValueString, QString maxValueString)
  {

//    std::cout<<"Clean Heat Map does not account for parents - so be warned!"<<std::endl;
    bool set_min = minValueString == "NA" ? false : true;
    bool set_max = maxValueString == "NA" ? false : true;
    //IGNORE PARENT LABEL PROBLEMS FOR NOW
  //  bool use_parents = selectedMap == "Parent labels" ? true : false;
    double min_value = 0;
    double max_value = 0;
    if(set_min)
      min_value = minValueString.toDouble();
    if(set_max)
      max_value = maxValueString.toDouble();

    std::set<int> labels;
    AttrMap<int, float>& heatData = m->labelHeat();
    AttrMap<int, float> cleanHeatData;
//    IntIntAttr & parents = mesh->parents();

    if(useSelection){
      labels = findAllSelectedLabels(m->graph());
    }
    else{
      labels = findAllLabelsSet(m->graph());
    }

    forall(const int &l, labels){
      if(set_min && heatData[l]<min_value)
        continue;
      if(set_max && heatData[l]>max_value)
        continue;
      cleanHeatData[l] = heatData[l];
    }
    m->labelHeat() = cleanHeatData;

    m->heatMapBounds() = m->calcHeatMapBounds();
    m->showHeat();
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(CleanHeatMap);


  bool CreateAttrMap::run(Mesh* m, Mesh* m2, QString prefix, QString name, QString key, QString value, QString meshNumber, bool cleanMesh)
  {
    Mesh* useMesh = m;
    if(meshNumber!="Active Mesh"){
      useMesh = m2;
    }

    AttrMap<int, double>& attrData = useMesh->attributes().attrMap<int, double>(prefix + " " + name);
    attrData.clear();

    std::vector<int> labels = findAllLabels(m->graph());

    if(key == "Parent" and value == "Parent Heat"){
      forall(IntFloatPair p, m->labelHeat())
          attrData[p.first] = p.second;
    } else if(key == "Parent" and value == "Label Heat"){
      forall(IntFloatPair p, m->labelHeat())
          attrData[m->parents()[p.first]] = p.second;
    } else if(key == "Label" and value == "Parent Heat"){
      // go through all labels, take parent value
      forall(const int& l, labels){
        if(l < 1) continue;
        attrData[l] = m->labelHeat()[m->parents()[l]];
      }
    } else { // Label, Label Heat
      std::cout<<"Saving the heat-map"<<std::endl;
      forall(IntFloatPair p, m->labelHeat()){
        bool add_value = !cleanMesh;

        if(cleanMesh){
          forall(const int& l, labels)
            if(l == p.first)
              add_value = true;
        }

        if(add_value){
          attrData[p.first] = p.second;
        }

      }
    }



    return true;
  }
  REGISTER_PROCESS(CreateAttrMap);

  QStringList getAllMeasures(Mesh* m, QString subfolder)
  {

    // find measure processes
    QStringList procs = mgx::listProcesses();

    QStringList addProcs;
    std::set<QString> procsSet;
    QStringList addAttr;
    forall(const QString& name, procs) {
      mgx::ProcessDefinition* def = mgx::getProcessDefinition(name);
      QStringList list = def->name.split("/");

      if(list[0] == "Mesh" and list[1] == "Heat Map" and list[2] == subfolder){
        QString newProc = list[3] + "/" + list[4];
        addProcs << newProc;
        procsSet.insert(newProc);
        //qDebug() << newProc << "\n";
      }

    }

    // now the attr maps
    Attributes *attributes;
    attributes = &m->attributes();
    QStringList attr = attributes->getAttrList();
    QStringList measureAttr;

    // forall(const QString &text, attr) {
    //   qDebug() << "s " << text << "\n";
    // }

    attr.sort(Qt::CaseSensitive);

    forall(const QString &text, attr) {
      //qDebug() << "a" << text << "\n";
      QStringList list = text.split(" ");

      if(list.size() < 4) continue;
      if(list[0] != "Measure") continue;
      if(list[2] != "Double") continue; // only "real" heat maps
      QString name;
      for(uint i = 3; i < (uint)list.size(); i++){
        name = name + list[i] + " ";
      }
      name = name.trimmed();
      if(procsSet.find(name) == procsSet.end()){
        addProcs << name;
      }
    }

    return addProcs;
  }

  void fillTreeWidgetWithMeasures(Mesh* m, QTreeWidget* tree, QString subfolder)
  {
    QStringList allProcs = getAllMeasures(m, subfolder);

    tree->clear();

    forall(QString s, allProcs){
      QTreeWidgetItem *itemNew = new QTreeWidgetItem(QStringList() << s);
      itemNew->setFlags(Qt::ItemIsEnabled | Qt::ItemIsSelectable | Qt::ItemIsUserCheckable);
      tree->addTopLevelItem(itemNew);
    }

  }

  bool CombineAttrMaps::initialize(QWidget *parent)
  {

    Ui_HeatMapManipulateDialog ui;
    QDialog dlg(parent);
    ui.setupUi(&dlg);
    Mesh* m = currentMesh();

    fillTreeWidgetWithMeasures(m, ui.mapTree, "Measures");
    fillTreeWidgetWithMeasures(m, ui.mapTree2, "Measures");

    // fill gui with current parameter
    QTreeWidgetItemIterator it(ui.mapTree);
    while (*it) {
      if ((*it)->text(0) == parm("Attr Map 1"))
        (*it)->setSelected(true);
      ++it;
    }

    QTreeWidgetItemIterator it2(ui.mapTree2);
    while (*it2) {
      if ((*it2)->text(0) == parm("Attr Map 2"))
        (*it2)->setSelected(true);
      ++it2;
    }
    if(dlg.exec() == QDialog::Accepted){

      if(ui.mapTree->currentItem()) setParm("Attr Map 1", ui.mapTree->currentItem()->text(0));
      if(ui.mapTree2->currentItem() and ui.combinationGroup->isChecked()) setParm("Attr Map 2", ui.mapTree2->currentItem()->text(0));
      setParm ("Combination", ui.combinationGroup->isChecked() ? "Yes" : "No");
      setParm ("Combination Type", ui.combinationType->currentText());
      setParm ("Transformation", ui.transformationGroup->isChecked() ? "Yes" : "No");
      setParm ("Transformation Type", ui.transformationType->currentText());
      setParm ("Transformation Type 2", ui.transformationType2->currentText());
      setParm ("Ignore Parent", ui.transformationIgnoreParent->text());
      setParm ("Lower Threshold", QString::number(ui.transformationLower->value()));
      setParm ("Upper Threshold", QString::number(ui.transformationUpper->value()));
      setParm ("New Attr Map", ui.combinedName->text());
      return true;

    } else {
      return false;
    }

  }

  double minVectorValue(std::vector<double>& values)
  {

    if(values.size() == 0) return 0;

    double res = HUGE_VAL;

    forall(double v, values){
      if(v < res) res = v;
    }

    return res;
  }

  double maxVectorValue(std::vector<double>& values)
  {

    if(values.size() == 0) return 0;

    double res = -HUGE_VAL;

    forall(double v, values){
      if(v > res) res = v;
    }

    return res;
  }

  double medianVectorValue(std::vector<double> values)
  {

    if(values.size() == 0) return 0;

    double res = 0.;
    sort(values.begin(), values.end());

    if (values.size() % 2 == 0){
      res = (values[values.size() / 2 - 1] + values[values.size() / 2]) / 2;
    } else {
      res = values[values.size() / 2];
    }

    return res;
  }

  double meanVectorValue(std::vector<double>& values)
  {

    if(values.size() == 0) return 0;

    double res = 0.;

    forall(double v, values){
      res += v;
    }

    return res/values.size();
  }

  double meanWeightedVectorValue(std::vector<double>& values, std::vector<double>& weights)
  {

    if(values.size() == 0) return 0;
    if(values.size() != weights.size()) return 0;

    double res = 0.;
    double sumWeights = 0.;

    //forall(double v, values){
    for(uint i = 0; i < values.size(); i++){
      res += values[i] * weights[i];
      sumWeights += weights[i];
    }
    if(sumWeights != 0)
      return res/sumWeights;
    else
      return 0;
  }


  double varianceVectorValue(std::vector<double>& values)
  {

    if(values.size() == 0) return 0;

    double mean = 0.;

    forall(double v, values){
      mean += v;
    }
    mean /= values.size();

    double res = 0;

    forall(double v, values){
      res += (mean - v) * (mean - v);
    }

    return res;
  }


  double geoMeanVectorValue(std::vector<double>& values){
  //  std::cout<<"In geo mean"<<std::endl;
    if(values.size() == 0) return 0;

    double res = 1.;
    double size = values.size();
    forall(double v, values){
      std::cout<<v<<std::endl;
      if(v<=0){
        std::cout<<"Geometric Mean: invalid value (<=0), ignored"<<std::endl;
      }

      res *= pow(v,1.0f/size);
    }
    std::cout<<res<<std::endl;
    return res;//pow(res,1.0f/(double)values.size());
  }


  double extractValue(QString type, std::vector<double>& values, std::vector<double>& weights)
  {

    if(type == "Mean"){
      return meanVectorValue(values);
    } else if(type == "Mean (Area Weighted)"){ // TODO!
      return meanWeightedVectorValue(values, weights);
    } else if(type == "Median"){
      return medianVectorValue(values);
    } else if(type == "Variance"){
      return varianceVectorValue(values);
    } else if(type == "Geometric Mean"){
      return geoMeanVectorValue(values);
    } else if(type == "Minimum"){
      return minVectorValue(values);
    } else if(type == "Maximum"){
      return maxVectorValue(values);
    } else if(type == "Max-Min"){
      return maxVectorValue(values) - minVectorValue(values);
    } else if(type == "Max/Min"){
      return maxVectorValue(values) / minVectorValue(values);
    }
    return 0;
  }

  bool CombineAttrMaps::run(Mesh* m, QString name1, QString name2,
      bool combination, QString combinationType,
      bool transformation, QString transformationType, QString transformationType2, int ignoreParent, double lowerT, double upperT,
      QString combinedAttrMap)
  {

    if(combination and transformation){
      return setErrorMessage("Please select either Combination or Transformation, not both!");
    }

    MeasureArea *Area1;
    if(!getProcess("Mesh/Heat Map/Measures/Geometry/Area", Area1))
      throw(QString("MeshProcessHeatMap:: Unable to make MeasureArea process"));
    //std::cout<<"got parms"<<std::endl;
    Area1->run();
    AttrMap<int, double>& attrDataArea = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Area");

    QString prefix = "Measure Label Double";

    vvGraph& S = m->graph();

    AttrMap<int, double>& attrData1 = m->attributes().attrMap<int, double>(prefix + " " + name1);

    QString measureFolder = "Mesh/Heat Map/Measures/";

    QString proName = measureFolder + name1;

    //qDebug() << "attr map \t" << proAttr << "  size  " << measureMap.size() << "\n";
    if(attrData1.empty()){
      Process *pro = makeProcess(proName);
      if(!pro) return setErrorMessage("Measure doesn't exist!");
      qDebug() << "Process call\t" << proName << "\n";
      pro->run();
    }

    if(attrData1.empty()) return setErrorMessage("Attr Map 1 doesn't exist!");

    AttrMap<int, double>& attrDataComb = m->attributes().attrMap<int, double>(prefix + " " + combinedAttrMap);
    attrDataComb.clear();
    if(combination){

      AttrMap<int, double>& attrData2 = m->attributes().attrMap<int, double>(prefix + " " + name2);

      if(attrData2.empty()){
        proName = measureFolder + name2;
        Process *pro = makeProcess(proName);
        if(!pro) return setErrorMessage("Measure doesn't exist!");
        qDebug() << "Process call\t" << proName << "\n";
        pro->run();
      }

      if(attrData2.empty()) return setErrorMessage("Attr Map 2 doesn't exist!");

      forall(IntDoublePair p, attrData1){
        if(combinationType == "Addition (A + B)"){
          attrDataComb[p.first] = p.second + attrData2[p.first];
        } else if(combinationType == "Multiplication (A * B)"){
          attrDataComb[p.first] = p.second * attrData2[p.first];
        } else if(combinationType == "Ratio (A/B)"){
          attrDataComb[p.first] = p.second / attrData2[p.first];
        } else {
          attrDataComb[p.first] = pow(p.second,attrData2[p.first]);
        }
      }
    } else { // transformation

      if(transformationType == "Parent Label"){
        std::map<int, std::set<int> > parentLabelMap;
        typedef std::pair<int, std::set<int> > IntSetIntP;
        std::map<int, double> parentValues;

        forall(const vertex& v, S){
          if(m->parents()[v->label] < 1) continue;
          parentLabelMap[m->parents()[v->label]].insert(v->label);
        }

        forall(IntSetIntP p, parentLabelMap){
          std::vector<double> values, weights;
          forall(int l, p.second){
            if(attrData1[l] >= lowerT and attrData1[l] <= upperT and m->parents()[p.first] != ignoreParent and m->parents()[l] != ignoreParent){
              values.push_back(attrData1[l]);
              weights.push_back(attrDataArea[l]);
            }
          }
          parentValues[p.first] = extractValue(transformationType2, values, weights);

          forall(int l, p.second){
            attrDataComb[l] = parentValues[p.first];
          }
        }

      } else {
        // create neighbor graph
        std::map<int, std::set<int> > cellNeighborMap;
        std::map<IntInt, double> neighMap;
        typedef std::pair<IntInt, double> IntIntDoubleP;
        typedef std::pair<int, std::set<int> > IntSetIntP;
        neighborhood2D(S, neighMap, "");

        forall(IntIntDoubleP p, neighMap){
          cellNeighborMap[p.first.first].insert(p.first.second);
        }

        forall(IntSetIntP p, cellNeighborMap){
          std::vector<double> values, weights;
          forall(int l, p.second){
            if(attrData1[l] >= lowerT and attrData1[l] <= upperT and m->parents()[p.first] != ignoreParent and m->parents()[l] != ignoreParent){
              values.push_back(attrData1[l]);
              weights.push_back(attrDataArea[l]);
            }
          }
          attrDataComb[p.first] = extractValue(transformationType2, values, weights);
        }


      }

    }

    // now visualize the new heat map
    forall(IntDoublePair p, attrData1){
      m->labelHeat()[p.first] = p.second;
    }

    forall(IntDoublePair p, attrDataComb){
      m->labelHeat()[p.first] = p.second;
    }

    m->heatMapBounds() = m->calcHeatMapBounds();
    m->showHeat();
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(CombineAttrMaps);


  bool SelectByHeat::run(Mesh* m, double lowerThreshold, double upperThreshold, bool resetSelection)
  {

    vvGraph& S = m->graph();

    if(resetSelection){
      forall(const vertex& v, S){
        v->selected = false;
      }
    }

    if(m->useParents()){

      forall(const vertex& v, S){
        int pLabel = m->parents()[v->label];
        if(pLabel < 1 or m->labelHeat().find(pLabel) == m->labelHeat().end()) continue;
        if(m->labelHeat()[pLabel] >= lowerThreshold and m->labelHeat()[pLabel] <= upperThreshold){
          v->selected = true;
        }
      }


    } else {

      forall(const vertex& v, S){
        if(v->label < 1 or m->labelHeat().find(v->label) == m->labelHeat().end()) continue;
        if(m->labelHeat()[v->label] >= lowerThreshold and m->labelHeat()[v->label] <= upperThreshold){
          v->selected = true;
        }

      }

    }

    // take care of borders in 2D meshes
    forall(const vertex& v, S){
      if(v->label > 0) continue;
      forall(const vertex& n, S.neighbors(v)){
        if(n->selected and n->label > 0) v->selected = true;
      }
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(SelectByHeat);


  bool SwitchHeatParent::run(Mesh* m, QString direction)
  {

    if(direction == "Heat -> Parent"){
      m->parents().clear();
      forall(IntFloatPair p, m->labelHeat()){
        m->parents()[p.first] = int(p.second);
      }
    } else if(direction == "Parent -> Heat"){
      m->labelHeat().clear();
      forall(IntIntPair p, m->parents()){
        m->labelHeat()[p.first] = p.second;
      }
    } else if(direction == "Parent <-> Heat"){
      IntIntAttr backupHeat;
      forall(IntFloatPair p, m->labelHeat()){
        backupHeat[p.first] = int(p.second);
      }
      IntFloatAttr backupParents;
      forall(IntIntPair p, m->parents()){
        backupParents[p.first] = p.second;
      }
      m->labelHeat() = backupParents;
      m->parents() = backupHeat;
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(SwitchHeatParent);


  void fillTreeWidgetWithMeasureProcesses(Mesh* m, QTreeWidget* tree)
  {
    QStringList procs = mgx::listProcesses();

    //ui.mapTree->header()->resizeSection(1, 50);
    //ui.mapTree->header()->resizeSection(0, 300);

    //populateMapTree(mesh);

    forall(const QString& name, procs) {
      mgx::ProcessDefinition* def = mgx::getProcessDefinition(name);
      QStringList list = def->name.split("/");
      QString measureFolder = "Mesh/Heat Map/Measures";
      QStringList folder = measureFolder.split("/");

      QString subfolder = folder[2];

      if(list[0] == "Mesh" and list[1] == "Heat Map" and list[2] == subfolder){
        QTreeWidgetItem *item2 = new QTreeWidgetItem(QStringList() << list[3] + "/" + list[4]);
        item2->setFlags(Qt::ItemIsEnabled | Qt::ItemIsSelectable | Qt::ItemIsUserCheckable);
        tree->addTopLevelItem(item2);
      }
    }

    //ui.mapTree->sortItems(0, Qt::AscendingOrder);
    tree->sortItems(0, Qt::AscendingOrder);
    //for(int i = 0; i < 1; i++)
    //  ui.mapTree->resizeColumnToContents(i);

  }

  void HeatAndChangeMap::changeMapType(const QString& type)
  {
    Ui_HeatMapDialogNew& ui = this->ui;

    if(ui.HeatMapType->currentText() == "Area 2D"){
      fillTreeWidgetWithMeasures(m, ui.measureProcessTree, "Measures");
    } else if(ui.HeatMapType->currentText() == "Volume 3D"){
      fillTreeWidgetWithMeasures(m, ui.measureProcessTree, "---"); // Measures 3D
    } else {

    }

  }

  bool HeatAndChangeMap::initialize(QWidget *parent)
  {
    if(!checkState().mesh(MESH_ANY))
      return false;

    if(parent and parm("Measure").isEmpty()) {
      //qDebug() << "init " << "\n";
  
      Ui_HeatMapDialogNew& ui = this->ui;
      QDialog dlg(parent);
      ui.setupUi(&dlg);
  
      Mesh* m = currentMesh();
      this->m = m;
      Mesh* m2;
      if(currentMesh() == mesh(0)){
        m2 = mesh(1);
        ui.label_active_mesh->setText("Mesh 1");
        ui.label_other_mesh->setText("Mesh 2");
      } else{
        m2 = mesh(0);
        ui.label_active_mesh->setText("Mesh 2");
        ui.label_other_mesh->setText("Mesh 1");
      }
  
      if(m->useParents()){
        ui.label_active_labels->setText("Parents");
      } else {
        ui.label_active_labels->setText("Labels");
      }
  
      if(m2->useParents()){
        ui.label_other_labels->setText("Parents");
      } else {
        ui.label_other_labels->setText("Labels");
      }
  
  //    if(ui.HeatMapType->findText(parms[12], Qt::MatchExactly) >= 0)
  //      ui.HeatMapType->setCurrentIndex(ui.HeatMapType->findText(parms[12], Qt::MatchExactly));
  
      if(ui.HeatMapType->currentText() == "Measures 2D"){
        fillTreeWidgetWithMeasures(m, ui.measureProcessTree, "Measures");
      } else if(ui.HeatMapType->currentText() == "Measures 3D"){
        fillTreeWidgetWithMeasures(m, ui.measureProcessTree, "---"); // Measures 3D
      } else {
  
      }
  
      // fill gui with current parameter
      QTreeWidgetItemIterator it(ui.measureProcessTree);
      while (*it) {
        if ((*it)->text(0) == parm("Measure"))
          (*it)->setSelected(true);
        ++it;
      }
      connect(ui.HeatMapType, SIGNAL(currentIndexChanged(const QString &)), this, SLOT(changeMapType(const QString &)));
  
      ui.forceRecalcCheckBox->setChecked(stringToBool(parm("Force Recalc")));
  
      //ui.parentBasedMap->setChecked(stringToBool(parm("Parent based")));
      if(ui.IdenticalLabelsComboBox->findText(parm("Identical Label Behavior"), Qt::MatchExactly) >= 0)
        ui.IdenticalLabelsComboBox->setCurrentIndex(ui.IdenticalLabelsComboBox->findText(parm("Identical Label Behavior"), Qt::MatchExactly));
  
      ui.HeatMapChange->setChecked(stringToBool(parm("Change Map")));
  
      if(ui.HeatMapIncr->findText(parm("Change Map Dir"), Qt::MatchExactly) >= 0)
        ui.HeatMapIncr->setCurrentIndex(ui.HeatMapIncr->findText(parm("Change Map Dir"), Qt::MatchExactly));
      if(ui.HeatMapDiffType->findText(parm("Change Map Type"), Qt::MatchExactly) >= 0)
        ui.HeatMapDiffType->setCurrentIndex(ui.HeatMapDiffType->findText(parm("Change Map Type"), Qt::MatchExactly));
  
      //qDebug() << "init " << "\n";
      QTreeWidgetItem* it0 = ui.measureProcessTree->topLevelItem(0);
      if(it0) 
        ui.measureProcessTree->setCurrentItem(it0);
  
  
      //qDebug() << "init " << "\n";
  
      if(dlg.exec() == QDialog::Accepted){
        // qDebug() << "acc " << "\n";
        // fill parameters
        if(ui.measureProcessTree->currentItem())
          setParm("Measure", ui.measureProcessTree->currentItem()->text(0));
  
        setParm("Force Recalc", ui.forceRecalcCheckBox->isChecked() ? "Yes" : "No");
  
        //if(ui.parentBasedMap->isChecked()){
          //setParm("Parent based", "Yes");
          setParm("Identical Label Behavior", ui.IdenticalLabelsComboBox->currentText());
        //} else {
          //setParm("Parent based", "No");
          //setParm("Identical Label Behavior", "");
        //}
        //   qDebug() << "4 " << "\n";
        if(ui.HeatMapChange->isChecked()){
          setParm("Change Map", "Yes");
          setParm("Change Map Dir", ui.HeatMapIncr->currentText());
          setParm("Change Map Type", ui.HeatMapDiffType->currentText());
          //parms[7] = ui.changeMapAssoc->currentText();
        } else {
          setParm("Change Map", "No");
          setParm("Change Map Dir", "");
          setParm("Change Map Type", "");
          setParm("Use Labels Active Mesh", "");
        }
        setParm("Use Labels Active Mesh", m->useParents()? "No" : "Yes");
        setParm("Use Labels Inactive Mesh", m2->useParents()? "No" : "Yes");
        //  qDebug() << "acc done " << "\n";
      } else
        return false;
//      {
//        setParm("Measure", "Aborted");
//      }

      if(parm("Measure").isEmpty())
        return false;
    }
    //  qDebug() << "gui done " << "\n";
    return true;
  }

  bool HeatAndChangeMap::getData(Mesh* m, QString measureFolder, QString measure, bool forceRecalc, AttrMap<int, double>& data)
  {
    int currentMeshId = this->currentMeshId();
    int otherMeshId = 1;
    if(currentMeshId == 1) otherMeshId=0;


    QString proName = measureFolder + measure;
    QString proAttr = "Measure Label Double " + measure;
    AttrMap<int, double>& measureMap = m->attributes().attrMap<int, double>(proAttr);

    Process *pro = makeProcess(proName);

    //qDebug() << "attr map \t" << proAttr << "  size  " << measureMap.size() << "\n";
    if((measureMap.empty() or forceRecalc) and pro){
      Process *pro = makeProcess(proName);
      qDebug() << "Process call\t" << proName << "\n";
      if(m != currentMesh())
        pro->setCurrentMeshId(otherMeshId);
      pro->run();
    } else {
      if(measureMap.empty()) return false;
    }
    //qDebug() << "attr map filled! \t" << proAttr << "  size  " << measureMap.size() << "\n";
    data = measureMap;
    //qDebug() << "attr map filled! \t" << proAttr << "  size  " << data.size() << "\n";
    return true;

  }

  double getChangeValue(bool dirIncreasing, QString type, double value1, double value2)
  {
    if(type == "Ratio"){
      if(dirIncreasing){
        double ratio = value2==0 ? -1E20 : value1/value2;
        return ratio;
      } else {
        double ratio = value1==0 ? -1E20 : value2/value1;
        return ratio;
      }
    } else if(type == "Difference"){
      if(dirIncreasing){
        double diff = value1-value2;
        return diff;
      } else {
        double diff = value2-value1;
        return diff;
      }
    }

    return -1E20;
  }

  bool HeatAndChangeMap::run(Mesh* m, Mesh* m2, QString measure, bool forceRecalc, QString identicalLabel,
                             bool changeMap, QString changeDir, QString changeType, bool m1Labels, bool m2Labels)
  {
    m->labelHeat().clear();

    if(measure == "Aborted"){
      m->updateTriangles();
      return setErrorMessage("Aborted");
    }
    if(measure == ""){
      m->updateTriangles();
      return setErrorMessage("No measure selected");
    } 

    qDebug() << "run " << "\n";
    QString measureFolder = "Mesh/Heat Map/Measures/";

    bool dirIncreasing = (changeDir == "Increasing");

    AttrMap<int, double> attr1, attr2;
    if(!getData(m, measureFolder, measure, forceRecalc, attr1)){
      m->updateTriangles();
      return setErrorMessage("Could not fetch the data");
    } 
    
    qDebug() << "attr Map size " << attr1.size() << "\n";

    if(changeMap){
      if(!getData(m2, measureFolder, measure, forceRecalc, attr2)) return setErrorMessage("Attribute Map does not exist for Mesh 2 and is not a valid Measure Process!");

      m2->labelHeat().clear();

      if(m1Labels and m2Labels){ // both using labels, just compare the labels

        forall(IntDoublePair p, attr1){
          double changeValue = getChangeValue(dirIncreasing,changeType,attr2[p.first], p.second);
          if(changeValue != -1E20 and attr2.find(p.first) != attr2.end()) m->labelHeat()[p.first] = changeValue;
        }

      } else if((!m1Labels and m2Labels) or (m1Labels and !m2Labels)){ // one mesh parents, the other labels

        std::map<int, double> parentValueMap;
        std::map<int, double> parentMaxMap;
        std::map<int, double> parentMinMap;
        std::map<int,int> parentCount;

        AttrMap<int, double> attr1P, attr2P;

        if(!m1Labels){
          attr1P = attr1;
          attr2P = attr2;
        } else {
          attr2P = attr1;
          attr1P = attr2;
        }

	      if(!m1Labels)
          dirIncreasing = !dirIncreasing;

        forall(IntDoublePair p, attr1P){
	        if( (!m1Labels && m->parents().find(p.first)==m->parents().end()) || (m1Labels && m2->parents().find(p.first)==m2->parents().end()))
		        continue;

	        int pLabel = !m1Labels ? m->parents()[p.first] : m2->parents()[p.first];

	        if(pLabel == 0)
		        continue;

	        parentCount[pLabel]++;

          if(identicalLabel == "Max" || identicalLabel == "Min" || identicalLabel == "Max-Min"){
	          if(parentMaxMap.find(pLabel) == parentMaxMap.end()){
		          parentMaxMap[pLabel] = p.second;
		          parentMinMap[pLabel] = p.second;
	          } else{
		          parentMaxMap[pLabel] = p.second > parentMaxMap[pLabel] ? p.second : parentMaxMap[pLabel];
		          parentMinMap[pLabel] = p.second < parentMinMap[pLabel] ? p.second : parentMinMap[pLabel];
            }
	        } else {//  if(identicalLabel == "Average" || identicalLabel == "Sum")
	          parentValueMap[pLabel] += p.second;
          }

        }

        forall(IntDoublePair p, attr2P){
	        double parentValue;
          if(parentCount[p.first] < 1) continue;
          parentValue = parentValueMap[p.first];
          if(identicalLabel == "Average"){
            parentValue /= (double)parentCount[p.first];
	        } else if(identicalLabel == "Max"){
		        parentValue = parentMaxMap[p.first];
	        } else if(identicalLabel == "Min"){
		        parentValue = parentMinMap[p.first];
	        } else if(identicalLabel == "Max-Min"){
		        parentValue = parentMaxMap[p.first]-parentMinMap[p.first];
          }



          double changeValue = getChangeValue(dirIncreasing,changeType,parentValue,p.second);

          if(changeValue != -1E20) m->labelHeat()[p.first] = changeValue;
        }

      } else { // parents of both, should that even be allowed?
        return setErrorMessage("Parents can't be applied for both meshes");
      }

    } else { //Parent based heat-map for a single mesh

      if(!m1Labels){
        forall(IntDoublePair p, attr1){
	        if(m->parents().find(p.first)==m->parents().end() || m->parents()[p.first] == 0) continue;
          m->labelHeat()[m->parents()[p.first]] = 0;
        }

	    if(identicalLabel == "Average" || identicalLabel == "Sum"){
        std::map<int,int> parentCount;
        forall(IntDoublePair p, attr1){ // Sum up values
  	      if(m->parents().find(p.first)==m->parents().end() || m->parents()[p.first] == 0) continue;
            m->labelHeat()[m->parents()[p.first]] += p.second;
            parentCount[m->parents()[p.first]]++;
          }

          if(identicalLabel == "Average"){
            forall(IntInt p, parentCount){
              m->labelHeat()[p.first]/=p.second;
            }
          }
	    } else if(identicalLabel == "Max" || identicalLabel == "Min" || identicalLabel == "Max-Min"){
        std::map<int, double> parentMaxMap;
        std::map<int, double> parentMinMap;
        forall(IntDoublePair p, attr1){ //Calc max/min
  	      if(m->parents().find(p.first)==m->parents().end() || m->parents()[p.first] == 0) continue;

	        int pLabel = m->parents()[p.first];
  	      if(parentMaxMap.find(pLabel) == parentMaxMap.end()){
	          parentMaxMap[pLabel] = p.second;
	          parentMinMap[pLabel] = p.second;
	        } else{
	          parentMaxMap[pLabel] = p.second > parentMaxMap[pLabel] ? p.second : parentMaxMap[pLabel];
	          parentMinMap[pLabel] = p.second < parentMinMap[pLabel] ? p.second : parentMinMap[pLabel];
          }
	      }
        
        forall(IntDoublePair p, attr1){ //Calc max/min
  	    if(m->parents().find(p.first)==m->parents().end() || m->parents()[p.first] == 0) continue;

	      int pLabel = m->parents()[p.first];
//	    std::cout<<p.first<<" "<<parentMaxMap[p.first]<<std::endl;
  	    if(identicalLabel == "Max")
		      m->labelHeat()[pLabel] = parentMaxMap[pLabel];
	      else if(identicalLabel == "Min")
		      m->labelHeat()[pLabel] = parentMinMap[pLabel];
	      else if(identicalLabel == "Max-Min")
		      m->labelHeat()[pLabel] = parentMaxMap[pLabel]-parentMinMap[pLabel];
        }
	    } else if(identicalLabel == "None"){ // load the attribute map as a plain parent attribute map
        forall(IntDoublePair p, attr1){
          m->labelHeat()[p.first] = p.second;
        }
      }



      } else {
        forall(IntDoublePair p, attr1){
          m->labelHeat()[p.first] = p.second;
        }
      }

    }

    //if((!m1Labels and changeMap) or (parentMap and m2Labels)) m->setUseParents(true);
    //else m->setUseParents(false);

    // forall(auto p, m->labelHeat()){
    //   std::cout << "heat " << p.first << "/" << p.second << std::endl;
    // }


    m->heatMapBounds() = m->calcHeatMapBounds();
    // if(useRange) m->heatMapBounds() = Point2f(rangeMin, rangeMax);
    m->setShowLabel("Label Heat");
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(HeatAndChangeMap);

  bool RunAllParentHeat::run(){
     if(!checkState().mesh(MESH_ANY))
   	return false;
     Mesh *m = currentMesh();
     Mesh* m2;
     if(currentMesh() == mesh(0))
        	m2 = mesh(1);
     else
        	m2 = mesh(0);
     QStringList identicalLabel = (QStringList() << "Sum" << "Average" << "Max"<<"Min"<<"Max-Min");
     for(int i=0;i<identicalLabel.size();i++){
	 HeatAndChangeMap *makeHeat;
	 if(!getProcess("Mesh/Heat Map/Heat Map", makeHeat))
     		throw(QString("MeshProcessHeatMap:: Unable to make Heat Map process: %1"));
	 makeHeat->run(m,m2,parm("Attr map"), false, identicalLabel[i],false, "NA", "NA", true, true);
	 QString AttrName = parm("Prefix")+identicalLabel[i];
	 CreateAttrMap *save_attr;
	 if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
     		throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process: %1"));
	 save_attr->setParm("Name",AttrName);
	 save_attr->run();
     }
     return true;
  }
  REGISTER_PROCESS(RunAllParentHeat);

  bool HeatMapCellType::run(Mesh* m, QString mode, int bins)
  {

    if(mode == "Bins Cells"){
      if(bins < 1) return setErrorMessage("Error: The bin number hast to be larger than 1");


      std::vector<double> heatValues;
      std::map<double, std::set<int> > heatLabelMap;

      forall(IntFloatPair p, m->labelHeat()){
        heatLabelMap[p.second].insert(p.first);
        heatValues.push_back(p.second);
      }

      double binSize = (double)heatValues.size() / (double)bins;


      // sort the heat values
      std::sort(heatValues.begin(), heatValues.end());

      int counter = 0;
      int parentLabel = 1;
      int heatLast = *(heatValues.begin()) - 1;

      forall(double v, heatValues){
        counter++;
        if(heatLast != v){
          forall(int l, heatLabelMap[v]){
            m->parents()[l] = parentLabel;
          }
        }

        if(counter > binSize and parentLabel < bins){
          parentLabel++;
          counter -= binSize;
        }
      }

    } else if(mode == "Bins Heat"){
      if(bins < 1) return setErrorMessage("Error: The bin number hast to be larger than 1");

      double minHeat = 1E20, maxHeat = -1E20;

      forall(IntFloatPair p, m->labelHeat()){
        if(p.second < minHeat) minHeat = p.second;
        if(p.second > maxHeat) maxHeat = p.second;
      }

      double binSize = (maxHeat-minHeat) / (double)bins;

      std::cout << "p " <<  binSize << std::endl;

      forall(IntFloatPair p, m->labelHeat()){
        int parent = ((double)((p.second - minHeat) / binSize)) + 1;
        if(parent > bins) parent = bins;
        m->parents()[p.first] = parent;
      }

    } else { // Threshold mode

      SelectByHeat *sbh;
      if(!getProcess("Mesh/Heat Map/Heat Map Select", sbh))
        throw(QString("MeshProcessHeatMap:: Unable to make SelectByHeat process"));
      sbh->run();

      forall(const vertex& v, m->graph()){
        if(v->selected){
          m->parents()[v->label] = 1;
          v->selected = false;
        } else {
          m->parents()[v->label] = 2;
        }
      }

    }

    m->setUseParents(true);
    m->setShowLabel("Label");
    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(HeatMapCellType);

  bool HeatMapAllMeasures2D::run(bool excludeSlow)
  {

    QStringList procs = mgx::listProcesses();
    QString folderName = "Mesh/Heat Map/Measures";

    QStringList slowProcesses;
    slowProcesses << "Mesh/Heat Map/Measures/Lobeyness/Visibility Stomata";
    slowProcesses << "Mesh/Heat Map/Measures/Lobeyness/Visibility Pavement";

    forall(const QString& proName, procs) {
      bool slow = false;
      forall(QString p, slowProcesses){
        if(p == proName){
          slow = true;
          qDebug() << "Exclude Process\t" << proName << "\n";
        }
      }
      if(slow) continue;

      mgx::ProcessDefinition* def = mgx::getProcessDefinition(proName);
      QStringList list = def->name.split("/");
      QStringList folder = folderName.split("/");

      QString subfolder = folder[2];

      if(list[0] == "Mesh" and list[1] == "Heat Map" and list[2] == subfolder){
        if(list[3] == "Location") continue; // location heat maps have to be run manually as they mostly require additional inputs

        Process *pro = makeProcess(proName);
        qDebug() << "Process call\t" << proName << "\n";
        pro->run();
      }
    }
    return true;
  }
  REGISTER_PROCESS(HeatMapAllMeasures2D);


  bool HeatMapSmooth::run(Mesh* m, int passes, QString mode, bool rescale, QString attrName)
  {

  if(passes < 1)
    return setErrorMessage("Please enter a number of Passes > 0");

  if(passes > 1000)
    return setErrorMessage("Please enter a number of Passes < 1000");

  IntFloatAttr& heatmap = m->labelHeat();
  IntFloatAttr newHeat;
  IntFloatAttr area;

  vvGraph& S = m->graph();

  MeasureArea ma(*this);
  ma.calculateArea(S, area);

  // create neighbor graph
  std::map<int, std::set<int> > cellNeighborMap;
  std::map<IntInt, double> neighMap;
  typedef std::pair<IntInt, double> IntIntDoubleP;
  typedef std::pair<int, std::set<int> > IntSetIntP;
  neighborhood2D(S, neighMap, "");

  forall(IntIntDoubleP p, neighMap){
    cellNeighborMap[p.first.first].insert(p.first.second);
  }

  while(passes > 0){
    passes--;
    std::map<int, double> labelHeatSmoothed;
    std::map<int, double> labelWeight;

    forall(IntSetIntP p, cellNeighborMap){
      p.second.insert(p.first);
      forall(int l, p.second){
        double weight = 1.;
        if(mode == "Area"){
          weight = area[l];
        }

        labelHeatSmoothed[p.first] += heatmap[l] * weight;
        labelWeight[p.first] += weight;
      }
      
    }

    forall(auto p, labelHeatSmoothed){
      if(labelWeight[p.first] > 0)
        labelHeatSmoothed[p.first] /= labelWeight[p.first];
      else
        continue;

      heatmap[p.first] = labelHeatSmoothed[p.first];
    }
  }



  if(attrName != ""){
    QString fullAttrName = "Measure Label Double " + attrName; 
    AttrMap<int, double>& attr = m->attributes().attrMap<int, double>(fullAttrName);
    attr.clear();

    forall(auto p, heatmap){
      attr[p.first] = p.second;
    }
  }

  if(rescale){
    m->heatMapBounds() = m->calcHeatMapBounds();
  }

  m->setShowLabel("Label Heat");
  m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(HeatMapSmooth);

  bool ApplyHeat::run(Mesh* mesh, const QString& heatUnit, double heat, double minHeat, double maxHeat)
  {
    const vvGraph& vs = mesh->selectedVertices();
  
    std::set<int> labels;
    forall(const vertex& v, vs)
      labels.insert(v->label);
  
    IntFloatAttr &labelHeat = mesh->labelHeat();
    forall(int label, labels)
      labelHeat[label] = heat;

    mesh->heatMapBounds() = Point2f(minHeat, maxHeat);
    mesh->heatMapUnit() = heatUnit;
    mesh->setShowLabel("Label Heat");
    mesh->updateTriangles();
  
    return(true);
  }
  REGISTER_PROCESS(ApplyHeat);

    // bool HeatMapCellDistanceSignal::run(Mesh* m){

    //   vvGraph& S = m->graph();


    //   // forall(vertex v, S){
    //   //   v->signal = 0;
    //   //   if(v->label == -1){
    //   //     if(cellBorderHandling == "Ignore") continue;
    //   //     double avg = 0;
    //   //     int counter = 0;
    //   //     forall(vertex n, S.neighbors(v)){
    //   //       if(m->labelHeat().find(n->label) == m->labelHeat().end()) continue;
    //   //       counter++;
    //   //       avg += m->labelHeat()[n->label];
    //   //     }
    //   //     if(counter > 0) v->signal = avg/(double)counter;
    //   //   }
    //   //   if(m->labelHeat().find(v->label) == m->labelHeat().end()) continue;
    //   //   v->signal = m->labelHeat()[v->label];
    //   // }


    //   m->signalBounds() = calcSignalBounds2(S);
    //   m->updateTriangles();


    //   return true;
    // }
  // REGISTER_PROCESS(HeatMapCellDistanceSignal);

}
