//
// 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 "MeshProcessMeasures.hpp"

//#include <Triangulate.hpp>
//#include <ThirdParty/triangle.h>
//#include <GraphUtils.hpp>
//#include <limits>
#include "Progress.hpp"
#include "GraphUtils.hpp" // getLabel & neighborhood2D
#include "PCAnalysis.hpp" // aspect ratio
#include "Triangulate.hpp" // calcNearestPointOnBezierGrid
#include "MeshProcessHeatMap.hpp" // CombineAttrMaps

//#include <QFileDialog>
//#include <QMessageBox>

//#include <SystemProcessSave.hpp>

//#include "ui_HeatMap.h"
//#include "ui_LoadHeatMap.h"
#include <cmath>


namespace mgx 
{

  void calcSignal(Mesh* mesh, IntFloatAttr &data, QString signalType, bool signalAvg, double borderSize)
  {

  vvGraph& S = mesh->graph();

  IntFloatAttr signalBorder, signalInterior, signalTotal;
  IntFloatAttr labelBorderArea, labelInteriorArea, labelTotalArea;

  bool doborder = true;
  if(signalType == "Cell Total") doborder = false;

  if(doborder)
    mesh->updateWallGeometry(borderSize);

      // 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;
          int label = mesh->getLabel(v, n, m, mesh->parents());
          if(label <= 0)
            continue;
          //if(DoSelected and SelectedLabels.count(label) == 0)
          //  continue;

          // Total area, signal, and label center
          float area = triangleArea(v->pos, n->pos, m->pos);
          signalTotal[label] += area * (v->signal + n->signal + m->signal) / 3.0;
          labelTotalArea[label] += area;

          // 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;
            signalBorder[label] += area * bamt * (v->signal + n->signal + m->signal) / 3.0;
            labelBorderArea[label] += area * bamt;

            // Interior area and signal
            float iamt = float(3 - border) / 3.0;
            signalInterior[label] += area * iamt * (v->signal + n->signal + m->signal) / 3.0;
            labelInteriorArea[label] += area * iamt;
        }
      }
    }

      forall(const IntFloatPair& p, labelTotalArea) {
        int label = p.first;

        // Calculate label center
        float totalSignal, bordSignal = 0, intSignal = 0;

        // Set signal amounts, average over area if required
        totalSignal = signalTotal[label];
        if(doborder) {
          bordSignal = signalBorder[label];
          intSignal = signalInterior[label];
        }
        if(signalAvg) {
          signalTotal[label] = (labelTotalArea[label] == 0 ? 0 : totalSignal / labelTotalArea[label]);
          if(doborder) {
            signalBorder[label] = (labelBorderArea[label] == 0 ? 0 : bordSignal / labelBorderArea[label]);
            signalInterior[label] = (labelInteriorArea[label] == 0 ? 0 : intSignal / labelInteriorArea[label]);
          }
        }
        if(signalType == "Cell Border/Total") signalBorder[label] = signalBorder[label] / signalTotal[label];
        if(signalType == "Cell Interior/Total") signalInterior[label] = signalInterior[label] / signalTotal[label];
        if(signalType == "Cell Border/Interior") signalBorder[label] = signalBorder[label] / signalInterior[label];
      }

    if(signalType == "Cell Total"){
      data = signalTotal;
    } else if(signalType == "Cell Border" or signalType == "Cell Border/Total" or signalType == "Cell Border/Interior"){
      data = signalBorder;
    } else {
      data = signalInterior;
    }
}

bool MeasureSignal::run(Mesh* mesh, IntFloatAttr &data, QString type, bool signalAvg, double borderSize)
  {
      calcSignal(mesh, data, type, signalAvg, borderSize);
      mesh->labelHeat() = data;
      mesh->setShowLabel("Label Heat");
      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
      mesh->updateTriangles();

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Signal/Signal");
      attrData.clear();

      forall(IntFloatPair p, data)
        attrData[p.first] = p.second;

        return true;

  }
  REGISTER_PROCESS(MeasureSignal);

bool MeasureSignalBorder::run(Mesh* mesh, IntFloatAttr &data)
  {

      MeasureSignal sig(*this);
      calcSignal(mesh, data, "Cell Border", stringToBool(sig.parm("Average")), sig.parm("Border Size").toDouble());
      mesh->labelHeat() = data;
      mesh->setShowLabel("Label Heat");
      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
      mesh->updateTriangles();

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Signal/Signal Border");
      attrData.clear();

      forall(IntFloatPair p, data)
        attrData[p.first] = p.second;

        return true;

  }
  REGISTER_PROCESS(MeasureSignalBorder);

bool MeasureSignalInterior::run(Mesh* mesh, IntFloatAttr &data)
  {
      
      MeasureSignal sig(*this);
      calcSignal(mesh, data, "Cell Interior", stringToBool(sig.parm("Average")), sig.parm("Border Size").toDouble());
      mesh->labelHeat() = data;
      mesh->setShowLabel("Label Heat");
      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
      mesh->updateTriangles();

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Signal/Signal Interior");
      attrData.clear();

      forall(IntFloatPair p, data)
        attrData[p.first] = p.second;

        return true;

  }
  REGISTER_PROCESS(MeasureSignalInterior);

bool MeasureSignalTotal::run(Mesh* mesh, IntFloatAttr &data)
  {
      MeasureSignal sig(*this);
      calcSignal(mesh, data, "Cell Total", stringToBool(sig.parm("Average")), sig.parm("Border Size").toDouble());
      mesh->labelHeat() = data;
      mesh->setShowLabel("Label Heat");
      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
      mesh->updateTriangles();

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Signal/Signal Total");
      attrData.clear();

      forall(IntFloatPair p, data)
        attrData[p.first] = p.second;

        return true;

  }
  REGISTER_PROCESS(MeasureSignalTotal);

  void MeasureArea::calculateArea(vvGraph &S, IntFloatAttr &data)
  {
    data.clear();

    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;
        int label = getLabel(v,n,m);
        if(label > 0)
          data[label] += triangleArea(v->pos, n->pos, m->pos);  
      }
    }
  }

  bool MeasureArea::run(Mesh* mesh, IntFloatAttr &data)
  {

    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) 
      return false;

    vvGraph &S = mesh->graph();
    calculateArea(S, data);
    mesh->labelHeat() = data;
    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM2;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    // copy to attribute
    AttrMap<int, double> &attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Area");
    attrData.clear();

    forall(IntFloatPair p, data)
      attrData[p.first] = p.second;

    return true;

  }
  REGISTER_PROCESS(MeasureArea);

  // find cell center, find connected tris, add edges between other two vertices
  void MeasurePerimeter::calculatePerimeter(vvGraph &S, IntFloatAttr &data)
  {
    data.clear();

    forall(const vertex &v, S){
      if(v->label == -1) continue;
      forall(const vertex &n, S.neighbors(v)){
          vertex m = S.nextTo(v, n);
          if(S.edge(m,n) or S.edge(m, n))
              data[v->label] += norm(n->pos-m->pos);
      }
    }
  }
  bool MeasurePerimeter::run(Mesh* mesh, IntFloatAttr &data)
  {

    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();

    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    calculatePerimeter(T, data);
    mesh->labelHeat() = data;
    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Perimeter");
    attrData.clear();

    forall(IntFloatPair p, data)
      attrData[p.first] = p.second;

    return true;

  }
  REGISTER_PROCESS(MeasurePerimeter);

  // find consecutive border segments nm, ml; add up norm(nm ^ ml)
  void MeasureBending::calculateBending(vvGraph &S, IntFloatAttr &data){
      data.clear();

      forall(const vertex &v, S){
          if(v->label == -1) continue;
          forall(const vertex &n, S.neighbors(v)){
              vertex m = S.nextTo(v, n);
              vertex l = S.nextTo(v, m);
              Point3d nm =  m->pos - n->pos;
              nm /= norm(nm);
              Point3d ml =  l->pos - m->pos;
              ml /= norm(ml);
              float s = norm(nm ^ ml);
              float c = (nm * ml);
              if(c < 0) s = 1 + (1-s);
              data[v->label] += s;
          }
      }
  }
  bool MeasureBending::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();

    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    calculateBending(T, data);
    mesh->labelHeat() = data;

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

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/Bending");
    attrData.clear();

    forall(IntFloatPair p, data)
      attrData[p.first] = p.second;

    return true;

  }
  REGISTER_PROCESS(MeasureBending);

  // calculate bending of cell together with its closest neighbor
  void MeasureStomatalBending::calculateStomatalBending(vvGraph &T, IntFloatAttr &stomataBending){ 
      stomataBending.clear();
      //vvGraph &T = mesh->graph();
      //AttrMap<int, double>& ar = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/Common Bending");
      std::set<VertexPr> NeighbourMap;
      std::map<VertexPr, Set> NhbrElem;
      IntFloatAttr Bending;

      MeasureCommonNhbrs cmnNhbr(*this);
      cmnNhbr.calculateClosestNhbr(T, NeighbourMap); 
     
      MeasureVariabilityRadius vari(*this);
      vari.findCommonNeighbors(T, NhbrElem);

      MeasureBending bend(*this);
      bend.calculateBending(T, Bending); 
      
      forall(VertexPr Pr, NeighbourMap) {  
          double stomaBending = 0;
          if(NhbrElem[Pr].size() <= 1)
              continue;
          forall(const vertex &n, T.neighbors(Pr.first)){
              vertex m = T.nextTo(Pr.first, n);
              vertex l = T.nextTo(Pr.first, m);
              if(NhbrElem[Pr].find(m) != NhbrElem[Pr].end() and NhbrElem[Pr].find(l) != NhbrElem[Pr].end()){   
                  Point3d nm =  m->pos - n->pos;
                  nm /= norm(nm);
                  Point3d ml =  l->pos - m->pos;
                  ml /= norm(ml);
                  float s = norm(nm ^ ml);
                  float c = (nm * ml);
                  if(c < 0)
                    s = 1 + (1-s);
                  stomaBending += s;
              }
          }  
          forall(const vertex &n, T.neighbors(Pr.second)){
              vertex m = T.nextTo(Pr.second, n);
              vertex l = T.nextTo(Pr.second, m);
              if(NhbrElem[Pr].find(m) != NhbrElem[Pr].end() and NhbrElem[Pr].find(l) != NhbrElem[Pr].end()){   
                  Point3d nm =  m->pos - n->pos;
                  nm /= norm(nm);
                  Point3d ml =  l->pos - m->pos;
                  ml /= norm(ml);
                  float s = norm(nm ^ ml);
                  float c = (nm * ml);
                  if(c < 0)
                    s = 1 + (1-s);
                  stomaBending += s;
              }
          }
          stomataBending[Pr.first->label] = Bending[Pr.first->label] + Bending[Pr.second->label] - stomaBending;
          stomataBending[Pr.second->label] = stomataBending[Pr.first->label]; 
      }
  }

  bool MeasureStomatalBending::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

      calculateStomatalBending(T, data);
      mesh->labelHeat() = data;

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

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/Common Bending");
      attrData.clear();

      forall(IntFloatPair p, data)
        attrData[p.first] = p.second;

        return true;

  }
  REGISTER_PROCESS(MeasureStomatalBending);

  bool MeasureAspectRatio::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    PCAnalysis2D pc(*this);
    Point3f cf;
    cf = sqrt(3);
    pc.run(mesh, "", 0.0, 0.0, "ratio", cf, data, false);
  
    mesh->labelHeat() = data;

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

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Aspect Ratio");
    attrData.clear();

    forall(IntFloatPair p, data)
      attrData[p.first] = p.second;

      return true;

  }
  REGISTER_PROCESS(MeasureAspectRatio);
  

  bool MeasureMajorAxis::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    PCAnalysis2D pc(*this);
    Point3f cf;
    cf = sqrt(3);
    pc.run(mesh, "", 0.0, 0.0, "max", cf, data, false);
  
    mesh->labelHeat() = data;

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

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Length Major Axis");
    attrData.clear();

    forall(IntFloatPair p, data)
      attrData[p.first] = p.second;

    return true;

  }
  REGISTER_PROCESS(MeasureMajorAxis);
  

  bool MeasureMinorAxis::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    PCAnalysis2D pc(*this);
    Point3f cf;
    cf = sqrt(3);
    pc.run(mesh, "", 0.0, 0.0, "min", cf, data, false);
  
    mesh->labelHeat() = data;

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

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Length Minor Axis");
    attrData.clear();

    forall(IntFloatPair p, data)
      attrData[p.first] = p.second;


    return true;

  }
  REGISTER_PROCESS(MeasureMinorAxis);
  


  // // go through neighbors of neighbors of vertices, add them to a set if they have label >0
  void MeasureNeighbors::findNeighbourVertices(vvGraph &S, std::map<int, std::set<vertex> > &Nhbr){
      Nhbr.clear();
      //vvGraph &S = mesh->graph();
      
      forall(const vertex &v, S){
        if(v->label == -1) continue;
        forall(const vertex &n, S.neighbors(v)){
          if(n->label == -1){
            forall(const vertex &q, S.neighbors(n)){
              if(q->label != -1 and q->label != v->label and Nhbr[v->label].find(q) == Nhbr[v->label].end())
                Nhbr[v->label].insert(q);
            }
          }
        }
      }
  }

  void MeasureNeighbors::calculateNeighbors(vvGraph &S, IntFloatAttr &neighbors){
        neighbors.clear();
        //vvGraph &S = mesh->graph();
        std::map<int, std::set<vertex> > Nhbr;
        findNeighbourVertices(S, Nhbr);

        forall(const vertex &v, S){
          if(v->label == -1) continue;
          neighbors[v->label] = Nhbr[v->label].size();
        }
  }
  bool MeasureNeighbors::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();

    std::map<IntInt, double> neighMap;
    neighborhood2D(S, neighMap);


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

      // calculateNeighbors(T, data);

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Network/Neighbors");
      attrData.clear();

      std::map<int, std::set<int> > neighborNumberMap;
      forall(auto p, neighMap){
        if(p.second > 0.0001){
          neighborNumberMap[p.first.first].insert(p.first.second);
          neighborNumberMap[p.first.second].insert(p.first.first);
        }
      }

      forall(auto p, neighborNumberMap){
        attrData[p.first] = p.second.size();
        data[p.first] = attrData[p.first];
      }

      mesh->labelHeat() = data;
      mesh->setShowLabel("Label Heat");
      mesh->heatMapUnit() = "#";
      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
      mesh->updateTriangles();

      return true;
  }
  REGISTER_PROCESS(MeasureNeighbors);

  void measureRadius(vvGraph& S, QString option, IntFloatAttr &data)
  {
    data.clear();
    forall(const vertex& v, S){
      if(v->label == -1) continue;

      double maxVal = -1E20;
      double minVal = 1E20;
      double avgVal = 0.;
      double counter = 0.;
      forall(const vertex &n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
         if(m->label > -1) std::cout << "sth wrong " << m->label << std::endl;
        double dis = distLinePoint(n->pos, m->pos, v->pos, true); // min distance
        double vn = norm(n->pos-v->pos);
        double vm = norm(m->pos-v->pos);
        if(dis < minVal) minVal = dis;
        if(vn > maxVal) maxVal = vn;
        if(vm > maxVal) maxVal = vm;

        avgVal += (vn+vm)/2. * norm(n->pos-m->pos);
        counter+= norm(n->pos-m->pos);
      }
      if(option == "max"){
        data[v->label] = maxVal;
      } else if(option == "min"){
        data[v->label] = minVal;
      } else {
        data[v->label] = avgVal/counter;
      }
    }


  }

  bool MeasureMaxRadius::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    measureRadius(T, "max", data);

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Maximum Radius");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();
    return true;

  }
  REGISTER_PROCESS(MeasureMaxRadius);


  bool MeasureMinRadius::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    measureRadius(T, "min", data);

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Minimum Radius");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();
    return true;

  }
  REGISTER_PROCESS(MeasureMinRadius);

  bool MeasureAvgRadius::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();

    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    measureRadius(T, "avg", data);

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Average Radius");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();
    return true;

  }
  REGISTER_PROCESS(MeasureAvgRadius);

  void disToBezier(vvGraph& T, Bezier& b, Matrix4d mGLTot, IntFloatAttr &data, std::map<int,Point3d> &dirs, 
    std::map<int, Point2i> &idxs, int dataPointsBezier){

    std::vector<std::vector<Point3d> > bezGrid;

    b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);

    data.clear();
    forall(const vertex& v, T){
      if(v->label == -1) continue;
      Point2i idx;
      Point3d pGrid = calcNearestPointOnBezierGrid(v->pos, bezGrid, idx);
      data[v->label] = norm(pGrid - v->pos);
      dirs[v->label] = pGrid - v->pos;
      dirs[v->label] /= norm(dirs[v->label]);
      idxs[v->label] = idx;

    }

  }

  void disToBezier(AttrMap<int, Point3d>& cellCenters, Bezier& b, Matrix4d mGLTot, IntFloatAttr &data, 
    std::map<int,Point3d> &dirs, std::map<int, Point2i> &idxs, int dataPointsBezier){

    std::vector<std::vector<Point3d> > bezGrid;

    b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);

    data.clear();
    forall(auto p, cellCenters){
      int label = p.first;
      Point3d pos = p.second;
      Point2i idx;
      Point3d pGrid = calcNearestPointOnBezierGrid(pos, bezGrid, idx);
      data[label] = norm(pGrid - pos);
      dirs[label] = pGrid - pos;
      dirs[label] /= norm(dirs[label]);
      idxs[label] = idx;

    }

  }

  bool MeasureMajorAxisTheta::run(const Stack *s1, Mesh* mesh, IntFloatAttr &data, QString refDir)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    AttrMap<int, SymmetricTensor>& cusDir = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    if(refDir == "Custom X" or refDir == "Custom Y" or refDir == "Custom Z"){
      if(cusDir.empty()) return setErrorMessage("No Custom Directions found!");
    }

    Matrix4d rotMatrixS1;
    s1->getFrame().getMatrix(rotMatrixS1.data());

    PCAnalysis2D pc(*this);
    Point3f cf;
    cf = sqrt(3);
    pc.run(mesh, "", 0.0, 0.0, "", cf, data, false);

    data.clear();

    Point3d axisX(1, 0, 0);
    if(refDir == "Y-Axis"){
      axisX = Point3d(0,1,0);
    } else if(refDir == "Z-Axis"){
      axisX = Point3d(0,0,1);
    }

    forall(const vertex& v, T){
      if(v->label == -1) continue;
      axisX = multMatrix4Point3(transpose(inverse(rotMatrixS1)),axisX);

      if(refDir == "Custom X" or refDir == "Custom Y" or refDir == "Custom Z"){
        if(cusDir.find(v->label) == cusDir.end()) continue;
        if(refDir == "Custom X"){
          axisX = Point3d(cusDir[v->label].ev1());
        } else if(refDir == "Custom Y"){
          axisX = Point3d(cusDir[v->label].ev2());
        } else if(refDir == "Custom Z"){
          axisX = Point3d(cusDir[v->label].evals());
        }
      }

      SymmetricTensor& tensor = mesh->cellAxis()[v->label];
      Point3d tensorMajor = Point3d(tensor.ev1().x(), tensor.ev1().y(), tensor.ev1().z());

      data[v->label] = angleVectors(axisX, tensorMajor);
    }


    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Major Axis Theta");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = "°";
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();
    return true;

  }
  REGISTER_PROCESS(MeasureMajorAxisTheta);

  bool MeasureDistanceToBezier::run(const Stack *s1, Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph &S = mesh->graph();

    std::map<int, triVector> cellTriangles;
    generateLabelTriangleMap(S, cellTriangles);

    std::map<int, double> labArea;
    AttrMap<int, Point3d>& cellCenters = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCenters");

    mesh->updateCentersNormals();

    CuttingSurface* cutSurf = cuttingSurface();
    Matrix4d rotMatrixS1, rotMatrixCS;
    s1->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());

    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    std::map<int,Point3d> dirs;
    std::map<int,Point2i> idxs;

    Bezier b = cutSurf->bezier();

    disToBezier(cellCenters, b, mGLTot, data, dirs, idxs, 200);

    if(data.empty()) return false;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Distance To Bezier");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;

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

    return true;

  }
  REGISTER_PROCESS(MeasureDistanceToBezier);


  bool MeasureBezierCoord::run(const Stack *s1, Mesh* mesh, QString bezCoord, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

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

    std::map<int, triVector> cellTriangles;
    generateLabelTriangleMap(S, cellTriangles);

    std::map<int, double> labArea;
    AttrMap<int, Point3d>& cellCenters = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCenters");

    mesh->updateCentersNormals();

    CuttingSurface* cutSurf = cuttingSurface();
    Matrix4d rotMatrixS1, rotMatrixCS;

    s1->getFrame().getMatrix(rotMatrixS1.data());

    cutSurf->frame().getMatrix(rotMatrixCS.data());

    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    std::map<int,Point3d> dirs;
    std::map<int,Point2i> idxs;

    Bezier b = cutSurf->bezier();
    int dataPointsBezier = 1000.0;
    disToBezier(cellCenters, b, mGLTot, data, dirs, idxs, dataPointsBezier);

    if(data.empty()) return false;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Bezier Grid Coord");
    attrData.clear();

    forall(auto p, cellCenters){
      int label = p.first;
      if(label < 1) continue;
      if(idxs.find(label) == idxs.end()) continue;
      if(bezCoord == "X") attrData[label] = idxs[label].x() / (double)dataPointsBezier;
      else attrData[label] = idxs[label].y() / (double)dataPointsBezier;

      mesh->labelHeat()[label] = attrData[label];
    }

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = "Bezier";
    mesh->heatMapBounds() = Point2f(0.0,1.0);//mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    return true;

  }
  REGISTER_PROCESS(MeasureBezierCoord);


  bool MeasureBezierLineCoord::run(const Stack *s1, Mesh* mesh, bool reverse, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph &S = mesh->graph();

    std::map<int, triVector> cellTriangles;
    generateLabelTriangleMap(S, cellTriangles);

    std::map<int, double> labArea;
    AttrMap<int, Point3d>& cellCenters = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCenters");

    mesh->updateCentersNormals();

    CuttingSurface* cutSurf = cuttingSurface();
    Matrix4d rotMatrixS1, rotMatrixCS;

    s1->getFrame().getMatrix(rotMatrixS1.data());

    cutSurf->frame().getMatrix(rotMatrixCS.data());

    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    std::map<int,Point3d> dirs;
    std::map<int,Point2i> idxs;

    Bezier b = cutSurf->bezier();

    std::vector<std::vector<Point3d> > bezGrid;
    std::vector<Point3d> discretizedBez, differentialBez;
    double totalLength;

    int dataPointsBezier = 200;
    b.discretizeLineEqualWeight(dataPointsBezier, 0.5, mGLTot, discretizedBez, differentialBez, totalLength);

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Bezier Line Coord");
    attrData.clear();

    forall(auto p, cellCenters){
      int label = p.first;
      if(label < 1) continue;

      Point3d pos = p.second;
      
      int bezIdx;
      double bezWeight;
      Point3d bezInterpPoint = calcNearestPointOnBezierLine(pos, discretizedBez, bezIdx, bezWeight);

      attrData[label] = (bezIdx + bezWeight) * totalLength/dataPointsBezier;
      if(reverse) attrData[label] = totalLength - attrData[label];

      mesh->labelHeat()[label] = attrData[label];
    }

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

    return true;

  }
  REGISTER_PROCESS(MeasureBezierLineCoord);



  bool MeasurePolarCoord::run(const Stack *s1, Mesh* mesh, QString axis, QString coord, IntFloatAttr& heatMap)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    Matrix4d rotMatrixS1;
    s1->getFrame().getMatrix(rotMatrixS1.data());

    // fix the rotations
    Matrix4d mGLRotMat = transpose(rotMatrixS1);

    AttrMap<int, double>& attrLong = mesh->attributes().attrMap<int, double>("Measure Label Double Location/PolarCoordLong");
    attrLong.clear();
    AttrMap<int, double>& attrRad = mesh->attributes().attrMap<int, double>("Measure Label Double Location/PolarCoordRad");
    attrRad.clear();
    AttrMap<int, double>& attrCirc = mesh->attributes().attrMap<int, double>("Measure Label Double Location/PolarCoordCirc");
    attrCirc.clear();

    Point3d centralAxis(0,0,1);
    Point3d refDir(1,0,0);
    bool zAxis = true, yAxis = false, xAxis = false;
    if(axis == "X"){
      centralAxis = Point3d(1,0,0);
      refDir = Point3d(0,0,1);
      zAxis = false;
      yAxis = false;
      xAxis = true;
    } else if(axis == "Y"){
      centralAxis = Point3d(0,1,0);
      refDir = Point3d(1,0,0);
      zAxis = false;
      yAxis = true;
      xAxis = false;
    }

    mesh->updateCentersNormals();

    double radConversion = 180./3.14159;
    Point3d refDir2 = centralAxis ^ refDir;

    // cell centers
    IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());
    Point3d origin(0,0,0);

    forall(auto p, centers){

      int cellLabel = p.first;
      Point3d cellCenter = Point3d(p.second);

      cellCenter = multMatrix4Point3(mGLRotMat,cellCenter); // correct rotations
      // calculate angle
      Point3d pointPlane = projectPointOnPlane(cellCenter, origin, centralAxis);

      double dis = 0;
      if(zAxis) dis = pointPlane.x()*pointPlane.x() + pointPlane.y()*pointPlane.y();
      else if(yAxis) dis = pointPlane.x()*pointPlane.x() + pointPlane.z()*pointPlane.z();
      else if(xAxis) dis = pointPlane.z()*pointPlane.z() + pointPlane.y()*pointPlane.y();
      dis = sqrt(dis);

      attrRad[cellLabel] = dis;

      Point3d v1 = pointPlane;
      double scalar = v1 * refDir;
      double angle = acos(scalar / norm(v1))*radConversion;
      if(v1 * refDir2 < 0) angle = 360-angle;

      attrCirc[cellLabel] = angle;

      attrLong[cellLabel] = cellCenter.z();
      if(axis == "X"){
        attrLong[cellLabel] = cellCenter.x();
      } else if(axis == "Y"){
        attrLong[cellLabel] = cellCenter.y();
      }

      if(coord == "Longitudinal")
        mesh->labelHeat()[cellLabel] = attrLong[cellLabel];
      else if(coord == "Radial")
        mesh->labelHeat()[cellLabel] = attrRad[cellLabel];
      else
        mesh->labelHeat()[cellLabel] = attrCirc[cellLabel];

    }
    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = "";
    if(coord == "Longitudinal" or coord == "Radial")
      mesh->heatMapUnit() = UM;
    else
      mesh->heatMapUnit() = "°";
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    return true;

  }
  REGISTER_PROCESS(MeasurePolarCoord);



  bool MeasureDistanceToMesh::run(const Stack *s1, const Stack *s2, Mesh* mesh, Mesh* m2, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph &S = mesh->graph();

    std::map<int, triVector> cellTriangles;
    generateLabelTriangleMap(S, cellTriangles);

    std::map<int, double> labArea;
    AttrMap<int, Point3d>& cellCenters = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCenters");

    mesh->updateCentersNormals();
    Matrix4d rotMatrixS1, rotMatrixS2;
    s1->getFrame().getMatrix(rotMatrixS1.data());
    s2->getFrame().getMatrix(rotMatrixS2.data());

    Matrix4d mGLTot = transpose(inverse(rotMatrixS2)) * transpose(rotMatrixS1);

    if(m2->graph().empty()) return false;

    forall(auto p, cellCenters){
      int label = p.first;
      Point3d pos = p.second;
      Point3d vCorrected = multMatrix4Point3(mGLTot,pos);
      Point3d pM = nearestPointOnMesh(vCorrected, m2->graph());
      data[label] = norm(pM - vCorrected);
    }

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Distance To Mesh");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();
    return true;

  }
  REGISTER_PROCESS(MeasureDistanceToMesh);


  // create a map of common neighbors of two neighboring cells
  void MeasureVariabilityRadius::findCommonNeighbors(vvGraph &T, std::map<VertexPr, Set> &NhbrElem){
    NhbrElem.clear();
    //vvGraph &T = mesh->graph();
    std::map<int, std::set<vertex> > Nhbr;
    MeasureNeighbors neighbor(*this);
    neighbor.findNeighbourVertices(T, Nhbr);
    forall(const vertex &v, T) {
      if(v->label == -1) continue;
      forall(const vertex &n, Nhbr[v->label]){
        if(n->label == -1) continue;
        if(NhbrElem.find(std::make_pair(v, n)) != NhbrElem.end() or NhbrElem.find(std::make_pair(n, v)) != NhbrElem.end()) continue;

        forall(const vertex &x, T.neighbors(v)){    //Iterate through the neighbors of v
          forall(const vertex &y, T.neighbors(n)){    //Iterate through the neighbors of n
            if(x == y and x != n and y != v){           
              NhbrElem[std::make_pair(v, n)].insert(x);
              NhbrElem[std::make_pair(v, n)].insert(y);
            }
          }
        }
        if(NhbrElem[std::make_pair(v, n)].size() > 1){
          NhbrElem[std::make_pair(v, n)] = NhbrElem[std::make_pair(v, n)];
        }
      } 
    } 
  }

  // calculate the variability radius for each of the cell pairs, pick the lowest value
  void MeasureVariabilityRadius::calculateVariabilityRadius(vvGraph &T, IntFloatAttr &variabilityRadius){  
      variabilityRadius.clear();
      //vvGraph &T = mesh->graph();
      //AttrMap<int, double>& ar = mesh->attributes().attrMap<int, double>("Measure Label Double /Shape/Variability Radius");
      IntFloatAttr Nhbr;
      MeasureNeighbors neighbor(*this);
      neighbor.calculateNeighbors(T, Nhbr);
      std::map<VertexPr, Set> NhbrElem;
      std::map<VertexPr, Point3d> meanNhbr;
      findCommonNeighbors(T, NhbrElem);
      forall(VertexPrSet CNhbr, NhbrElem){   //Calculate center of the border of neighboring vertices
        meanNhbr[CNhbr.first] = 0;   //cc.NhbrElem contains common neighbors (border vertices)
        forall(const vertex &elem, CNhbr.second)  //Iterate through neighbors
            meanNhbr[CNhbr.first] += elem->pos;  //Add positions
        meanNhbr[CNhbr.first] /= CNhbr.second.size();  //Take average of positions
        variabilityRadius[CNhbr.first.first->label] = HUGE_VAL;  
        variabilityRadius[CNhbr.first.second->label] = HUGE_VAL;
      }
      forall(VertexPrSet CNhbr, NhbrElem){
        double stdev = 0, countele = 0;  
        forall(const vertex &elem, T.neighbors(CNhbr.first.first)){
          if(CNhbr.second.find(elem) == CNhbr.second.end()){
            stdev += pow((norm(meanNhbr[CNhbr.first] - elem->pos)),2);
            countele++; 
          }
        }
        forall(const vertex &elem, T.neighbors(CNhbr.first.second)){
          if(CNhbr.second.find(elem) == CNhbr.second.end()){
            stdev += pow((norm(meanNhbr[CNhbr.first] - elem->pos)),2);
            countele++; 
          }
        }
        stdev /= countele;
        stdev = pow(stdev, 0.5); 
        if(stdev < variabilityRadius[CNhbr.first.first->label])
            variabilityRadius[CNhbr.first.first->label] = stdev;
        if(stdev < variabilityRadius[CNhbr.first.second->label])
            variabilityRadius[CNhbr.first.second->label] = stdev;
      }
  }

  bool MeasureVariabilityRadius::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    calculateVariabilityRadius(T, data);

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/Variability Radius");
    attrData.clear();

    forall(IntFloatPair p, data){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = data;
    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = UM;
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeasureVariabilityRadius);

  // 
  void MeasureCommonNhbrs::calculateClosestNhbr(vvGraph &T, std::set<VertexPr> &NeighbourMap){
      NeighbourMap.clear();
      //vvGraph &T = mesh->graph();
      //Compute closest neighbor map
      std::set<vertex> computed;
      std::map<vertex, vertex> NhbrPr;
      
      std::map<VertexPr, Set> NhbrElem;
      std::map<int, std::set<vertex> > Nhbr;
      MeasureNeighbors neighbor(*this);
      neighbor.findNeighbourVertices(T, Nhbr);

      IntFloatAttr variabilityRadius;
      MeasureVariabilityRadius vari(*this);

      vari.calculateVariabilityRadius(T, variabilityRadius); 
      vari.findCommonNeighbors(T, NhbrElem);    
      forall(const vertex &v, T){
        forall(const vertex &n, Nhbr[v->label]){  
          if(computed.find(v) == computed.end() and variabilityRadius[v->label] == variabilityRadius[n->label] and NhbrElem[std::make_pair(v, n)].size() > 1 ){
            computed.insert(v);
            computed.insert(n);
            NhbrPr[v] = n;
            NhbrPr[n] = v;
          } 
        }
      }
      forall(VertexPr Pr, NhbrPr)
          NeighbourMap.insert(std::make_pair(Pr.first, Pr.second));
  }

  void MeasureCommonNhbrs::calculateCommonNhbrs(vvGraph &T, IntFloatAttr &commonNeighbors){
      commonNeighbors.clear();
      //vvGraph &T = mesh->graph();
      //AttrMap<int, double>& ar = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/Common Neighbors");
      std::map<int, std::set<vertex> > Nhbr;
      IntFloatAttr Neighbors;
      MeasureNeighbors neighbor(*this);
      neighbor.findNeighbourVertices(T, Nhbr);
      neighbor.calculateNeighbors(T, Neighbors);  

      std::set<VertexPr> NeighbourMap;
      calculateClosestNhbr(T, NeighbourMap);
      forall(VertexPr Pr, NeighbourMap){
          int RepeatedNeighbors = 0;
          forall(const vertex &elem, Nhbr[Pr.first->label]){
              if(Nhbr[Pr.second->label].find(elem) != Nhbr[Pr.second->label].end()){
                  RepeatedNeighbors++; 
              }
          }
          commonNeighbors[Pr.first->label] = Neighbors[Pr.first->label] + Neighbors[Pr.second->label] - RepeatedNeighbors - 2;
          commonNeighbors[Pr.second->label] = Neighbors[Pr.first->label] + Neighbors[Pr.second->label] - RepeatedNeighbors- 2;
      }
  }

  bool MeasureCommonNhbrs::run(Mesh* mesh, IntFloatAttr &data)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

      calculateCommonNhbrs(T, data);

      AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/CommonNeighbors");
      attrData.clear();

      forall(IntFloatPair p, data){
        attrData[p.first] = p.second;
      }

      mesh->labelHeat() = data;

      mesh->setShowLabel("Label Heat");
      mesh->heatMapUnit() = "";
      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
      mesh->updateTriangles();
      return true;
  }
  REGISTER_PROCESS(MeasureCommonNhbrs);

  void MeasureStomataArea::calculateStomataArea(vvGraph &T, IntFloatAttr &stomataArea){
      stomataArea.clear();
      //vvGraph &T = mesh->graph();
      
      IntFloatAttr area, perimeter;
      MeasureArea ar(*this);
      ar.calculateArea(T, area);
      MeasurePerimeter peri(*this);
      peri.calculatePerimeter(T, perimeter);
      std::map<int, std::set<vertex> > Nhbr;
      MeasureNeighbors neighbor(*this);
      neighbor.findNeighbourVertices(T, Nhbr);
      forall(const vertex &v, T){
        int flag = 0, max = 0;
        vertex x;
        float totArea = 0, totPeri = 0;
        if(stomataArea[v->label] != 0)
            continue;
        forall(const vertex &k, Nhbr[v->label]){ // set flag if area difference is <20  
          if(abs(area[v->label] - area[k->label]) < 20){
            flag = 1;
            x = k;
            break;
          } 
        }
        if(flag == 1){
          totArea = area[v->label] + area[x->label];
          std::set <vertex> vSet;
          forall(const vertex &m, T.neighbors(v)){ 
            forall(const vertex &n, T.neighbors(x)){
              if(n == m){
                vSet.insert(n); // common border
              }
            }
          }
          forall(const vertex &p, vSet){
            forall(const vertex &q, vSet){
              if(norm(p->pos - q->pos) > max)
                max = norm(p->pos - q->pos); // max distance between two vertices of the common border
            }
          }
          totPeri = pow(((perimeter[v->label]-max) + (perimeter[x->label] - max)), 2);
          if(1 - (((totArea/totPeri) * 4 * 3.14)) < 0.25){
            stomataArea[v->label] = totArea/totPeri;
            stomataArea[x->label] = totArea/totPeri;
          }
        }
      }
  }

  bool MeasureStomataArea::run(Mesh* mesh, IntFloatAttr &stomataArea)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    calculateStomataArea(T, stomataArea);

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Shape/Common Bending");
    attrData.clear();

    forall(IntFloatPair p, stomataArea){
      attrData[p.first] = p.second;
    }

    mesh->labelHeat() = stomataArea;
    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = "";
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(MeasureStomataArea);

  bool MeasureNeighborhoodArea::run(Mesh* mesh, IntFloatAttr &avgArea)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Neighborhood/Area");
    attrData.clear();

    CombineAttrMaps manipulateHeat(*this);
    manipulateHeat.run(mesh, "/Geometry/Area", "", stringToBool("No"), "", stringToBool("Yes"), "Neighborhood", "Mean", -1, 0, 100000, "Neighborhood/Area");
    
    forall(IntDoublePair p, attrData){
      mesh->labelHeat()[p.first] = p.second;
    }

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

    return true;
  }
  REGISTER_PROCESS(MeasureNeighborhoodArea);

  bool MeasureNeighborhoodPerimeter::run(Mesh* mesh, IntFloatAttr &avgPerimeter)
  {

    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Neighborhood/Perimeter");
    attrData.clear();

    CombineAttrMaps manipulateHeat(*this);
    manipulateHeat.run(mesh, "/Geometry/Area", "", stringToBool("No"), "", stringToBool("Yes"), "Neighborhood", "Mean", -1, 0, 100000, "Neighborhood/Perimeter");

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

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

    return true;
  }
  REGISTER_PROCESS(MeasureNeighborhoodPerimeter);

  bool MeasureNeighborhoodAspectRatio::run(Mesh* mesh, IntFloatAttr &avgLenBr)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Neighborhood/Aspect Ratio");
    attrData.clear();

    CombineAttrMaps manipulateHeat(*this);
    manipulateHeat.run(mesh, "/Geometry/Area", "", stringToBool("No"), "", stringToBool("Yes"), "Neighborhood", "Mean", -1, 0, 100000, "Neighborhood/Aspect Ratio");

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

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

    return true;
  }
  REGISTER_PROCESS(MeasureNeighborhoodAspectRatio);

  bool MeasureNeighborhoodNeighbors::run(Mesh* mesh, IntFloatAttr &avgNeighbors)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Neighborhood/Neighbors");
    attrData.clear();

    CombineAttrMaps manipulateHeat(*this);
    manipulateHeat.run(mesh, "Network/Neighbors", "", stringToBool("No"), "", stringToBool("Yes"), "Neighborhood", "Mean", -1, 0, 100000, "Neighborhood/Neighbors");

    forall(IntDoublePair p, attrData){
      mesh->labelHeat()[p.first] = p.second;
    }
    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = "";
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(MeasureNeighborhoodNeighbors);

  bool MeasureNeighborhoodVariabilityRadius::run(Mesh* mesh, IntFloatAttr &avgVariabilityRadius)
  {
    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph T;
    vvGraph &S = mesh->graph();


    if(mesh->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Neighborhood/Variability Radius");
    attrData.clear();

    CombineAttrMaps manipulateHeat(*this);
    manipulateHeat.run(mesh, "Geometry/Area", "", stringToBool("No"), "", stringToBool("Yes"), "Neighborhood", "Mean", -1, 0, 100000, "Neighborhood/Variability Radius");

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

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

    return true;
  }
  REGISTER_PROCESS(MeasureNeighborhoodVariabilityRadius);


  bool MeasureCellDistance::run(Mesh* mesh, QString weight, bool cellTypes, IntFloatAttr &heatmap)
  {

    if(mesh->meshType() != "MGXC" and mesh->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(mesh->graph().empty()) return false;

    vvGraph &S = mesh->graph();

    std::set<int> selectedCells = findAllSelectedLabels(S);
    std::set<int> allCells = findAllLabelsSet(S);

    AttrMap<int, double>& attr = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Cell Distance");
    
    std::map<IntInt, double> neighMap;
    neighborhood2D(S, neighMap, weight);

    typedef std::pair<IntInt, double> IntIntDoubleP;
    
    if(cellTypes){
      std::map<IntInt, double> neighMapCellTypes;

      forall(IntIntDoubleP p, neighMap){
        if(mesh->parents()[p.first.first] != mesh->parents()[p.first.second]) continue;
        neighMapCellTypes[p.first] = p.second;
      }

      neighMap = neighMapCellTypes;

    } else if(mesh->useParents()){
      std::set<int> selectedCellsParents;

      forall(int l, selectedCells){
        if(mesh->parents()[l] > 0)
          selectedCellsParents.insert(mesh->parents()[l]);
      }
      selectedCells = selectedCellsParents;

      std::set<int> allCellsParents;

      forall(int l, allCells){
        allCellsParents.insert(mesh->parents()[l]);
      }
      allCells = allCellsParents;

      // convert to a parent neighbor map
      std::map<IntInt, double> neighMapParents;
      forall(IntIntDoubleP p, neighMap){
        if(mesh->parents()[p.first.first] == mesh->parents()[p.first.second]) continue;
        std::pair<int,int> newP = std::make_pair(mesh->parents()[p.first.first], mesh->parents()[p.first.second]);
        neighMapParents[newP] += p.second;
      }
      neighMap = neighMapParents;
    }

    if(selectedCells.empty()) return setErrorMessage("No cells selected!");

    bool equalWeights = weight == "1" ? true : false;

    std::map<int, double> heatDijk = dijkstra(allCells, neighMap, selectedCells, equalWeights, HUGE_VAL, 0.0001);

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    mesh->labelHeat().clear();

    forall(auto p, heatDijk){
    //forall(const vertex& v, S){
      int currentLabel = p.first;//v->label;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents() and !cellTypes) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      if(!std::isfinite(heatDijk[currentLabel])) continue; // check if NAN or INF
      mesh->labelHeat()[currentLabel] = p.second;//heatDijk[currentLabel];
      attr[currentLabel] = p.second;// heatDijk[currentLabel];
      if(maxValue < p.second) maxValue = p.second;//heatDijk[currentLabel];
      if(minValue > p.second) minValue = p.second;//heatDijk[currentLabel];
    }

    mesh->setShowLabel("Label Heat");
    mesh->heatMapUnit() = "";
    if(weight == "Euclidean") mesh->heatMapUnit() = UM;
    if(weight == "1") mesh->heatMapUnit() = "cells";
    if(weight == "1 / Wall Length") mesh->heatMapUnit() = "";
    mesh->heatMapBounds() = mesh->calcHeatMapBounds();
    mesh->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(MeasureCellDistance);

  bool MeasureCellCoord::run(Stack* s, Mesh* mesh, QString dim, IntFloatAttr &heatmap)
  {

    if(mesh->graph().empty()) return false;

    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Location/Cell Coordinate");
    attrData.clear();
    mesh->labelHeat().clear();

    mesh->updateCentersNormals();
    IntPoint3fAttr &centers = (mesh->useParents() ? mesh->parentCenterVis() : mesh->labelCenterVis());

    Matrix4d rotMatrixS1;
    s->getFrame().getMatrix(rotMatrixS1.data());

    // fix the rotations
    Matrix4d mGLTot = transpose(rotMatrixS1);

    forall(auto p, centers){
      int l = p.first;
      Point3d c = Point3d(p.second);
      Point3d centroidRotated = multMatrix4Point3(mGLTot,c); // correct rotations
      double dis = 0;
      if(dim == "X"){
        dis = centroidRotated.x();
      } else if(dim == "Y"){
        dis = centroidRotated.y();
      } else if(dim == "Z"){
        dis = centroidRotated.z();
      } else {
        dis = norm(centroidRotated);
      }
      attrData[l] = dis;
      mesh->labelHeat()[l] = dis;
    }

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

    return true;
  }
  REGISTER_PROCESS(MeasureCellCoord);


 // returns the labels of all neighbors of a vertex
 std::set<int> getNeighborLabels(vvGraph& S, const vertex& v)
 {
    std::set<int> neighbors;
    forall(const vertex& n, S.neighbors(v)){
      if(n->label > -1) neighbors.insert(n->label);
    }
    return neighbors;
 }


 // // counts the number of non-border vertices of a vertex
 // // in a cell mesh this is the number of different neighbors
 // int countNeighbors(vvGraph& S, const vertex& v)
 // {
 //    int neighborCount = 0;
 //    forall(const vertex& n, S.neighbors(v)){
 //      if(n->label > -1) neighborCount++;
 //    }
 //    return neighborCount;
 // }

  // // returns the labels of two common neighbors of two sets of neighbor labels
  // // input:  two sets of int labels
  // // output: pair of int labels that occur in both sets, return (-1,-1) if there are more or less than 2 common neighbors
  // std::pair<int,int> findCommonNeighbors(std::set<int> v1, std::set<int> v2)
  // {
  //   std::pair<int,int> common = std::make_pair(-1,-1);
    
  //   std::vector<int> commonV;

  //   forall(const int& n, v1){
  //     if(v2.find(n) == v2.end()) continue;
  //     commonV.push_back(n);
  //   }

  //   if(commonV.size() != 2) return common;

  //   return std::make_pair(commonV[0], commonV[1]);
    
  // }


 bool JunctionDistance::run(Mesh *m, QString mode, bool onlyDirect, bool ignoreBorderCells, IntFloatAttr& heatMap)
  {

    const QString meshType = m->tissue().meshType();
    if(meshType != "MGX2D" and meshType != "MGXM" and meshType != "MGXC") return setErrorMessage(QString("Mesh Type must be MGX2D, MGXM or MGXC"));

    vvGraph& S = m->graph();

    AttrMap<int, double>& attr = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Junction Distance");
    attr.clear();

    std::map<IntInt, double> neighMap;

    neighborhood2D(S, neighMap, "");

    std::map<int, double> labelHeatValueMap;

    std::set<int> allLabels = findAllLabelsSet(S);

    forall(const int& l, allLabels){
      if(mode == "Min") labelHeatValueMap[l] = 1E20;
        else labelHeatValueMap[l] = -1E20;
    }

    std::set<int> borderCells;
    forall(const vertex& v, S){
      if(!v->margin) continue;
      forall(const vertex& n, S.neighbors(v)){
        if(n->label == -1) continue;
        borderCells.insert(n->label);
      }
    }

    forall(auto p, neighMap){

      //std::cout << "c " << p.first.first << "/" << p.first.second << "//" << p.second << "//" << neighMap[p.first] << std::endl;

      if(mode == "Min"){
        labelHeatValueMap[p.first.first] = min(labelHeatValueMap[p.first.first],p.second);
        labelHeatValueMap[p.first.second] = min(labelHeatValueMap[p.first.second],p.second);
      } else if(mode == "Max"){
        labelHeatValueMap[p.first.first] = max(labelHeatValueMap[p.first.first],p.second);
        labelHeatValueMap[p.first.second] = max(labelHeatValueMap[p.first.second],p.second);
      }

    }

    if(!onlyDirect){
      std::map<int, std::set<int> > labelNeighborsMap;
      forall(auto p, neighMap){
        labelNeighborsMap[p.first.first].insert(p.first.second);
        labelNeighborsMap[p.first.second].insert(p.first.first);
      }

      forall(auto p, labelNeighborsMap){
        std::vector<std::pair<int, int> > cellNeighPairs;
        for(auto it1 = p.second.begin(); it1 != p.second.end(); it1++){
          for(auto it2 = it1; it2 != p.second.end(); it2++){
            if(it2 == it1) continue;
            std::pair<int, int> neighP = std::make_pair(*it1,*it2);
            cellNeighPairs.push_back(neighP);
          }
        }

        forall(auto nP, cellNeighPairs){
          if(neighMap.find(nP) == neighMap.end()) continue;
          if(neighMap[nP] < labelHeatValueMap[p.first] and mode == "Min") labelHeatValueMap[p.first] = neighMap[nP];
          if(neighMap[nP] > labelHeatValueMap[p.first] and mode == "Max") labelHeatValueMap[p.first] = neighMap[nP];
        }
      }

    }

    m->labelHeat().clear();
    forall(auto p, labelHeatValueMap){
      if(p.second == 1E20 or p.second == -1E20) continue;
      if(ignoreBorderCells and borderCells.find(p.first) != borderCells.end()) continue;
      m->labelHeat()[p.first] = p.second;
    }

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

    m->setShowLabel("Label Heat");
    m->heatMapUnit() = UM;
    m->heatMapBounds() = m->calcHeatMapBounds();
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(JunctionDistance);


//  bool CellAxisAspectRatio::run(Mesh* mesh, IntFloatAttr &heatmap)
//  {
//      vvGraph &S = mesh->graph();
//      
//
//      forall(const IntSymTensorPair &p, mesh->cellAxis()) {
//      
//        const SymmetricTensor& tensor = p.second;//mesh->cellAxis()[p.first];
//
//        double ratio = 0;
//        if(tensor.evals()[1] > 0) ratio = tensor.evals()[0]/tensor.evals()[1];
//        if(ratio < 1) ratio = 1;
//
//        mesh->labelHeat()[p.first] = ratio;
//
//      }
//
//      //mesh->labelHeat() = avgCommonNeighbors;
//
//      mesh->setShowLabel("Label Heat");
//      mesh->heatMapUnit() = "ratio";
//      mesh->heatMapBounds() = mesh->calcHeatMapBounds();
//      mesh->updateTriangles();
//
//      heatmap = mesh->labelHeat();
//
//      return true;
//  }
//  REGISTER_PROCESS(CellAxisAspectRatio);

}
