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


#include <MeshProcessMeasures.hpp>
#include <Triangulate.hpp>

//#include <RootCellProcessing.hpp> // 3D cell properties from Cell atlas
//#include <RootCellAnalyzing.hpp> // 3D cell properties from Cell atlas

#include <MeshProcessCellAnalysis3D.hpp>

//#include <CGALUtils.hpp> // pointInPolygon
#include <ThirdParty/triangle.h> // delaunay2D
#include <MeshBuilder.hpp>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>


#include <QFileDialog>

#include <algorithm>
#include <chrono>
#include <thread>

#include "Dir.hpp" // stripCurrentDir


namespace mgx
{

  typedef std::pair<int, set<int> > IntIntSetP;
  typedef std::pair<int, Point3f> IntPoint3fP;
  typedef std::pair<vertex, Point3d> VtxP3dP;


  std::map<int, std::set<vertex> > getCellJunctions(vvGraph& S)
  {
    std::map<int, std::set<vertex> > junctions;

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

      if(neighborLabels.size() > 2 or (neighborLabels.size() == 2 and v->margin)){
        forall(int l, neighborLabels){
          junctions[l].insert(v);
        }
      }
    }


    return junctions;
  }

  std::map<int, std::set<vertex> > getCellContour(vvGraph& S)
  {
    std::map<int, std::set<vertex> > contour;

    forall(const vertex& v, S){
      if(v->label > -1) continue;
      std::set<int> neighborLabels;
      forall(const vertex& n, S.neighbors(v)){
        if(n->label < 0) continue;
        contour[n->label].insert(v);
      }
    }


    return contour;
  }

   void createJunctionAssociation(Mesh* m, Mesh* m2, vvGraph& T1, vvGraph T2, JunctionInformation& jInfo)
   {

   jInfo.J1J2.clear();
   jInfo.J2J1.clear();

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


   // get the cell junctions
   JunctionCellsMap junctionLabels1, junctionLabels2;
   CellJunctionsMap cellJunctions1, cellJunctions2;
   CellJunctionsInOrderMap cellJunctionInOrder1;

   findCellJunctions(m, T1, jInfo);
   findCellJunctionsParents(m2, T2, jInfo);
   findCellJunctionsInOrder(m2, T1, jInfo);

   forall(const VertLabelsP& j1, jInfo.junctionLabels1){
     forall(const VertLabelsP& j2, jInfo.junctionLabels2){
       if(j1.second == j2.second){
       	std::cout << "j test " << j1.second.size() << "/";
       	  forall(int l, j1.second){
            std::cout << l << "/";
       	  }
       	  std::cout << std::endl;

         jInfo.J1J2[j1.first] = j2.first;
         jInfo.J2J1[j2.first] = j1.first;
       }
     }
   }

   // special treatment for the junctions at the border of the known region:
   // those are not complete as one parent label is still missing, but can be found as two cells are known
   std::set<VertLabelsP> borderCellJunctions;

   forall(const VertLabelsP& j2, jInfo.junctionLabels2){
   	 if(j2.second.size() == 2){
   	 	borderCellJunctions.insert(j2);
   	 }
   }
//std::cout << "bla2 " << J1J2.size() << "/" << J2J1.size() << std::endl;
  forall(VertLabelsP p, borderCellJunctions){
  	std::set<VertLabelsP> foundJ;
  	forall(const VertLabelsP& j1, jInfo.junctionLabels1){
  	  std::vector<int> sharedN;
      std::set_intersection (p.second.begin(), p.second.end(), j1.second.begin(), j1.second.end(), std::back_inserter(sharedN));
      if(sharedN.size() > 1){
       	foundJ.insert(j1);
      }
  	}

    forall(const VertLabelsP& j, foundJ){
      if(jInfo.J1J2.find(j.first) == jInfo.J1J2.end()){
        jInfo.J1J2[j.first] = p.first;
        jInfo.J2J1[p.first] = j.first;
      }
    }

  }

   }

   void writeInputVectors(QString inputType, std::vector<Point3d>& data, std::vector<Point3d>& values, JunctionInformation& jInfo, double distanceThreshold, std::set<int>& ignoreLabels)
   {

    data.clear();
    values.clear();

    std::map<vertex, Point3d> replaceVtxsMap1, replaceVtxsMap2;
    std::set<Point3d> addedAvgPoints;


    if(inputType == "Junctions" or inputType == "Both"){


      forall(CellJunctionsInOrderP p, jInfo.cellJunctionInOrder1){

        std::vector<vertex> borderInOrder = p.second;

        for(uint i = 0; i < borderInOrder.size()-1; i++){
      	  vertex v1,v2;
      	  v1 = borderInOrder[i];
      	  v2 = borderInOrder[borderInOrder.size() - 1];
      	  if(i > 0) v2 = borderInOrder[i + 1];

          vertex v21, v22;
          if(jInfo.J1J2.find(v1) == jInfo.J1J2.end() or jInfo.J1J2.find(v2) == jInfo.J1J2.end()) continue;

          v21 = jInfo.J1J2[v1];
          v22 = jInfo.J1J2[v2];

          if(norm(v1->pos - v2->pos) < distanceThreshold or norm(v21->pos - v22->pos) < distanceThreshold){

            // only 1 merging allowed currently
            if(replaceVtxsMap1.find(v1) != replaceVtxsMap1.end() or
               replaceVtxsMap1.find(v2) != replaceVtxsMap1.end() or
               replaceVtxsMap2.find(v21) != replaceVtxsMap2.end() or
               replaceVtxsMap2.find(v22) != replaceVtxsMap2.end())
            	continue;

            Point3d avgPos = (v1->pos + v2->pos)/2.;
            replaceVtxsMap1[v1] = avgPos;
            replaceVtxsMap1[v2] = avgPos;

            Point3d avgPos2 = (v21->pos + v22->pos)/2.;
            replaceVtxsMap2[v21] = avgPos2;
            replaceVtxsMap2[v22] = avgPos2;

            // debug output
            // v1->selected = true;
            // v2->selected = true;

            // v21->selected = true;
            // v22->selected = true;
          }

        }

      }

    std::cout << "merged vertices: " << replaceVtxsMap1.size() << "/" << replaceVtxsMap2.size() << std::endl;
    std::cout << "size of correspondance map: " << jInfo.J1J2.size() << std::endl;

      forall(VtxVtxP p, jInfo.J1J2){
      	if(replaceVtxsMap1.find(p.first) != replaceVtxsMap1.end() and replaceVtxsMap2.find(p.second) != replaceVtxsMap2.end()){
          if(addedAvgPoints.find(replaceVtxsMap1[p.first]) != addedAvgPoints.end()) continue;

          data.push_back(replaceVtxsMap1[p.first]);
          values.push_back(replaceVtxsMap2[p.second]);
          addedAvgPoints.insert(replaceVtxsMap1[p.first]);

          std::cout << "data/values (J/R) " << p.first->pos << "/" << p.second->pos << std::endl;
          p.first->selected = true;
          p.second->selected = true;
      	} else {
          data.push_back(p.first->pos);
          values.push_back(p.second->pos);

          std::cout << "data/values (J) " << p.first->pos << "/" << p.second->pos << std::endl;
      	}
      }
    }
    if(inputType == "Centers" or inputType == "Both"){
      forall(IntPoint3dP p, jInfo.labelDaughtersCenterMap){

        if(ignoreLabels.find(p.first) != ignoreLabels.end()){
          std::cout << "ignored label " << p.first << std::endl;
          continue;
        }

        if(jInfo.labelParentCenterMap.find(p.first) == jInfo.labelParentCenterMap.end()) continue;
      //forall(const vertex& v, S){
      //  if(labelPosMap.find(v->label) != labelPosMap.end()){

      	//TODO for some reason there are both points with 0-values in there?
        if(jInfo.labelParentCenterMap[p.first] == Point3d(0,0,0) and p.second == Point3d(0,0,0)) continue;

            data.push_back(jInfo.labelParentCenterMap[p.first]);
            values.push_back(p.second);

            std::cout << "data/values (C) " << p.first << "/" << jInfo.labelParentCenterMap[p.first] << "/" << p.second << std::endl;
         // }
        }
    }

   }


  bool DeformMesh::run(Mesh* m, Mesh* m2)
  {

   vvGraph& S = m->graph();
   vvGraph& S2 = m2->graph();

   AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");


   origMeshPos.clear();
   forall(const vertex& v, S){
     origMeshPos[v] = v->pos;
   }

   vvGraph T1, T2;
   mgxmToMgxc(S, T1, 0);
   mgxmToMgxc(S2, T2, 0);

   CreateImplicitDeformation2D *createImpDef;
    if(!getProcess("Mesh/Deformation/Mesh 2D/Create Deformation Cell Surface", createImpDef))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateImplicitDeformation2D process"));

   AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
   bool forceRecalc = stringToBool(createImpDef->parm("Force Recalculation"));

   if(impDefAttr.empty() or forceRecalc)
     createImpDef->run(m, m2, T1, T2, createImpDef->parm("Direction"), createImpDef->parm("CG Steps").toInt(), createImpDef->parm("Junction Merging Threshold (um)"),  createImpDef->parm("Junction Merging Threshold (um)").toDouble());

   ImplicitDeformation& impDef = impDefAttr[0];//createImpDef->impDef;

  if(createImpDef->parm("Direction") == "active -> other"){
    forall(const vertex& v, S){
      Point3d deformedPos = impDef.RBFdeformation(v->pos);
      v->pos = deformedPos;
    }
    m->updateAll();
  } else {
    forall(const vertex& v, S2){
      Point3d deformedPos = impDef.RBFdeformation(v->pos);
      v->pos = deformedPos;
    }
    m2->updateAll();
  }

    return true;
  }
  REGISTER_PROCESS(DeformMesh);

   bool DeformationProjection::run(Mesh* m, Mesh* m2, QString projVtxs)
  {

   vvGraph& S = m->graph();
   vvGraph& S2 = m2->graph();

   vvGraph T1, T2;
   mgxmToMgxc(S, T1, 0);
   mgxmToMgxc(S2, T2, 0);

   // create deformation
   CreateImplicitDeformation2D *createImpDef;
    if(!getProcess("Mesh/Deformation/Mesh 2D/Create Deformation Cell Surface", createImpDef))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateImplicitDeformation2D process"));

   AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
   bool forceRecalc = stringToBool(createImpDef->parm("Force Recalculation"));

   if(impDefAttr.empty() or forceRecalc)
     createImpDef->run(m, m2, T1, T2, createImpDef->parm("Direction"), createImpDef->parm("CG Steps").toInt(), createImpDef->parm("Junction Merging Threshold (um)"),  createImpDef->parm("Junction Merging Threshold (um)").toDouble());

   ImplicitDeformation& impDef = impDefAttr[0];//createImpDef->impDef;
   JunctionInformation& jInfo = createImpDef->jInfo;

   IntPoint3fAttr centers = m->labelCenter();

   std::set<int> cellsM1, cellsM2, newCells;


  forall(const vertex& v, S){
    if(v->label < 1) continue;
    cellsM1.insert(v->label);
    //if(labelPosMap.find(v->label) == labelPosMap.end()){
    //if(newCells.find(v->label) != newCells.end()){
      //v->selected = true;
    //}
  }

  std::map<int,Point3d> centers2;// = m2->labelCenter();
  forall(const vertex& v, S2){
    if(v->label < 1) continue;
    cellsM2.insert(v->label);
    centers2[v->label] = v->pos;
  }

  std::vector<vertex> newVtxs;
  std::vector<Point3d> pointsToProject;

  if(projVtxs == "All"){
  	if(createImpDef->parm("Direction") == "active -> other"){
  	  forall(const vertex& v, S){
  	    pointsToProject.push_back(v->pos);
  	  }
  	} else if(createImpDef->parm("Direction") == "other -> active"){
  	  forall(const vertex& v, S2){
  	    pointsToProject.push_back(v->pos);
  	  }
    }

  } else {
    if(projVtxs == "Centers" or projVtxs == "Centers + Junctions"){

  	  if(createImpDef->parm("Direction") == "active -> other"){
  	    newCells = cellsM1;
        forall(int l, newCells){
      	  pointsToProject.push_back(Point3d(centers[l]));
        }
  	  } else if(createImpDef->parm("Direction") == "other -> active"){
  	    newCells = cellsM2;
        forall(int l, newCells){
      	  pointsToProject.push_back(Point3d(centers2[l]));
        }

  	  }


    }

    if(projVtxs == "Junctions" or projVtxs == "Centers + Junctions"){

      // get the cell junctions
      JunctionCellsMap junctionLabels;

  	  if(createImpDef->parm("Direction") == "active -> other"){
  	    junctionLabels = jInfo.junctionLabels1;
  	  } else if(createImpDef->parm("Direction") == "other -> active"){
  	    junctionLabels = jInfo.junctionLabels2;
  	  }

      forall(VertLabelsP p, junctionLabels){
      	pointsToProject.push_back(Point3d(p.first->pos));
      }

    }
  }

  forall(Point3d p, pointsToProject){
  	Point3d deformedPos = impDef.RBFdeformation(p);
  	//if(projectOnMesh) deformedPos = nearestPointOnMesh(deformedPos, m2->graph());
  	vertex vNew;
    vNew->pos = deformedPos;
    vNew->selected = true;
    newVtxs.push_back(vNew);
  }

    forall(vertex v, newVtxs){
      if(createImpDef->parm("Direction") == "active -> other") S2.insert(v);
      else S.insert(v);
    }

    if(createImpDef->parm("Direction") == "active -> other") m2->updateAll();
    else m->updateAll();

    return true;
  }
  REGISTER_PROCESS(DeformationProjection);

  bool MorphMesh::run(Mesh* m, Mesh* m2, int steps, int msStep)
  {

    using namespace std::this_thread; // sleep_for
    using namespace std::chrono; // nanoseconds

    vvGraph& S = m->graph();
    vvGraph& S2 = m2->graph();

    vvGraph T1, T2;
    mgxmToMgxc(S, T1, 0);
    mgxmToMgxc(S2, T2, 0);

    // save the original vertex positions of the mesh that is to be changed
    std::map<vertex, Point3d> origPos;

    vvGraph meshToMorph;

    // create deformation
   CreateImplicitDeformation2D *createImpDef;
   if(!getProcess("Mesh/Deformation/Mesh 2D/Create Deformation Cell Surface", createImpDef))
     throw(QString("MeshProcessHeatMap:: Unable to make CreateImplicitDeformation2D process"));



   AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
   bool forceRecalc = stringToBool(createImpDef->parm("Force Recalculation"));

   if(impDefAttr.empty() or forceRecalc)
     createImpDef->run(m, m2, T1, T2, createImpDef->parm("Direction"), createImpDef->parm("CG Steps").toInt(), createImpDef->parm("Junction Merging Threshold (um)"),  createImpDef->parm("Junction Merging Threshold (um)").toDouble());

   ImplicitDeformation& impDef = impDefAttr[0];//createImpDef->impDef;

    if(createImpDef->parm("Force Recalculation") == "m1 -> m2"){
      meshToMorph = S;
    } else if(createImpDef->parm("Force Recalculation") == "m2 -> m1"){
  	  meshToMorph = S2;
    } else {
  	  return setErrorMessage("Direction not defined properly");
    }

    forall(const vertex& v, meshToMorph){
      origPos[v] = v->pos;
    }

    progressStart(QString("Morphing mesh-%1").arg(m->userId()), steps);

    // morph mesh
    for(int i = 0; i < steps; i++){
      if(!progressAdvance(i))
        userCancel();

      std::cout << "step " << i+1 << " running ....";
      auto startTime = std::chrono::system_clock::now();

      forall(const vertex& v, meshToMorph){
        Point3d deformedPos = impDef.RBFdeformation(origPos[v]);

        Point3d newPos = ((i+1) * deformedPos + (steps-i-1) * origPos[v])/((double)steps);
        //std::cout << "new pos" << v->pos << "/" << origPos[v] << "/" << newPos << std::endl;
        v->pos = newPos;
      }
      if(createImpDef->parm("Force Recalculation") == "m1 -> m2") m->updateAll();
      else m2->updateAll();

      updateState();
      updateViewer();

      std::cout << " done" << std::endl;

      sleep_until(startTime + nanoseconds(msStep*1000));
    }

    return true;
  }
  REGISTER_PROCESS(MorphMesh);

  void morphMeshStep(vvGraph& meshToMorph, std::vector<AttrMap<vertex, Point3d> >& deformedPosVector, int deformationStep, int i, int steps)
  {
    forall(const vertex& v, meshToMorph){
      Point3d currentOrigPos = deformedPosVector[deformationStep-1][v];
      Point3d deformedPos = deformedPosVector[deformationStep][v];
      Point3d newPos = ((i+1) * deformedPos + (steps-i-1) * currentOrigPos)/((double)steps);
    //   //std::cout << "new pos" << v->pos << "/" << origPos[v] << "//" << deformedPos << "/" << newPos << std::endl;
      v->pos = newPos;
    }
  }

 void morphMeshStepFix(vvGraph& meshToMorph, std::vector<AttrMap<vertex, Point3d> >& deformedPosVector, int deformationStep, int i, int steps,
  vertex fixV, Point3d fixDispl, Point3d& newDispl)
  {
    std::map<vertex, Point3d> vtxDisplacementMap;

    forall(const vertex& v, meshToMorph){
      Point3d origPosD = deformedPosVector[deformationStep-1][v];

      Point3d deformedPos = deformedPosVector[deformationStep][v];
      Point3d newPos = ((i+1) * deformedPos + (steps-i-1) * origPosD)/((double)steps);
      //std::cout << "new pos" << v->pos << "/" << origPos[v] << "//" << deformedPos << "/" << newPos << std::endl;
      vtxDisplacementMap[v] = newPos - deformedPosVector[0][v];
      //v->pos = newPos;
    }


    //fixDispl = fixDispl + vtxDisplacementMap[fixV];
    newDispl = vtxDisplacementMap[fixV];
    //std::cout << "fix " << fixDispl << "/" << newDispl << std::endl;

    forall(const vertex& v, meshToMorph){
      v->pos = deformedPosVector[0][v] + vtxDisplacementMap[v] - newDispl;
      //v->pos += vtxDisplacementMap[v] - fixDispl - newDispl;
    }

  }

 void calcDeformedPoints(vvGraph& meshToMorph, std::vector<ImplicitDeformation>& impDefVec, std::vector<AttrMap<vertex, Point3d> >& deformedPosVector)
 {
  deformedPosVector.clear();

  // save the initial points
  AttrMap<vertex, Point3d> vtxPMap;
  forall(const vertex& v, meshToMorph){
    vtxPMap[v] = v->pos;
  }
  deformedPosVector.push_back(vtxPMap);

  for(uint i=0; i < impDefVec.size(); i++){
    AttrMap<vertex, Point3d> vtxPMap;
    ImplicitDeformation& impDef = impDefVec[i];


    //forall(const vertex& v, meshToMorph){
    #pragma omp parallel for
    for(int j = 0; j < meshToMorph.size(); j++){
      vertex v = meshToMorph[j];
      Point3d origPos = deformedPosVector[i][v];
      Point3d deformedPos = impDef.RBFdeformation(origPos);
      vtxPMap[v] = deformedPos;
    }
    deformedPosVector.push_back(vtxPMap);
  }
 }

 // https://en.wikipedia.org/wiki/Catmull-Rom_spline
 Point3d calcCatmullRomSpline(Point3d p0, Point3d p1, double t, Point3d pPrev, Point3d pNext)
 {

  double h00 = 2*t*t*t - 3*t*t + 1;
  double h10 = t*t*t - 2*t*t + t;
  double h01 = -2*t*t*t + 3*t*t;
  double h11 = t*t*t - t*t;

  Point3d m0 = (p1 - pPrev)/2.;
  Point3d m1 = (pNext - p0)/2.;

  Point3d result = h00*p0 + h10*m0 + h01*p1 + h11*m1;

  return result;
 }


 void morphMeshStepFixSpline(vvGraph& meshToMorph, std::vector<AttrMap<vertex, Point3d> >& deformedPosVector, int deformationStep, int i, int steps,
  vertex fixV, vertex fixAxis, std::map<vertex, Point3d>& displacementMap)
  {

    Point3d newDispl;
    double t = (double)(i+1) / (double)(steps);

    std::cout << "spline " << i << "/" << steps << "/" << t << "/" << deformationStep << std::endl;
    std::cout << "spline " << deformedPosVector[deformationStep].size() << std::endl;

    std::map<vertex, Point3d> vtxDisplacementMap;

    Point3d zeroP(0,0,0);
    forall(const vertex& v, meshToMorph){
      vtxDisplacementMap[v] = zeroP;
    }

    #pragma omp parallel for
    for(int i = 0; i < meshToMorph.size(); i++){
      const vertex& v = meshToMorph[i];
      Point3d origPosD = deformedPosVector[deformationStep-1][v];

      Point3d deformedPos = deformedPosVector[deformationStep][v];

      Point3d newPosOld = ((i+1) * deformedPos + (steps-i-1) * origPosD)/((double)steps);
      Point3d newPos;
      if(deformationStep == 1 and deformedPosVector.size() == 2){
        newPos = calcCatmullRomSpline(deformedPosVector[deformationStep-1][v], deformedPosVector[deformationStep][v], t, deformedPosVector[deformationStep-1][v], deformedPosVector[deformationStep][v]);
      } else if(deformationStep == deformedPosVector.size()-1){
        newPos = calcCatmullRomSpline(deformedPosVector[deformationStep-1][v], deformedPosVector[deformationStep][v], t, deformedPosVector[deformationStep-2][v], deformedPosVector[deformationStep][v]);
      } else if(deformationStep == 1) {
        newPos = calcCatmullRomSpline(deformedPosVector[deformationStep-1][v], deformedPosVector[deformationStep][v], t, deformedPosVector[deformationStep-1][v], deformedPosVector[deformationStep+1][v]);
      } else {
        newPos = calcCatmullRomSpline(deformedPosVector[deformationStep-1][v], deformedPosVector[deformationStep][v], t, deformedPosVector[deformationStep-2][v], deformedPosVector[deformationStep+1][v]);
      }

      //std::cout << "new pos " << newPos << "/" << newPosOld << "//" << deformedPos << "//" << origPosD << std::endl;
      vtxDisplacementMap[v] = newPos - deformedPosVector[0][v];
      //v->pos = newPos;
    }


    forall(const vertex& v, meshToMorph){
      // translate
      v->pos = deformedPosVector[0][v] + vtxDisplacementMap[v];// - newDispl;
      //v->pos += vtxDisplacementMap[v] - fixDispl - newDispl;
    }
    displacementMap = vtxDisplacementMap;

  }


   bool MorphMeshCreatePositions::run(Mesh* m, QString def1, QString def2, QString def3, QString def4, QString def5, QString def6)
  {

    vvGraph& S = m->graph();

    AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");
    origMeshPos.clear();

    AttrMap<int, ImplicitDeformation>& impDefAttr1 = m->attributes().attrMap<int, ImplicitDeformation>(def1);
    if(impDefAttr1.empty()) std::cout << "Implicit Deformation 1 doesn't exist!" << std::endl;

    AttrMap<int, ImplicitDeformation>& impDefAttr2 = m->attributes().attrMap<int, ImplicitDeformation>(def2);
    if(impDefAttr2.empty()) std::cout << "Implicit Deformation 2 doesn't exist!" << std::endl;

    AttrMap<int, ImplicitDeformation>& impDefAttr3 = m->attributes().attrMap<int, ImplicitDeformation>(def3);
    if(impDefAttr3.empty()) std::cout << "Implicit Deformation 3 doesn't exist!" << std::endl;

    AttrMap<int, ImplicitDeformation>& impDefAttr4 = m->attributes().attrMap<int, ImplicitDeformation>(def4);
    if(impDefAttr4.empty()) std::cout << "Implicit Deformation 4 doesn't exist!" << std::endl;

    AttrMap<int, ImplicitDeformation>& impDefAttr5 = m->attributes().attrMap<int, ImplicitDeformation>(def5);
    if(impDefAttr5.empty()) std::cout << "Implicit Deformation 5 doesn't exist!" << std::endl;

    AttrMap<int, ImplicitDeformation>& impDefAttr6 = m->attributes().attrMap<int, ImplicitDeformation>(def6);
    if(impDefAttr6.empty()) std::cout << "Implicit Deformation 6 doesn't exist!" << std::endl;

    std::map<vertex, Point3d> origPos;
    std::vector<ImplicitDeformation> morphingDeformations;
    std::vector<QString> impDefAttr;

    forall(const vertex& v, S){
      origMeshPos[v] = v->pos;
    }

    if(!impDefAttr1[0].data.empty()){
      morphingDeformations.push_back(impDefAttr1[0]);
      impDefAttr.push_back(def1);
    }
    if(!impDefAttr2[0].data.empty()){
      morphingDeformations.push_back(impDefAttr2[0]);
      impDefAttr.push_back(def2);
    }
    if(!impDefAttr3[0].data.empty()){
      morphingDeformations.push_back(impDefAttr3[0]);
      impDefAttr.push_back(def3);
    }
    if(!impDefAttr4[0].data.empty()){
      morphingDeformations.push_back(impDefAttr4[0]);
      impDefAttr.push_back(def4);
    }
    if(!impDefAttr5[0].data.empty()){
      morphingDeformations.push_back(impDefAttr5[0]);
      impDefAttr.push_back(def5);
    }
    if(!impDefAttr6[0].data.empty()){
      morphingDeformations.push_back(impDefAttr6[0]);
      impDefAttr.push_back(def6);
    }


    std::cout << "Calculate Deformed Positions" << std::endl;
    std::vector<AttrMap<vertex, Point3d> > deformedPosVector;
    calcDeformedPoints(S, morphingDeformations, deformedPosVector);

    for(uint i = 1; i < deformedPosVector.size(); i++){
      AttrMap<vertex, Point3d>& impDefPos = m->attributes().attrMap<vertex, Point3d>(impDefAttr[i-1] + "_positions");
      impDefPos.clear();
      forall(VtxP3dP p, deformedPosVector[i]){
        impDefPos[p.first] = p.second;
      }
    }


    return true;

  }
    REGISTER_PROCESS(MorphMeshCreatePositions);

  double interpolatedSignal(Point3d pos, const std::vector<vertex>& vtxVec)
  {
    double result = 0;

    // find nearest vertex
    double minDis = HUGE_VAL;
    vertex minV;
    forall(vertex v, vtxVec){
      double dis = norm(pos - v->pos);
      if(dis < minDis){
        minDis = dis;
        minV = v;
      }
    }

    return minV->signal;
  }


   bool CreateSignalProjection::run(Mesh* m, Mesh* m2, int useDef, bool useParentLabels)
  {

    if(useDef < 1 or useDef > 6) return setErrorMessage("Use Def Parameter must be between 1 and 6");

    //QStringList createMorphingParms;
    //if(!getProcessParms<MorphMeshCreatePositions>(this, createMorphingParms)) throw(QString("Unable to get Create Morphing Parameters"));

    MorphMeshCreatePositions *createMorphingParms;
    if(!getProcess("Mesh/Deformation/Morphing/Create Morphing", createMorphingParms))
      throw(QString("MeshProcessHeatMap:: Unable to make createMorphingParms process"));

    QString def1 = createMorphingParms->parm("Impl Def Attr Map 1");
    QString def2 = createMorphingParms->parm("Impl Def Attr Map 2");
    QString def3 = createMorphingParms->parm("Impl Def Attr Map 3");
    QString def4 = createMorphingParms->parm("Impl Def Attr Map 4");
    QString def5 = createMorphingParms->parm("Impl Def Attr Map 5");
    QString def6 = createMorphingParms->parm("Impl Def Attr Map 6");

    QString attrName = "";

    if(useDef == 6){
      attrName = def6;
    } else if(useDef == 5){
      attrName = def5;
    } else if(useDef == 4){
      attrName = def4;
    } else if(useDef == 3){
      attrName = def3;
    } else if(useDef == 2){
      attrName = def2;
    } else if(useDef == 1){
      attrName = def1;
    }

    AttrMap<vertex, double>& signalMap = m->attributes().attrMap<vertex, double>(attrName + "_signal");
    signalMap.clear();

    AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");

    AttrMap<vertex, Point3d>& impDefPos = m->attributes().attrMap<vertex, Point3d>(attrName + "_positions");
    if(impDefPos.empty()){
      std::cout << "Attr Map doesn't exist!" << std::endl;
      qDebug() << attrName + "_positions\n";
    } else {
      qDebug() << "Loaded Attr Map" << attrName + "\n";
    }


    // project the mesh points
    vvGraph &S = m->graph();
    vvGraph &S2 = m2->graph();

    forall(const vertex& v, S){
      signalMap[v] = 0;
    }

    std::map<int, std::vector<vertex> > labelVtxVec;
    std::vector<vertex> allVtxs;

    forall(const vertex& v, S2){
      if(v->label > 0){
        labelVtxVec[v->label].push_back(v);
      } else {
        forall(const vertex& n, S2.neighbors(v)){
          if(v->label > 0){
            labelVtxVec[v->label].push_back(n);
          }
        }
      }
      allVtxs.push_back(v);
    }

    forall(auto p, labelVtxVec){
      std::cout << "labelVtxVec " << p.second.size() << std::endl;
    }

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

    std::map<int, std::set<int> > labelNeighborMap;

    forall(auto p, neighMapS1){
      if(p.second < 0.001) continue;
      labelNeighborMap[p.first.first].insert(p.first.second);
      labelNeighborMap[p.first.second].insert(p.first.first);
    }

    // extend the labelVtxMap by the neighbor vtxs
    std::map<int, std::vector<vertex> > labelVtxVecExtended;

    forall(auto p, labelVtxVec){
      forall(vertex v, p.second){
        labelVtxVecExtended[p.first].push_back(v);
      }
      forall(int l, labelNeighborMap[p.first]){
        forall(vertex v, labelVtxVec[l]){
          if(v->label == 0) continue;
          labelVtxVecExtended[p.first].push_back(v);
        }
      }
      std::cout << "label " << p.first << "/" << p.second.size() << "/" <<  labelVtxVecExtended[p.first].size() << std::endl;
    }

    if(useParentLabels){
      #pragma omp parallel for
      for(int i = 0; i < S.size(); i++){
        vertex v = S[i];

        int l = m->parents()[v->label];

        if(l < 1){
          forall(vertex n, S.neighbors(v)){
            if(m->parents()[n->label] > 0) l = m->parents()[n->label];
          }
        }

        if(l < 1){
          signalMap[v] = v->signal;
          continue;
        }

        // use correct attrMap
        Point3d projectedP = impDefPos[v];

        signalMap[v] = interpolatedSignal(projectedP, labelVtxVecExtended[l]);
      }
    } else{
      #pragma omp parallel for
      for(int i = 0; i < S.size(); i++){
        vertex v = S[i];
        int l = v->label;
        if(l == 0 or v->margin){
          signalMap[v] = 0;//v->signal;
          continue;
        }
        Point3d projectedP = impDefPos[v];
        signalMap[v] = interpolatedSignal(projectedP, allVtxs);
      }
    }

    return true;

  }
    REGISTER_PROCESS(CreateSignalProjection);


   bool CreateHeatProjection::run(Mesh* m, int useDef)
  {

    if(useDef < 0 or useDef > 6) return setErrorMessage("Use Def Parameter must be between 0 and 6");

    //QStringList createMorphingParms;
    //if(!getProcessParms<MorphMeshCreatePositions>(this, createMorphingParms)) throw(QString("Unable to get Create Morphing Parameters"));

    MorphMeshCreatePositions *createMorphingParms;
    if(!getProcess("Mesh/Deformation/Morphing/Create Morphing", createMorphingParms))
      throw(QString("MeshProcessHeatMap:: Unable to make createMorphingParms process"));

    QString def1 = createMorphingParms->parm("Impl Def Attr Map 1");
    QString def2 = createMorphingParms->parm("Impl Def Attr Map 2");
    QString def3 = createMorphingParms->parm("Impl Def Attr Map 3");
    QString def4 = createMorphingParms->parm("Impl Def Attr Map 4");
    QString def5 = createMorphingParms->parm("Impl Def Attr Map 5");
    QString def6 = createMorphingParms->parm("Impl Def Attr Map 6");

    QString attrName = "";

    if(useDef == 6){
      attrName = def6;
    } else if(useDef == 5){
      attrName = def5;
    } else if(useDef == 4){
      attrName = def4;
    } else if(useDef == 3){
      attrName = def3;
    } else if(useDef == 2){
      attrName = def2;
    } else if(useDef == 1){
      attrName = def1;
    }
    QString attrNameHeat  = attrName + "_heat";

    if(useDef == 0){
      attrNameHeat = "OrigHeat";
    }

    AttrMap<int, double>& defHeatMap = m->attributes().attrMap<int, double>(attrNameHeat);
    defHeatMap.clear();

    // get the current heat map and save it to the attribute map
    forall(auto p, m->labelHeat()){
      defHeatMap[p.first] = p.second;
    }

    return true;

  }
    REGISTER_PROCESS(CreateHeatProjection);


   bool LoadSignalProjection::run(Mesh* m, QString attrName)
  {

    AttrMap<vertex, double>& attrSignal = m->attributes().attrMap<vertex, double>(attrName);
    if(attrSignal.empty()) std::cout << "Attr Map doesn't exist!" << std::endl;

    forall(const vertex &v, m->graph()){
      v->signal = attrSignal[v];
    }

    m->updateAll();


    return true;

  }
    REGISTER_PROCESS(LoadSignalProjection);


   void createDeformedPosVec(Mesh* m, std::vector<AttrMap<vertex, Point3d> >& deformedPosVector, QString def1, QString def2, QString def3, QString def4, QString def5){

    deformedPosVector.clear();

    AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");

    AttrMap<vertex, Point3d>& impDefPos1 = m->attributes().attrMap<vertex, Point3d>(def1 + "_positions");
    if(impDefPos1.empty()) std::cout << "Implicit Deformation 1 doesn't exist!" << std::endl;

    AttrMap<vertex, Point3d>& impDefPos2 = m->attributes().attrMap<vertex, Point3d>(def2 + "_positions");
    if(impDefPos2.empty()) std::cout << "Implicit Deformation 2 doesn't exist!" << std::endl;

    AttrMap<vertex, Point3d>& impDefPos3 = m->attributes().attrMap<vertex, Point3d>(def3 + "_positions");
    if(impDefPos3.empty()) std::cout << "Implicit Deformation 3 doesn't exist!" << std::endl;

    AttrMap<vertex, Point3d>& impDefPos4 = m->attributes().attrMap<vertex, Point3d>(def4 + "_positions");
    if(impDefPos4.empty()) std::cout << "Implicit Deformation 4 doesn't exist!" << std::endl;

    AttrMap<vertex, Point3d>& impDefPos5 = m->attributes().attrMap<vertex, Point3d>(def5 + "_positions");
    if(impDefPos5.empty()) std::cout << "Implicit Deformation 5 doesn't exist!" << std::endl;

    deformedPosVector.push_back(origMeshPos);

    if(!impDefPos1.empty()){
      deformedPosVector.push_back(impDefPos1);
    }
    if(!impDefPos2.empty()){
      deformedPosVector.push_back(impDefPos2);
    }
    if(!impDefPos3.empty()){
      deformedPosVector.push_back(impDefPos3);
    }
    if(!impDefPos4.empty()){
      deformedPosVector.push_back(impDefPos4);
    }
    if(!impDefPos5.empty()){
      deformedPosVector.push_back(impDefPos5);
    }

   }


   void createDeformedPosVec(Mesh* m, std::vector<AttrMap<vertex, Point3d> >& deformedPosVector, QStringList defList){

    deformedPosVector.clear();

    AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");

    deformedPosVector.push_back(origMeshPos);

    forall(QString def, defList){
      qDebug() << "createDefPosVec " << def << "\n";
      AttrMap<vertex, Point3d>& impDefPos = m->attributes().attrMap<vertex, Point3d>(def + "_positions");
      if(impDefPos.empty()) std::cout << "Implicit Deformation doesn't exist!" << std::endl;

      if(!impDefPos.empty()){
        deformedPosVector.push_back(impDefPos);
      }

    }

   }

 /*  void createDeformedSignalVec(Mesh* m, std::vector<AttrMap<vertex, double> >& deformedSignalVector, QString def1, QString def2, QString def3, QString def4, QString def5){

    deformedSignalVector.clear();

    AttrMap<vertex, double>& origMeshPos = m->attributes().attrMap<vertex, double>("OrigSignal");
    origMeshPos.clear();

    forall(const vertex& v, m->graph()){
      origMeshPos[v] = v->signal;
    }

    AttrMap<vertex, double>& impDefPos1 = m->attributes().attrMap<vertex, double>(def1 + "_signal");
    if(impDefPos1.empty()) std::cout << "Implicit Deformation 1 doesn't exist!" << std::endl;

    AttrMap<vertex, double>& impDefPos2 = m->attributes().attrMap<vertex, double>(def2 + "_signal");
    if(impDefPos2.empty()) std::cout << "Implicit Deformation 2 doesn't exist!" << std::endl;

    AttrMap<vertex, double>& impDefPos3 = m->attributes().attrMap<vertex, double>(def3 + "_signal");
    if(impDefPos3.empty()) std::cout << "Implicit Deformation 3 doesn't exist!" << std::endl;

    AttrMap<vertex, double>& impDefPos4 = m->attributes().attrMap<vertex, double>(def4 + "_signal");
    if(impDefPos4.empty()) std::cout << "Implicit Deformation 4 doesn't exist!" << std::endl;

    AttrMap<vertex, double>& impDefPos5 = m->attributes().attrMap<vertex, double>(def5 + "_signal");
    if(impDefPos5.empty()) std::cout << "Implicit Deformation 5 doesn't exist!" << std::endl;

    deformedSignalVector.push_back(origMeshPos);

    if(!impDefPos1.empty()){
      deformedSignalVector.push_back(impDefPos1);
    }
    if(!impDefPos2.empty()){
      deformedSignalVector.push_back(impDefPos2);
    }
    if(!impDefPos3.empty()){
      deformedSignalVector.push_back(impDefPos3);
    }
    if(!impDefPos4.empty()){
      deformedSignalVector.push_back(impDefPos4);
    }
    if(!impDefPos5.empty()){
      deformedSignalVector.push_back(impDefPos5);
    }

   }
*/

   void createDeformedHeatVec(Mesh* m, std::vector<AttrMap<int, double> >& deformedHeatVector, QStringList defList){

    deformedHeatVector.clear();

    forall(QString def, defList){
      AttrMap<int, double>& impDefPos = m->attributes().attrMap<int, double>(def + "_heat");
      if(impDefPos.empty()) qDebug() << "Implicit Deformation doesn't exist: " << def << "_heat" << "\n";
      if(!impDefPos.empty()){
        deformedHeatVector.push_back(impDefPos);
      }
    }


   }

   void createDeformedSignalVec(Mesh* m, std::vector<AttrMap<vertex, double> >& deformedSignalVector, QStringList defList){

    deformedSignalVector.clear();

    AttrMap<vertex, double>& origMeshPos = m->attributes().attrMap<vertex, double>("OrigSignal");
    origMeshPos.clear();

    forall(const vertex& v, m->graph()){
      origMeshPos[v] = v->signal;
    }
    deformedSignalVector.push_back(origMeshPos);


    forall(QString def, defList){
      AttrMap<vertex, double>& impDefPos = m->attributes().attrMap<vertex, double>(def + "_signal");
      if(impDefPos.empty()) std::cout << "Implicit Deformation doesn't exist!" << std::endl;
      if(!impDefPos.empty()){
        deformedSignalVector.push_back(impDefPos);
      }
    }


   }

   void readFixVertices(Mesh *m, vertex& vOrigin, vertex& vAxis, vertex& vAxis2){

    AttrMap<vertex,Point3d>& morphAxis1 = m->attributes().attrMap<vertex, Point3d>("Morphing Axis1");
    AttrMap<vertex,Point3d>& morphAxis2 = m->attributes().attrMap<vertex, Point3d>("Morphing Axis2");
    AttrMap<vertex,Point3d>& morphOrig = m->attributes().attrMap<vertex, Point3d>("Morphing Origin");

    if(!morphOrig.empty()){
      auto vv =  morphOrig.begin();
      vOrigin = vv->first;
    } else {
      std::cout << "No Origin set!" << std::endl;
      vOrigin = NULL;
    }

    if(!morphAxis1.empty()){
      auto vv =  morphAxis1.begin();
      vAxis = vv->first;
    } else {
      std::cout << "No Axis 1 set!" << std::endl;
      vAxis = NULL;
    }


    if(!morphAxis2.empty()){
      auto vv =  morphAxis2.begin();
      vAxis2 = vv->first;
    } else {
      std::cout << "No Axis 2 set!" << std::endl;
      vAxis2 = NULL;
    }

   }

   bool MorphMeshMulti::run(Stack *s, Stack* s2, Mesh* m, Mesh* m2, int steps, int msStep, QString heatmap, bool useSpline, double heatMin,
    double heatMax, bool morphOther, bool signalBlending, double rotX, double rotY, double rotZ)
  {

    vvGraph& S = m->graph();
    vvGraph& S2 = m2->graph();

    using namespace std::this_thread; // sleep_for

    MorphMeshCreatePositions *createMorphingParms;
    if(!getProcess("Mesh/Deformation/Morphing/Create Morphing", createMorphingParms))
      throw(QString("MeshProcessHeatMap:: Unable to make createMorphingParms process"));

    QString def1 = createMorphingParms->parm("Impl Def Attr Map 1");
    QString def2 = createMorphingParms->parm("Impl Def Attr Map 2");
    QString def3 = createMorphingParms->parm("Impl Def Attr Map 3");
    QString def4 = createMorphingParms->parm("Impl Def Attr Map 4");
    QString def5 = createMorphingParms->parm("Impl Def Attr Map 5");
    QString def6 = createMorphingParms->parm("Impl Def Attr Map 6");

    QStringList defList = QStringList() << def1 << def2 << def3 << def4 << def5 << def6;
    // for(int i = 0; i < createMorphingParms.size(); i++){
    //   defList << createMorphingParms[i];
    //   qDebug() << "Morph " << createMorphingParms[i] << "\n";
    // }



    std::vector<AttrMap<vertex, Point3d> > deformedPosVector, deformedPosVector2;

    std::vector<AttrMap<vertex, double> > deformedPosVectorSignal;

    std::vector<AttrMap<int, double> > deformedPosVectorHeat;

    //createDeformedPosVec(m, deformedPosVector, def1, def2, def3, def4, def5);
    //if(morphOther) createDeformedPosVec(m2, deformedPosVector2, def1, def2, def3, def4, def5);

    createDeformedPosVec(m, deformedPosVector, defList);
    if(morphOther) createDeformedPosVec(m2, deformedPosVector2, defList);

    //createDeformedSignalVec(m, deformedPosVectorSignal, def1, def2, def3, def4, def5);

    createDeformedSignalVec(m, deformedPosVectorSignal, defList);

    createDeformedHeatVec(m, deformedPosVectorHeat, defList);

    vertex vFix, vAxis, vFix2, vAxis2;
    Point3d normalPos, normalAxis, normalPos2, normalAxis2;
    //Point3d

    std::set<int> allLabels;
    forall(vertex v, S){
      allLabels.insert(v->label);
    }

    readFixVertices(m, vFix, vAxis, vAxis2);
    if(vAxis){
      normalAxis = vAxis->nrml;
      normalPos = vAxis->pos;
    }
    if(vAxis2 and vFix){
      normalAxis2 = vAxis2->pos - vFix->pos;
      Point3d pNrml = vFix->nrml;
        if(vAxis){
          pNrml = vAxis->nrml;
        }
      Point3d proj = projectPointOnPlane(vAxis2->pos, vFix->pos, pNrml);
      //Point3d orth = normalAxis % normalAxis2;
      //normalAxis2 = orth % normalAxis;
      normalAxis2 = proj - vFix->pos;
      normalAxis2 /= norm(normalAxis2);
    }



    //normalPos2 = vAxis2->pos;
    // if(morphOther){
    //   readFixVertices(m2, vFix2, vAxis2);
    //   normalAxis2 = vAxis2->nrml;
    //   normalPos2 = vAxis2->pos;
    // }

    int deformationStep = 0;

    IntFloatAttr heatMapT0;

    if(heatmap == "Growth" or heatmap == "Growth Backwards" or heatmap == "Cell Area"){
      // calc initial cell sizes
      MeasureArea ma(*this);

      ma.run(m, heatMapT0);
    }



    Matrix4d rotMatrixS1, rotMatrixStep, updatedMat,rotMatrixS2, updatedMat2;
    s->getFrame().getMatrix(rotMatrixS1.data());
    s2->getFrame().getMatrix(rotMatrixS2.data());
    updatedMat = rotMatrixS1;
    updatedMat2 = rotMatrixS2;


    for(uint j = 1; j < deformedPosVector.size(); j++){
    //forall(ImplicitDeformation id, morphingDeformations){

      std::cout << "s " << deformedPosVector.size() << "/" << deformedPosVector[deformationStep].size() << std::endl;

      std::cout << "Implicit Deformation " << ++deformationStep << std::endl;

      // morph mesh
      for(int i = 0; i < steps; i++){

        std::cout << "step " << i+1 << " running ...." << std::endl;
        //auto startTime = std::chrono::system_clock::now();

        if(!progressAdvance(1))
          userCancel();

        std::map<vertex, Point3d> displacementMap;

        if(useSpline){
          morphMeshStepFixSpline(S, deformedPosVector, deformationStep, i, steps, vFix, vAxis, displacementMap);
          if(morphOther) morphMeshStepFixSpline(S2, deformedPosVector2, deformationStep, i, steps, vFix, vAxis, displacementMap);
        } //else if(vFix){
          //morphMeshStepFix(S, deformedPosVector, deformationStep, i, steps, vFix, fixDispl, newDispl);
        //}
        else{
          morphMeshStep(S, deformedPosVector, deformationStep, i, steps);
        }

        // if
        if(vAxis or vFix){
          // rotate
          m->updateCentersNormals();

          Point3d targetNormal(1,0,0);
          Point3d actualNormal(1,0,0);

          Point3d targetNormal2(0,1,0);
          Point3d actualNormal2(0,1,0);

          Point3d dispVAxis(0,0,0);


          if(vAxis){
            calcNormal(S, vAxis, vAxis->nrml);
            targetNormal = normalAxis;
            actualNormal = vAxis->nrml;
            dispVAxis = displacementMap[vAxis];

          std::cout << "normals " << targetNormal << "/" << actualNormal << std::endl;

          Matrix3d rotMat1 = calcRotMatrix(actualNormal, targetNormal);

          forall(const vertex& v, S){
            v->pos = v->pos - dispVAxis - normalPos;
            v->pos = rotMat1 * v->pos;
            v->pos = v->pos + dispVAxis + normalPos;
          }

          }

       if(vAxis2 and vFix){
        Point3d pNrml = vFix->nrml;
        if(vAxis){
          calcNormal(S, vAxis, vAxis->nrml);
          pNrml = vAxis->nrml;
        }
            Point3d proj = projectPointOnPlane(vAxis2->pos, vFix->pos, pNrml);
            //Point3d orthNew = dir % normalAxis;
            //dir = orthNew % normalAxis;
            targetNormal2 = normalAxis2;
            actualNormal2 = proj - vFix->pos;
            actualNormal2 /= norm(actualNormal2);

          std::cout << "normals2 " << targetNormal2 << "/" << actualNormal2 << std::endl;

          Matrix3d rotMat2 = calcRotMatrix(actualNormal2, targetNormal2);
          Point3d vFixPos = vFix->pos;

          forall(const vertex& v, S){
            v->pos = v->pos - vFixPos;
            v->pos = rotMat2 * v->pos;
            v->pos = v->pos + vFixPos;
          }

        pNrml = vFix->nrml;
        if(vAxis){
          calcNormal(S, vAxis, vAxis->nrml);
          pNrml = vAxis->nrml;
        }
            proj = projectPointOnPlane(vAxis2->pos, vFix->pos, pNrml);
            //orthNew = dir % normalAxis;
            //dir = orthNew % normalAxis;
            targetNormal2 = normalAxis2;
            actualNormal2 = proj - vFix->pos;
            actualNormal2 /= norm(actualNormal2);

          }

          // if(vAxis){
          //   calcNormal(S, vAxis, vAxis->nrml);
          //   std::cout << "applied " << vAxis->nrml << "/" << actualNormal2 << "/" << (vAxis2->pos - vFix->pos)/norm(vAxis2->pos - vFix->pos) << std::endl;
          // }
          if(vFix){
            Point3d vFixPos = vFix->pos;
            forall(const vertex& v, S){
              v->pos = v->pos - vFixPos + deformedPosVector[0][vFix];
            }
          }
        }

        if(signalBlending){

          forall(const vertex& v, S){
            double signalBefore = 0;
            if(deformedPosVectorSignal.size() > deformationStep-1)
              signalBefore = deformedPosVectorSignal[deformationStep-1][v];
            double signalAfter = 0;
            if(deformedPosVectorSignal.size() > deformationStep)
              signalAfter = deformedPosVectorSignal[deformationStep][v];
            v->signal = ((i+1) * signalAfter + (steps-i-1) * signalBefore)/((double)steps);
          }

        }

        if(heatmap == "Growth" or heatmap == "Growth Backwards" or heatmap == "Cell Area"){
          // calc new cell sizes & update heat
          MeasureArea ma(*this);
          IntFloatAttr heatMapT1, change;
          ma.run(m, heatMapT1);

          forall(IntFloatPair p, heatMapT0){
            if(heatMapT0[p.first] == 0 or heatMapT1[p.first] == 0) continue;

            change[p.first] = heatMapT1[p.first]/heatMapT0[p.first];
            if(heatmap == "Growth Backwards") change[p.first] = 1/change[p.first];
            if(heatmap == "Cell Area") change[p.first] = heatMapT1[p.first];
          }
          heatMapT0 = heatMapT1;

          m->labelHeat() = change;

          m->showHeat();
          m->heatMapBounds() = m->calcHeatMapBounds();
          if(heatMin != 0 and heatMax != 0)
            m->heatMapBounds() = Point2f(std::min(heatMin, heatMax), std::max(heatMin, heatMax));
          m->setUseParents(false);

        } else if(heatmap == "Attribute (TimeLapse)"){
          // load heat from attr map

          m->labelHeat().clear();

          forall(int l, allLabels){

            int mapCounter = deformationStep-1;
            if(i>steps/2){
              mapCounter++;
            }

            double heatBefore = 0;
            if(deformedPosVectorHeat.size() > mapCounter-1)
              heatBefore = deformedPosVectorHeat[mapCounter-1][l];
            double heatAfter = 0;
            if(deformedPosVectorHeat.size() > mapCounter)
              heatAfter = deformedPosVectorHeat[mapCounter][l];

            if(j==1 and i<=steps/2){
              m->labelHeat()[l] = heatAfter;
              if(heatAfter == 0) m->labelHeat().erase(l);
            } else if(j==deformedPosVector.size()-1 and i>steps/2){
              m->labelHeat()[l] = heatBefore;
              if(heatBefore == 0) m->labelHeat().erase(l);
            } else if(i<=steps/2){
              m->labelHeat()[l] = ((i+steps/2) * heatAfter + (steps-i-steps/2) * heatBefore)/((double)steps);
              if(heatBefore == 0 or heatAfter == 0) m->labelHeat().erase(l);
            } else if(i>steps/2){
              m->labelHeat()[l] = ((i-steps/2) * heatAfter + (steps-i+steps/2) * heatBefore)/((double)steps);
              if(heatBefore == 0 or heatAfter == 0) m->labelHeat().erase(l);
            }

          }

          m->showHeat();
          m->heatMapBounds() = m->calcHeatMapBounds();
          if(heatMin != 0 and heatMax != 0)
            m->heatMapBounds() = Point2f(std::min(heatMin, heatMax), std::max(heatMin, heatMax));
          m->setUseParents(false);
        } else if(heatmap == "Attribute (Standard)"){
          // load heat from attr map

          m->labelHeat().clear();

          forall(int l, allLabels){

            int mapCounter = deformationStep-1;

            double heatBefore = 0;
            //std::cout << "Heat Blending " << deformedPosVectorHeat.size() << "/" << mapCounter << "/" << mapCounter << "/" << (deformedPosVectorHeat.size() >= mapCounter) << std::endl;
            if(deformedPosVectorHeat.size() >= mapCounter){
              //std::cout << "here1 " << std::endl;
              if(mapCounter == 0){
                AttrMap<int, double>& origHeat = m->attributes().attrMap<int, double>("OrigHeat");
                heatBefore = origHeat[l];
                //std::cout << "here " << l << "/" << heatBefore << std::endl;
              } else {
                heatBefore = deformedPosVectorHeat[mapCounter-1][l];
                //std::cout << "here2 " << std::endl;
              }
            }

            double heatAfter = 0;
            if(deformedPosVectorHeat.size() > mapCounter)
              heatAfter = deformedPosVectorHeat[mapCounter][l];

            if(heatBefore != 0 and heatAfter != 0)
              m->labelHeat()[l] = ((i) * heatAfter + (steps-i) * heatBefore)/((double)steps);
            std::cout << "Heat Blending " << l << "/" << deformedPosVectorHeat.size() << "/" << mapCounter 
            << "/" << i << "/" << steps << "/" << heatAfter << "/" << heatBefore << std::endl;

          }

          m->showHeat();
          m->heatMapBounds() = m->calcHeatMapBounds();
          if(heatMin != 0 or heatMax != 0)
            m->heatMapBounds() = Point2f(std::min(heatMin, heatMax), std::max(heatMin, heatMax));
          m->setUseParents(false);
        }


        m->updateAll();
        if(morphOther) m2->updateAll();

        // rotation
    double xRotStep = rotX * 3.14159 / 180. / (double)(steps*deformedPosVector.size());
    double yRotStep = rotY * 3.14159 / 180. / (double)(steps*deformedPosVector.size());
    double zRotStep = rotZ * 3.14159 / 180. / (double)(steps*deformedPosVector.size());

    Matrix4d matX = Matrix4d::rotation(Point4d(1,0,0,1), xRotStep);
    Matrix4d matY = Matrix4d::rotation(Point4d(0,1,0,1), yRotStep);
    Matrix4d matZ = Matrix4d::rotation(Point4d(0,0,1,1), zRotStep);

        rotMatrixStep = matX * matY * matZ;
        updatedMat = updatedMat * rotMatrixStep;
        updatedMat2 = updatedMat2 * rotMatrixStep;
        s->getFrame().setFromMatrix(updatedMat.data());
        s2->getFrame().setFromMatrix(updatedMat2.data());

        updateState();
        updateViewer();

        if(j == deformedPosVector.size()-1 and (int)i == steps-1) break;
        std::this_thread::sleep_for(std::chrono::milliseconds(msStep));
      }
      //fixDispl += newDispl;

    }

    s->getFrame().setFromMatrix(rotMatrixS1.data());
    s2->getFrame().setFromMatrix(rotMatrixS2.data());

    return true;
  }
  REGISTER_PROCESS(MorphMeshMulti);



  bool ResetMeshMorphing::run(Mesh* m, bool resetMesh, bool resetTransf)
  {

    vvGraph& S = m->graph();

    if(resetMesh){
      AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");

      AttrMap<vertex, double>& origSignal = m->attributes().attrMap<vertex, double>("OrigSignal");

      AttrMap<int, double>& origHeat = m->attributes().attrMap<int, double>("OrigHeat");

      forall(const vertex& v, S){
        v->pos = origMeshPos[v];
        v->signal = origSignal[v];
        m->labelHeat()[v->label] = origHeat[v->label];
      }
    } else if(resetTransf){
      // TODO
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(ResetMeshMorphing);

  bool CreateImplicitDeformation2D::run(Mesh* m, Mesh* m2, vvGraph& T1, vvGraph& T2, QString direction, int cgSteps, QString projInput, double distanceThreshold)
  {
    std::set<int> ignoreLabels;
    return run(m, m2, T1, T2, direction, cgSteps, projInput, distanceThreshold, ignoreLabels);
  }

  bool CreateImplicitDeformation2D::run(Mesh* m, Mesh* m2, vvGraph& T1, vvGraph& T2, QString direction, int cgSteps, QString projInput, double distanceThreshold, std::set<int> ignoreLabels,bool updateNormals)
  {

   std::cout << "start" << std::endl;
    QTime currentProcessTime;
    currentProcessTime.start();
    int time0 = currentProcessTime.elapsed();
    int time1,time2,time3,time4,time5,time6=0;
    int time1a,time1b,time1c,time1d,time1e,time1f,time2a,time2b,time2c=0;

   // get labels and parents
   //vvGraph& S = m->graph();
   vvGraph& S2 = m2->graph();

   if(updateNormals){
    m->updateCentersNormals();
    m2->setUseParents(true);
    m2->updateCentersNormals();
   }

    std::cout<<"Deformation Time 1: "<< currentProcessTime.elapsed() - time0<<std::endl;
    time1 = currentProcessTime.elapsed();


   JunctionInformation jInfoNew;
   impDef.clear();
   jInfo = jInfoNew;

   createJunctionAssociation(m, m2, T1, T2, jInfo);
   std::cout << "found corresponding junctions" << std::endl;

    std::cout<<"Deformation Time 2: "<< currentProcessTime.elapsed() - time1<<std::endl;
    time2 = currentProcessTime.elapsed();

   // get cell centers of existing parents in mesh 2

   std::map<int, Point3d> labelPosMap;

   IntPoint3fAttr centers = m->labelCenter();
   IntPoint3fAttr centers2 = m2->labelCenter();
   IntPoint3fAttr centersP2 = m2->parentCenter();

   std::vector<Point3d> data, values;

   std::map<int, Point3d> cellCentersM1, cellCentersM2;

   forall(const vertex& v, S2){
     if(v->label < 1) continue;
     cellCentersM2[v->label] = Point3d(centers2[v->label]);

     labelPosMap[m2->parents()[v->label]] = Point3d(centersP2[m2->parents()[v->label]]);
   }

   forall(const vertex& v, T1){
     if(v->label < 1) continue;
  	 cellCentersM1[v->label] = Point3d(centers[v->label]);
  	 v->selected = true;
   }

   jInfo.labelDaughtersCenterMap = labelPosMap;
   jInfo.labelParentCenterMap = cellCentersM1;

   //qDebug() << "bal " << projInput << "\n";
   //std::cout << "bla2 "<< labelPosMap.size() << "/" << cellCentersM1.size() << "/" << ignoreLabels.size() << std::endl;

   writeInputVectors(projInput, data, values, jInfo, distanceThreshold, ignoreLabels);
   if(stringToBool(parm("Use Manual Landmarks"))){
     AttrMap<int, Point3d>& landmarks_m1 = m->attributes().attrMap<int, Point3d>("Manual Landmarks");
     AttrMap<int, Point3d>& landmarks_m2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

     if(landmarks_m1.empty())
       return setErrorMessage("No manual landmarks found!");
     
     forall(auto p, landmarks_m1){
       data.push_back(p.second);
       values.push_back(landmarks_m2[p.first]);
     }
   }

    std::cout<<"Deformation Time 3: "<< currentProcessTime.elapsed() - time2<<std::endl;
    time3 = currentProcessTime.elapsed();

   std::cout << "size of data/value vector: " << data.size() << "/" << values.size() << std::endl;

   if(data.empty() or values.empty()){
     std::cout << "Data Vector emtpy! Process aborted, no deformation function generated!" << std::endl;
     return setErrorMessage("Data Vector emtpy! Process aborted, no deformation function generated!");
   } 

   if(direction == "active -> other"){
     impDef.initialize(data,values);
   } else if(direction == "other -> active"){
     impDef.initialize(values,data);
   } else {
     return setErrorMessage("Direction not defined properly");
   }


    std::cout<<"Deformation Time 4: "<< currentProcessTime.elapsed() - time3<<std::endl;
    time4 = currentProcessTime.elapsed();

   if(cgSteps > 0)
     impDef.solve(cgSteps);
   else
     impDef.solveInv();


   if(direction == "active -> other"){
      AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
      impDefAttr.clear();
      impDefAttr[0] = impDef;
      impDef.initialize(data,values);
   } else if(direction == "other -> active"){
     AttrMap<int, ImplicitDeformation>& impDefAttrOther = m2->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
     impDefAttrOther.clear();
     impDefAttrOther[0] = impDef;
   }

    std::cout<<"Deformation Time 5: "<< currentProcessTime.elapsed() - time4<<std::endl;
    time5 = currentProcessTime.elapsed();


   return true;
  }
  REGISTER_PROCESS(CreateImplicitDeformation2D);


  bool CreateImplicitDeformation3D::run(Mesh* m, Mesh* m2, QString direction, int cgSteps, bool inputCenters, bool inputWalls,
   bool inputWallsOut, bool inputJunctions, bool useManualLandmarks, bool considerExplosion)
  {

   if(!inputCenters and !inputWalls and !inputWallsOut and !inputJunctions and !useManualLandmarks) return setErrorMessage("At least one input type needed");

   bool forceRecalcAnalysis = false;

   // get labels and parents
   //vvGraph& S = m->graph();
   vvGraph& S2 = m2->graph();

  AttrMap<int, Point3d>& explDispl = m->attributes().attrMap<int, Point3d>("Explosion Displacement");
  AttrMap<int, Point3d>& centroidsM1 = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
  AttrMap<IntInt, double>& wallArea = m->attributes().attrMap<IntInt, double>("Shared Wall Areas");

  AttrMap<int, Point3d>& centroidsM2 = m2->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
  AttrMap<int, Point3d>& centroidsM2P = m2->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroidsParents");

  // if(centroidsM1.empty() or forceRecalcAnalysis){

  // }

  // if(centroidsM2.empty() or forceRecalcAnalysis){

  // }

  // if(wallArea.empty() or forceRecalcAnalysis){

  // }


  // check parents of m2
  std::map<int,std::set<int> > parentLabelsMap;
  forall(const vertex& v, S2){
    if(m2->parents()[v->label] < 1) continue;
    parentLabelsMap[m2->parents()[v->label]].insert(v->label);
  }

  forall(const vertex& v, S2){
    if(m2->parents()[v->label] < 1) continue;
    parentLabelsMap[m2->parents()[v->label]].insert(v->label);
  }

  //std::cout << "bla " << centroidsM1.size() << "/" << centroidsM2.size() << std::endl;

  std::map<int, set<int> > parentDaughterMap;


  forall(IntInt p, m2->parents()){
    if(p.second < 1) continue;

    //std::cout << "ppp " << p.first << "/" << p.second << std::endl;
    parentDaughterMap[p.second].insert(p.first);
  }

  //std::cout << "parent " << parentDaughterMap.size() << std::endl;

  std::vector<Point3d> data, values;


  // fill data vectors
  if(inputCenters){

  forall(IntIntSetP p, parentDaughterMap){
    if(p.second.empty()) continue;

    Point3d avgPos(0,0,0);
    forall(int l, p.second){
      avgPos += centroidsM2[l];
      //std::cout << "parent l " << p.first << "/" << l << "//" << centroidsM2[l] << std::endl;
    }
    avgPos /= p.second.size();

    if(norm(centroidsM1[p.first]) < 1E-3 and norm(avgPos)< 1E-3) continue;

    data.push_back(centroidsM1[p.first]);
    values.push_back(avgPos);
    //std::cout << "CENTROIDS " << centroidsM1[p.first] << " / " << avgPos << std::endl;
  }

  }


  AttrMap<IntInt, Point3d>& wallCenters = m->attributes().attrMap<IntInt, Point3d>("Cell Wall Centers");
  AttrMap<IntInt, Point3d>& wallCenters2 = m2->attributes().attrMap<IntInt, Point3d>("Cell Wall Centers");

  if(inputWalls){ // inside walls

    forall(auto p, wallCenters){
      IntInt walls = p.first;
      Point3d wallCenter = p.second;



      if(walls.first == 0 or walls.second == 0) continue;
      if(walls.first > walls.second) continue; // TODO check if needed

      std::set<int> daughters1 = parentDaughterMap[walls.first];
      std::set<int> daughters2 = parentDaughterMap[walls.second];

      Point3d avgPos(0,0,0);
      int counter = 0;
      forall(int l, daughters1){
        forall(int l2, daughters2){
          IntInt wall2 = std::make_pair(l,l2);
          if(l2 < l) wall2 = std::make_pair(l2,l);
          if(wallCenters2.find(wall2) == wallCenters2.end()) continue;
          counter++;
          avgPos += wallCenters2[wall2];
        }
      }

      if(counter == 0) continue;
      avgPos /= (double)counter;

      data.push_back(wallCenter);
      values.push_back(avgPos);
      std::cout << "INSIDE WALL "<< walls.first << "/" << walls.second << "//" << wallCenter << "/" << avgPos << std::endl;
    }

  }

  if(inputWallsOut){ // outside walls

    forall(auto p, wallCenters){
      IntInt walls = p.first;
      Point3d wallCenter = p.second;

      if(walls.first == 0) continue;
      if(walls.second != 0) continue;

      // std::cout << "OUTWALL "<< walls.first << "/" << walls.second << std::endl;

      // walls.first should now be the label and walls.second == 0

      std::set<int> daughters = parentDaughterMap[walls.first];

      // std::cout << "daughters "<< daughters.size() << std::endl;

      Point3d avgPos(0,0,0);
      int counter = 0;
      forall(int l, daughters){
        IntInt wall = std::make_pair(l,0);
        if(wallCenters2.find(wall) == wallCenters2.end()) continue;
        counter++;
        avgPos += wallCenters2[wall];
      }

      // std::cout << "counter "<< counter << std::endl;

      if(counter == 0) continue;
      avgPos /= (double)counter;

      data.push_back(wallCenter);
      values.push_back(avgPos);
      std::cout << "OUTSIDE WALL "<< walls.first << "/" << walls.second << "//" << wallCenter << "/" << avgPos << std::endl;
    }

  }

  AttrMap<Point4i, Point3d>& junctionsPos = m->attributes().attrMap<Point4i, Point3d>("Junction Positions");
  AttrMap<Point4i, Point3d>& junctionsPos2 = m2->attributes().attrMap<Point4i, Point3d>("Junction Positions");
  typedef std::pair<Point4i, Point3d> P4iP3dP;
  // if(inputJunctions and !junctionsPos.empty() and !junctionsPos2.empty()){

  //   // forall(P4iP3dP p, junctionsPos){

  //   // }

  // }

  if(useManualLandmarks){
    AttrMap<int, Point3d>& landmarks_m1 = m->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& landmarks_m2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    if(landmarks_m1.empty())
      return setErrorMessage("No manual landmarks found!");
     
    forall(auto p, landmarks_m1){
      data.push_back(p.second);
      values.push_back(landmarks_m2[p.first]);
    }
  }

  double threshold = parm("Point Distance Threshold").toDouble();
  if(threshold > 0){
    std::vector<Point3d> dataNew, valuesNew;

    // check for points that are too close

    // this is n squared and should be optimized!
    for(int i = 0; i < data.size(); i++){
      bool add = true;
      for(int j = i; j < data.size(); j++){
        if(j==i) continue;
        double dis = norm(data[i] - data[j]);
        if(dis < threshold){
          add = false;
        }
      }
      if(add){
        dataNew.push_back(data[i]);
        valuesNew.push_back(values[i]);
      }
    }

    std::cout << "reduced data vector due to threshold from: " << data.size() << " to " << dataNew.size() << std::endl;

    data = dataNew;
    values = valuesNew;


  }

   // create projection
   if(direction == "active -> other"){
     impDef.initialize(data,values);
   } else if(direction == "other -> active"){
     impDef.initialize(values,data);
   } else {
     return setErrorMessage("Direction not defined properly");
   }

   if(cgSteps > 0)
     impDef.solve(cgSteps);
   else
     impDef.solveInv();


   AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
   impDefAttr.clear();

   impDefAttr[0] = impDef;

   return true;
  }
  REGISTER_PROCESS(CreateImplicitDeformation3D);

   bool SetCorrectParents::run(Mesh* m, Mesh* m2)
  {

    AttrMap<int,int>& fixedParents = m2->attributes().attrMap<int, int>("Fixed Correct Parents");

    fixedParents.clear();

    forall(const vertex& v, m2->graph()){
      if(m2->parents()[v->label] > 0) fixedParents[v->label] = m2->parents()[v->label];
    }

    return true;
  }
  REGISTER_PROCESS(SetCorrectParents);

   bool RestoreCorrectParents::run(Mesh* m, Mesh* m2)
  {
    AttrMap<int,int>& fixedParents = m2->attributes().attrMap<int, int>("Fixed Correct Parents");

    m2->parents().clear();

    forall(IntInt p, fixedParents){
      m2->parents()[p.first] = p.second;
    }

    m2->updateAll();

    return true;
  }
  REGISTER_PROCESS(RestoreCorrectParents);

   bool CompareParents::run(Mesh* m, Mesh* m2, bool selectErrors)
  {

    int total = 0, missingIn1 = 0, missingIn2 = 0, errors = 0, correct = 0;

    std::set<int> errorLabel;

    forall(IntInt p, m->parents()){
      if(p.second == 0) continue;

      total++;
      if(m2->parents().find(p.first) == m2->parents().end()){
        missingIn1++;
        continue;
      }

      if(m2->parents()[p.first] == 0){
        missingIn1++;
        continue;
      }

      if(p.second != m2->parents()[p.first]){
        errorLabel.insert(p.first);
        errors++;
      } else {
        correct++;
      }

    }

    forall(IntInt p, m2->parents()){
      if(p.second == 0) continue;

      if(m->parents().find(p.first) == m->parents().end()){
        missingIn2++;
        total++;
        continue;
      }

      if(m->parents()[p.first] == 0){
        missingIn2++;
        total++;
        continue;
      }
    }

    if(selectErrors){
      forall(const vertex& v, m->graph()){
        if(errorLabel.find(v->label) != errorLabel.end()) v->selected = true;
        else v->selected = false;
      }
      forall(const vertex& v, m2->graph()){
        if(errorLabel.find(v->label) != errorLabel.end()) v->selected = true;
        else v->selected = false;
      }
      m->updateAll();
      m2->updateAll();
    }

    std::cout << "Total Parent Labeled Cells: " << total << std::endl;
    std::cout << "Identical in both meshes: " << correct << std::endl;
    std::cout << "Different in the meshes: " << errors << std::endl;
    std::cout << "Missing in Mesh 1: " << missingIn1 << std::endl;
    std::cout << "Missing in Mesh 2: " << missingIn2 << std::endl;

    return true;
  }
  REGISTER_PROCESS(CompareParents);

   bool SetMorphOrigin::run(Mesh* m)
  {

    m->updateCentersNormals();
    IntPoint3fAttr centers = m->labelCenter();

    AttrMap<vertex,Point3d>& morphOrig = m->attributes().attrMap<vertex, Point3d>("Morphing Origin");
    morphOrig.clear();

    forall(const vertex& v, m->graph()){
      if(!v->selected or v->label < 1) continue;
      morphOrig[v] = v->pos;
      std::cout << "Origin set: " << v->label << "/" << v->pos << std::endl;
      return true;
    }

    return setErrorMessage("No cell selected!");
  }
  REGISTER_PROCESS(SetMorphOrigin);

   bool SetMorphAxis::run(Mesh* m)
  {

    m->updateCentersNormals();

    AttrMap<vertex,Point3d>& morphAxis = m->attributes().attrMap<vertex, Point3d>("Morphing Axis1");
    morphAxis.clear();

    forall(const vertex& v, m->graph()){
      if(!v->selected or v->label < 1) continue;
      morphAxis[v] = v->nrml;
      std::cout << "Axis set: " << v->label << "/" << v->nrml << std::endl;
      return true;
    }

    return true;
  }
  REGISTER_PROCESS(SetMorphAxis);

   bool SetMorphAxis2::run(Mesh* m)
  {

    m->updateCentersNormals();

    AttrMap<vertex,Point3d>& morphAxis = m->attributes().attrMap<vertex, Point3d>("Morphing Axis2");
    morphAxis.clear();

    forall(const vertex& v, m->graph()){
      if(!v->selected or v->label < 1) continue;
      morphAxis[v] = v->nrml;
      std::cout << "Axis set: " << v->label << "/" << v->nrml << std::endl;
      return true;
    }

    return true;
  }
  REGISTER_PROCESS(SetMorphAxis2);


   bool DeformationParentLabeling::run(Mesh* m, Mesh* m2, vvGraph& T1, vvGraph& T2, bool useFixedParents, QString findCells, double omitCells, bool correctErrors,std::map<IntInt, double> * _neighMapS1)
  {
    QTime currentProcessTime;
    currentProcessTime.start();
    int time0 = currentProcessTime.elapsed();
    int time1,time2,time3,time4,time5,time6=0;
    int time1a,time1b,time1c,time1d,time1e,time1f,time2a,time2b,time2c=0;


   std::cout<<"Time 1a: "<< currentProcessTime.elapsed() - time0<<std::endl;
   time1a = currentProcessTime.elapsed();


   if(_neighMapS1!=NULL)//SHOULD BE INDICATED AS A PARAMETER - CHANGING TO INCREASE SPEED
       m->updateCentersNormals();

   m2->setUseParents(true);
   m2->updateCentersNormals();

    std::cout<<"Time 1b: "<< currentProcessTime.elapsed() - time1a<<std::endl;
    time1b = currentProcessTime.elapsed();


   vvGraph& S = m->graph();
   vvGraph& S2 = m2->graph();

    std::cout<<"Time 1c: "<< currentProcessTime.elapsed() - time1b<<std::endl;
    time1c = currentProcessTime.elapsed();


   std::map<IntInt, double> neighMapS1;
   if(_neighMapS1 == NULL)
       neighborhood2D(S, neighMapS1);
   else
       neighMapS1 = (*_neighMapS1);


    std::cout<<"Time 1d: "<< currentProcessTime.elapsed() - time1c<<std::endl;
    time1d = currentProcessTime.elapsed();


   DistanceWaterfront dw(*this);
   IntFloatAttr heatMap;
   bool updateNormals = false;
   dw.run(m, m2, neighMapS1, heatMap, updateNormals);

    std::cout<<"Time 1e: "<< currentProcessTime.elapsed() - time1e<<std::endl;
    time1e = currentProcessTime.elapsed();


   std::set<int> omitLabels;

   std::cout << "omits" << std::endl;

   forall(IntFloatPair p, heatMap){
     if(p.second > 2){
       if(std::rand() % 100 < omitCells){
         omitLabels.insert(p.first);
         std::cout << "omit " << p.first << std::endl;
       }
     }
   }
    std::cout<<"Time 1f: "<< currentProcessTime.elapsed() - time1e<<std::endl;
    time1f = currentProcessTime.elapsed();

    std::cout<<"Time 1: "<< currentProcessTime.elapsed() - time0<<std::endl;
    time1 = currentProcessTime.elapsed();

   // create deformation
   CreateImplicitDeformation2D *createImpDef;

    if(!getProcess("Mesh/Deformation/Mesh 2D/Create Deformation Cell Surface", createImpDef))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateImplicitDeformation2D process"));

   AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
   //bool forceRecalc = stringToBool(parms[4]);

   //if(impDefAttr.empty() or forceRecalc)
   bool res = createImpDef->run(m, m2, T1, T2, createImpDef->parm("Direction"), createImpDef->parm("CG Steps").toInt(), createImpDef->parm("Projection Input"),  createImpDef->parm("Junction Merging Threshold (um)").toDouble());

   if(!res){
     return setErrorMessage("Error when creating the deformation!");
   }

    std::cout<<"Time 2a: "<< currentProcessTime.elapsed() - time1<<std::endl;
    time2a = currentProcessTime.elapsed();

   ImplicitDeformation& impDef = impDefAttr[0];//createImpDef->impDef;
   JunctionInformation& jInfo = createImpDef->jInfo;

   AttrMap<int, int> fixedParents = m2->attributes().attrMap<int, int>("Fixed Correct Parents");

   bool missing = false;

   std::map<int, Point3d> cellCentersM1, cellCentersM2;
   std::map<int, Point3d> cellNormalsM2;

  forall(const vertex& v, T2){
    if(v->label < 1) continue;
    cellCentersM2[v->label] = v->pos;
    cellNormalsM2[v->label] = v->nrml;
    if(m2->parents()[v->label] < 1){
      missing = true;
      continue;
    }
  }

   if(!missing) return true;

  // go through the neighbors of already parent labelled cells in mesh 1 and find their deforemd pos
  // then in mesh 2 find the two closest cells to this position
  //forall neighbors

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

  std::set<int> newCells;

  forall(IntPoint3dP p, jInfo.labelDaughtersCenterMap){
    int parentLabel = p.first;
    forall(IntIntDoubleP np, neighMapS1){
      if(np.first.first == parentLabel){
        if(jInfo.labelDaughtersCenterMap.find(np.first.second) == jInfo.labelDaughtersCenterMap.end()){
          newCells.insert(np.first.second);
        }
      }
    }
  }

  std::cout<<"Time 2b: "<< currentProcessTime.elapsed() - time2a<<std::endl;
  time2b = currentProcessTime.elapsed();

  std::cout<<"Time 2: "<< currentProcessTime.elapsed() - time1<<std::endl;
  time2 = currentProcessTime.elapsed();

  std::set<vertex> junctionsToBeFound;


  if(findCells == "Nearest Center" or findCells == "Nearest Center Threshold"){

  std::map<int,Point3d> centers2;// = m2->labelCenter();
  forall(const vertex& v, S2){
    if(v->label < 1) continue;
    centers2[v->label] = v->pos;
  }



  // then in mesh 2 find the two closest cells to this position
  forall(int l, newCells){
    Point3d deformedPos = impDef.RBFdeformation(Point3d(jInfo.labelParentCenterMap[l]));

    // project the deformed pos onto the mesh for better accuracy
    Point3d deformedPosMesh = nearestPointOnMesh(deformedPos, m2->graph());

    double minDis1 = HUGE_VAL, minDis2 = HUGE_VAL;
    int label1 = -1, label2 = -1;
    forall(IntPoint3dP p, centers2){
      double dis = norm(deformedPosMesh - Point3d(p.second));
      if(dis < minDis1){
        minDis2 = minDis1;
        minDis1 = dis;
        label2 = label1;
        label1 = p.first;
      } else if(dis < minDis2){
        minDis2 = dis;
        label2 = p.first;
      }
    }

    std::cout << "r " << l << "/" << label1 << "/" << label2 << "/" << minDis1 << "//" << minDis2 << "//" << minDis2/minDis1 << std::endl;

    if(findCells == "Nearest Center"){
      m2->parents()[label1] = l;
      continue;
    }

    if(minDis2/minDis1 > 5){
      m2->parents()[label1] = l;
    } else {
      m2->parents()[label1] = l;
      m2->parents()[label2] = l;
    }

  }

  } else { // findCells == Junction Projection


  forall(int l, newCells){
  	// get junctions in order; project them into m2-space and then onto a plane
    std::vector<vertex> junctionsInOrder = jInfo.cellJunctionInOrder1[l];
    std::vector<vertex> junctionsInOrderM2(junctionsInOrder.size());
    std::vector<Point3d> junctionPointsInOrder;

    for(uint i = 0; i < junctionsInOrder.size(); i++){
      vertex v;
      v->selected = true;
      v->pos = impDef.RBFdeformation(junctionsInOrder[i]->pos);
      junctionPointsInOrder.push_back(v->pos);
      junctionsInOrderM2[i] = v;
      //if(l == 277) S2.insert(v);
    }

    Point3d planePos, planeNrml;
    findPolygonPlane(junctionPointsInOrder, planePos, planeNrml);

    // forall(Point3d &p, junctionPointsInOrder){
    //   vertex v2;
    //   v2->pos = p;
    //   v2->selected = true;
    //   if(l == 266) S2.insert(v2);
    // }

    Matrix3d rotMat;
    rotatePointsIntoXY(planeNrml, junctionPointsInOrder, rotMat);

    forall(Point3d &p, junctionPointsInOrder){
      p.z() = 0;
      // vertex v2;
      // v2->pos = p;
      // v2->selected = true;
      // if(l == 266) S2.insert(v2);
    }

    // check which cell centers are inside
    forall(IntPoint3dP p, cellCentersM2){
      if(m2->parents()[p.first] > 0) continue;
      Point3d cellCenter = Point3d(p.second);
    //forall(IntPoint3dP p, labelPosMap){
      Point3d projectedCellCenter;
      double sN;
      planeLineIntersect(planePos, planeNrml, cellCenter, cellCenter-cellNormalsM2[p.first], sN, projectedCellCenter);

      // if(l == 266){
      //   vertex v2;
      //   v2->pos = projectedCellCenter;
      //   v2->selected = false;
      //   v2->label = p.first;
      //   v2->type = 'l';
      //   S2.insert(v2);
      // }

      std::vector<Point3d> projectedCellCenterV;
      projectedCellCenterV.push_back(projectedCellCenter);
      rotatePointsIntoXY(planeNrml, projectedCellCenterV, rotMat);

      // if(l == 266){
      //   projectedCellCenterV[0].z() = 0;
      //   // vertex v2;
      //   // v2->pos = projectedCellCenterV[0];
      //   // v2->selected = false;
      //   // v2->label = p.first;
      //   // v2->type = 'l';
      //   // S2.insert(v2);
      // }

      // label accordingly


  	  if(pointInPolygon(projectedCellCenterV[0], junctionPointsInOrder)){
  	  	std::cout << "new cell found " << l << "/" << p.first << std::endl;
  	    m2->parents()[p.first] = l;
  	  }
  	}

  }

  }

    std::cout<<"Time 3: "<< currentProcessTime.elapsed() - time2<<std::endl;
    time3 = currentProcessTime.elapsed();

    CorrespondenceJunctions cc(*this);
    cc.checkNhbdNew(m, m2, T1, T2, false);


    std::cout<<"Time 4: "<< currentProcessTime.elapsed() - time3<<std::endl;
    time4 = currentProcessTime.elapsed();

    forall(vertex v, S2){
      if(v->label < 1) continue;
      if(m2->labelHeat()[m2->parents()[v->label]] == 1){
        if(correctErrors)
          m2->parents()[v->label] = 0;
      }
    }

    if(useFixedParents){
      forall(IntInt p, fixedParents){
        m2->parents()[p.first] = p.second;
      }
    }

    m->updateAll();
    m2->updateAll();

    std::cout<<"Time 5: "<< currentProcessTime.elapsed() - time4<<std::endl;
    time5 = currentProcessTime.elapsed();

    return true;
  }
  REGISTER_PROCESS(DeformationParentLabeling);


   bool AutoFillParents::run(Mesh* m, Mesh* m2, int steps, int fixDistance, double retreat)
  {

    if(retreat < 0) retreat = 0;
    else if(retreat > 50) retreat = 50;

    AttrMap<int, int>& fixedParents = m2->attributes().attrMap<int, int>("Fixed Correct Parents");


    vvGraph& S = m->graph();
    vvGraph& S2 = m2->graph();
    m->updateCentersNormals();

    vvGraph T1, T2;
    mgxmToMgxc(S, T1, 0);
    mgxmToMgxc(S2, T2, 0);
    std::map<IntInt, double> neighMapS1;
    neighborhood2D(S, neighMapS1);


    std::map<int, int> parentsBefore, parentsAfter;


    progressStart(QString("Morphing mesh-%1").arg(m->userId()), steps);

  for(int i = 0; i < steps; i++){
    if(!progressAdvance(i))
      userCancel();

    DeformationParentLabeling *proc;

    if(!getProcess("Mesh/Deformation/Mesh 2D/Auto Parent Labeling Step 2D", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DeformationParentLabeling process"));

    proc->run(m, m2, T1, T2, stringToBool(proc->parm("Use Fixed Parents")), proc->parm("Find Cells"), proc->parm("Omit Cells %").toDouble(), stringToBool(proc->parm("Correct Errors")),&neighMapS1);

    // fix parents that are deep enough in the tissue, they can assumed to be correct
    if(fixDistance > 0 or retreat > 0){

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

      DistanceWaterfront dw(*this);
      IntFloatAttr heatMap;
      dw.run(m, m2, neighMapS1, heatMap);

      forall(IntFloatPair p, heatMap){
        if(retreat > 0 and retreat >= std::rand() % 100){
          forall(const vertex& v, S2){
            if(m2->parents()[v->label] == 1){
              m2->parents()[v->label] = 0;
            }
          }
        } else if(p.second > fixDistance and fixDistance > 0){
          forall(const vertex& v, S2){
          	if(m2->parents()[v->label] == p.first){
          	  fixedParents[v->label] = p.first;
          	}
          }
          //fixedParents[p.first];
        }
      }

    }
    m2->setShowLabel("Label");
    m->updateTriangles();
    m2->updateTriangles();

    updateState();
    updateViewer();

  }

    return true;
  }
  REGISTER_PROCESS(AutoFillParents);

   bool DeformationParentLabeling3D::run(Mesh* m, Mesh* m2, int cgStep, bool autoCorrect, bool considerExplosion)
  {
    bool useFixedParents = true;
    bool forceRecalcAnalysis = false;


  	AttrMap<int, int>& fixedParents = m2->attributes().attrMap<int, int>("Fixed Correct Parents");

    AttrMap<int, Point3d>& explDispl = m->attributes().attrMap<int, Point3d>("Explosion Displacement");
    AttrMap<int, Point3d>& centroidsM1 = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    AttrMap<IntInt, double>& wallArea = m->attributes().attrMap<IntInt, double>("Shared Wall Areas");

    AttrMap<int, Point3d>& centroidsM2 = m2->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    AttrMap<int, Point3d>& centroidsM2P = m2->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroidsParents");

    if(centroidsM1.empty() or forceRecalcAnalysis){

    }

    if(centroidsM2.empty() or forceRecalcAnalysis){

    }

    if(wallArea.empty() or forceRecalcAnalysis){

    }

    // if(cellAtlasCoords and cellAtlasAttr1->empty()){
    //   return setErrorMessage("No Cell Atlas Data Present");
    // }

    vvGraph& S = m->graph();
    vvGraph& S2 = m2->graph();

    int selectedLabel = -1;
    forall(const vertex &v, S) {
      //lvMap1[v->label].push_back(v);
      if(v->selected) selectedLabel = v->label;
    }

    // check parents of m2
    std::map<int,std::set<int> > parentLabelsMap;
    forall(const vertex& v, S2){
      if(m2->parents()[v->label] < 1) continue;
      parentLabelsMap[m2->parents()[v->label]].insert(v->label);
    }

    std::cout << "bla " << centroidsM1.size() << "/" << centroidsM2.size() << std::endl;

    // create projection

    std::map<int, set<int> > parentDaughterMap;


    forall(IntInt p, m2->parents()){
      if(p.second < 1) continue;

      std::cout << "ppp " << p.first << "/" << p.second << std::endl;
      parentDaughterMap[p.second].insert(p.first);
    }

    std::set<int> cellsMissingParents, cellsMissingDaughters;

    forall(IntPoint3dP p, centroidsM2){
      if(m2->parents()[p.first] < 1){
        cellsMissingParents.insert(p.first);
      }
    }

    forall(IntPoint3dP p, centroidsM1){
      if(!parentDaughterMap[p.first].empty()) continue;
      cellsMissingDaughters.insert(p.first);
    }

    std::set<int> cellsToCheck;

    forall(IntIntDoubleP p, wallArea){
      if(!parentDaughterMap[p.first.first].empty() and cellsMissingDaughters.find(p.first.second) != cellsMissingDaughters.end()){
        cellsToCheck.insert(p.first.second);
      }
    }

    // forall(const vertex& v, S){
    //   if(cellsToCheck.find(v->label) != cellsToCheck.end()){
    //   	v->selected = true;
    //   }
    // }

    CreateImplicitDeformation3D *createImpDef;

    if(!getProcess("Mesh/Deformation/Mesh 3D/Create Deformation Cell Volumes", createImpDef))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateImplicitDeformation3D process"));

    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
    bool forceRecalc = stringToBool(createImpDef->parm("Force Recalculation"));


    if(impDefAttr.empty() or forceRecalc)
      createImpDef->run(m, m2, createImpDef->parm("Direction"), createImpDef->parm("CG Steps").toInt(), 
        stringToBool(createImpDef->parm("Input Centers")), stringToBool(createImpDef->parm("Input Inside Walls")), 
        stringToBool(createImpDef->parm("Input Outside Walls")), stringToBool(createImpDef->parm("Input Junctions")), 
        stringToBool(createImpDef->parm("Use Manual Landmarks")), stringToBool(createImpDef->parm("Force Recalculation")));

    ImplicitDeformation impDef = impDefAttr[0];

    std::cout << "sel Cell " << selectedLabel << "/" << cellsToCheck.size() << std::endl;
    m2->labelHeat().clear();

    std::map<int, std::map<int, double> > cellLabelParentsCandidatesMaps;

    forall(int l, cellsToCheck){
      std::cout << "check Cell " << l << std::endl;
      Point3d centerProjected = impDef.RBFdeformation(centroidsM1[l]);

      double minDis1 = HUGE_VAL, minDis2 = HUGE_VAL;
      int label1 = -1, label2 = -1;
      forall(IntPoint3dP p, centroidsM2){
        if(m2->parents()[p.first] > 0) continue;
        double dis = norm(centerProjected - Point3d(p.second));
        if(l == selectedLabel){
          m2->labelHeat()[p.first] = dis;
          std::cout << "c " << p.first << "/" << dis << "/" << l << std::endl;
        }
        if(dis < minDis1){
          minDis2 = minDis1;
          minDis1 = dis;
          label2 = label1;
          label1 = p.first;
        } else if(dis < minDis2){
          minDis2 = dis;
          label2 = p.first;
        }
      }

      if(minDis2/minDis1 > 5 and minDis1 < 2){
        //m2->parents()[label1] = l;
        cellLabelParentsCandidatesMaps[label1][l] = minDis1;
        std::cout << "new daughters " << label1 << std::endl;
      } else if(minDis2/minDis1 <= 5 and minDis1 < 10 and minDis2 < 10){
      	//m2->parents()[label1] = l;
      	cellLabelParentsCandidatesMaps[label1][l] = minDis1;
      	//m2->parents()[label2] = l;
      	cellLabelParentsCandidatesMaps[label2][l] = minDis2;
      	std::cout << "new daughters " << label1 << "/" << label2 << std::endl;
      }

    }

    typedef std::pair<int, std::map<int, double> > IntIntDMapP;
    typedef std::pair<int, double> IntDoubleP;
    forall(IntIntDMapP p, cellLabelParentsCandidatesMaps){
      double minDis = 1E20;
      int parentLabel = -1;
      forall(IntDoubleP dp, p.second){
        if(minDis > dp.second){
          minDis = dp.second;
          parentLabel = dp.first;
        }
      }
      m2->parents()[p.first] = parentLabel;
    }

    // forall(IntPoint3dP p, rcp1.rootData.cellCentroids){
    //   if(parentDaughterMap[p.first].size()
    //   //vertex n;
    //   Point3d centerProjected = impDef.RBFdeformation(p.second);

    //   S2.insert(n);
    // }

    m->updateAll();

    m2->updateAll();

    if(!autoCorrect) return true;

    // QStringList parms;
    // CheckCorrespondence3D *checkCor = getProcessParms<CheckCorrespondence3D>(this, parms);
    // if(!checkCor)
    //   throw(QString("Unable to find process"));

    // CheckCorrespondence3D cp(*this);
    // cp.run(m,m2,info1,info2, parms[0].toDouble(), parms[1], parms[2].toDouble(), parms[3].toDouble());

    // double errorThreshold = parms[3].toDouble();

    // forall(vertex v, S2){
    //   if(v->label < 1) continue;
    //   if(m2->labelHeat()[m2->parents()[v->label]] >= errorThreshold){
    //     m2->parents()[v->label] = 0;
    //   }
    // }

    // if(useFixedParents){
    //   forall(IntInt p, fixedParents){
    //     m2->parents()[p.first] = p.second;
    //   }
    // }

    m->updateAll();

    m2->updateAll();


    return true;
  }
  REGISTER_PROCESS(DeformationParentLabeling3D);

   bool DeformMesh3D::run(Mesh* m, Mesh* m2, bool forceRecalc)
  {

    vvGraph& S = m->graph();
    vvGraph& S2 = m2->graph();

    AttrMap<vertex, Point3d>& origMeshPos = m->attributes().attrMap<vertex, Point3d>("OrigMeshPos");


    origMeshPos.clear();
    forall(const vertex& v, S){
      origMeshPos[v] = v->pos;
    }

    CreateImplicitDeformation3D *createImpDef;

    if(!getProcess("Mesh/Deformation/Mesh 3D/Create Deformation Cell Volumes", createImpDef))
      throw(QString("MeshProcessHeatMap:: Unable to make CreateImplicitDeformation3D process"));

    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");

    if(impDefAttr.empty() or forceRecalc)
      createImpDef->run(m, m2, createImpDef->parm("Direction"), createImpDef->parm("CG Steps").toInt(), stringToBool(createImpDef->parm("Input Centers")), stringToBool(createImpDef->parm("Input Inside Walls")), stringToBool(createImpDef->parm("Input Outside Walls")), stringToBool(createImpDef->parm("Input Junctions")), stringToBool(createImpDef->parm("Use Manual Landmarks")), stringToBool(createImpDef->parm("Force Recalculation")));

    ImplicitDeformation impDef = impDefAttr[0];

    if(createImpDef->parm("Direction") == "active -> other"){
      #pragma omp parallel for
      for(uint i = 0; i < S.size(); i++){
        const vertex& v = S[i];
        v->pos = impDef.RBFdeformation(v->pos);
      }
      m->updateAll();
    } else {
      #pragma omp parallel for
      for(uint i = 0; i < S2.size(); i++){
        const vertex& v = S2[i];
        v->pos = impDef.RBFdeformation(v->pos);
      }
      m2->updateAll();
    }

    return true;
  }
  REGISTER_PROCESS(DeformMesh3D);

  // compares the neighborhoods of cells in m and the the parents of m2
  // creates a heat map with the numbers of wrong neighbors
  bool CheckCorrespondence3D::run(Mesh* m, Mesh* m2, double minSharedWall, QString type, double factor, double threshold)
  {

    AttrMap<IntInt, double>& wallAreas1 = m->attributes().attrMap<IntInt, double>("Shared Wall Areas");
    AttrMap<int, double>& attrDataWallArea1 = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");

    if(wallAreas1.empty()){
      double vertexTolerance = 0.001;
      NhbdGraphInfo info1;
      neighborhoodGraph(m->graph(), vertexTolerance, info1);

      forall(IntIntDoubleP p, info1.sharedWallArea){
        wallAreas1[p.first] = p.second;
      }
    }

    AttrMap<IntInt, double>& wallAreas2 = m2->attributes().attrMap<IntInt, double>("Shared Wall Areas");
    AttrMap<int, double>& attrDataWallArea2 = m2->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");

    if(wallAreas2.empty()){
      double vertexTolerance = 0.001;
      NhbdGraphInfo info2;
      neighborhoodGraph(m2->graph(), vertexTolerance, info2);

      forall(IntIntDoubleP p, info2.sharedWallArea){
        wallAreas2[p.first] = p.second;
      }
    }


    std::map<int, std::set<int> > neighborsM1, neighborM2Parent;
    std::map<int, int> parentErrorCounter;
    std::map<int, std::set<int> > parentDaughterMap;

    if(factor < 0) factor = 1;

    // create map from parent to labels (m2)
    forall(IntInt p, m2->parents()){
      if(p.second < 1) continue;
      parentDaughterMap[p.second].insert(p.first);
    }


    if(type == "Neighbors"){

      // create map from label to neighbor labels (m1)
      forall(IntIntDoubleP p, wallAreas1){
        if(p.second < minSharedWall) continue;
        neighborsM1[p.first.first].insert(p.first.second);
      }

      // create map from parents to neighbor parent labels (m2)
      forall(IntIntDoubleP p, wallAreas2){
        if(p.second < minSharedWall) continue;
        if(m2->parents()[p.first.second] > 0)
          neighborM2Parent[m2->parents()[p.first.first]].insert(m2->parents()[p.first.second]);
      }

      // count neighborhood errors by checking if all of the neighbors of mesh1 are also in mesh2
      forall(IntIntSetP p, neighborsM1){
        forall(int l, neighborM2Parent[p.first]){
        	if(l == p.first) continue; // ignore same parent (cells divide and have same parents, shouldn't count as error)
        	if(p.second.find(l) == p.second.end()){
        	  parentErrorCounter[l]++;
        	  parentErrorCounter[p.first]++;
            //std::cout << "invalid neighbors " << p.first << "/" << l << std::endl;
        	}
        }
      }

      // vice versa
      forall(IntIntSetP p, neighborM2Parent){
        forall(int l, neighborsM1[p.first]){
          if(l == p.first) continue; // ignore same parent (cells divide and have same parents, shouldn't count as error)
          if(p.second.find(l) == p.second.end()){
            parentErrorCounter[l]++;
            parentErrorCounter[p.first]++;
            //std::cout << "invalid neighbors2 " << p.first << "/" << l << std::endl;
          }
        }
      }

      // fill heat maps
      forall(const vertex& v, m->graph()){
        if(v->label < 1) continue;
        if(parentDaughterMap.find(v->label) == parentDaughterMap.end()) continue;
        // divide by 2 as we counted all errors twice twice
        m->labelHeat()[v->label] = parentErrorCounter[v->label]/2;
        m2->labelHeat()[v->label] = parentErrorCounter[v->label]/2;
      }

    } else {

      std::map<int, double> labelErrorMap;

      std::map<IntInt, double> sharedWallArea2Parents;

      std::map<int, double> cellArea2Parents;

      // first convert the labels of mesh2 into parents for the sharedWallAreas
      forall(IntIntDoubleP p, wallAreas2){
        IntInt parentL = std::make_pair(m2->parents()[p.first.first], m2->parents()[p.first.second]);
        sharedWallArea2Parents[parentL] = p.second;

        if(parentL.first == parentL.second){
          cellArea2Parents[parentL.first] -= 2*p.second;
        }
      }

      // same for the total cell areas
      typedef std::pair<int, double> IntDoubleP;
      forall(IntDoubleP p, attrDataWallArea2){
        cellArea2Parents[m2->parents()[p.first]] += p.second;
      }

      // now compare the two maps and add up the differences
      forall(IntIntDoubleP p, wallAreas1){
        if(p.first.first < p.first.second) continue;
        double wall2 = 0;
        if(sharedWallArea2Parents.find(p.first) != sharedWallArea2Parents.end()){
          wall2 = sharedWallArea2Parents[p.first];
        }
        double cellGrowth = cellArea2Parents[p.first.first]/attrDataWallArea1[p.first.first];
        double dif1 =  std::abs(wall2-p.second*cellGrowth)/cellArea2Parents[p.first.first];
        if(wall2 > minSharedWall and p.second > minSharedWall) dif1 /= factor;
        labelErrorMap[p.first.first] += dif1;

        double cellGrowth2 = cellArea2Parents[p.first.second]/attrDataWallArea1[p.first.second];
        double dif2 = std::abs(wall2-p.second*cellGrowth2)/cellArea2Parents[p.first.second];
        if(wall2 > minSharedWall and p.second > minSharedWall) dif2 /= factor;
        labelErrorMap[p.first.second] += dif2;

      }

      // fill heat maps
      forall(const vertex& v, m->graph()){
        if(v->label < 1) continue;
        if(parentDaughterMap.find(v->label) == parentDaughterMap.end()) continue;
        m->labelHeat()[v->label] = labelErrorMap[v->label];
        m2->labelHeat()[v->label] = labelErrorMap[v->label];
      }

    }

    // visualize them
    m->heatMapBounds() = Point2f(0.,threshold);
    m2->heatMapBounds() = Point2f(0.,threshold);

    m->showHeat();
    m2->showHeat();

    m->updateTriangles();
    m2->updateTriangles();

    updateViewer();

  	return true;
  }
  REGISTER_PROCESS(CheckCorrespondence3D);


   bool DistanceWaterfront::run(Mesh* m, Mesh* m2, std::map<IntInt, double>& neighMapS1, IntFloatAttr& heatMap,bool updateNormals)
  {

   if(updateNormals){
     m->updateCentersNormals();
     m2->setUseParents(true);
     m2->updateCentersNormals();
   }

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

   // get cell centers of existing parents in mesh 2

   std::set<int> cellsM1, cellsM2, alreadyParents, notParents;

   std::map<int, int> cellsDistanceM2;

  std::map<int, Point3d> labelPosMap;

  IntPoint3fAttr centers = m->labelCenter();
  IntPoint3fAttr centersP2 = m2->parentCenter();

  std::vector<Point3d> data, values;

  std::map<int, Point3d> cellCentersM1, cellCentersM2;

  forall(const vertex& v, S2){
    if(v->label < 1) continue;
    cellCentersM2[v->label] = v->pos;
    cellsM2.insert(v->label);
    //if(m2->parents()[v->label] < 1){
    //  missing = true;
    //  continue;
    //}
    cellsM1.insert(m2->parents()[v->label]);
    alreadyParents.insert(m2->parents()[v->label]);

    // calc cell center
    //Point3d cellCenter = v->pos;
    //initialParents[v->label] = m2->parents()[v->label];
    labelPosMap[m2->parents()[v->label]] = Point3d(centersP2[m2->parents()[v->label]]);
  }


  // find the cells directly at the waterfront
  //std::cout << "cells " << cellsM1.size() << std::endl;
  std::set<int> toDel;
  forall(int l, cellsM1){
  	forall(IntIntDoubleP p, neighMapS1){
  	  if(p.first.first == l and alreadyParents.find(p.first.second) == alreadyParents.end()){
  	  	//std::cout << "neigh " << l << "/" << p.first.second << "/" << m2->parents()[p.first.second] << std::endl;
  	  	cellsDistanceM2[l] = 1;
  	  	toDel.insert(l);
  	  }
  	}
  }
  forall(int l, toDel){
    cellsM1.erase(l);
  }

  //std::cout << "cells " << cellsM1.size() << std::endl;
  // now label the inside
  int d = 1;
  while(!cellsM1.empty()){

    std::set<int> toDel;
    bool changed = false;

    forall(int l, cellsM1){
      //if(cellsDistanceM2[l] == d){
        // go through neighbors
      forall(IntIntDoubleP p, neighMapS1){
        if(p.first.first == l and cellsDistanceM2[p.first.second] == d){
          cellsDistanceM2[l] = d+1;
          toDel.insert(l);
          changed = true;
        }
      }
      //}
    }

    forall(int l, toDel){
      cellsM1.erase(l);
    }

    d++;

    if(!changed) cellsM1.clear();

  }

    forall(IntInt p, cellsDistanceM2){
      m2->labelHeat()[p.first] = p.second;
    }

    m2->heatMapBounds() = m2->calcHeatMapBounds();
    m2->updateTriangles();

    heatMap = m2->labelHeat();

    return true;
  }
  REGISTER_PROCESS(DistanceWaterfront);


   bool FindIsolatedParents3D::run(Mesh* m2, NhbdGraphInfo& info2, std::set<int>& isolatedParents, bool selectP)
  {

    std::map<int, std::set<int> > neighborsParent;

    // create map from label to neighbor labels (m1)
    forall(IntIntDoubleP p, info2.sharedWallArea){
      if(m2->parents()[p.first.first] < 1) continue;

      neighborsParent[p.first.first].insert(p.first.second);
    }

    std::set<int> isolatedParentsMap;

    forall(IntIntSetP p, neighborsParent){
      forall(int l, p.second){
      	if(l > 0) continue;
      	isolatedParentsMap.insert(p.first);
      }
    }

    if(selectP){
      forall(const vertex& v, m2->graph()){
        if(isolatedParentsMap.find(v->label) == isolatedParentsMap.end()) continue;
        v->selected = true;
      }
      m2->updateAll();
    }
    isolatedParents = isolatedParentsMap;

    return true;
  }
  REGISTER_PROCESS(FindIsolatedParents3D);


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

   bool ExportDeformation::run(Mesh* m, Mesh* m2, QString exportTo, QString filename)
  {

    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");

    if(impDefAttr.empty()) return setErrorMessage("No Implicit Deformation present. Create one first!");

    ImplicitDeformation& impDef = impDefAttr[0];

    std::vector<Point3d>& data = impDef.data;
    std::vector<Point3d>& values = impDef.values;
    std::vector<std::vector<double> >& rb_coeffs = impDef.rb_coeffs;
    std::vector<std::vector<double> >& poly_coeffs = impDef.poly_coeffs;

    if(data.empty() or values.empty() or rb_coeffs.empty() or poly_coeffs.empty()) return setErrorMessage("The deformation attribute map is invalid or empty. Create a new one first!");

    int length = data.size();

    if(exportTo == "CSV File"){
      // file for root cell data
      QFile file(filename);
      if(!file.open(QIODevice::WriteOnly))
      {
        std::cout << "File cannot be opened for writing";
        return false;
      }
      QTextStream out(&file);

      out.setRealNumberPrecision(32);

      out << "data X,data Y,data Z,values X,values Y,values Z,rb coeffs X,rb coeffs Y,rb coeffs Z,poly coeffs X,poly coeffs Y,poly coeffs Z";
      out << endl;

      for(int i = 0; i<std::max(length,4); i++){
        if(length > 4 or i < 4){
          out << data[i].x() << "," << data[i].y() << "," << data[i].z() << ",";
          out << values[i].x() << "," << values[i].y() << "," << values[i].z() << ",";
          out << rb_coeffs[0][i] << "," << rb_coeffs[1][i] << "," << rb_coeffs[2][i] << ",";
          //std::cout << "d " << i << "/" << data[i].x() << "/" << values[i].x() << "/" << rb_coeffs[0][i] << std::endl;
        }
        if(i < 4){
          out << poly_coeffs[0][i] << "," << poly_coeffs[1][i] << "," << poly_coeffs[2][i] << ",";
          //std::cout << "p " << i << "/" << poly_coeffs[0][i] << std::endl;
        }

        out << endl;
      }

    } else if(exportTo == "Attribute Map"){
      AttrMap<int, ImplicitDeformation>& impDefExportAttr = m->attributes().attrMap<int, ImplicitDeformation>(filename);

      impDefExportAttr.clear();

      impDefExportAttr = impDefAttr;
    }

    std::cout << "Exported Deformation with " << impDef.data.size() << " data points" << std::endl;

    return true;
  }
  REGISTER_PROCESS(ExportDeformation);


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

   bool ImportDeformation::run(Mesh* m, Mesh* m2, QString importFrom, QString filename)
  {
    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");
    ImplicitDeformation impDef;

    impDef.initCoeffs();

    std::vector<Point3d>& data = impDef.data;
    std::vector<Point3d>& values = impDef.values;
    std::vector<std::vector<double> >& rb_coeffs = impDef.rb_coeffs;
    std::vector<std::vector<double> >& poly_coeffs = impDef.poly_coeffs;

    if(importFrom == "CSV File"){

      // open file
      QFile file(filename);
      if(!file.open(QIODevice::ReadOnly))
      {
        setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
        return false;
      }
      QTextStream ss(&file);
      QString line = ss.readLine();
      QStringList fields = line.split(",");

      bool ok;
      int lineCounter = 0;

      while(ss.status() == QTextStream::Ok)
      {
        fields = ss.readLine().split(",");
        if(fields[0].isEmpty()) break;

        std::cout << "l " << lineCounter << std::endl;

        double dataX = fields[0].toDouble(&ok);
        double dataY = fields[1].toDouble(&ok);
        double dataZ = fields[2].toDouble(&ok);

        double valuesX = fields[3].toDouble(&ok);
        double valuesY = fields[4].toDouble(&ok);
        double valuesZ = fields[5].toDouble(&ok);

        double rbX = fields[6].toDouble(&ok);
        double rbY = fields[7].toDouble(&ok);
        double rbZ = fields[8].toDouble(&ok);

        data.push_back(Point3d(dataX, dataY, dataZ));
        values.push_back(Point3d(valuesX, valuesY, valuesZ));

        std::cout << "d " << dataX << std::endl;

        rb_coeffs[0].push_back(rbX);
        rb_coeffs[1].push_back(rbY);
        rb_coeffs[2].push_back(rbZ);

        std::cout << "r " << rbX << "/" << rbY << "/" << rbZ << std::endl;

        if(lineCounter < 4){
          double polyX = fields[9].toDouble(&ok);
          double polyY = fields[10].toDouble(&ok);
          double polyZ = fields[11].toDouble(&ok);

          poly_coeffs[0][lineCounter] = polyX;
          poly_coeffs[1][lineCounter] = polyY;
          poly_coeffs[2][lineCounter] = polyZ;

          std::cout << "p " << polyX << "/" << polyY << "/" << polyZ << std::endl;
        }

        lineCounter++;
      }

    impDefAttr[0] = impDef;

    } else if(importFrom == "Attribute Map"){
      AttrMap<int, ImplicitDeformation>& impDefExportAttr = m->attributes().attrMap<int, ImplicitDeformation>(filename);

      if(impDefExportAttr.empty()) return setErrorMessage("Attribute Map doesn't exist!");

      for(uint i = 0; i < impDefExportAttr[0].data.size(); i++){
        std::cout << "l " << impDefExportAttr[0].data[i] << std::endl;
        std::cout << "d " << impDefExportAttr[0].values[i] << std::endl;
        std::cout << "r " << impDefExportAttr[0].rb_coeffs[0][i] << std::endl;
        if(i < 4) std::cout << "p " << impDefExportAttr[0].poly_coeffs[0][i] << std::endl;
      }

      impDefAttr = impDefExportAttr;
    }

    std::cout << "Imported Deformation with " << impDef.data.size() << " data points" << std::endl;

    return true;
  }
  REGISTER_PROCESS(ImportDeformation);

/*
  void calcEigenVecsValues(Matrix3d& M, Point3d& p1, Point3d& p2, Point3d& p3, Point3d& ev){

    gsl_matrix* mat = gsl_matrix_alloc(3, 3);
    gsl_vector* eval = gsl_vector_alloc(3);
    gsl_matrix* evec = gsl_matrix_alloc(3, 3);
    gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc(3);


    for(int i = 0; i < 3; ++i)
      for(int j = 0; j < 3; ++j)
        gsl_matrix_set(mat, i, j, M[i][j]);

    gsl_eigen_symmv(mat, eval, evec, w);
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

    p1 = Point3d(gsl_matrix_get(evec, 0, 0), gsl_matrix_get(evec, 1, 0), gsl_matrix_get(evec, 2, 0));
    p2 = Point3d(gsl_matrix_get(evec, 0, 1), gsl_matrix_get(evec, 1, 1), gsl_matrix_get(evec, 2, 1));
    p3 = Point3d(gsl_matrix_get(evec, 0, 2), gsl_matrix_get(evec, 1, 2), gsl_matrix_get(evec, 2, 2));
    ev = Point3d(gsl_vector_get(eval, 0), gsl_vector_get(eval, 1), gsl_vector_get(eval, 2));

  }

*/
   SymmetricTensor convertToSymTensor(Matrix3d t)
   {
    SymmetricTensor st;

    st.ev1() = Point3f(t[0][0], t[0][1], t[0][2]);
    st.ev2() = Point3f(t[1][0], t[1][1], t[1][2]);
    st.evals() = Point3f(t[2][0], t[2][1], t[2][2]);

    return st;
   }

   Matrix3d convertToMatTensor(SymmetricTensor st)
   {
    Matrix3d tensor;

    tensor[0][0] = st.ev1().x();
    tensor[0][1] = st.ev1().y();
    tensor[0][2] = st.ev1().z();

    tensor[1][0] = st.ev2().x();
    tensor[1][1] = st.ev2().y();
    tensor[1][2] = st.ev2().z();

    tensor[2][0] = st.evals().x();
    tensor[2][1] = st.evals().y();
    tensor[2][2] = st.evals().z();

    return tensor;
   }


   SymmetricTensor deformationTensorOfPosition(ImplicitDeformation& impDef, Point3d& pos, Point3d nrml = Point3d(0,0,0))
   {
    SymmetricTensor st;

    Matrix3d t;

    if(nrml != Point3d(0,0,0)){
      t = impDef.calcProjectedGrowthTensorOnPos(pos, nrml);
    } else {
      t = impDef.calcGrowthTensor(pos);
    }

    st.ev1() = Point3f(t[0][0], t[0][1], t[0][2]);
    st.ev2() = Point3f(t[1][0], t[1][1], t[1][2]);
    st.evals() = Point3f(t[2][0], t[2][1], t[2][2]);

    return st;
   }


   SymmetricTensor convertTensorToPDG(SymmetricTensor tensor)
   {
    SymmetricTensor pdg;

    Matrix3d t = convertToMatTensor(tensor);

    // calc Eigenvalues
    Point3d p1, p2, p3, ev;
    calcEigenVecsValues(t, p1, p2, p3, ev);

    pdg.ev1() = Point3f(p1.x(), p1.y(), p1.z());
    pdg.ev2() = Point3f(p2.x(), p2.y(), p2.z());
    pdg.evals() = Point3f(sqrt(ev.x()), sqrt(ev.y()), sqrt(ev.z()));

    return pdg;
   }

   SymmetricTensor convertTensorToPDG3D(SymmetricTensor tensor, Point3d& evals)
   {
    SymmetricTensor pdg;

    Matrix3d t = convertToMatTensor(tensor);

    // calc Eigenvalues
    Point3d p1, p2, p3, ev;
    calcEigenVecsValues(t, p1, p2, p3, ev);

    evals = Point3d(sqrt(ev.x()), sqrt(ev.y()), sqrt(ev.z()));

    pdg.ev1() = Point3f(p1.x(), p1.y(), p1.z()) * evals.x();
    pdg.ev2() = Point3f(p2.x(), p2.y(), p2.z()) * evals.y();
    pdg.evals() = Point3f(p3.x(), p3.y(), p3.z()) * evals.z();

    return pdg;
   }


   bool ComputeDeformationGradient::run(Mesh* m, Mesh* m2, QString mode, QString modeCell, QString modeVtx, bool projectOnCell, bool otherMesh, bool validP)
  {

    vvGraph& S = m->graph();

    IntSymTensorAttr& cellAxis = m->cellAxis();
    IntSymTensorAttr& cellAxis2 = m2->cellAxis();

    VtxSymTensorAttr& vertexAxis = m->vertexAxis();
    VtxSymTensorAttr& vertexAxis2 = m2->vertexAxis();

    // get the current deformation
    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");

    AttrMap<int, SymmetricTensor>& IDAttr = m->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor DeformationGradient");
    IDAttr.clear();

    AttrMap<int, SymmetricTensor>& defGrad = m->attributes().attrMap<int, SymmetricTensor>("Deformation Gradient");
    AttrMap<vertex, SymmetricTensor>& defGradVtx = m->attributes().attrMap<vertex, SymmetricTensor>("Deformation Gradient Vertex");

    //AttrMap<int, Matrix3d>& labelTensorMap = m->attributes().attrMap<int, Matrix3d>("Deformation Gradient");
    //labelTensorMap.clear();

    if(impDefAttr.empty())
      return setErrorMessage("No Deformation Attribute Map present. Please create or import one.");

    ImplicitDeformation& impDef = impDefAttr[0];

    std::map<int, std::set<int> > parentLabelMap;

    forall(IntInt p, m2->parents()){
      parentLabelMap[p.second].insert(p.first);
    }

    m->setCellAxisType("PDG");


    if(mode == "Cell"){
      cellAxis.clear();
      cellAxis2.clear();
      defGrad.clear();

      std::map<int, std::set<vertex> > cellContour = getCellContour(S);
      std::map<int, Point3d> labelCellNormalMap;


      // do pca on junctions to find cell plane
      forall(auto p, cellContour){
        std::vector<Point3d> pointsPolygon;
        forall(vertex v, p.second){
          pointsPolygon.push_back(v->pos);
        }
        Point3d planePos, planeNrml;

        findPolygonPlane(pointsPolygon, planePos, planeNrml);

        labelCellNormalMap[p.first] = planeNrml;
      }


      if(modeCell == "Center"){
        m->updateCentersNormals();
        IntPoint3fAttr& centers = (m->useParents() ? m->parentCenterVis() : m->labelCenterVis());

        forall(IntPoint3fP p, centers){
          int cellLabel = p.first;
          if(cellLabel < 1) continue;
          Point3d cellCenter = Point3d(p.second);
          //Matrix3d t = impDef.calcGrowthTensor(cellCenter);

          if(!validP or parentLabelMap[cellLabel].size() > 0){
            if(projectOnCell){
              defGrad[cellLabel] = deformationTensorOfPosition(impDef, cellCenter, labelCellNormalMap[cellLabel]);
            } else {
              defGrad[cellLabel] = deformationTensorOfPosition(impDef, cellCenter);
            }
          }


        }
      } else if(modeCell == "Average"){
        // average the tensor based on triangle size and then project it onto the cell

        std::map<int, double> labelAreaCounter;

        std::map<int, SymmetricTensor> labelSymTensorMap;

        forall(const vertex& v, S){
          forall(const vertex& n, S.neighbors(v)){
            vertex n2 = S.nextTo(v,n);
            if(!S.uniqueTri(v,n,n2)) continue;

            int cellLabel = m->getLabel(v,n,n2, m->parents());

            Point3d triCenter = (v->pos + n->pos + n2->pos)/3.;
            //Point3d triNrml = (v->nrml + n->nrml + n2->nrml)/3.;
            double triArea = triangleArea(v->pos, n->pos, n2->pos);
            labelAreaCounter[cellLabel] += triArea;

            SymmetricTensor st = deformationTensorOfPosition(impDef, triCenter);

            if(!validP or parentLabelMap[cellLabel].size() > 0){
              labelSymTensorMap[cellLabel].ev1() += triArea * st.ev1();
              labelSymTensorMap[cellLabel].ev2() += triArea * st.ev2();
              labelSymTensorMap[cellLabel].evals() += triArea * st.evals();
            }

          }
        }

        typedef std::pair<int, double> IntDoubleP;
        forall(IntDoubleP p, labelAreaCounter){
          int cellLabel = p.first;
          if(!validP or parentLabelMap[cellLabel].size() > 0){
            defGrad[cellLabel].ev1() = labelSymTensorMap[cellLabel].ev1()/p.second;
            defGrad[cellLabel].ev2() = labelSymTensorMap[cellLabel].ev2()/p.second;
            defGrad[cellLabel].evals() = labelSymTensorMap[cellLabel].evals()/p.second;
            Point3d cellNormal = labelCellNormalMap[cellLabel];
            if(projectOnCell){
              Matrix3d projTensor = calcProjectedGrowthTensor(convertToMatTensor(defGrad[cellLabel]), cellNormal);
              defGrad[cellLabel] = convertToSymTensor(projTensor);
            }
          }

        }
      }

      // convert the 3D tensor into the PDG format
      forall(auto p, defGrad){
        int cellLabel = p.first;
        cellAxis[cellLabel] = convertTensorToPDG(p.second);
        if(otherMesh) cellAxis2[cellLabel] = convertTensorToPDG(p.second); // TODO parent

          //std::cout << "tensor " << cellLabel << "/" << defGrad[cellLabel].ev1() << "/" << defGrad[cellLabel].ev2() << "/" << defGrad[cellLabel].evals() << std::endl;
         //std::cout << "evecs " << cellLabel << "/" << cellAxis[cellLabel].ev1() << "/" << cellAxis[cellLabel].ev2() << "/" << cellAxis[cellLabel].evals() << std::endl;
        IDAttr[p.first]=p.second;
      }

      m->setCellAxisType("PDG");
      if(otherMesh) m2->setCellAxisType("PDG");

      DisplayPDGs *proc;
      if(!getProcess("Mesh/Cell Axis/Deformation Gradient/Display Cell Gradient", proc))
        throw(QString("ImplicitDeformationProcesses:: Unable to make DisplayPDGs process"));

      proc->run(m);
      if(otherMesh) proc->run(m2);

    } else if(mode == "Vertex"){ // Vertex based PDGs (for subcellular)
      vertexAxis.clear();
      defGradVtx.clear();

      forall(const vertex& v, S){
        if(projectOnCell){
          defGradVtx[v] = deformationTensorOfPosition(impDef, v->pos, v->nrml);
        } else {
          defGradVtx[v] = deformationTensorOfPosition(impDef, v->pos);
        }
        vertexAxis[v] = convertTensorToPDG(defGradVtx[v]);
      }

      DisplayPDGsVertex* displ;
      if(!getProcess("Mesh/Cell Axis/Deformation Gradient/Display Vertex Gradient", displ))
        throw(QString("ImplicitDeformationProcesses:: Unable to make DisplayPDGsVertex process"));
      displ->run(m);

    } else {
      return setErrorMessage("Invalid Mode!");
    }

    return true;
  }
  REGISTER_PROCESS(ComputeDeformationGradient);


   bool ComputeDeformationGradient3D::run(Mesh* m, Mesh* m2, Stack* s, Store* s1, QString mode, bool projectOnCell, bool otherMesh, bool validP)
  {

    vvGraph& S = m->graph();

    // get the current deformation
    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");

    AttrMap<int, SymmetricTensor>& defGrad = m->attributes().attrMap<int, SymmetricTensor>("Deformation Gradient");
    AttrMap<vertex, SymmetricTensor>& defGradVtx = m->attributes().attrMap<vertex, SymmetricTensor>("Deformation Gradient Vertex");

    IntSymTensorAttr &cellAxis = m->cellAxis();
    VtxSymTensorAttr &vertexAxis = m->vertexAxis();

    if(impDefAttr.empty())
      return setErrorMessage("No Deformation Attribute Map present. Please create or import one.");

    ImplicitDeformation& impDef = impDefAttr[0];

    m->updateCentersNormals();
    IntPoint3fAttr centers;// = (m->useParents() ? m->parentCenterVis() : m->labelCenterVis());

    std::map<int, Point3d> labelCellNormalMap;

    // get the centroids from the attr map
    AttrMap<int, Point3d>& centroidsM1 = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    centers.clear();
    forall(IntPoint3dP p, centroidsM1){
      centers[p.first] = Point3f(p.second);
    }

    std::map<int, std::set<int> > parentLabelMap;

    forall(IntInt p, m2->parents()){
      parentLabelMap[p.second].insert(p.first);
    }

    if(mode == "Cell"){
      cellAxis.clear();
      defGrad.clear();

      forall(IntPoint3fP p, centers){
        int cellLabel = p.first;
        if(cellLabel < 1) continue;
        Point3d cellCenter = Point3d(p.second);
        //Matrix3d t = impDef.calcGrowthTensor(cellCenter);
        //defGrad[cellLabel] = deformationTensorOfPosition(impDef, cellCenter);

        if(!validP or parentLabelMap[cellLabel].size() > 0){
          defGrad[cellLabel] = deformationTensorOfPosition(impDef, cellCenter);
          Point3d evals;
          //cellAxis[cellLabel] = convertTensorToPDG(defGrad[cellLabel]);
          cellAxis[cellLabel] = convertTensorToPDG3D(defGrad[cellLabel], evals);
          std::cout << "tensor " << cellLabel << "/" << defGrad[cellLabel].ev1() << "/" << defGrad[cellLabel].ev2() << "/" << defGrad[cellLabel].evals() << std::endl;

          std::cout << "evecs " << cellLabel << "/" << evals << "/" << cellAxis[cellLabel].ev1() << "/" << cellAxis[cellLabel].ev2() << "/" << cellAxis[cellLabel].evals() << std::endl;
          cellAxis[cellLabel].ev1() /= norm(cellAxis[cellLabel].ev1());
          cellAxis[cellLabel].ev2() /= norm(cellAxis[cellLabel].ev2());
          cellAxis[cellLabel].evals() = Point3f(evals);

        }

      }

      DisplayDeformationGradient3D* displ;
      if(!getProcess("Mesh/Cell Axis 3D/Deformation Gradient/Display Cell Gradient", displ))
        throw(QString("MeshProcessHeatMap:: Unable to make DisplayDeformationGradient3D process"));

      displ->run(m);
      if(otherMesh) displ->run(m2);

    } else if(mode == "Stack"){

      const HVecUS& data = s1->data();
      Point3i imgSize(s->size());

      std::map<int, double> labelVoxelCounter;
      std::map<int, SymmetricTensor> labelSymTensorMap;

      #pragma omp parallel for schedule(guided)
      for(int z=0; z<imgSize.z(); z++){
        for(int x=0; x<imgSize.y(); x++){
          for(int y=0; y<imgSize.x(); y++){
            Point3i currentPoint(x,y,z);
            int label = data[s->offset(currentPoint)];
            if(label == 0) continue;
            Point3d pWorld(s->imageToWorld(currentPoint));
            SymmetricTensor st = deformationTensorOfPosition(impDef, pWorld);

            #pragma omp critical
            {
              labelSymTensorMap[label].ev1() += st.ev1();
              labelSymTensorMap[label].ev2() += st.ev2();
              labelSymTensorMap[label].evals() += st.evals();
              labelVoxelCounter[label]++;
            }
          }
        }
      }

      forall(auto p, labelVoxelCounter){
        int cellLabel = p.first;
        labelSymTensorMap[cellLabel].ev1() /= p.second;
        labelSymTensorMap[cellLabel].ev2() /= p.second;
        labelSymTensorMap[cellLabel].evals() /= p.second;

        Point3d evals;

        cellAxis[cellLabel] = convertTensorToPDG3D(labelSymTensorMap[cellLabel], evals);

          std::cout << "tensorS " << cellLabel << "/" << labelSymTensorMap[cellLabel].ev1() << "/" << labelSymTensorMap[cellLabel].ev2() << "/" << labelSymTensorMap[cellLabel].evals() << std::endl;

          std::cout << "evecsS " << cellLabel << "/" << evals << "/" << cellAxis[cellLabel].ev1() << "/" << cellAxis[cellLabel].ev2() << "/" << cellAxis[cellLabel].evals() << std::endl;
          cellAxis[cellLabel].ev1() /= norm(cellAxis[cellLabel].ev1());
          cellAxis[cellLabel].ev2() /= norm(cellAxis[cellLabel].ev2());
          cellAxis[cellLabel].evals() = Point3f(evals);

      }

      DisplayDeformationGradient3D* displ;
      if(!getProcess("Mesh/Cell Axis 3D/Deformation Gradient/Display Cell Gradient", displ))
        throw(QString("MeshProcessHeatMap:: Unable to make DisplayDeformationGradient3D process"));

      displ->run(m);
      if(otherMesh) displ->run(m2);

    } else if(mode == "Vertex"){
      vertexAxis.clear();
      defGradVtx.clear();

      forall(const vertex& v, S){
        if(projectOnCell){
          defGradVtx[v] = deformationTensorOfPosition(impDef, v->pos, v->nrml);
        } else {
          defGradVtx[v] = deformationTensorOfPosition(impDef, v->pos);
        }
        if(!validP or parentLabelMap[v->label].size() > 0){
          vertexAxis[v] = convertTensorToPDG(defGradVtx[v]);
        }
      }

      DisplayPDGsVertex* displ;
      if(!getProcess("Mesh/Cell Axis/Deformation Gradient/Display Vertex Gradient", displ))
        throw(QString("MeshProcessHeatMap:: Unable to make DisplayPDGsVertex process"));
      displ->run(m);
    } else if(mode == "Cell Wall"){
      // run "Label Cell Walls" for mesh1
      QStringList parms;
      LabelCellWalls* lcw1;
      if(!getProcess("Mesh/Cell Mesh/Tools/Label Cell Walls", lcw1))
        throw(QString("MeshProcessHeatMap:: Unable to make LabelCellWalls process"));
      lcw1->run(m, 0, false);

      //LabelCellWalls* lcw2 = getProcessParms<LabelCellWalls>(this, parms);
      //lcw2->run(m2, 0, false);

      RelabelFragments* rf;
      if(!getProcess("Mesh/Cell Mesh/Tools/Relabel Fragments", rf))
        throw(QString("MeshProcessHeatMap:: Unable to make RelabelFragments process"));
      rf->run(m, 0);

      // run 2D Compute & Display Process for cell average

      ComputeDeformationGradient* cdg;
      if(!getProcess("Mesh/Cell Axis/Deformation Gradient/Compute", cdg))
        throw(QString("MeshProcessHeatMap:: Unable to make ComputeDeformationGradient process"));

      cdg->run(m, m2, "Cell", "Average", "", projectOnCell, false, false);
    } else {
      return setErrorMessage("Invalid Mode");
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(ComputeDeformationGradient3D);


  REGISTER_PROCESS(DisplayDeformationGradient);
  REGISTER_PROCESS(DisplayDeformationGradientCellWall);


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

  bool SaveCellAxis3D::run(Mesh* mesh, const QString& filename, const QString& type)
  {
    QFile outputFile;
    QTextStream ts;
    if(filename.isEmpty()) {
      setErrorMessage(QString("Unable to write to output file"));
      return false;
    }

    outputFile.setFileName(filename);
    if(!outputFile.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("Unable to write to output file"));
      return false;
    }

    ts.setDevice(&outputFile);
    if(mesh->cellAxisType() == "PDG"){
      if(type == "Original Tensor")
        ts << QString("Label,xMax,yMax,zMax,xMid,yMid,zMid,xMin,yMin,zMin,lengthMax,lengthMid,lengthMin") << endl;
      else if(type == "Custom Tensor")
        ts << QString("Label,xMax,yMax,zMax,xMid,yMid,zMid,xMin,yMin,zMin,lengthMax,lengthMid,lengthMin") << endl;
      else
        ts << QString("Label,xMax,yMax,zMax,xMid,yMid,zMid,xMin,yMin,zMin,lengthMax,lengthMid,lengthMin") << endl;
    }
    else if(mesh->cellAxisType() == "PCA"){
      if(type == "Original Tensor")
        ts << QString("Label,xMax,yMax,zMax,xMid,yMid,zMid,xMin,yMin,zMin,lengthMax,lengthMid,lengthMin") << endl;
      else if(type == "Custom Tensor")
        ts << QString("Label,xMax,yMax,zMax,xMid,yMid,zMid,xMin,yMin,zMin,lengthMax,lengthMid,lengthMin") << endl;
      else
        ts << QString("Label,xMax,yMax,zMax,xMid,yMid,zMid,xMin,yMin,zMin,lengthMax,lengthMid,lengthMin") << endl;
      }
    else {
      //setErrorMessage(QString("Unknown cell axis type"));
      std::cout << "Unknwon cell axis type! Saving Default file" << std::endl;
      return false;
    }

    if(type == "Original Tensor"){

      forall(const IntSymTensorPair& p, mesh->cellAxis()) {
        int label = p.first;
        const SymmetricTensor& tensor = p.second;
        const Point3f& vMax = tensor.ev1();
        const Point3f& vMin = tensor.ev2();
        const Point3f& eig = tensor.evals();
        ts << label << "," << vMax.x() << "," << vMax.y() << "," << vMax.z() << "," << vMin.x() << "," << vMin.y()
           << "," << vMin.z() << "," << eig.x() << "," << eig.y() << "," << eig.z() << "," << norm(vMax) << "," << norm(vMin) << "," << norm(eig) << endl;
      }

    } else if(type == "Custom Tensor"){

      const IntSymTensorAttr& cellAxisCustom = currentMesh()->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
      if(cellAxisCustom.empty()) return setErrorMessage("No Custom Directions found.");

      forall(const IntSymTensorPair& p, mesh->cellAxis()) {
        int label = p.first;
        Point3f customX = cellAxisCustom[label].ev1();
        Point3f customY = cellAxisCustom[label].ev2();
        const SymmetricTensor& tensor = p.second;

        SymmetricTensor tensorCustom = customDirectionTensor(tensor, customX, customY, false);

        const Point3f& vMax = tensorCustom.ev1();
        const Point3f& vMin = tensorCustom.ev2();
        const Point3f& eig = tensorCustom.evals();
        ts << label << "," << vMax.x() << "," << vMax.y() << "," << vMax.z() << "," << vMin.x() << "," << vMin.y()
           << "," << vMin.z() << "," << eig.x() << "," << eig.y() << "," << eig.z() << endl;
      }

    } else { // save visualized tensor

      forall(const auto& p, mesh->cellAxisVis()) {
        int label = p.first;
        const Point3f& vMax = Point3f(p.second[0][0], p.second[0][1], p.second[0][2]);
        const Point3f& vMin = Point3f(p.second[1][0], p.second[1][1], p.second[1][2]);
        const Point3f& eig = Point3f(p.second[2][0], p.second[2][1], p.second[2][2]);
        ts << label << "," << vMax.x() << "," << vMax.y() << "," << vMax.z() << "," << vMin.x() << "," << vMin.y()
             << "," << vMin.z() << "," << eig.x() << "," << eig.y() << "," << eig.z() << endl;
      }
    }

    return true;
  }
  REGISTER_PROCESS(SaveCellAxis3D);



   bool DisplayDeformationGradient3D::run(Mesh* m, QString heatmapOption, QString axisOption, QColor colorPos, QColor colorNeg, double axisWidth,
    double axisScale, double axisOffset, bool projectOnSurface)
  {

    vvGraph& S = m->graph();
    QString displayAxis = axisOption;
    QString displayHeatMap = heatmapOption;
    bool attrMaps = false;

    AttrMap<int, SymmetricTensor>& customDirections = m->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    IntPoint3fAttr &centers = m->labelCenterVis();
    centers.clear();

    forall(auto p, centroids){
      centers[p.first] = Point3f(p.second);
    }

    typedef std::pair<int,Matrix3d> IntM3dP;

    IntMatrix3fAttr &cellAxisVis = m->cellAxisVis();
    IntSymTensorAttr &cellAxis = m->cellAxis();
    cellAxisVis.clear();

    std::map<int, int> parentLabelMap;

    IntVec3ColorbAttr &cellAxisColor = m->cellAxisColor();

    // Scale PDGs for visualization
    m->setAxisWidth(axisWidth);
    m->setAxisOffset(axisOffset);
    Colorb colorExpansion(colorPos);
    Colorb colorShrinkage(colorNeg);

    IntFloatAttr& labelHeatMap = m->labelHeat();
    if(displayHeatMap != "None") {
      labelHeatMap.clear();
    }

    forall(IntSymTensorPair p, cellAxis){

      int cellLabel = p.first;
      SymmetricTensor tensor = p.second;
      // convert SymmetricTensor into 3D tensor
      //Matrix3d t = convertToMatTensor(p.second);

      Point3f axisDirX, axisDirY, axisDirZ;
      float axisLengthX, axisLengthY, axisLengthZ;
      Point3f axisDirCustomX, axisDirCustomY, axisDirCustomZ;
      float axisLengthCustomX = -1, axisLengthCustomY = -1, axisLengthCustomZ = -1;

       // if any "custom" option is selected, calculate the custom directions
      if(displayAxis == "Custom X" or displayAxis == "Custom Y" or displayAxis == "Custom Z" or displayAxis == "Custom All" or
        displayHeatMap == "Custom X" or displayHeatMap == "Custom Y" or displayHeatMap == "Custom Z" or (attrMaps and !customDirections.empty())){

        if(customDirections.empty())
          throw(QString("No custom cell axis found. Run 'Cell Axis Custom Directions' first."));
        if(customDirections.find(cellLabel) == customDirections.end())
          continue;

        Point3f customX = customDirections[cellLabel].ev1();
        Point3f customY = customDirections[cellLabel].ev2();
        Point3f customZ = customDirections[cellLabel].evals();

        // length of vectors
        axisLengthCustomX = lengthCustomDirTensor3D(tensor, customX);
        axisLengthCustomY = lengthCustomDirTensor3D(tensor, customY);
        axisLengthCustomZ = lengthCustomDirTensor3D(tensor, customZ);

        // direction of vectors
        axisDirCustomX = customX;
        axisDirCustomX /= (norm(axisDirCustomX));
        axisDirCustomY = customY;
        axisDirCustomY /= (norm(axisDirCustomY));
        axisDirCustomZ = customZ;
        axisDirCustomZ /= (norm(axisDirCustomZ));
      }

      // now the non-custom cell axis
      // length of vectors
      if(displayAxis == "Custom X" or displayAxis == "Custom Y" or displayAxis == "Custom Z"  or displayAxis == "Custom All"){
        axisLengthX = axisLengthCustomX -1;
        axisLengthY = axisLengthCustomY -1;
        axisLengthZ = axisLengthCustomZ -1;

        axisDirX = axisDirCustomX;
        axisDirY = axisDirCustomY;
        axisDirZ = axisDirCustomZ;


      } else {
        axisLengthX = tensor.evals()[0] -1;
        axisLengthY = tensor.evals()[1] -1;
        axisLengthZ = tensor.evals()[2] -1;

        axisDirX = tensor.ev1();
        axisDirY = tensor.ev2();
        axisDirZ = axisDirX ^ axisDirY;
        axisDirZ /= (norm(axisDirZ));

      }

      // direction of vectors
      std::cout << "tensor disp " << cellLabel << "/" << tensor.ev1() << "/" << tensor.ev2() << "/" << tensor.evals() << std::endl;


      // set the axes
      Point3b showAxis;
      showAxis[0] = (displayAxis == "Stretch Max" or displayAxis == "Custom X" or displayAxis == "All" or displayAxis == "Custom All");
      showAxis[1] = (displayAxis == "Stretch Mid" or displayAxis == "Custom Y" or displayAxis == "All" or displayAxis == "Custom All");
      showAxis[2] = (displayAxis == "Stretch Min" or displayAxis == "Custom Z" or displayAxis == "All" or displayAxis == "Custom All");

      cellAxisVis.erase(cellLabel);

      cellAxisVis[cellLabel][0] = (showAxis[0] ? axisLengthX : 0.f) * axisDirX * axisScale;
      cellAxisVis[cellLabel][1] = (showAxis[1] ? axisLengthY : 0.f) * axisDirY * axisScale;
      cellAxisVis[cellLabel][2] = (showAxis[2] ? axisLengthZ : 0.f) * axisDirZ * axisScale;

      cellAxisColor[cellLabel][0] = axisLengthX < 0 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
      cellAxisColor[cellLabel][1] = axisLengthY < 0 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
      cellAxisColor[cellLabel][2] = axisLengthZ < 0 ? Colorb(colorShrinkage) : Colorb(colorExpansion);

      // set the heat map

      if(displayHeatMap != "None") {

        int cell = p.first;
        float value = 0;

        axisLengthX += 1;
        axisLengthY += 1;
        axisLengthZ += 1;

        if(displayHeatMap == "Stretch Max")
          value = axisLengthX; // times 2 to capture the full length of the cell (not 1/2 or the radius of it)
        else if(displayHeatMap == "Stretch Mid")
          value = axisLengthY;
        else if(displayHeatMap == "Stretch Min")
          value = axisLengthZ;
        else if(displayHeatMap == "Ratio Max/Mid")
          value = axisLengthX/axisLengthY;
        else if(displayHeatMap == "Ratio Max/Min")
          value = axisLengthX/axisLengthZ;
        else if(displayHeatMap == "Ratio Mid/Min")
          value = axisLengthY/axisLengthZ;
        else if(displayHeatMap == "Product")
          value = axisLengthX * axisLengthY * axisLengthZ;
        else if(displayHeatMap == "Anisotropy")
          value = (axisLengthX-0.5*axisLengthY-0.5*axisLengthZ)/(axisLengthX+axisLengthY+axisLengthZ);
        else if(displayHeatMap == "Custom X")
          value = axisLengthX;
        else if(displayHeatMap == "Custom Y")
          value = axisLengthY;
        else if(displayHeatMap == "Custom Z")
          value = axisLengthZ;

        labelHeatMap[cell] = value;
      }

    //   if(attrMaps){

    //     shape3dProduct[cell] = axisLengthX * axisLengthY * axisLengthZ *8;
    //     shape3dMax[cell] = axisLengthX *2;
    //     shape3dMid[cell] = axisLengthY * 2;
    //     shape3dMin[cell] = axisLengthZ *2;
    //     shape3dMaxMid[cell] = axisLengthX/axisLengthY;
    //     shape3dMaxMin[cell] = axisLengthX/axisLengthZ;
    //     shape3dMidMin[cell] = axisLengthY/axisLengthZ;
    //     shape3dAniso[cell] = (axisLengthX-0.5*axisLengthY-0.5*axisLengthZ)/(axisLengthX+axisLengthY+axisLengthZ);
    //     shape3dCustomX[cell] = axisLengthCustomX *2;
    //     shape3dCustomY[cell] = axisLengthCustomY *2;
    //     shape3dCustomZ[cell] = axisLengthCustomZ *2;

    //   }
    // }

    }

    // display it
    m->setAxisView("Cell Axis");
    m->setCellAxisType("PDG");

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

    return true;
  }
  REGISTER_PROCESS(DisplayDeformationGradient3D);

   bool DisplayDeformationGradientVertex3D::run(Mesh* m, QString heatmapOption, QString axisOption, QColor colorPos, QColor colorNeg, double axisWidth,
    double axisScale, double axisOffset, double minDis, bool projectOnSurface)
  {

    //AttrMap<vertex, Matrix3d>& vertexTensorMap = m->attributes().attrMap<vertex, Matrix3d>("Deformation Gradient Vertex");

    VtxSymTensorAttr& vertexAxis = m->vertexAxis();

    vvGraph& S = m->graph();

    // get the current deformation
    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");

    if(impDefAttr.empty())
      return setErrorMessage("No Deformation Attribute Map present. Please create or import one.");

    ImplicitDeformation& impDef = impDefAttr[0];

    // forall(const vertex& v, S){

    //   if(projectOnSurface){
    //     vertexTensorMap[v] = impDef.calcProjectedGrowthTensorOnPos(v->pos, v->nrml);
    //   } else {
    //     vertexTensorMap[v] = impDef.calcGrowthTensor(v->pos);
    //   }
    // }

    std::map<int, std::set<vertex> > labelVertexVisMap;

    //AttrMap<vertex, SymmetricTensor>& customDirAttr = m->attributes().attrMap<vertex, SymmetricTensor>("Custom Directions");

    //if(vertexTensorMap.empty())
    //  return setErrorMessage("No Gradient found. Please compute one first!");

    typedef std::pair<int,Matrix3d> IntM3dP;

    // if there is no cell axis stored, display error message
    const VtxSymTensorAttr& cellAxisV = m->vertexAxis();
    // if(cellAxis.size() == 0 or m->cellAxisType() != "fibril"){
    //   setErrorMessage(QString("No fibril axis stored in active mesh!"));
    //   return false;
    // }

    // Check cell axis and prepare them for display. Populate cellAxisScaled in Mesh.
    VtxMatrix3fAttr &vtxAxisVis = m->vertexAxisVis();
    vtxAxisVis.clear();
    VtxVec3ColorbAttr &vtxAxisColor = m->vertexAxisColor();
    vtxAxisColor.clear();

    std::map<int, int> parentLabelMap;

    if(m->useParents()){

    }

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

    // QStringList parms;
    // DeformationGradient *dg = getProcessParms<DeformationGradient>(this, parms);


    IntVec3ColorbAttr &cellAxisColor = m->cellAxisColor();

    // Scale PDGs for visualization
    m->setAxisWidth(axisWidth);
    m->setAxisOffset(axisOffset);
    Colorb colorExpansion(colorPos);
    Colorb colorShrinkage(colorNeg);

    m->labelHeat().clear();

    typedef std::pair<vertex, Matrix3d> VtxM3dP;

    double minSig = 1E20, maxSig = -1E20;
    std::unordered_map<vertex, Point3d> vtxEigMap;

    forall(auto p, vertexAxis){

      //int cellLabel = p.first;
      vertex v = p.first;

      int minDisVtx = 1E20;
      if(minDis > 0){
        forall(vertex w, labelVertexVisMap[v->label]){
          double dis = norm(w->pos - v->pos);
          if(dis <  minDisVtx) minDisVtx = dis;
        }
        //if(minDisVtx < minDis) continue;
      }


      Point3d p1, p2, p3, ev;

      p1 = Point3d(p.second.ev1());
      p2 = Point3d(p.second.ev2());
      ev = Point3d(p.second.evals());

      Point3f zeroP(0,0,0);
      double evCustom = 0;

      vtxEigMap[v] = ev;


      if(minDisVtx > minDis){
        vtxAxisVis[v][0] = axisOption == "Max" or axisOption == "All" or axisOption == "Cell Plane" ? Point3f(p1.x(), p1.y(), p1.z()) * axisScale * (ev.x()-1) : zeroP;
        vtxAxisVis[v][1] = axisOption == "Mid" or axisOption == "All" or axisOption == "Cell Plane" ? Point3f(p2.x(), p2.y(), p2.z()) * axisScale * (ev.y()-1) : zeroP;
        vtxAxisVis[v][2] = axisOption == "Min" or axisOption == "All" ? Point3f(p3.x(), p3.y(), p3.z()) * axisScale * (ev.z()-1) : zeroP;

        vtxAxisColor[v][0] = ev.x() < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
        vtxAxisColor[v][1] = ev.y() < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
        vtxAxisColor[v][2] = ev.z() < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);

        labelVertexVisMap[v->label].insert(v);
      }

      // if(axisOption == "Custom"){

      //   Point3d customDir = Point3d(customDirAttr[v].ev1());
      //   evCustom = (t * customDir) * customDir;

      //   vtxAxisVis[v][0] = Point3f(customDir) * axisScale * (evCustom-1);
      //   vtxAxisColor[v][0] = evCustom < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
      // }


      //std::cout << "cell " << cellLabel << "/" << norm(t[0]) << "/" << norm(t[1]) << "/" << norm(t[2]) << "//" << t[0] << "/" << t[1] << "/" << t[2] << std::endl;

      //std::cout << "ev " << cellLabel << "/" << ev << "//" << p1 << "/" << p2 << "/" << p3 << std::endl;



      if(heatmapOption == "Stretch Max"){
        v->signal = ev.x();
      } else if(heatmapOption == "Stretch Mid"){
        v->signal = ev.y();
      } else if(heatmapOption == "Stretch Min"){
        v->signal = ev.z();
      } else if(heatmapOption == "Aniso Max/Mid"){
        v->signal = ev.y() > 1e-3 ? ev.x() / ev.y() : 0.0;
      } else if(heatmapOption == "Aniso Max/Min"){
        v->signal = ev.z() > 1e-3 ? ev.x() / ev.z() : 0.0;
      } else if(heatmapOption == "Aniso Mid/Min"){
        v->signal = ev.z() > 1e-3 ? ev.y() / ev.z() : 0.0;
      } else if(heatmapOption == "Stretch Product"){
        v->signal = ev.x() * ev.y() * ev.z();
      } else if(heatmapOption == "Stretch Custom"){
        v->signal = evCustom;
      }

      if(v->signal < minSig) minSig = v->signal;
      if(v->signal > maxSig) maxSig = v->signal;
    }

    // forall(auto p, vtxAxisVis){

    //   double colorV = (vtxEigMap[p.first].x() - minV)/(maxV - minV) * 255;

    //   int colorI = (int)colorV;

    //   QColor col1(colorI,colorI,colorI,0);

    //   vtxAxisColor[p.first][0] = Colorb(col1);
    //   vtxAxisColor[p.first][1] = Colorb(col1);
    //   vtxAxisColor[p.first][2] = Colorb(col1);

    // }

    // display it
    m->setAxisView("Vertex Axis");
    m->setCellAxisType("PDG");
    m->signalBounds() = Point2f(minSig, maxSig);
    m->setShowVertex("Signal");
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(DisplayDeformationGradientVertex3D);


   bool DisplayDeformationGradientOtherMesh::run(Mesh* m, QString heatmapOption, QString axisOption, QColor colorPos, QColor colorNeg, double axisWidth,
    double axisScale, double axisOffset, double minDis, bool projectOnSurface)
  {

    AttrMap<vertex, Matrix3d>& vertexTensorMap = m->attributes().attrMap<vertex, Matrix3d>("Deformation Gradient Vertex");
    vertexTensorMap.clear();

    vvGraph& S = m->graph();

    // get the current deformation
    AttrMap<int, ImplicitDeformation>& impDefAttr = m->attributes().attrMap<int, ImplicitDeformation>("Implicit Deformation");

    if(impDefAttr.empty())
      return setErrorMessage("No Deformation Attribute Map present. Please create or import one.");

    ImplicitDeformation& impDef = impDefAttr[0];

    forall(const vertex& v, S){

      if(projectOnSurface){
        vertexTensorMap[v] = impDef.calcProjectedGrowthTensorOnPos(v->pos, v->nrml);
      } else {
        vertexTensorMap[v] = impDef.calcGrowthTensor(v->pos);
      }
    }

    std::map<int, std::set<vertex> > labelVertexVisMap;

    //AttrMap<vertex, SymmetricTensor>& customDirAttr = m->attributes().attrMap<vertex, SymmetricTensor>("Custom Directions");

    if(vertexTensorMap.empty())
      return setErrorMessage("No Gradient found. Please compute one first!");

    // if there is no cell axis stored, display error message
    const VtxSymTensorAttr& cellAxisV = m->vertexAxis();
    // if(cellAxis.size() == 0 or m->cellAxisType() != "fibril"){
    //   setErrorMessage(QString("No fibril axis stored in active mesh!"));
    //   return false;
    // }

    // Check cell axis and prepare them for display. Populate cellAxisScaled in Mesh.
    VtxMatrix3fAttr &vtxAxisVis = m->vertexAxisVis();
    vtxAxisVis.clear();
    VtxVec3ColorbAttr &vtxAxisColor = m->vertexAxisColor();
    vtxAxisColor.clear();

    std::map<int, int> parentLabelMap;

    if(m->useParents()){

    }

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

    // QStringList parms;
    // DeformationGradient *dg = getProcessParms<DeformationGradient>(this, parms);


    IntVec3ColorbAttr &cellAxisColor = m->cellAxisColor();

    // Scale PDGs for visualization
    m->setAxisWidth(axisWidth);
    m->setAxisOffset(axisOffset);
    Colorb colorExpansion(colorPos);
    Colorb colorShrinkage(colorNeg);

    m->labelHeat().clear();

    typedef std::pair<vertex, Matrix3d> VtxM3dP;

    double minSig = 1E20, maxSig = -1E20;
    std::unordered_map<vertex, Point3d> vtxEigMap;

    forall(VtxM3dP p, vertexTensorMap){



      //int cellLabel = p.first;
      vertex v = p.first;
      Matrix3d t = p.second;

      int minDisVtx = 1E20;
      if(minDis > 0){
        forall(vertex w, labelVertexVisMap[v->label]){
          if(norm(w->pos - v->pos) <  minDisVtx) minDisVtx = norm(w->pos - v->pos);
        }
        //if(minDisVtx < minDis) continue;
      }



      Point3d p1, p2, p3, ev;
      calcEigenVecsValues(t, p1, p2, p3, ev);

      Point3f zeroP(0,0,0);
      double evCustom = 0;

      vtxEigMap[v] = ev;


      if(minDisVtx > minDis){
        vtxAxisVis[v][0] = axisOption == "Max" or axisOption == "All" or axisOption == "Cell Plane" ? Point3f(p1.x(), p1.y(), p1.z()) * axisScale * (ev.x()-1) : zeroP;
        vtxAxisVis[v][1] = axisOption == "Mid" or axisOption == "All" or axisOption == "Cell Plane" ? Point3f(p2.x(), p2.y(), p2.z()) * axisScale * (ev.y()-1) : zeroP;
        vtxAxisVis[v][2] = axisOption == "Min" or axisOption == "All" ? Point3f(p3.x(), p3.y(), p3.z()) * axisScale * (ev.z()-1) : zeroP;

        vtxAxisColor[v][0] = ev.x() < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
        vtxAxisColor[v][1] = ev.y() < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
        vtxAxisColor[v][2] = ev.z() < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);

        labelVertexVisMap[v->label].insert(v);
      }

      // if(axisOption == "Custom"){

      //   Point3d customDir = Point3d(customDirAttr[v].ev1());
      //   evCustom = (t * customDir) * customDir;

      //   vtxAxisVis[v][0] = Point3f(customDir) * axisScale * (evCustom-1);
      //   vtxAxisColor[v][0] = evCustom < 1 ? Colorb(colorShrinkage) : Colorb(colorExpansion);
      // }


      //std::cout << "cell " << cellLabel << "/" << norm(t[0]) << "/" << norm(t[1]) << "/" << norm(t[2]) << "//" << t[0] << "/" << t[1] << "/" << t[2] << std::endl;

      //std::cout << "ev " << cellLabel << "/" << ev << "//" << p1 << "/" << p2 << "/" << p3 << std::endl;



      if(heatmapOption == "Stretch Max"){
        v->signal = ev.x();
      } else if(heatmapOption == "Stretch Mid"){
        v->signal = ev.y();
      } else if(heatmapOption == "Stretch Min"){
        v->signal = ev.z();
      } else if(heatmapOption == "Aniso Max/Mid"){
        v->signal = ev.y() > 1e-3 ? ev.x() / ev.y() : 0.0;
      } else if(heatmapOption == "Aniso Max/Min"){
        v->signal = ev.z() > 1e-3 ? ev.x() / ev.z() : 0.0;
      } else if(heatmapOption == "Aniso Mid/Min"){
        v->signal = ev.z() > 1e-3 ? ev.y() / ev.z() : 0.0;
      } else if(heatmapOption == "Stretch Product 3D"){
        v->signal = ev.x() * ev.y() * ev.z();
      } else if(heatmapOption == "Stretch Product 2D"){
        v->signal = ev.x() * ev.y();
      } else if(heatmapOption == "Stretch Custom"){
        v->signal = evCustom;
      }

      if(v->signal < minSig) minSig = v->signal;
      if(v->signal > maxSig) maxSig = v->signal;
    }

    // forall(auto p, vtxAxisVis){

    //   double colorV = (vtxEigMap[p.first].x() - minV)/(maxV - minV) * 255;

    //   int colorI = (int)colorV;

    //   QColor col1(colorI,colorI,colorI,0);

    //   vtxAxisColor[p.first][0] = Colorb(col1);
    //   vtxAxisColor[p.first][1] = Colorb(col1);
    //   vtxAxisColor[p.first][2] = Colorb(col1);

    // }

    // display it
    m->setAxisView("Vertex Axis");
    m->setCellAxisType("PDG");

    //m->heatMapBounds() = m->calcHeatMapBounds();

    m->signalBounds() = Point2f(minSig, maxSig);

    m->setShowVertex("Signal");
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(DisplayDeformationGradientOtherMesh);


   bool Junctions3D::run(Mesh* m, Mesh* m2, bool junctionMode)
  {
    vvGraph& S = m->graph();
    vvGraph& S2 = m2->graph();

    // std::map<vertex, int> vtxWallMap;
    // std::map<vertex, int> vtxCellCounterMap;

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

    // vtxCellCounterMap = info.vtxCellCounterMap;
    // vtxWallMap = info.vtxNewLabelMap;

    // std::map<IntInt, int> allWalls = info.allWalls;
    // typedef std::pair<IntInt, int> IntIntPIntP;

    // // find all 3-way junctions in wall mesh

    // std::map<

    // forall(IntIntPIntP p, allWalls){

    // }

    // find all 3-way junctions on wall mesh
if(junctionMode){
    std::map<std::set<int>, std::set<vertex> > junctionMap; // map from set of labels to vertices
    typedef std::pair<std::set<int>, std::set<vertex> > IntSetVtxSetP;

    AttrMap<int, IntInt>& wallLabelCellLabelsMap = m->attributes().attrMap<int, IntInt>("Wall Label Cell Labels");
    //typedef std::pair<IntInt, int> IntIntPIntP;

    AttrMap<Point4i, Point3d>& junctionsPos = m->attributes().attrMap<Point4i, Point3d>("Junction Positions");

    forall(const vertex& v, S){
      if(v->label > -1) continue;

      std::set<int> neighborsV, parentsV;

      forall(const vertex& n, S.neighbors(v)){
        if(n->label == -1) continue;
        neighborsV.insert(n->label);
      }

      if(neighborsV.size() == 3){ // find 3-way junctions on the individual cells
        v->selected = true;
        forall(int l, neighborsV){ // add the pairs of cells that belong to the 3 cell walls (should be 4 in total, 0 for outside)
          parentsV.insert(wallLabelCellLabelsMap[l].first);
          parentsV.insert(wallLabelCellLabelsMap[l].second);
        }
        junctionMap[parentsV].insert(v);
      }

    }
    forall(IntSetVtxSetP p, junctionMap){
      if(p.second.empty()) continue;
      Point3d avgP(0,0,0);

      forall(int l, p.first){
        std::cout << "a " << l << std::endl;
      }

      forall(vertex v, p.second){
        avgP += v->pos;
      }
      avgP /= p.second.size();

      std::cout << "s " << avgP << "/" << p.second.size() << std::endl;

      vertex v;
      v->pos = avgP;
      S2.insert(v);
    }
} else {


}
    m->updateAll();
    m2->updateAll();

    return true;
  }
  REGISTER_PROCESS(Junctions3D);

  // void findNearestVtx(const vertex& refVtx, const std::set<vertex>& vtxSet, vertex& nearestVtx, vertex& runnerUpVtx, vertex& vertex3rd)
  // {

  //   double minDis1 = 1E18;
  //   double minDis2 = 1E19;
  //   double minDis3 = 1E20;

  //   vertex n1, n2, n3;

  //   forall(const vertex& v, vtxSet){
  //     double dis = norm(v->pos - refVtx->pos);
  //     if(dis < minDis1){
  //       n3 = n2;
  //       minDis3 = minDis2;

  //       n2 = n1;
  //       minDis2 = minDis1;

  //       n1 = v;
  //       minDis1 = dis;
  //     } else if(dis < minDis2){
  //       n3 = n2;
  //       minDis3 = minDis2;

  //       n2 = v;
  //       minDis2 = dis;
  //     } else if(dis < minDis3){
  //       n3 = v;
  //       minDis3 = dis;
  //     }
  //   }

  //   nearestVtx = n1;
  //   runnerUpVtx = n2;
  //   vertex3rd = n3;

  //   return;
  // }

   bool FlattenMesh::run(Mesh* m)
  {
    vvGraph& S = m->graph();

    // get junctions of all cells
    //std::map<int, std::set<vertex> > cellJunctions = getCellJunctions(S);

    std::map<int, std::set<vertex> > cellJunctions = getCellContour(S);

    std::cout << "s " << cellJunctions.size() << std::endl;

    std::map<int, Point3d> cellNormals;
    std::map<int, Point3d> cellPoints;

    std::unordered_map<int, std::vector<Point3d> > cellJunctionPosMap;
    std::unordered_map<int, std::vector<vertex> > cellJunctionsOrdered;

    std::unordered_map<int, std::set<vertex> > cellVerticesMap;

    forall(const vertex& v, S){
      if(v->label < 0) continue;

      cellVerticesMap[v->label].insert(v);
    }

    // do pca on junctions to find cell plane
    forall(auto p, cellJunctions){

      std::vector<Point3d> pointsPolygon;

      forall(vertex v, p.second){
        pointsPolygon.push_back(v->pos);
        cellJunctionsOrdered[p.first].push_back(v);
      }
      cellJunctionPosMap[p.first] = pointsPolygon;

      Point3d planePos, planeNrml;

      findPolygonPlane(cellJunctionPosMap[p.first], planePos, planeNrml);

      cellNormals[p.first] = planeNrml;
      cellPoints[p.first] = planePos;
    }

    // project junctions onto plane, rotate into XY
    forall(auto p, cellJunctions){
      if(p.second.size() < 3) continue;
      Matrix3d m;
      rotatePointsIntoXY(cellNormals[p.first], cellJunctionPosMap[p.first], m);

      std::vector<Point2f> pList;
      std::vector<Point2f> ptList;
      std::vector<Point3i> triList;

      forall(Point3d pos, cellJunctionPosMap[p.first]){
        pList.push_back(Point2f(pos.x(), pos.y()));
      }

      // delaunay triangulation 2D on the projected junctions
      triangulateDelaunay2D(pList, ptList, triList);

      // create list of triangles and original vtxs
      // for(int i = 0; i<pList.size(); i++){
      //   std::cout << "pl " << pList[i] << "/" << ptList[i] << std::endl;
      // }

      std::map<vertex, int> vertexIndexMap;

      // for(int i = 0; i<cellJunctionsOrdered[p.first].size(); i++){
      //   vertexIndexMap[cellJunctionsOrdered[p.first][i]] = i;
      // }

      // foreach vertex: find the nearest 2 junctions, find their triangle & normal, project points onto it
      forall(const vertex& v, cellVerticesMap[p.first]){
        //if(!v->selected) continue;
        //vertex n1, n2, n3;
        //findNearestVtx(v, cellJunctions[p.first], n1, n2, n3);

        Point3i foundTri;


        forall(Point3i tri, triList){

          std::vector<Point3d> pos;
          pos.push_back(v->pos);
          Matrix3d m;
          rotatePointsIntoXY(cellNormals[p.first], pos, m);

          std::vector<Point3d> triP;

          triP.push_back(Point3d(ptList[tri.x()].x(), ptList[tri.x()].y(), 0));
          triP.push_back(Point3d(ptList[tri.y()].x(), ptList[tri.y()].y(), 0));
          triP.push_back(Point3d(ptList[tri.z()].x(), ptList[tri.z()].y(), 0));

          // if(v->selected and pointInPolygon(pos[0], triP) != pointInPolygonCGAL(pos[0], triP)){
          //   pointInPolygon(pos[0], triP,true);
          //   std::cout << "p " << pointInPolygon(pos[0], triP) << "/" << pointInPolygonCGAL(pos[0], triP) << std::endl;
          //   std::cout << "t " << pos[0] << "/" << triP[0] << "/" << triP[1] << "/" << triP[2] << std::endl;
          // }

          // if(pointInPolygonNew(pos[0], triP) and v->selected){
          //   std::cout << "found " << pos[0] << "/" << triP[0] << "/" << triP[1] << "/" << triP[2] << std::endl;
          // }

          if(pointInPolygon(pos[0], triP)){
            foundTri = tri;

            // get tri normal
            Point3d triNormal = (cellJunctionsOrdered[p.first][foundTri.z()]->pos - cellJunctionsOrdered[p.first][foundTri.x()]->pos) %
                                (cellJunctionsOrdered[p.first][foundTri.y()]->pos - cellJunctionsOrdered[p.first][foundTri.x()]->pos);

            triNormal /= norm(triNormal);

            //cellJunctionsOrdered[p.first][foundTri.z()]->selected = true;
            //cellJunctionsOrdered[p.first][foundTri.y()]->selected = true;
            //cellJunctionsOrdered[p.first][foundTri.x()]->selected = true;

            //Point3d triNormal = (n3->pos - n1->pos) % (n2->pos - n1->pos);

            //Point3d triNormal = n1->nrml;

            // project the vertex onto the tri
            double sN;
            Point3d intersectN;
            planeLineIntersect(cellJunctionsOrdered[p.first][foundTri.x()]->pos, triNormal, v->pos, v->pos-triNormal, sN, intersectN);

            //std::cout << "p " << v->pos << "/" << triNormal << "/" << intersectN << std::endl;

            v->pos = intersectN;

            break;
          }

        }

      }
    }

    // foreach vertex: project onto

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

    // forall(auto p, cellJunctions){
    //   forall(const vertex& v, p.second){
    //     v->selected = true;
    //   }
    // }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(FlattenMesh);


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

   bool LoadPointsCSV::run(Mesh* m, Mesh* m2, QString filename)
  {
    AttrMap<int, Point3d>& pc1 = m->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& pc2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    if(pc1.size() != pc2.size()){
      return setErrorMessage("the Manual Landmarks attribute maps have a different length! Clear them first before running this process.");
    }

      // open file
      QFile file(filename);
      if(!file.open(QIODevice::ReadOnly))
      {
        setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
        return false;
      }
      QTextStream ss(&file);
      QString line = ss.readLine();
      QStringList fields = line.split(",");

      bool ok;
      int lineCounter = 0;

      while(ss.status() == QTextStream::Ok)
      {
        fields = ss.readLine().split(",");
        if(fields[0].isEmpty()) break;

        std::cout << "l " << lineCounter << std::endl;

        double dataP1X = fields[1].toDouble(&ok);
        double dataP1Y = fields[2].toDouble(&ok);
        double dataP1Z = fields[3].toDouble(&ok);
        double dataP2X = fields[4].toDouble(&ok);
        double dataP2Y = fields[5].toDouble(&ok);
        double dataP2Z = fields[6].toDouble(&ok);

        Point3d p1(dataP1X, dataP1Y, dataP1Z);
        Point3d p2(dataP2X, dataP2Y, dataP2Z);

        std::cout << "d " << p1 << "/" << p2 << std::endl;

        int l = pc1.size();

        pc1[l] = p1;
        pc2[l] = p2;

        lineCounter++;
      }


    return true;
  }
  REGISTER_PROCESS(LoadPointsCSV);


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

   bool SavePointsCSV::run(Mesh* m, Mesh* m2, QString filename, bool saveCenters, bool saveJunctions, bool saveCentroids, bool saveManual)
  {
    AttrMap<int, Point3d>& pc1 = m->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& pc2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    vvGraph& S1 = m->graph();
    vvGraph& S2 = m2->graph();

    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream outS(&file);

    outS << "Point, mesh1X, mesh1Y, mesh1Z, mesh2X, mesh2Y, mesh2Z" << endl;


    if(saveCenters){
      m->updateCentersNormals();
      m2->updateCentersNormals();
      IntPoint3fAttr centers = m->labelCenter();
      IntPoint3fAttr centersP2 = m2->parentCenter();
      //std::cout << "saveC" << std::endl;
      forall(auto p, centers){
        if(centersP2.find(p.first) != centersP2.end()){
          outS << "C_" << p.first << "," << centers[p.first].x() << "," << centers[p.first].y() << ","<< centers[p.first].z() << ","<< centersP2[p.first].x() << ","<< centersP2[p.first].y() << ","<< centersP2[p.first].z() << endl;
      //std::cout << "saveC" << p.first << std::endl;
        }
      }
    }

    if(saveJunctions){ // this code is copied from the CreateDeformation2D process

      vvGraph& S = m->graph();
      vvGraph& S2 = m2->graph();

      vvGraph T1,T2;

      mgxmToMgxc(S, T1, 0);
      mgxmToMgxc(S2, T2, 0);

      JunctionInformation jInfo;

      createJunctionAssociation(m, m2, T1, T2, jInfo);
      std::cout << "found corresponding junctions" << std::endl;

      double distanceThreshold = -1;

      std::map<int, Point3d> labelPosMap;

      IntPoint3fAttr centers = m->labelCenter();
      IntPoint3fAttr centers2 = m2->labelCenter();
      IntPoint3fAttr centersP2 = m2->parentCenter();

      std::map<int, Point3d> cellCentersM1, cellCentersM2;

      forall(const vertex& v, S2){
        if(v->label < 1) continue;
        cellCentersM2[v->label] = Point3d(centers2[v->label]);
        labelPosMap[m2->parents()[v->label]] = Point3d(centersP2[m2->parents()[v->label]]);
      }

      forall(const vertex& v, T1){
        if(v->label < 1) continue;
        cellCentersM1[v->label] = Point3d(centers[v->label]);
        v->selected = true;
      }

      jInfo.labelDaughtersCenterMap = labelPosMap;
      jInfo.labelParentCenterMap = cellCentersM1;


    std::map<vertex, Point3d> replaceVtxsMap1, replaceVtxsMap2;
    std::set<Point3d> addedAvgPoints;


      forall(CellJunctionsInOrderP p, jInfo.cellJunctionInOrder1){

        std::vector<vertex> borderInOrder = p.second;

        for(uint i = 0; i < borderInOrder.size()-1; i++){
          vertex v1,v2;
          v1 = borderInOrder[i];
          v2 = borderInOrder[borderInOrder.size() - 1];
          if(i > 0) v2 = borderInOrder[i + 1];

          vertex v21, v22;
          if(jInfo.J1J2.find(v1) == jInfo.J1J2.end() or jInfo.J1J2.find(v2) == jInfo.J1J2.end()) continue;

          v21 = jInfo.J1J2[v1];
          v22 = jInfo.J1J2[v2];

          if(norm(v1->pos - v2->pos) < distanceThreshold or norm(v21->pos - v22->pos) < distanceThreshold){

            // only 1 merging allowed currently
            if(replaceVtxsMap1.find(v1) != replaceVtxsMap1.end() or
               replaceVtxsMap1.find(v2) != replaceVtxsMap1.end() or
               replaceVtxsMap2.find(v21) != replaceVtxsMap2.end() or
               replaceVtxsMap2.find(v22) != replaceVtxsMap2.end())
              continue;

            Point3d avgPos = (v1->pos + v2->pos)/2.;
            replaceVtxsMap1[v1] = avgPos;
            replaceVtxsMap1[v2] = avgPos;

            Point3d avgPos2 = (v21->pos + v22->pos)/2.;
            replaceVtxsMap2[v21] = avgPos2;
            replaceVtxsMap2[v22] = avgPos2;

          }

        }

      }

      forall(VtxVtxP p, jInfo.J1J2){
        if(replaceVtxsMap1.find(p.first) != replaceVtxsMap1.end() and replaceVtxsMap2.find(p.second) != replaceVtxsMap2.end()){
          if(addedAvgPoints.find(replaceVtxsMap1[p.first]) != addedAvgPoints.end()) continue;

          addedAvgPoints.insert(replaceVtxsMap1[p.first]);

          outS << "J" << "," << replaceVtxsMap1[p.first].x() << "," << replaceVtxsMap1[p.first].y() << ","
          << replaceVtxsMap1[p.first].z() << ","<< replaceVtxsMap2[p.second].x() << ","<< replaceVtxsMap2[p.second].y() 
          << ","<< replaceVtxsMap2[p.second].z() << endl;
      

          //std::cout << "data/values (J/R) " << p.first->pos << "/" << p.second->pos << std::endl;
          p.first->selected = true;
          p.second->selected = true;
        } else {
          outS << "J" << "," << p.first->pos.x() << "," << p.first->pos.y() << ","<< p.first->pos.z() << ","<< p.second->pos.x() << ","<< p.second->pos.y() << ","<< p.second->pos.z() << endl;
      

          //std::cout << "data/values (J) " << p.first->pos << "/" << p.second->pos << std::endl;
        }
      }
   
    }

    if(saveCentroids){
      AttrMap<int, Point3d>& centroidsM1 = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
      AttrMap<int, Point3d>& centroidsM2 = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

      forall(auto p, centroidsM1){
        if(centroidsM2.find(p.first) != centroidsM2.end()){
          outS << "C3_" << p.first << "," << centroidsM1[p.first].x() << "," << centroidsM1[p.first].y() << ","<< centroidsM1[p.first].z() << ","<< centroidsM2[p.first].x() << ","<< centroidsM2[p.first].y() << ","<< centroidsM2[p.first].z() << endl;
      
        }
      }
    }

    if(saveManual){
      int l = pc1.size();
      for(int i = 1; i<l+1 ; i++){
        outS << "M_" << i << "," << pc1[i].x() << "," << pc1[i].y() << ","<< pc1[i].z() << ","<< pc2[i].x() << ","<< pc2[i].y() << ","<< pc2[i].z() << endl;
      }
    }

    return true;
  }
  REGISTER_PROCESS(SavePointsCSV);


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

  bool RootJunctionData::run(const QString& filename)
  {
	Mesh* mesh1 = currentMesh(), *mesh2 =  otherMesh();
   vvGraph& S1 = mesh1->graph();
   vvGraph& S2 = mesh2->graph();
   vvGraph T1, T2;
   mgxmToMgxc(S1, T1, 0);
   mgxmToMgxc(S2, T2, 0);
  JunctionInformation jInfo;
  createJunctionAssociation(mesh1, mesh2, T1,T2, jInfo);


    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg("junctionData"));
      return false;
    }


    QTextStream out(&file);
	forall(const VertLabelsP& j1, jInfo.junctionLabels1){
     forall(const VertLabelsP& j2, jInfo.junctionLabels2){
       if(j1.second == j2.second){
	     out << j1.first->pos.x() << "," << j2.first->pos.x()<< endl;

       }
     }
   }

  }
  REGISTER_PROCESS(RootJunctionData);

  bool AddManualLandmarks::run(Mesh* m1, Mesh* m2)
  {
    AttrMap<int, Point3d>& landmarks_m1 = m1->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& landmarks_m2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    std::cout << "current number of manual landmarks: " << landmarks_m1.size() << std::endl;

    if(landmarks_m1.size() != landmarks_m2.size()){
      return setErrorMessage("Landmark Attr Maps have different sizes! Correct or clear them first!");
    }

    vvGraph& S1 = m1->graph();
    vvGraph& S2 = m2->graph();

    std::vector<vertex> sel1 = m1->selectedVertices();
    std::vector<vertex> sel2 = m2->selectedVertices();

    if(sel1.empty() or sel2.empty()){
      return setErrorMessage("There must be at least 1 vertex selected in BOTH meshes.");
    }

    Point3d p1(0,0,0);
    for(vertex v : sel1){
      p1 += v->pos;
    }
    p1/=sel1.size();

    Point3d p2(0,0,0);
    for(vertex v : sel2){
      p2 += v->pos;
    }
    p2/=sel2.size();

    int newIdx = landmarks_m1.size() + 1;

    landmarks_m1[newIdx] = p1;
    landmarks_m2[newIdx] = p2;

    std::cout << "List of all manual landmarks: " << std::endl;
    for(auto p : landmarks_m1){
      std::cout << p.first << "/" << p.second << "/" << landmarks_m2[p.first] << std::endl;
    }

    return true;
  }
  REGISTER_PROCESS(AddManualLandmarks);

  bool ClearManualLandmarks::run(Mesh* m1, Mesh* m2)
  {

    AttrMap<int, Point3d>& landmarks_m1 = m1->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& landmarks_m2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    landmarks_m1.clear();
    landmarks_m2.clear();
    return true;
  }
  REGISTER_PROCESS(ClearManualLandmarks);

 // return true if pointToTest is inside a defined cone
 bool pointInOpenCone(Point3d axisVec, Point3d& coneTop, double coneAngle, Point3d& pointToTest)
 {
    Point3d apexToX = coneTop - pointToTest;

    double dp = apexToX * axisVec;
    double denom = norm(axisVec) * norm(apexToX);

    //double dir = (pointToTest - coneBase) * axisVec;

    if(dp/denom > cos(coneAngle)) return true;
    else return false;
 }



  bool SurfaceExtensionLandmarks::run(const Stack *s1, const Stack *s2, Mesh* m1, Mesh* m2)
  {

    AttrMap<int, Point3d>& landmarks_m1 = m1->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& landmarks_m2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    vvGraph& S1 = m1->graph();
    vvGraph& S2 = m2->graph();

    std::vector<vertex> selV = m1->selectedVertices();

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

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

    std::map<Point3d, Point3d> landmarks_map;

    forall(vertex v, selV){
      double minDis = HUGE_VAL;
      vertex minV;
      Point3d vCorrected = multMatrix4Point3(mGLTot,v->pos);

      std::map<double, Point3d> disPointMap;

      forall(vertex n, S2){
        double dis = norm(vCorrected - n->pos);
        while(disPointMap.find(dis) != disPointMap.end()){
          dis = dis + 1E-6;
        }
        disPointMap[dis] = n->pos;
      }

      Point3d nearest(0,0,0);

      forall(auto p, disPointMap){
        if(!pointInOpenCone(v->nrml,vCorrected,parm("Cone Angle").toDouble(),p.second) and !pointInOpenCone(-v->nrml,vCorrected,parm("Cone Angle").toDouble(),p.second)) continue;
        nearest = p.second;
        break;
      }


      //Point3d pM = nearestPointOnMesh(vCorrected, S2);
      if(norm(nearest) > 1E-3){
        landmarks_map[nearest] = v->pos;
      }
      
    }

    int c = landmarks_m1.size();

    forall(auto p, landmarks_map){
      landmarks_m1[c] = p.second;
      landmarks_m2[c] = p.first;
      c++;
    }

    return true;
  }
  REGISTER_PROCESS(SurfaceExtensionLandmarks);


  bool DrawManualLandmarks::run(const Stack *s1, const Stack *s2,Mesh* m1, Mesh* m2)
  {

    AttrMap<int, Point3d>& landmarks_m1 = m1->attributes().attrMap<int, Point3d>("Manual Landmarks");
    AttrMap<int, Point3d>& landmarks_m2 = m2->attributes().attrMap<int, Point3d>("Manual Landmarks");

    vvGraph& S1 = m1->graph();
    vvGraph& S2 = m2->graph();

    forall(const vertex& v, S1){
      v->selected = false;
    }
    forall(const vertex& v, S2){
      v->selected = false;
    }

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

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


    forall(auto p, landmarks_m1){
      vertex v;
      v->pos = p.second;
      v->selected = true;

      Point3d vCorrected = multMatrix4Point3(inverse(mGLTot),landmarks_m2[p.first]);

      vertex n;
      n->pos = vCorrected;
      n->selected = true;
      S1.insert(v);
      S1.insert(n);
      S1.insertEdge(v,n);
      S1.insertEdge(n,v);
    }

    // forall(auto p, landmarks_m2){
    //   vertex v;
    //   v->pos = p.second;
    //   v->selected = true;
    //   S2.insert(v);
    // }

    m1->updateAll();
    m2->updateAll();

    return true;
  }
  REGISTER_PROCESS(DrawManualLandmarks);





}
