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

#include "Progress.hpp"

#include <Triangulate.hpp>
#include <GraphUtils.hpp>

namespace mgx
{

  bool Measure3DVolume::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Volume");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DVolume);


  bool Measure3DNeighbors::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Network/Neighbors");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DNeighbors);


  bool Measure3DWallArea::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DWallArea);

  bool Measure3DCoordLong::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Longitudinal");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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


    return true;
  }
  REGISTER_PROCESS(Measure3DCoordLong);

  bool Measure3DCoordCirc::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Circumferential");
    if(attrData.size() == 0)
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DCoordCirc);

  bool Measure3DCoordRad::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Radial");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DCoordRad);

  bool Measure3DSizeLong::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Longitudinal");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeLong);

  bool Measure3DSizeCirc::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Circumferential");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeCirc);

  bool Measure3DSizeRad::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Radial");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeRad);

  bool Measure3DSizeX::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Length X");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeX);

  bool Measure3DSizeY::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Length Y");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeY);

  bool Measure3DSizeZ::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Length Z");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeZ);

  bool Measure3DSizeCustom::run(Mesh* mesh, QString dim, QString attr1, QString attr2, QString attr3, IntFloatAttr &data)
  {

    AttrMap<int, double>& attrData1 = mesh->attributes().attrMap<int, double>("Measure Label Double " + attr1);
    attrData1.clear();
    AttrMap<int, double>& attrData2 = mesh->attributes().attrMap<int, double>("Measure Label Double " + attr2);
    attrData2.clear();
    AttrMap<int, double>& attrData3 = mesh->attributes().attrMap<int, double>("Measure Label Double " + attr3);
    attrData3.clear();

    std::map<int, triVector> cellTriangles;

    vvGraph& S = mesh->graph();
    generateLabelTriangleMap(S, cellTriangles);

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

    if(customDirAttr.empty()){
      return setErrorMessage("No valid custom directions found!");
    }

    forall(auto p, customDirAttr){
      std::cout << "customDir  " << p.first << "/" << p.second.ev1() << std::endl;
    }

    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    if(centroids.empty()){
      return setErrorMessage("No valid cell centroids found!");
    }

    forall(auto p, centroids){
      std::cout << "centr  " << p.first << "/" << p.second << std::endl;
    }

    typedef std::pair<int, triVector> labelTriP;

    data.clear();

    forall(labelTriP p, cellTriangles){
      int l = p.first; // label
      if(centroids.find(l) == centroids.end() or customDirAttr.find(l) == customDirAttr.end()) continue;

      Point3d customDirX = Point3d(customDirAttr[l].ev1());
      Point3d customDirY = Point3d(customDirAttr[l].ev2());
      Point3d customDirZ = Point3d(customDirAttr[l].evals());

      Point3d customSel;

      double sizeX, sizeY, sizeZ;

      sizeX = estimateCellLength(l, centroids[l], customDirX, p.second);
      sizeY = estimateCellLength(l, centroids[l], customDirY, p.second);
      sizeZ = estimateCellLength(l, centroids[l], customDirZ, p.second);

      if(dim == "X"){
        customSel = customDirX;
        data[l] = sizeX;
      }
      if(dim == "Y"){
        customSel = customDirY;
        data[l] = sizeY;
      }
      if(dim == "Z"){
        customSel = customDirZ;
        data[l] = sizeZ;
      }

      if(norm(customSel - Point3d(0,0,0)) < 0.001) data.erase(l);

      attrData1[l] = sizeX;
      attrData2[l] = sizeY;
      attrData3[l] = sizeZ;

    }

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

    return true;
  }
  REGISTER_PROCESS(Measure3DSizeCustom);

  bool Measure3DVolumeSurfaceRatio::run(Mesh* mesh, IntFloatAttr &data, double powerVol, double powerWall)
  {
    AttrMap<int, double>& attrDataVol = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Volume");
    if(attrDataVol.size() == 0)
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    AttrMap<int, double>& attrDataWall = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");
    if(attrDataWall.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

     data.clear();

    forall(IntDouble p, attrDataVol){
      data[p.first] = std::pow(p.second,1./powerVol)/std::pow(attrDataWall[p.first], 1./powerWall);
    }

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

    return true;
  }
  REGISTER_PROCESS(Measure3DVolumeSurfaceRatio);

  bool Measure3DOutsideWallArea::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Outside Wall Area");
    if(attrData.size() == 0) 
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    typedef std::pair<int,double> IntDouble;

     data.clear();

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

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

    return true;
  }
  REGISTER_PROCESS(Measure3DOutsideWallArea);

  bool Measure3DOutsideWallAreaPercent::run(Mesh* mesh, IntFloatAttr &data)
  {
    AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Outside Wall Area");
    if(attrData.size() == 0)
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    AttrMap<int, double>& attrDataCW = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");
    if(attrDataCW.size() == 0)
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    AttrMap<int, double>& attrDataRatio = mesh->attributes().attrMap<int, double>("Measure Label Double Geometry/Outside Wall Area Ratio");

    typedef std::pair<int,double> IntDouble;

    data.clear();
    attrDataRatio.clear();

    forall(IntDouble p, attrData){
      data[p.first] = p.second/attrDataCW[p.first];
      attrDataRatio[p.first] = data[p.first];
    }

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

    return true;
  }
  REGISTER_PROCESS(Measure3DOutsideWallAreaPercent);


  bool Measure3DCellDistance::run(Mesh* mesh, QString weight, bool cellTypes, IntFloatAttr &data)
  {
    vvGraph& S = mesh->graph();

    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

    // create cell graph
    NhbdGraphInfo info;
    neighborhoodGraph(S, 0.0001, info); // now in graphutils

    std::map<IntInt, double> neighMap = info.sharedWallArea;

    // find selected cells
    std::set<int> selectedCells = findAllSelectedLabels(S);
    std::set<int> allCells = findAllLabelsSet(S);

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

      forall(auto 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(auto 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;
    }

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

    if(weight == "Euclidean"){
      std::map<IntInt, double> neighMapNew;
      forall(auto p, neighMap){
        neighMapNew[p.first] = 1./norm(centroids[p.first.first] - centroids[p.first.second]);
      }
      neighMap = neighMapNew;
    }


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

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    // update heat map
    mesh->labelHeat().clear();
    forall(const vertex& v, S){
      int currentLabel = v->label;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents()) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      mesh->labelHeat()[currentLabel] = heatDijk[currentLabel];
      attrData[currentLabel] = heatDijk[currentLabel];
      if(maxValue < heatDijk[currentLabel]) maxValue = heatDijk[currentLabel];
      if(minValue > heatDijk[currentLabel]) minValue = heatDijk[currentLabel];
    }

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

    return true;
  }
  REGISTER_PROCESS(Measure3DCellDistance);


  bool Measure3DCellDistanceBezierRing::run(Mesh* mesh, QString weight, bool cellTypes, double thresholdDis, bool selectedAsDirect, 
      bool considerOrientation, bool selectCells, bool keep, IntFloatAttr& heatMap)
  {
    vvGraph& S = mesh->graph();

    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

    MeasureDistanceToBezier3D proc(*this);

    IntFloatAttr heat;
    Stack* s1 = currentStack();
    proc.run(s1, mesh, false, false, heat);

    std::set<int> alreadySelectedCells = findAllSelectedLabels(S);
    std::set<int> selectedCells;
    std::map<int,double> initialCells;

    forall(auto p, heat){
      if(p.second < thresholdDis or (selectedAsDirect and alreadySelectedCells.find(p.first) != alreadySelectedCells.end())){
        selectedCells.insert(p.first);
        if(weight == "Euclidean") initialCells[p.first] = p.second;
        else initialCells[p.first] = 0;
      }
    }

    if(selectCells){
      forall(const vertex& v, S) {
        if(!keep) 
          v->selected = false;
        if(selectedCells.find(v->label) == selectedCells.end()) continue;
          v->selected = true;
      }
      mesh->updateAll();
    }

    //if(considerOrientation){
    IntFloatAttr heat2;
    // run Distance to Plane from Bezier
    MeasureDistanceToBezierPlane3D* proc2;
    if(!getProcess("Mesh/Heat Map/Measures 3D/Location/Distance to Plane from Bezier", proc2))
      throw(QString("MeshProcessHeatMap:: Unable to make SelectByHeat process"));
    proc2->run();

    forall(auto p, mesh->labelHeat())
      heat2[p.first] = p.second;

    // create cell graph
    NhbdGraphInfo info;
    neighborhoodGraph(S, 0.0001, info); // now in graphutils

    std::map<IntInt, double> neighMap = info.sharedWallArea;

    // find selected cells
    // std::set<int> selectedCells = findAllSelectedLabels(S);
    std::set<int> allCells = findAllLabelsSet(S);

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

      forall(auto 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(auto 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;
    }

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

    if(weight == "Euclidean"){
      std::map<IntInt, double> neighMapNew;
      forall(auto p, neighMap){
        neighMapNew[p.first] = 1./norm(centroids[p.first.first] - centroids[p.first.second]);
      }
      neighMap = neighMapNew;
    }


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

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

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    // update heat map
    mesh->labelHeat().clear();
    forall(const vertex& v, S){
      int currentLabel = v->label;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents()) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      if(heatDijk.find(currentLabel) == heatDijk.end()) continue;
      mesh->labelHeat()[currentLabel] = heatDijk[currentLabel];
      if(considerOrientation){
        if(heat2[currentLabel] < 0) mesh->labelHeat()[currentLabel] *= -1;
      }
      attrData[currentLabel] = heatDijk[currentLabel];
      if(maxValue < heatDijk[currentLabel]) maxValue = heatDijk[currentLabel];
      if(minValue > heatDijk[currentLabel]) minValue = heatDijk[currentLabel];
    }

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

    return true;
  }
  REGISTER_PROCESS(Measure3DCellDistanceBezierRing);



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

    vvGraph& S = mesh->graph();

    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

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

    if(centroids.empty()) return setErrorMessage("No centroids found. Run Cell Analysis 3D first");

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

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

    forall(int l, allCells){
      if(centroids.find(l) == centroids.end()) continue;
      Point3d centroidRotated = multMatrix4Point3(mGLTot,centroids[l]); // 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(Measure3DCellCoord);


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

    vvGraph& S = mesh->graph();
    vvGraph& S2 = mesh2->graph();

    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

    if(centroids.empty()) return setErrorMessage("No centroids in the active mesh.");

    std::map<int, double> newHeat;

    forall(auto p, centroids){

      double minDis = HUGE_VAL;

      forall(vertex v, S2){
        double dis = norm(v->pos - p.second);
        if(dis < minDis){
          minDis = dis;
        }
      }
      newHeat[p.first] = minDis;

    }

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    // update heat map
    mesh->labelHeat().clear();
    forall(const vertex& v, S){
      int currentLabel = v->label;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents()) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      mesh->labelHeat()[currentLabel] = newHeat[currentLabel];
      attrData[currentLabel] = newHeat[currentLabel];
      if(maxValue < newHeat[currentLabel]) maxValue = newHeat[currentLabel];
      if(minValue > newHeat[currentLabel]) minValue = newHeat[currentLabel];
    }

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

    std::cout << "new size " << attrData.size() << std::endl;

    return true;
  }
  REGISTER_PROCESS(Measure3DDistanceToMesh);


  bool MeasureDistanceToBezier3D::run(const Stack *s1, Mesh* mesh, bool orient, bool reverse, IntFloatAttr& heatMap)
  {
    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

    CuttingSurface* cutSurf = cuttingSurface();

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

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

    Bezier b = cutSurf->bezier();

    int dataPointsBezier = 200;
    std::vector<std::vector<Point3d> > bezGrid;

    b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);

    std::map<int, double> newHeat;

    forall(auto p, centroids){
      Point2i idx;
      Point3d pGrid = calcNearestPointOnBezierGrid(p.second, bezGrid, idx);

      if(idx.x() == dataPointsBezier-1) idx.x()--;
      if(idx.x() == 0) idx.x()++;
      if(idx.y() == dataPointsBezier-1) idx.y()--;
      if(idx.y() == 0) idx.y()++;

      Point3d bezDerivX = bezGrid[idx.x()+1][idx.y()] - bezGrid[idx.x()-1][idx.y()];
      Point3d bezDerivY = bezGrid[idx.x()][idx.y()+1] - bezGrid[idx.x()][idx.y()-1];

      Point3d bezNrml = bezDerivX ^ bezDerivY;

      newHeat[p.first] = norm(pGrid-p.second);

      double orientation = bezNrml * (p.second - pGrid);

      if(orientation < 0 and !reverse and orient) newHeat[p.first] = -newHeat[p.first];
      else if(orientation > 0 and reverse and orient) newHeat[p.first] = -newHeat[p.first];
    }

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    // update heat map
    mesh->labelHeat().clear();
    forall(auto p, centroids){
      int currentLabel = p.first;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents()) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      mesh->labelHeat()[currentLabel] = newHeat[currentLabel];
      attrData[currentLabel] = newHeat[currentLabel];
      if(maxValue < newHeat[currentLabel]) maxValue = newHeat[currentLabel];
      if(minValue > newHeat[currentLabel]) minValue = newHeat[currentLabel];
    }

    heatMap = mesh->labelHeat();

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

    return true;
  }
  REGISTER_PROCESS(MeasureDistanceToBezier3D);


  bool MeasureBezierLineCoord3D::run(const Stack *s1, Mesh* mesh, IntFloatAttr& heatMap)
  {
    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

    CuttingSurface* cutSurf = cuttingSurface();

    Matrix4d rotMatrixS1, rotMatrixCS;
    double totalBezLength;
    int dataPointsBezier = 200;
    s1->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());

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

    Bezier b = cutSurf->bezier();

    std::vector<Point3d> bezierVec, bDiff;

    std::map<int, double> newHeat;

    b.discretizeLineEqualWeight(dataPointsBezier, 0.0, mGLTot, bezierVec, bDiff, totalBezLength);

    double bezPointDis = totalBezLength / (double)dataPointsBezier;

    forall(auto p, centroids){
      int idx;
      double idxWeight;
      //Point3d pGrid = calcNearestPointOnBezierGrid(p.second, bezGrid, idx);
      Point3d closest = calcNearestPointOnBezierLine(p.second, bezierVec, idx, idxWeight);

      double bezDis = ((idx-1) + idxWeight) * bezPointDis;

      newHeat[p.first] = bezDis;
    }

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    // update heat map
    mesh->labelHeat().clear();
    forall(auto p, centroids){
      int currentLabel = p.first;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents()) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      mesh->labelHeat()[currentLabel] = newHeat[currentLabel];
      attrData[currentLabel] = newHeat[currentLabel];
      if(maxValue < newHeat[currentLabel]) maxValue = newHeat[currentLabel];
      if(minValue > newHeat[currentLabel]) minValue = newHeat[currentLabel];
    }

    heatMap = mesh->labelHeat();

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

    return true;
  }
  REGISTER_PROCESS(MeasureBezierLineCoord3D);



  bool MeasureDistanceToBezierPlane3D::run(const Stack *s1, Mesh* mesh, bool orient, bool reverse, IntFloatAttr& heatMap)
  {
    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

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

    CuttingSurface* cutSurf = cuttingSurface();

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

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

    Bezier b = cutSurf->bezier();

    int dataPointsBezier = 200;
    std::vector<std::vector<Point3d> > bezGrid;
    std::vector<Point3d> bezVec;

    // get all supporting points of Bezier
    b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);

    forall(auto p, bezGrid){
      forall(auto p2, p){
        bezVec.push_back(p2);
      }
    }

    // do a pca for plane fitting
    Point3d planePos, planeNrml;
    findPolygonPlane(bezVec, planePos, planeNrml);

    // calc distance plane to cell centroids
    std::map<int, double> newHeat;

    forall(auto p, centroids){
      Point2i idx;

      Point3d pPlane = projectPointOnPlane(p.second, planePos, planeNrml);

      newHeat[p.first] = norm(pPlane-p.second);

      double orientation = planeNrml * (p.second - pPlane);

      if(orientation < 0 and !reverse and orient) newHeat[p.first] = -newHeat[p.first];
      else if(orientation > 0 and reverse and orient) newHeat[p.first] = -newHeat[p.first];
    }

    double maxValue = -HUGE_VAL, minValue = HUGE_VAL;

    // update heat map
    mesh->labelHeat().clear();
    forall(auto p, centroids){
      int currentLabel = p.first;//(*cellAtlasAttr)[i].cellLabel;
      //if(mesh->useParents()) currentLabel = mesh->parents()[v->label];
      if(currentLabel < 1) continue;
      mesh->labelHeat()[currentLabel] = newHeat[currentLabel];
      attrData[currentLabel] = newHeat[currentLabel];
      if(maxValue < newHeat[currentLabel]) maxValue = newHeat[currentLabel];
      if(minValue > newHeat[currentLabel]) minValue = newHeat[currentLabel];
    }

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

    return true;
  }
  REGISTER_PROCESS(MeasureDistanceToBezierPlane3D);


}
