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

#include <MeshProcessCellAxis.hpp>
#include <Dir.hpp>
#include <Progress.hpp>
#include <SetVector.hpp>
#include <Triangulate.hpp>
#include <GraphUtils.hpp>
#include <MeshProcessMeasures.hpp> // MeasureArea
#include <MeshProcessHeatMap.hpp> // geoMeanVectorValue

namespace mgx
{
  using std::swap;

  namespace {

    typedef std::unordered_map<int, std::pair<VtxVec, VtxVec> > IntVtxVecPairMap;
    typedef std::pair<int, std::pair<VtxVec, VtxVec> > IntVtxVecPairPair;
    typedef std::unordered_map<int, set_vector<int> > IntLabelsMap;
    typedef std::unordered_map<vertex, set_vector<int> > VertLabelsMap;
    typedef std::unordered_map<int, vertex> VMap;
    typedef std::pair<int, vertex> VMPair;

    struct CorrespondenceMeshes {
      void clear()
      {
        J1J2.clear();
        J1J2forC1.clear();
        J1J2forC2.clear();
        centerNeighborsMapT1.clear();
        centerNeighborsMapT2.clear();
      }

      // To store the correspondence between junctions:
      // for the entire meshes:
      VtxVtxMap J1J2;
      // for each cell in mesh1:
      IntVtxVecPairMap J1J2forC1;
      IntVtxVecPairMap J1J2forC2;
      // Needed to check the connections between neighbors:
      IntLabelsMap centerNeighborsMapT1;
      IntLabelsMap centerNeighborsMapT2;
    };

    CorrespondenceMeshes cm;
  }

  SymmetricTensor customDirectionTensor(SymmetricTensor origTensor, Point3f customX, Point3f customY, bool shearB)
  {

    SymmetricTensor customTensor;

    Point3f customTensorEvals = Point3f(0,0,0);
    float cos, sin;

    Point3d xDir = Point3d(normalized(customX));
    Point3d yDir = Point3d(normalized(customY));

    cos = xDir * Point3d(normalized(origTensor.ev1()));
    sin = yDir * Point3d(normalized(origTensor.ev1()));

    float strainMax = origTensor.evals()[0] - 1;
    float strainMin = origTensor.evals()[1] - 1;

    float strainCustomX = cos * cos * strainMax + sin * sin * strainMin;
    float strainCustomY = sin * sin * strainMax + cos * cos * strainMin;

    float shear = (strainMin - strainMax) * sin * cos;

    if(shearB){
      strainCustomX += 2 * shear * sin * cos;
      strainCustomY += 2 * shear * sin * cos;
    }

    customTensorEvals = Point3f(strainCustomX+1, strainCustomY+1, shear);

    customTensor.ev1() = Point3f(xDir);
    customTensor.ev2() = Point3f(yDir);
    customTensor.evals() = customTensorEvals;

    return customTensor;
  }

  // TODO NEEDS TO BE TESTED
  double lengthCustomDirTensor3D(SymmetricTensor origTensor, Point3f customDir)
  {
    double result = -1;

    customDir /= norm(customDir);

    // assemble tensor matrix
    Point3f ev1 = origTensor.ev1();
    ev1 *= origTensor.evals().x();
    Point3f ev2 = origTensor.ev1();
    ev2 *= origTensor.evals().y();
    Point3f ev3 = ev1 ^ ev2;
    ev3 *= origTensor.evals().z();

    Matrix3d tensor = origTensor.toMatrix();

    // tensor[0][0] = ev1.x();
    // tensor[0][1] = ev1.y();
    // tensor[0][2] = ev1.z();
    //
    // tensor[1][0] = ev2.x();
    // tensor[1][1] = ev2.y();
    // tensor[1][2] = ev2.z();
    //
    // tensor[2][0] = ev3.x();
    // tensor[2][1] = ev3.y();
    // tensor[2][2] = ev3.z();

    Point3f TV;
    for(int i = 0; i < 3; i++){
      for(int j = 0; j < 3; j++){
        TV[i] += tensor[i][j] * customDir[j];
      }
    }

    //Point3f TV = tensor * customDir;
    result = customDir * TV;

    return result;
  }


  // take the correspondence between junctions for the whole mesh, compute the one for each cell
  bool J1J2forCell(Mesh *mesh, const vvGraph& T1, const IntIntMap& labelsT1, const IntLabelsMap& parentDaughtersMap,
                   set_vector<vertex>& problemJunctionT1)
  {
    // Get label count for each vertex
    VtxIntMap labCount;
    getLabelCount(T1, labCount);

    forall(const vertex &c1, T1) {
      std::vector<vertex> J1forC1;
      std::vector<vertex> J2forC2;

      // c1 has to be a center, connected to  at least 3 vertices
      {
        IntIntMap::const_iterator it = labelsT1.find(c1->label);
        if(c1->type != 'c' or T1.valence(c1) < 3 or it == labelsT1.end() or it->second == 0)
          continue;
      }
      IntLabelsMap::const_iterator it = parentDaughtersMap.find(c1->label);
      if(it == parentDaughtersMap.end())
        continue;
      const set_vector<int>& lDaughters = it->second;
      if(lDaughters.count(0) > 0)
        continue;

      // pick a vertex in neighborhood of c1 to start searching for correspondence.
      vertex v1 = T1.anyIn(c1);
      vertex start1 = v1;
      do
        v1 = T1.nextTo(c1, v1);
      while(labCount[v1] < 3 and v1 != start1);
      // if couldn't find any vertex connected to 3 neighbors, go to next cell in T1
      if(labCount[v1] < 3) {
        Information::out << "could not find a 3 way junction in cell: " << c1->label << " in T1" << endl;
        continue;
        //break;
      }

      vertex j1 = v1;     // starting vertex for the correspondance
      // check if there is a corresponding vertex in S2
      do {
        if(j1->type == 'j') {
          VtxVtxMap::iterator it = cm.J1J2.find(j1);
          if(it != cm.J1J2.end()) {
            vertex j2 = it->second;
            J1forC1.push_back(j1);
            J2forC2.push_back(j2);
          } else {
            problemJunctionT1.insert(j1);
          }
        }
        j1 = T1.nextTo(c1, j1);
      } while(j1 != v1);

      // if the correspondence was OK for the whole cell, store it into  J1J2forC1
      {
        IntLabelsMap::iterator it = cm.centerNeighborsMapT1.find(c1->label);
        if(it != cm.centerNeighborsMapT1.end() and it->second.size() != 0 and it->second.size() == J2forC2.size()) {
          cm.J1J2forC1[c1->label].first = J1forC1;
          cm.J1J2forC1[c1->label].second = J2forC2;
        } else
          Information::out << "Label " << c1->label << " has only " << J2forC2.size()
                           << " corresponding vertices in T1 but " << /*J1forC1.size()*/ it->second.size() << " in T0" << endl;
      }
    }
    return true;
  }

  // Check the correspondence between the junctions
  bool CorrespondenceJunctions::run(Mesh* mesh1, Mesh* mesh2, bool ShowVVCorrespondence, bool convert, bool newMethod)
  {
    // clean up the previous correspondence
    cm.clear();
    mesh1->setShowLabel("Label");
    mesh2->setShowLabel("Label");

    // Define temp vvgraphs
    vvGraph T1s, T2s;
    // Use the same graph if convert set
    vvGraph &T1 = convert ? mesh1->graph() : T1s;
    vvGraph &T2 = convert ? mesh2->graph() : T2s;

    // Create new graph containing only junctions
    if(convert)
      mesh1->updateAll();
    else
      mesh1->updateSelection(); // In case it fails
    if(!mgxmToMgxc(mesh1->graph(), T1, 0))
      return false;

    if(convert)
      mesh2->updateAll();
    else
      mesh2->updateSelection(); // In case it fails
    if(!mgxmToMgxc(mesh2->graph(), T2, 0))
      return false;
    if(!newMethod)
      return(checkNhbd(mesh1, mesh2, T1, T2, ShowVVCorrespondence));

    return(checkNhbdNew(mesh1, mesh2, T1, T2, ShowVVCorrespondence));
  }
  REGISTER_PROCESS(CorrespondenceJunctions);

  bool CorrespondenceJunctions::checkNhbd(Mesh* mesh1, Mesh* mesh2, vvGraph &T1, vvGraph &T2, bool ShowVVCorrespondence){
    // Create a map of labels in both meshes, to check parent labeling or create it
    // we will later add 0 values in the label map for cells which don't have daughters
    // (in T1) or parent (in T2)
    IntIntMap labelsT1, labelsT2;
    forall(const vertex &c, T1)
        if(c->type == 'c')
        labelsT1[c->label] = c->label;
    forall(const vertex &c, T2)
        if(c->type == 'c')
        labelsT2[c->label] = c->label;

    // daughterParentMap will be either labels or parents depending on what is selected
    // If cell does not have parent in T1, key = label in T2, value = 0
    IntIntAttr &parents = mesh2->parents();
    IntIntAttr daughterParentMap;
    if(mesh2->useParents()) {
      Information::out << "Found parent labels, using parent labeling to check correspondence between junctions." << endl;
      // Check the parents for cells which are not in the mesh1 or mesh2
      for(IntIntAttr::const_iterator it = parents.begin(); it != parents.end(); it++) {
        int lParent = it->second;
        int lDaughter = it->first;
        if(labelsT1.count(lParent) == 0 or labelsT2.count(lDaughter) == 0)
          // RSS If we don't find it just skip it.
          continue;
        //return setErrorMessage(QString("Problems with parent label: %1, daughter label: %2."
        //     " First run Correct Parents, then Save and Load parents.").arg(lParent).arg(lDaughter));
        daughterParentMap[lDaughter] = lParent;
      }
    } else {
      // If there is no parent map stored in mesh, the meshes should have the same labels
      // create a parent map based on the same cells in both meshes
      Information::out << "No parent labels found, using labels to check correspondence between junctions." << endl;
      for(IntIntMap::const_iterator l = labelsT1.begin(); l != labelsT1.end(); l++) {
        int label = l->first;
        if(labelsT2.count(label) != 0)
          daughterParentMap[label] = label;
      }
    }

    // Look at the daughters (in T2) for each cell in T1
    // parentDaughtersMap: key = parent in T1, value = labels of daughters in T2
    // If there is no parent, key = 0, values = labels of orphans (in T2)
    IntLabelsMap parentDaughtersMap;
    for(IntIntMap::const_iterator l = labelsT2.begin(); l != labelsT2.end(); l++) {
      int labelDaughter = l->first;
      int labelParent;
      set_vector<int> LDaughters;
      // add to daughterParentMap 0 values for the cells that do not have parent in T1.
      // map the label of cells in T2 without parent with 0
      {
        IntIntAttr::iterator it = daughterParentMap.find(labelDaughter);
        if(it == daughterParentMap.end() or it->second == 0 or labelDaughter == 0) {
          labelsT2[labelDaughter] = 0;
          daughterParentMap[labelDaughter] = 0;
          labelParent = 0;
        } else {
          labelParent = it->second;
        }
      }
      IntLabelsMap::iterator it = parentDaughtersMap.find(labelParent);
      if(it != parentDaughtersMap.end())
        LDaughters = it->second;
      LDaughters.insert(labelDaughter);
      parentDaughtersMap[labelParent] = LDaughters;
    }
    // add 0 values in parentDaughtersMap for cells in T1 that do not have daughters
    // map the label of cells in T1 without daughters with 0
    set_vector<int> noDaughter;
    noDaughter.insert(0);
    for(IntIntMap::const_iterator l = labelsT1.begin(); l != labelsT1.end(); l++) {
      int labelParent = l->first;
      if(parentDaughtersMap.count(labelParent) == 0) {
        parentDaughtersMap[labelParent] = noDaughter;
        labelsT1[labelParent] = 0;
      }
    }

    // Fill in a map of the labels associated with each junction in T1 and T2.
    // The list of labels is sorted (it is a set).
    // If there are only 2 neighbors, then the edge is at the mesh border and
    // we say it is in contact with the outside (label 0);
    // Labels of cells of T1 that do not have daughters are stored as 0.
    // This will be used to identify junctions at each time points.
    VertLabelsMap junctionLabelsT1;
    forall(const vertex &j, T1) {
      set_vector<int> labels;
      if(j->type == 'j') {
        // list the labels around junction
        forall(const vertex &n, T1.neighbors(j))
            if(n->label > 0) {
          IntIntMap::iterator it = labelsT1.find(n->label);
          if(it != labelsT1.end())
            labels.insert(it->second);
        }
        if(labels.size() == 2)
          labels.insert(0);
        junctionLabelsT1[j] = labels;
      }
    }
    VertLabelsMap junctionLabelsT2;
    forall(const vertex &j, T2) {
      set_vector<int> labels;
      if(j->type == 'j') {
        // list the labels around junction
        forall(const vertex &n, T2.neighbors(j))
            if(n->label > 0) {
          IntIntMap::iterator it = labelsT2.find(n->label);
          if(it != labelsT2.end())
            labels.insert(it->second);
        }
        if(labels.size() == 2)
          labels.insert(0);
        junctionLabelsT2[j] = labels;
      }
    }

    // Look for the correspondence between junctions in T1 and in T2
    forall(const vertex &j2, T2) {
      if(junctionLabelsT2.count(j2) > 0) {
        set_vector<int> labelsDaughters;
        VertLabelsMap::iterator it = junctionLabelsT2.find(j2);
        if(it != junctionLabelsT2.end())
          labelsDaughters = it->second;
        set_vector<int> labelsParents;
        forall(int lDaughter, labelsDaughters) {
          IntIntAttr::iterator it = daughterParentMap.find(lDaughter);
          int lParent = 0;
          if(it != daughterParentMap.end())
            lParent = it->second;
          labelsParents.insert(lParent);
        }
        // here compare the labels of parents of junction in T2 with junction in T1.
        // we considere only 3 ways junctions.
        if(labelsParents.size() == 3)
          for(VertLabelsMap::const_iterator it = junctionLabelsT1.begin();
              it != junctionLabelsT1.end(); it++) {
            const set_vector<int>& labelsJunction = it->second;
            if(labelsJunction == labelsParents)
              cm.J1J2[it->first] = j2;
          }
      }
    }

    // Show the correspondence between vertices.
    IntMatrix3fMap& J1J2vectors = mesh1->vvCorrespondence();
    J1J2vectors.clear();
    mesh1->showVVCorrespondence() = ShowVVCorrespondence;
    int i = 0;
    for(VtxVtxMap::const_iterator it = cm.J1J2.begin(); it != cm.J1J2.end(); it++, i++) {
      vertex j1 = it->first;
      vertex j2 = it->second;
      J1J2vectors[i][0] = Point3f(j1->pos);
      J1J2vectors[i][1] = Point3f(j2->pos);
    }

    // Create a map of neighbors (labels) for each cell in T1 and T2.
    // in T1, labels of neighbors without daughters are stored as 0
    forall(const vertex &c, T1) {
      set_vector<int> labels;
      if(c->type == 'c' and c->label > 0) {
        // list the labels at junctions around the center
        forall(const vertex &j, T1.neighbors(c)) {
          if(junctionLabelsT1.count(j) != 0) {
            VertLabelsMap::iterator it = junctionLabelsT1.find(j);
            if(it != junctionLabelsT1.end())
              labels.insert(it->second.begin(), it->second.end());
          }
        }
        labels.erase(c->label);
        cm.centerNeighborsMapT1[c->label] = labels;
      }
    }
    // in T2, labels of neighbors without parents are stored as 0
    forall(const vertex &c, T2) {
      set_vector<int> labels;
      if(c->type == 'c' and c->label > 0) {
        // list the labels at junctions around the center
        forall(const vertex &j, T2.neighbors(c)) {
          if(junctionLabelsT2.count(j) != 0) {
            VertLabelsMap::iterator it = junctionLabelsT2.find(j);
            if(it != junctionLabelsT2.end())
              labels.insert(it->second.begin(), it->second.end());
          }
        }
        labels.erase(c->label);
        cm.centerNeighborsMapT2[c->label] = labels;
      }
    }

    // Select the vertices in T1 which cause trouble
    set_vector<vertex> problemJunctionT1;
    Information::out << "Check problems with correspondence between junctions:" << endl;
    J1J2forCell(mesh1, T1, labelsT1, parentDaughtersMap, problemJunctionT1);
    Information::out << problemJunctionT1.size() << " problem junctions found in "
                     << labelsT1.size() << " cells" << endl;
    forall(const vertex &v, T1) {
      v->selected = false;
      if(problemJunctionT1.count(v) != 0)
        v->selected = true;
    }

    // Check that the neighbors around each cell are consistent in T1 and T2
    IntFloatMap heatMapDiffNeighborsT1;
    IntFloatMap heatMapDiffNeighborsT2;
    for(IntLabelsMap::const_iterator it = parentDaughtersMap.begin();
        it != parentDaughtersMap.end(); ++it) {
      IntIntMap::iterator ip = labelsT1.find(it->first);
      if(ip == labelsT1.end())
        continue;
      int labelParent = ip->second;
      if(labelParent == 0)
        continue;

      set_vector<int> LDaughters = it->second;
      // case of cells without daughters (giant cell + periphery)
      if(LDaughters.size() == 1 and LDaughters.count(0) == 1) {
        IntLabelsMap::iterator ip = parentDaughtersMap.find(0);
        if(ip != parentDaughtersMap.end())
          LDaughters = ip->second;
      }

      set_vector<int> neighborsT1;
      set_vector<int> neighborsT2;
      set_vector<int> parentsNeighborsT2;

      // list the neighbors around parent cell in T1
      {
        IntLabelsMap::iterator ip = cm.centerNeighborsMapT1.find(labelParent);
        if(ip != cm.centerNeighborsMapT1.end())
          neighborsT1 = ip->second;
        else
          neighborsT1.insert(0);
      }

      // list the neighbors around daughter cells in T2
      forall(int labelD, LDaughters) {
        if(cm.centerNeighborsMapT2.count(labelD) != 0) {
          IntLabelsMap::iterator ip = cm.centerNeighborsMapT2.find(labelD);
          if(ip != cm.centerNeighborsMapT2.end())
            neighborsT2.insert(ip->second.begin(), ip->second.end());
        }
      }
      // take out daughter cells from neighborsT2
      set_vector<int> neighborsT2temp;
      forall(int neighbor, neighborsT2)
          if(LDaughters.count(neighbor) == 0)
          neighborsT2temp.insert(neighbor);
      neighborsT2 = neighborsT2temp;

      // look at the parents of neighbors in T2
      forall(int neighbor, neighborsT2) {
        IntIntAttr::iterator ip = daughterParentMap.find(neighbor);
        if(ip != daughterParentMap.end())
          parentsNeighborsT2.insert(ip->second);
        else
          parentsNeighborsT2.insert(0);
      }

      // compare the neighbors around parent and daughter cells, fill in heat map with differences
      heatMapDiffNeighborsT1[labelParent] = heatMapDiffNeighborsT2[labelParent]
          = parentsNeighborsT2 != neighborsT1 ? 1 : 0;

      if(parentsNeighborsT2.size() != neighborsT1.size())
        Information::out << "problem with cell: " << labelParent
                         << ", the number of neighbors changes. " << parentsNeighborsT2.size()
                         << " parents of neighbors in T2 vs. " << neighborsT1.size()
                         << " neighbors in T1." << endl;
    }

    // Display the difference in neighbor number in a heatmap
    mesh1->labelHeat().clear();
    mesh2->labelHeat().clear();
    forall(const IntFloatPair& p, heatMapDiffNeighborsT1)
        mesh1->labelHeat()[p.first] = p.second;
    forall(const IntFloatPair& p, heatMapDiffNeighborsT2)
        mesh2->labelHeat()[p.first] = p.second;

    mesh1->heatMapBounds() = Point2f(0, 1);
    mesh1->setShowLabel("Label Heat");
    mesh1->updateTriangles();
    mesh1->updateSelection();
    mesh2->heatMapBounds() = Point2f(0, 1);
    mesh2->setShowLabel("Label Heat");
    mesh2->updateTriangles();

    SETSTATUS("Cells with change in number of neighbors appear in red, vertices which "
              "could not be identified are selected on the first mesh.");
    return true;
  }


 bool CorrespondenceJunctions::checkNhbdNew(Mesh* mesh1, Mesh* mesh2, vvGraph &T1, vvGraph &T2, bool ShowVVCorrespondence){
    // Create a map of labels in both meshes, to check parent labeling or create it
    // we will later add 0 values in the label map for cells which don't have daughters
    // (in T1) or parent (in T2)
   std::clock_t prev_time = std::clock();
   std::clock_t cur_time = std::clock();

    IntIntMap labelsT1, labelsT2;
    forall(const vertex &c, T1)
        if(c->type == 'c')
        labelsT1[c->label] = c->label;
    forall(const vertex &c, T2)
        if(c->type == 'c')
        labelsT2[c->label] = c->label;


    std::map<VMPair,VMPair> labelVertPairing;
    // daughterParentMap will be either labels or parents depending on what is selected
    // If cell does not have parent in T1, key = label in T2, value = 0
    IntIntAttr &parents = mesh2->parents();
    IntIntAttr daughterParentMap;
    if(mesh2->useParents()) {
      Information::out << "Found parent labels, using parent labeling to check correspondence between junctions." << endl;
      // Check the parents for cells which are not in the mesh1 or mesh2
      for(IntIntAttr::const_iterator it = parents.begin(); it != parents.end(); it++) {
        int lParent = it->second;
        int lDaughter = it->first;
        if(labelsT1.count(lParent) == 0 or labelsT2.count(lDaughter) == 0)
          // RSS If we don't find it just skip it.
          continue;
        //return setErrorMessage(QString("Problems with parent label: %1, daughter label: %2."
        //     " First run Correct Parents, then Save and Load parents.").arg(lParent).arg(lDaughter));
        daughterParentMap[lDaughter] = lParent;
      }
    } else {
      // If there is no parent map stored in mesh, the meshes should have the same labels
      // create a parent map based on the same cells in both meshes
      Information::out << "No parent labels found, using labels to check correspondence between junctions." << endl;
      for(IntIntMap::const_iterator l = labelsT1.begin(); l != labelsT1.end(); l++) {
        int label = l->first;
        if(labelsT2.count(label) != 0)
          daughterParentMap[label] = label;
      }
    }

std::cout<<std::endl;
   cur_time = std::clock();
   std::cout<<"Time at time point 1: "<<cur_time-prev_time;
   prev_time = cur_time;
std::cout<<std::endl;

    // Look at the daughters (in T2) for each cell in T1
    // parentDaughtersMap: key = parent in T1, value = labels of daughters in T2
    // If there is no parent, key = 0, values = labels of orphans (in T2)
    IntLabelsMap parentDaughtersMap;
    for(IntIntMap::const_iterator l = labelsT2.begin(); l != labelsT2.end(); l++) {
      int labelDaughter = l->first;
      int labelParent;
      set_vector<int> LDaughters;
      // add to daughterParentMap 0 values for the cells that do not have parent in T1.
      // map the label of cells in T2 without parent with 0
      {
        IntIntAttr::iterator it = daughterParentMap.find(labelDaughter);
        if(it == daughterParentMap.end() or it->second == 0 or labelDaughter == 0) {
          labelsT2[labelDaughter] = 0;
          daughterParentMap[labelDaughter] = 0;
          labelParent = 0;
        } else {
          labelParent = it->second;
        }
      }
      IntLabelsMap::iterator it = parentDaughtersMap.find(labelParent);
      if(it != parentDaughtersMap.end())
        LDaughters = it->second;
      LDaughters.insert(labelDaughter);
      parentDaughtersMap[labelParent] = LDaughters;
    }

std::cout<<std::endl;
   cur_time = std::clock();
   std::cout<<"Time at time point 2: "<<cur_time-prev_time;
   prev_time = cur_time;
std::cout<<std::endl;

    // add 0 values in parentDaughtersMap for cells in T1 that do not have daughters
    // map the label of cells in T1 without daughters with 0
    set_vector<int> noDaughter;
    noDaughter.insert(0);
    for(IntIntMap::const_iterator l = labelsT1.begin(); l != labelsT1.end(); l++) {
      int labelParent = l->first;
      if(parentDaughtersMap.count(labelParent) == 0) {
        parentDaughtersMap[labelParent] = noDaughter;
        labelsT1[labelParent] = 0;
      }
    }


std::cout<<std::endl;
   cur_time = std::clock();
   std::cout<<"Time at time point 3: "<<cur_time-prev_time;
   prev_time = cur_time;
std::cout<<std::endl;

    // Fill in a map of the labels associated with each junction in T1 and T2.
    // The list of labels is sorted (it is a set).
    // If there are only 2 neighbors, then the edge is at the mesh border and
    // we say it is in contact with the outside (label 0);
    // Labels of cells of T1 that do not have daughters are stored as 0.
    // This will be used to identify junctions at each time points.
    VertLabelsMap junctionLabelsT1;
    typedef std::map<int, std::vector<vertex> > LabelsVertMap;
    LabelsVertMap labelJunctionsT1;
    forall(const vertex &j, T1) {
      set_vector<int> labels;
      if(j->type == 'j') {
        // list the labels around junction
        forall(const vertex &n, T1.neighbors(j))
            if(n->label > 0) {
          IntIntMap::iterator it = labelsT1.find(n->label);
          if(it != labelsT1.end()){
            labels.insert(it->second);

//          VMPair cur_pair(it->second,j);
//          labelVertPairing[cur_pair]=cur_pair;
            (labelJunctionsT1[it->second]).push_back(j);
          }
        }
        if(labels.size() == 2){
          labels.insert(0);
          (labelJunctionsT1[0].push_back(j));
        }
        junctionLabelsT1[j] = labels;
      }
    }
    VertLabelsMap junctionLabelsT2;
//    LabelsVertMap labelJunctionsT2;
    forall(const vertex &j, T2) {
      set_vector<int> labels;
      if(j->type == 'j') {
        // list the labels around junction
        forall(const vertex &n, T2.neighbors(j))
            if(n->label > 0) {
          IntIntMap::iterator it = labelsT2.find(n->label);
          if(it != labelsT2.end()){
            labels.insert(it->second);
 //           (labelJunctionsT2[it->second]).push_back(j);
          }
        }
        if(labels.size() == 2){
          labels.insert(0);
//          (labelJunctionsT2[0]).push_back(j);
        }
        junctionLabelsT2[j] = labels;

      }
    }

std::cout<<std::endl;
   cur_time = std::clock();
   std::cout<<"Time at time point 4: "<<cur_time-prev_time;
   prev_time = cur_time;
std::cout<<std::endl;


   int loop1 = 0;
   int loop2 = 0;
   std::clock_t timer1;
   std::clock_t timer2;

    // Look for the correspondence between junctions in T1 and in T2
    set_vector<int> labelsDaughters;
    set_vector<int> labelsParents;
    IntIntAttr::iterator dPmapEnd = daughterParentMap.end();
    forall(const vertex &j2, T2) {
//      timer1 = std::clock();
      if(junctionLabelsT2.count(j2) > 0) {
        labelsDaughters.clear();
        VertLabelsMap::iterator it = junctionLabelsT2.find(j2);
        if(it != junctionLabelsT2.end())
          labelsDaughters = it->second;
        labelsParents.clear();
        forall(const int &lDaughter, labelsDaughters) {
          IntIntAttr::iterator it = daughterParentMap.find(lDaughter);
          int lParent = 0;
          if(it != dPmapEnd)
            lParent = it->second;
          labelsParents.insert(lParent);
        }
//        loop1+= std::clock()-timer1;
        // here compare the labels of parents of junction in T2 with junction in T1.
        // we considere only 3 ways junctions.
//        timer1 = std::clock();
        if(labelsParents.size() == 3){
          int cur_label = labelsParents[0] > labelsParents[1] ? labelsParents[0] : labelsParents[1]; //avoid boundary label
          std::vector<vertex> potenial_corespondences = labelJunctionsT1[cur_label];
          for(int i=0;i<potenial_corespondences.size();i++){
            const set_vector<int> & labelsJunction = junctionLabelsT1[potenial_corespondences[i]];
            if(labelsJunction == labelsParents){
                cm.J1J2[potenial_corespondences[i]] = j2;
                break;
            }

          }

/*
          VertLabelsMap::const_iterator jLT1end = junctionLabelsT1.end();
          for(VertLabelsMap::const_iterator it = junctionLabelsT1.begin();
              it != jLT1end; it++) {
            const set_vector<int>& labelsJunction = it->second;
            if(labelsJunction == labelsParents){
              cm.J1J2[it->first] = j2;
              break;
            }
          }
*/
        }
//        loop2+= std::clock()-timer1;
      }
    }

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

    cur_time = std::clock();
    std::cout<<"Time at time point 5: "<<cur_time-prev_time;
    prev_time = cur_time;
    std::cout<<std::endl;
    std::cout<<"loop1 is: "<<loop1<<std::endl;
    std::cout<<"loop2 is: "<<loop2<<std::endl;
std::cout<<std::endl;


    // Show the correspondence between vertices.
    IntMatrix3fMap& J1J2vectors = mesh1->vvCorrespondence();
    J1J2vectors.clear();
    mesh1->showVVCorrespondence() = ShowVVCorrespondence;
    int i = 0;
    for(VtxVtxMap::const_iterator it = cm.J1J2.begin(); it != cm.J1J2.end(); it++, i++) {
      vertex j1 = it->first;
      vertex j2 = it->second;
      J1J2vectors[i][0] = Point3f(j1->pos);
      J1J2vectors[i][1] = Point3f(j2->pos);
    }

    // Create a map of neighbors (labels) for each cell in T1 and T2.
    // in T1, labels of neighbors without daughters are stored as 0
    forall(const vertex &c, T1) {
      set_vector<int> labels;
      if(c->type == 'c' and c->label > 0) {
        // list the labels at junctions around the center
        forall(const vertex &j, T1.neighbors(c)) {
          if(junctionLabelsT1.count(j) != 0) {
            VertLabelsMap::iterator it = junctionLabelsT1.find(j);
            if(it != junctionLabelsT1.end())
              labels.insert(it->second.begin(), it->second.end());
          }
        }
        labels.erase(c->label);
        cm.centerNeighborsMapT1[c->label] = labels;
      }
    }
    // in T2, labels of neighbors without parents are stored as 0
    forall(const vertex &c, T2) {
      set_vector<int> labels;
      if(c->type == 'c' and c->label > 0) {
        // list the labels at junctions around the center
        forall(const vertex &j, T2.neighbors(c)) {
          if(junctionLabelsT2.count(j) != 0) {
            VertLabelsMap::iterator it = junctionLabelsT2.find(j);
            if(it != junctionLabelsT2.end())
              labels.insert(it->second.begin(), it->second.end());
          }
        }
        labels.erase(c->label);
        cm.centerNeighborsMapT2[c->label] = labels;
      }
    }

//std::cout<<std::endl;
//    cur_time = std::clock();
//    std::cout<<"Time at time point 6: "<<cur_time-prev_time;
//    prev_time = cur_time;
//std::cout<<std::endl;

    // Select the vertices in T1 which cause trouble
    set_vector<vertex> problemJunctionT1;
    Information::out << "Check problems with correspondence between junctions:" << endl;
    J1J2forCell(mesh1, T1, labelsT1, parentDaughtersMap, problemJunctionT1);
    Information::out << problemJunctionT1.size() << " problem junctions found in "
                     << labelsT1.size() << " cells" << endl;
    forall(const vertex &v, T1) {
      v->selected = false;
      if(problemJunctionT1.count(v) != 0)
        v->selected = true;
    }

    // Check that the neighbors around each cell are consistent in T1 and T2
    IntFloatMap heatMapDiffNeighborsT1;
    IntFloatMap heatMapDiffNeighborsT2;
    for(IntLabelsMap::const_iterator it = parentDaughtersMap.begin();
        it != parentDaughtersMap.end(); ++it) {
      IntIntMap::iterator ip = labelsT1.find(it->first);
      if(ip == labelsT1.end())
        continue;
      int labelParent = ip->second;
      if(labelParent == 0)
        continue;

      set_vector<int> LDaughters = it->second;
      // case of cells without daughters (giant cell + periphery)
      if(LDaughters.size() == 1 and LDaughters.count(0) == 1) {
        IntLabelsMap::iterator ip = parentDaughtersMap.find(0);
        if(ip != parentDaughtersMap.end())
          LDaughters = ip->second;
      }

      set_vector<int> neighborsT1;
      set_vector<int> neighborsT2;
      set_vector<int> parentsNeighborsT2;

      // list the neighbors around parent cell in T1
      {
        IntLabelsMap::iterator ip = cm.centerNeighborsMapT1.find(labelParent);
        if(ip != cm.centerNeighborsMapT1.end())
          neighborsT1 = ip->second;
        else
          neighborsT1.insert(0);
      }

      // list the neighbors around daughter cells in T2
      forall(int labelD, LDaughters) {
        if(cm.centerNeighborsMapT2.count(labelD) != 0) {
          IntLabelsMap::iterator ip = cm.centerNeighborsMapT2.find(labelD);
          if(ip != cm.centerNeighborsMapT2.end())
            neighborsT2.insert(ip->second.begin(), ip->second.end());
        }
      }
      // take out daughter cells from neighborsT2
      set_vector<int> neighborsT2temp;
      forall(int neighbor, neighborsT2)
          if(LDaughters.count(neighbor) == 0)
          neighborsT2temp.insert(neighbor);
      neighborsT2 = neighborsT2temp;

      // look at the parents of neighbors in T2
      forall(int neighbor, neighborsT2) {
        IntIntAttr::iterator ip = daughterParentMap.find(neighbor);
        if(ip != daughterParentMap.end())
          parentsNeighborsT2.insert(ip->second);
        else
          parentsNeighborsT2.insert(0);
      }

      // compare the neighbors around parent and daughter cells, fill in heat map with differences
      heatMapDiffNeighborsT1[labelParent] = heatMapDiffNeighborsT2[labelParent]
          = parentsNeighborsT2 != neighborsT1 ? 1 : 0;

      if(parentsNeighborsT2.size() != neighborsT1.size())
        Information::out << "problem with cell: " << labelParent
                         << ", the number of neighbors changes. " << parentsNeighborsT2.size()
                         << " parents of neighbors in T2 vs. " << neighborsT1.size()
                         << " neighbors in T1." << endl;
    }

//std::cout<<std::endl;
//    cur_time = std::clock();
//    std::cout<<"Time at time point 7: "<<cur_time-prev_time;
//    prev_time = cur_time;
//std::cout<<std::endl;


    // Display the difference in neighbor number in a heatmap
    mesh1->labelHeat().clear();
    mesh2->labelHeat().clear();
    forall(const IntFloatPair& p, heatMapDiffNeighborsT1)
        mesh1->labelHeat()[p.first] = p.second;
    forall(const IntFloatPair& p, heatMapDiffNeighborsT2)
        mesh2->labelHeat()[p.first] = p.second;

    mesh1->heatMapBounds() = Point2f(0, 1);
    mesh1->setShowLabel("Label Heat");
    mesh1->updateTriangles();
    mesh1->updateSelection();
    mesh2->heatMapBounds() = Point2f(0, 1);
    mesh2->setShowLabel("Label Heat");
    mesh2->updateTriangles();

    SETSTATUS("Cells with change in number of neighbors appear in red, vertices which "
              "could not be identified are selected on the first mesh.");


//std::cout<<std::endl;
//    cur_time = std::clock();
//    std::cout<<"Time at time point 8: "<<cur_time-prev_time;
//    prev_time = cur_time;




    return true;
  }


   // find cell junctions of a cell mesh in the correct order
  void findCellJunctionsInOrder(Mesh* mesh, vvGraph& S, JunctionInformation& jInfo)
  {

    findCellJunctions(mesh, S, jInfo);

    CellJunctionsInOrderMap cellJunctions;
    forall(const vertex &c, S) {
      if(c->type != 'c') continue;
      // cell
      forall(const vertex &j, S.neighbors(c)) {
        if(j->type != 'j') continue;
        if(jInfo.junctionLabels1.find(j) == jInfo.junctionLabels1.end()) continue;
        cellJunctions[c->label].push_back(j);
      }
    }
    jInfo.cellJunctionInOrder1 = cellJunctions;
  }

   // find cell junctions of a cell mesh
  void findCellJunctions(Mesh* mesh, vvGraph& S, JunctionInformation& jInfo)
  {

    bool parentLabels = false;

    // Fill in a map of the labels associated with each junction in T1 and T2.
    // The list of labels is sorted (it is a set).
    // If there are only 2 neighbors, then the edge is at the mesh border and
    // we say it is in contact with the outside (label 0);
    // Labels of cells of T1 that do not have daughters are stored as 0.
    // This will be used to identify junctions at each time points.
    JunctionCellsMap junctionLabels;
    CellJunctionsMap cellJunctions;
    forall(const vertex &j, S) {
      std::set<int> labels;
      if(j->type != 'j') continue;
      //j->selected = true;
      // list the labels around junction
      std::map<int, int> neighborLabelCount;
      forall(const vertex &n, S.neighbors(j)){
        if(n->label < 1) continue;
        if(parentLabels){
          neighborLabelCount[mesh->parents()[n->label]]++;
        } else {
          neighborLabelCount[n->label]++;
        }

      }

      bool newJ = false;
      forall(IntInt p, neighborLabelCount){
        if(p.second > 1 and p.first > 0){
          newJ = true;
        }
      }
      if(!newJ and neighborLabelCount.size() > 1){
        forall(const vertex &n, S.neighbors(j)){
          int label = parentLabels ? mesh->parents()[n->label] : n->label;
          if(label < 1) continue;
          cellJunctions[label].insert(j);
          labels.insert(label);
        }
        // std::cout << "findCellJ " << j->saveId << "/" << labels.size() << std::endl;
        if(labels.size() >= 3 or (labels.size()>=2 and j->margin))
          junctionLabels[j] = labels;

      }
    }

      jInfo.junctionLabels1 = junctionLabels;
      jInfo.cellJunctions1 = cellJunctions;
  }


   // find cell junctions of a cell mesh
  void findCellJunctionsParents(Mesh* mesh, vvGraph& S, JunctionInformation& jInfo)
  {

    bool parentLabels = true;

    // Fill in a map of the labels associated with each junction in T1 and T2.
    // The list of labels is sorted (it is a set).
    // If there are only 2 neighbors, then the edge is at the mesh border and
    // we say it is in contact with the outside (label 0);
    // Labels of cells of T1 that do not have daughters are stored as 0.
    // This will be used to identify junctions at each time points.
    JunctionCellsMap junctionLabels;
    CellJunctionsMap cellJunctions;
    forall(const vertex &j, S) {
      std::set<int> labels;
      if(j->type != 'j') continue;
      //j->selected = true;
      // list the labels around junction
      std::map<int, int> neighborLabelCount;
      forall(const vertex &n, S.neighbors(j)){
        if(n->label < 1) continue;
        if(parentLabels){
          neighborLabelCount[mesh->parents()[n->label]]++;
        } else {
          neighborLabelCount[n->label]++;
        }

      }

      bool newJ = false;
      forall(IntInt p, neighborLabelCount){
        if(p.second > 1 and p.first > 0){
          newJ = true;
        }
      }
      if(!newJ and neighborLabelCount.size() > 1){
        forall(const vertex &n, S.neighbors(j)){
          int label = parentLabels ? mesh->parents()[n->label] : n->label;
          if(label < 1) continue;
          cellJunctions[label].insert(j);
          labels.insert(label);
        }

        junctionLabels[j] = labels;

      }
    }

      jInfo.junctionLabels2 = junctionLabels;
      jInfo.cellJunctions2 = cellJunctions;
  }

  // PDGs computation

  // compute polygon center:
  Point3d cellCenter(const std::vector<Point3d>& polygon)
  {
    Point3d center = Point3d(0, 0, 0);
    int nbSides = polygon.size();
    for(int i = 0; i < nbSides; i++)
      center += polygon[i];

    center /= nbSides;
    return center;
  }

  // compute polygon normal:
  Point3d cellNormal(const std::vector<Point3d>& polygon, const Point3d& centerPolygon)
  {
    Point3d normalPolygon = Point3d(0, 0, 0);
    int nbSides = polygon.size();
    int prev = nbSides - 1;
    for(int i = 0; i < nbSides; i++) {
      normalPolygon += (polygon[prev] - centerPolygon) % (polygon[i] - centerPolygon);     // cross product
      prev = i;
    }
    normalPolygon.normalize();
    return normalPolygon;
  }

  // multiply vector of point3d by matrix:
  std::vector<Point3d> matrixMultiply(const std::vector<Point3d>& P, const Matrix3d& M)
  {
    int nbRows = P.size();
    std::vector<Point3d> result(nbRows);
    for(int i = 0; i < nbRows; ++i)
      result[i] = P[i] * M;
    return result;
  }

  // for two matrices P and M of size 2xnbRows, multiply: transpose(P) * M == the dot produc
  Matrix2d dotProduct(const std::vector<Point2d>& P, const std::vector<Point2d>& M)
  {
    Matrix2d prod;
    for(size_t i = 0; i < P.size(); ++i) {
      prod(0, 0) += P[i][0] * M[i][0];
      prod(0, 1) += P[i][0] * M[i][1];
      prod(1, 0) += P[i][1] * M[i][0];
      prod(1, 1) += P[i][1] * M[i][1];
    }
    return prod;
  }

  inline double square(double x) {
    return x * x;
  }

  // SVD decomposition of a 2D matrix M in: rotation(-theta)*[maxS,0; 0, minS]*rotation(psi)
  void SVDMatrix2x2(Matrix2d M, double& theta, double& psi, double& maxS, double& minS)
  {
    double c1, c2, c3, c4;
    c1 = sqrt(square(M(0, 0) + M(1, 1)) + square(M(0, 1) - M(1, 0)));
    c2 = sqrt(square(M(0, 0) - M(1, 1)) + square(M(0, 1) + M(1, 0)));
    c3 = atan2(M(0, 1) + M(1, 0), M(0, 0) - M(1, 1));
    c4 = atan2(M(0, 1) - M(1, 0), M(0, 0) + M(1, 1));

    maxS = (c1 + c2) / 2;
    minS = (c1 - c2) / 2;
    theta = (c3 - c4) / 2;
    psi = (c3 + c4) / 2;
  }

  // compute Principal Directions of Growth (PDGs) in 2D
  void computePDGs(const std::vector<vertex>& J1forC1, const std::vector<vertex>& J2forC2, SymmetricTensor& CellPDGs_1,
                   SymmetricTensor& CellPDGs_2)
  {
    int nbSides = J1forC1.size();
    std::vector<Point3d> polygon1(nbSides), polygon2(nbSides);

    for(int i = 0; i < nbSides; i++) {
      polygon1[i] = J1forC1[i]->pos;
      polygon2[i] = J2forC2[i]->pos;
    }

    // Compute transformation matrix for the two polygons
    // Compute polygon center (3D)
    const Point3d& centerPolygon1 = cellCenter(polygon1);
    const Point3d& centerPolygon2 = cellCenter(polygon2);

    // Compute polygon normal (3D)
    const Point3d& normalPolygon1 = cellNormal(polygon1, centerPolygon1);
    const Point3d& normalPolygon2 = cellNormal(polygon2, centerPolygon2);

    // Project polygon 3d points on the cell plane and center the polygon on (0,0,0)
    std::vector<Point3d> polygonProj1(nbSides), polygonProj2(nbSides);

    for(int i = 0; i < nbSides; i++) {
      // pos=Point to project (Point3d), p=Point of the plane(Point3d),n=Normal to the plane(point3d)
      polygonProj1[i] = projectPointOnPlane(polygon1[i] - centerPolygon1, Point3d(), normalPolygon1);
      polygonProj2[i] = projectPointOnPlane(polygon2[i] - centerPolygon2, Point3d(), normalPolygon2);
    }
    // Rotate polygon plane to align it with (x,y) plane (new code from triangulate.hpp/cpp, changed Dec 2015)
    Matrix3d invRotPolygonPlane1;
    rotatePointsIntoXY(normalPolygon1, polygonProj1, invRotPolygonPlane1);

    // same for the second polygon
    Matrix3d invRotPolygonPlane2;
    rotatePointsIntoXY(normalPolygon2, polygonProj2, invRotPolygonPlane2);

    // Take out 3rd dimension from projected polygons:
    std::vector<Point2d> polygon1_2D(nbSides), polygon2_2D(nbSides);
    for(int i = 0; i < nbSides; ++i) {
      polygon1_2D[i] = Point2d(polygonProj1[i].x(), polygonProj1[i].y());
      polygon2_2D[i] = Point2d(polygonProj2[i].x(), polygonProj2[i].y());
    }

    // Compute transformation matrix between c1 and c2(in 2D, on projected polygons):
    Matrix2d XTX = dotProduct(polygon1_2D, polygon1_2D);
    Matrix2d XTY = dotProduct(polygon1_2D, polygon2_2D);
    Matrix2d transfM = inverse(XTX) * XTY;

    // Decompose transformation matrix into rotations and growth:
    double theta, psi, maxStretchRatio, minStretchRatio;
    SVDMatrix2x2(transfM, theta, psi, maxStretchRatio, minStretchRatio);

    // Convert transformation matrix from 2D to 3D:
    Matrix3d transfM_3D;
    for(int i = 0; i < 2; ++i)
      for(int j = 0; j < 2; ++j)
        transfM_3D(i, j) = transfM(i, j);

    Matrix3d rotPolygonPlane2 = inverse(invRotPolygonPlane2);
    Matrix3d rotPolygonPlane1 = inverse(invRotPolygonPlane1);
    transfM_3D = invRotPolygonPlane1 * transfM_3D * rotPolygonPlane2;

    // Apply 3D transformation to polygon1 (cell before growth), in order to compute residuals
    // (difference between actual and computed transformation):
    std::vector<Point3d> polygon1Centered, polygon1Transf;
    std::vector<double> residuals;
    polygon1Centered.resize(nbSides);
    residuals.resize(nbSides);

    for(int i = 0; i < nbSides; ++i)
      polygon1Centered[i] = polygon1[i] - centerPolygon1;
    polygon1Transf = matrixMultiply(polygon1Centered, transfM_3D);

    for(int i = 0; i < nbSides; ++i) {
      polygon1Transf[i] = polygon1Transf[i] + centerPolygon1;
      residuals[i] = norm(polygon1Transf[i] - polygon2[i]);
    }

    // Return the Principal Directions for polygon2 (deformed configuration):
    CellPDGs_2.ev1() = Point3f(Point3d(cos(psi), sin(psi), 0) * rotPolygonPlane2);    // unit vector PDG max
    CellPDGs_2.ev2() = Point3f(Point3d(-sin(psi), cos(psi), 0) * rotPolygonPlane2);   // unit vector PDG min
    CellPDGs_2.evals() = Point3f(maxStretchRatio, minStretchRatio, 0);

    // Return the Principal Directions for polygon1 (original configuration):
    CellPDGs_1.ev1() = Point3f(Point3d(cos(theta), sin(theta), 0) * rotPolygonPlane1);
    CellPDGs_1.ev2() = Point3f(Point3d(-sin(theta), cos(theta), 0) * rotPolygonPlane1);
    CellPDGs_1.evals() = Point3f(maxStretchRatio, minStretchRatio, 0);
  }

  bool GrowthDirections::run(Mesh* mesh1, Mesh* mesh2)
  {
    mesh1->clearCellAxis();
    mesh2->clearCellAxis();
    // Get the meshes
    mesh1->updateCentersNormals();
    mesh2->updateCentersNormals();

    // use the method CellAxis to modify a map associating labels with a vector of
    // 3 Point3f (the cell axis)
    IntSymTensorAttr &CellPDGs_1 = mesh1->cellAxis();
    IntSymTensorAttr &CellPDGs_2 = mesh2->cellAxis();
    CellPDGs_1.clear();
    CellPDGs_2.clear();

    // Calculate the PDGs
    forall(const IntVtxVecPairPair &pr, cm.J1J2forC1)
        computePDGs(cm.J1J2forC1[pr.first].first, cm.J1J2forC1[pr.first].second,
        CellPDGs_1[pr.first], CellPDGs_2[pr.first]);

    mesh1->setCellAxisType("PDG");
    mesh2->setCellAxisType("PDG");
    mesh1->setCellAxisUnit("1");
    mesh2->setCellAxisUnit("1");

    Information::out << "Size of Cell PDG map in mesh1: " << CellPDGs_1.size() << endl;
    Information::out << "Size of Cell PDG map in mesh2: " << CellPDGs_2.size() << endl;

    AttrMap<int, SymmetricTensor>& pdgAttr1 = mesh1->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor PDGs");
    pdgAttr1.clear();
    AttrMap<int, SymmetricTensor>& pdgAttr2 = mesh2->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor PDGs");
    pdgAttr2.clear();

    forall(auto p, CellPDGs_1){
      pdgAttr1[p.first] = p.second;
    }
    forall(auto p, CellPDGs_2){
      pdgAttr2[p.first] = p.second;
    }

    // Run PDG display on both meshes, with the parameters from the GUI
    DisplayPDGs proc(*this);

    proc.run(mesh1);
    proc.run(mesh2);

    return true;
  }
  REGISTER_PROCESS(GrowthDirections);

  void DisplayPDGs::processParms()
  {
    DisplayHeatMap = parm("Heatmap");
    ScaleHeatMap = parm("ScaleHeat");
    RangeHeat = Point2d(parm("Heat min").toDouble(), parm("Heat max").toDouble());
    DisplayAxis = parm("Show Axis");
    ColorPos = parm("Color +");
    ColorNeg = parm("Color -");
    AxisWidth = parm("Line Width").toFloat();
    AxisScale = parm("Line Scale").toFloat();
    AxisOffset = parm("Line Offset").toFloat();
    AnisotropyThreshold = parm("Threshold").toFloat();
    CustomShear = stringToBool(parm("Custom Dir Shear"));
  }

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

    // Check values of parameters and change them in the GUI if needed.
    const IntSymTensorAttr& cellAxisCustom = currentMesh()->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    // If PDGs not projected on Bezier yet, make sure we won't try to display the projection:
    if(cellAxisCustom.size() == 0) {
      QString heatMap = parm("Heatmap");
      if(heatMap == "StretchCustomX" or heatMap == "StretchCustomY" or heatMap == "Shear"
        or heatMap == "AnisoCustom" or heatMap == "AnisoRatio"){
        setParm("Heatmap", "StretchMax");
        std::cout << "No Custom Directions available, switching to Default Heat Map Option." << std::endl;
      }

      QString axis = parm("Show Axis");
      if(axis == "StrainCustomX" or axis == "StrainCustomY" or axis == "BothCustom"){
        setParm("Show Axis", "StrainMax");
        std::cout << "No Custom Directions available, switching to Default Axis Option." << std::endl;
      }
    }

    return true;
  }

  // Display PDGs (separatly from each other) and heatmaps of anisotropy etc.
  // The stretch vectors are stored in cellAxis, so that the norm of the vector give us the stretch (always positive)
  // The strain (= stretch -1) vectors are visualized in cellAxisVis. No deformation ==> strain = 0 ==> cell axis null
  bool DisplayPDGs::run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
                        const QString displayAxis, const QColor& colorPos, const QColor& colorNeg, float axisWidth, float axisScale,
                        float axisOffset, float anisotropyThreshold, bool customShear)
  {
    // if there is no PDGs stored, display error message
    const IntSymTensorAttr& cellAxis = mesh->cellAxis();
    IntSymTensorAttr& cellAxisCustom = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    if(cellAxis.size() == 0)
      throw(QString("No cell axis stored in active mesh"));
    if(mesh->cellAxisType() != QString("PDG"))
      throw(QString("No PDGs stored in active mesh"));

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

    // Check if there is a cell center for each cell axis
    mesh->updateCentersNormals();
    IntPoint3fAttr &centers = (mesh->useParents() ? mesh->parentCenterVis() : mesh->labelCenterVis());
    int nProblems = 0;
    forall(const IntSymTensorPair &p, cellAxis) {
      int label = p.first;
      if(centers.count(label) == 0)
        nProblems++;
    }
    if(nProblems != 0)
      SETSTATUS("Warning, non-existing cell center found for " << nProblems << " cells.");

    // Scale the cell axis for display
    IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
    cellAxisVis.clear();
    IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();
    cellAxisColor.clear();

    forall(const IntSymTensorPair &p, cellAxis) {
      int cell = p.first;
      SymmetricTensor tensor = p.second;
      float strain1st, strain2nd, anisotropy;
      Point3f dir1st, dir2nd;

      // if any "custom" option is selected, calculate the custom directions
      if(displayAxis == "StrainCustomX" or displayAxis == "StrainCustomY" or displayAxis == "BothCustom" or
        displayHeatMap == "StretchCustomX" or displayHeatMap == "StretchCustomY" or displayHeatMap == "Shear" or displayHeatMap == "StretchCustomX-StretchCustomY"
        or displayHeatMap == "AnisoCustom" or displayHeatMap == "AnisoRatio" or displayHeatMap == "CustomX Growth Proportion" or displayHeatMap == "CustomY Growth Proportion"){

        if(cellAxisCustom.size() == 0)
          throw(QString("No custom cell axis. Create a custom cell axis first."));
        if(cellAxisCustom.find(cell) == cellAxisCustom.end())
          continue;

        Point3f customX = cellAxisCustom[cell].ev1();
        Point3f customY = cellAxisCustom[cell].ev2();
        cellAxisCustom[cell] = customDirectionTensor(tensor, customX, customY, customShear);
      }

      // if custom cell axis required, grab it
      if(displayAxis == "StrainCustomX" or displayAxis == "StrainCustomY" or displayAxis == "BothCustom" or displayHeatMap == "StretchCustomX-StretchCustomY"){
        tensor = cellAxisCustom[cell];
      }


      strain1st = tensor.evals()[0] - 1;
      strain2nd = tensor.evals()[1] - 1;
      dir1st = tensor.ev1();
      dir2nd = tensor.ev2();
      anisotropy = max(tensor.evals()[0],tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);

      if(displayAxis == "StretchMax-StretchMin"){
        strain1st = max(tensor.evals()[0], tensor.evals()[1]) - min(tensor.evals()[0], tensor.evals()[1]);
        strain2nd = 0;
      }

      Point3b showAxis;
      showAxis[0] = (displayAxis == "StrainMax" or displayAxis == "StrainCustomX" or displayAxis == "Both" or displayAxis == "BothCustom" or displayAxis == "StretchMax-StretchMin");
      showAxis[1] = (displayAxis == "StrainMin" or displayAxis == "StrainCustomY" or displayAxis == "Both" or displayAxis == "BothCustom");
      showAxis[2] = false;

      if(anisotropyThreshold != 0 and anisotropy < anisotropyThreshold)
        showAxis = Point3b(false, false, false);
      std::cout<<strain1st<<":"<<dir1st<<std::endl;
      cellAxisVis[cell][0] = (showAxis[0] ? strain1st : 0.f) * dir1st * axisScale;
      cellAxisVis[cell][1] = (showAxis[1] ? strain2nd : 0.f) * dir2nd * axisScale;
      cellAxisVis[cell][2] = Point3f(0, 0, 0);

      for(size_t i = 0; i < 3; ++i) {
        if(tensor.evals()[i] - 1 < 0)
          cellAxisColor[cell][i] = Colorb(colorShrinkage);
        else
          cellAxisColor[cell][i] = Colorb(colorExpansion);
      }
    }

    // Fill in heat map table
    std::vector<float> vecValues;
    if(displayHeatMap != "None") {
      mesh->labelHeat().clear();
      IntFloatAttr& labelHeatMap = mesh->labelHeat();
      forall(const IntSymTensorPair &p, cellAxis) {
        int cell = p.first;
        SymmetricTensor tensor = p.second;

        float value = 0;

        if(displayHeatMap == "StretchMax")
          value = max(tensor.evals()[0], tensor.evals()[1]);
        if(displayHeatMap == "StretchMin")
          value = min(tensor.evals()[0], tensor.evals()[1]);
        if(displayHeatMap == "Aniso=StretchMax/StretchMin")
          value = max(tensor.evals()[0], tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
        if(displayHeatMap == "StretchMax-StretchMin")
          value = max(tensor.evals()[0], tensor.evals()[1]) - min(tensor.evals()[0], tensor.evals()[1]);
        if(displayHeatMap == "ProductStretches")
          value = tensor.evals()[0] * tensor.evals()[1];

        if(displayHeatMap == "CustomX Growth Proportion" || displayHeatMap == "CustomY Growth Proportion"){
          SymmetricTensor cus_tensor = cellAxisCustom.find(cell)->second;
          if(displayHeatMap == "CustomX Growth Proportion")
              value = cus_tensor.evals()[0];
          if(displayHeatMap == "CustomY Growth Proportion")
              value = cus_tensor.evals()[1];
          value/=(cus_tensor.evals()[0]+cus_tensor.evals()[1]);
        }

        // if custom heat map is required, grab it
        if(displayHeatMap == "StretchCustomX" or displayHeatMap == "StretchCustomY" or displayHeatMap == "Shear" or
          displayHeatMap == "AnisoCustom" or displayHeatMap == "AnisoRatio" or displayHeatMap == "StretchCustomX-StretchCustomY"){
          if(cellAxisCustom.size() == 0)
            throw(QString("No custom cell axis. Run 'Cell Axis Custom Directions' first."));
          if(cellAxisCustom.find(cell) == cellAxisCustom.end())
            continue;
          tensor = cellAxisCustom.find(cell)->second;
          if(displayHeatMap == "StretchCustomX")
            value = tensor.evals()[0];
          if(displayHeatMap == "StretchCustomY")
            value = tensor.evals()[1];
          if(displayHeatMap == "StretchCustomX-StretchCustomY")
            value = tensor.evals()[0]-tensor.evals()[1];
          if(displayHeatMap == "Shear")
            value = fabs(tensor.evals()[2]);
          if(displayHeatMap == "AnisoCustom")
            value = max(tensor.evals()[0], tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
          if(displayHeatMap == "AnisoRatio"){
            double anisoPDG = max(p.second.evals()[0], p.second.evals()[1])/min(p.second.evals()[0], p.second.evals()[1]);
            double anisoCus = max(tensor.evals()[0], tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
            value = max(anisoPDG, anisoCus)/min(anisoPDG, anisoCus);
          }
        }

        labelHeatMap[cell] = value;
        vecValues.push_back(value);
      }
    }

    // Find bounds for heatmap scale.
    float lowBound = 0, highBound = 1;
    // Manual scaling by user input
    if(scaleHeatMap == "Manual"){
      lowBound  = rangeHeat.x();
      highBound = rangeHeat.y();
    }
    // Automatic scaling
    if(vecValues.size() != 0){
      sort(vecValues.begin(), vecValues.end());
      // Default scaling: take min and max values
      if(scaleHeatMap == "None"){
        lowBound = vecValues.front();
        highBound = vecValues.back();
      }
      // Auto scale: take low 5 and high 5 percentile
      else if(scaleHeatMap == "Auto"){
        lowBound = *(vecValues.begin()+int(vecValues.size()*.05));
        highBound = *(vecValues.begin()+int(vecValues.size()*.95));
      }
    }
    // Round the bounds values of heat map, to make sure we have the same bound in case
    // of loading PDGs on both meshes.
    lowBound = floor(lowBound * 1000) / 1000;
    highBound = ceil(highBound * 1000) / 1000;

    mesh->heatMapBounds() = Point2f(lowBound, highBound);
    mesh->heatMapUnit() = QString("ratio");
    mesh->setShowLabel("Label Heat");
    mesh->updateTriangles();

    if(displayAxis != "None")
      mesh->setAxisView("Cell Axis");

    return true;
  }
  REGISTER_PROCESS(DisplayPDGs);

  void DisplayPDGsVertex::processParms()
  {
    DisplayHeatMap = parm("Heatmap");
    ScaleHeatMap = parm("ScaleHeat");
    RangeHeat = Point2d(parm("Heat min").toDouble(), parm("Heat max").toDouble());
    DisplayAxis = parm("Show Axis");
    ColorPos = parm("Color +");
    ColorNeg = parm("Color -");
    AxisWidth = parm("Line Width").toFloat();
    AxisScale = parm("Line Scale").toFloat();
    AxisOffset = parm("Line Offset").toFloat();
    AnisotropyThreshold = parm("Threshold").toFloat();
    vtxMinDis = parm("Min Dis Vtx").toDouble();
  }

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

    // Check values of parameters and change them in the GUI if needed.
    const IntSymTensorAttr& cellAxisCustom = currentMesh()->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    // If PDGs not projected on Bezier yet, make sure we won't try to display the projection:
    if(cellAxisCustom.size() == 0) {
      QString heatMap = parm("Heatmap");
      if(heatMap == "StretchCustomX" or heatMap == "StretchCustomY" or heatMap == "Shear"
        or heatMap == "AnisoCustom" or heatMap == "AnisoRatio"){
        setParm("Heatmap", "StretchMax");
        std::cout << "No Custom Directions available, switching to Default Heat Map Option." << std::endl;
      }

      QString axis = parm("Show Axis");
      if(axis == "StrainCustomX" or axis == "StrainCustomY" or axis == "BothCustom"){
        setParm("Heatmap", "StrainMax");
        std::cout << "No Custom Directions available, switching to Default Axis Option." << std::endl;
      }
    }

    return true;
  }

  // Display PDGs (separatly from each other) and heatmaps of anisotropy etc.
  // The stretch vectors are stored in cellAxis, so that the norm of the vector give us the stretch (always positive)
  // The strain (= stretch -1) vectors are visualized in cellAxisVis. No deformation ==> strain = 0 ==> cell axis null
  bool DisplayPDGsVertex::run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
                        const QString displayAxis, const QColor& colorPos, const QColor& colorNeg, float axisWidth, float axisScale,
                        float axisOffset, float anisotropyThreshold, bool customShear, double vtxMinDis)
  {
    // if there is no PDGs stored, display error message
    const VtxSymTensorAttr& vertexAxis = mesh->vertexAxis();
    //IntSymTensorAttr& cellAxisCustom = mesh->attributes().attrMap<int, SymmetricTensor>("Custom Directions");

    if(vertexAxis.size() == 0)
      throw(QString("No vertex axis stored in active mesh"));
    //if(mesh->cellAxisType() != QString("PDG"))
    //  throw(QString("No PDGs stored in active mesh"));

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

    // Scale the cell axis for display
    VtxMatrix3fAttr &vertexAxisVis = mesh->vertexAxisVis();
    vertexAxisVis.clear();
    VtxVec3ColorbAttr &vtxAxisColor = mesh->vertexAxisColor();
    vtxAxisColor.clear();

    double minSig = 1E20;
    double maxSig = -1E20;
    std::map<int, std::set<vertex> > labelVertexVisMap;


    forall(const auto &p, vertexAxis) {
      vertex v = p.first;
      SymmetricTensor tensor = p.second;
      float strain1st, strain2nd, anisotropy;
      Point3f dir1st, dir2nd;

      // // if any "custom" option is selected, calculate the custom directions
      // if(displayAxis == "StrainCustomX" or displayAxis == "StrainCustomY" or displayAxis == "BothCustom" or
      //   displayHeatMap == "StretchCustomX" or displayHeatMap == "StretchCustomY" or displayHeatMap == "Shear"
      //   or displayHeatMap == "AnisoCustom" or displayHeatMap == "AnisoRatio"){

      //   if(cellAxisCustom.size() == 0)
      //     throw(QString("No custom cell axis. Create a custom cell axis first."));
      //   if(cellAxisCustom.find(cell) == cellAxisCustom.end())
      //     continue;

      //   Point3f customX = cellAxisCustom[cell].ev1();
      //   Point3f customY = cellAxisCustom[cell].ev2();
      //   cellAxisCustom[cell] = customDirectionTensor(tensor, customX, customY, customShear);
      // }

      // // if custom cell axis required, grab it
      // if(displayAxis == "StrainCustomX" or displayAxis == "StrainCustomY" or displayAxis == "BothCustom"){
      //   tensor = cellAxisCustom[cell];
      // }

      int minDisVtx = INT_MAX;
      if(vtxMinDis > 0){
        forall(vertex w, labelVertexVisMap[v->label]){
          double dis = norm(w->pos - v->pos);
          if(dis <  minDisVtx) minDisVtx = dis;
        }
      }

      strain1st = tensor.evals()[0] - 1;
      strain2nd = tensor.evals()[1] - 1;
      dir1st = tensor.ev1();
      dir2nd = tensor.ev2();
      anisotropy = max(tensor.evals()[0],tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);



      if(displayHeatMap == "StretchMax"){
        v->signal = tensor.evals()[0];
      } else if(displayHeatMap == "StretchMin"){
        v->signal = tensor.evals()[1];
      } else if(displayHeatMap == "Aniso"){
        v->signal = tensor.evals()[1] > 1e-3 ? tensor.evals()[0] / tensor.evals()[1] : 0.0;
      } else if(displayHeatMap == "ProductStretches"){
        v->signal = tensor.evals()[0] * tensor.evals()[1];
      } else if(displayHeatMap == "StretchMax-StretchMin"){
        v->signal = tensor.evals()[0] - tensor.evals()[1];
      } else {
        v->signal = 0;
      }

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

      if(minDisVtx > vtxMinDis){
        Point3b showAxis;
        showAxis[0] = (displayAxis == "StrainMax" or displayAxis == "StrainCustomX" or displayAxis == "Both" or displayAxis == "BothCustom");
        showAxis[1] = (displayAxis == "StrainMin" or displayAxis == "StrainCustomY" or displayAxis == "Both" or displayAxis == "BothCustom");
        showAxis[2] = false;

        if(anisotropyThreshold != 0 and anisotropy < anisotropyThreshold)
          showAxis = Point3b(false, false, false);

        vertexAxisVis[v][0] = (showAxis[0] ? strain1st : 0.f) * dir1st * axisScale;
        vertexAxisVis[v][1] = (showAxis[1] ? strain2nd : 0.f) * dir2nd * axisScale;
        vertexAxisVis[v][2] = Point3f(0, 0, 0);

        for(size_t i = 0; i < 3; ++i) {
          if(tensor.evals()[i] - 1 < 0)
            vtxAxisColor[v][i] = Colorb(colorShrinkage);
          else
            vtxAxisColor[v][i] = Colorb(colorExpansion);
        }

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

    }

    mesh->signalBounds() = Point2f(minSig, maxSig);
    mesh->setShowVertex("Signal");
    mesh->setAxisView("Vertex Axis");
    mesh->updateTriangles();


    return true;
  }
  REGISTER_PROCESS(DisplayPDGsVertex);


  bool AveragePDGs::run(Mesh *mesh,bool projSurf,QString avgType){
    //PROJ ON SURF AND AREA WEIGHTING NOT IMPLEMENTED AT PRESENT
    //NEED TO TEST THE TOO/FROM MATRIX CONVERSION
  vvGraph& S = mesh->graph();

  IntSymTensorAttr& cellAxis = mesh->cellAxis();
  if(mesh->cellAxisType() != QString("PDG"))
      throw(QString("No PDGs stored in active mesh"));

  // create MGXC if not already
  // vvGraph T;
  // if(mesh->tissue().meshType() == "MGXM")
  //   mgxmToMgxc(S, T, 1);
  // else
  //   T = S;

//  void MeasureArea::calculateArea(vvGraph &S, IntFloatAttr &data){

  IntFloatAttr areaHeat;
  if(avgType == "area" || avgType == "area weighted geometric mean"){
    MeasureArea *ma;
    if(!getProcess("Mesh/Heat Map/Measures/Geometry/Area", ma))
      throw(QString("MeshProcessHeatMap:: Unable to make MeasureArea process"));
    ma->calculateArea(S, areaHeat);

  }

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

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

  // Update the normals and center of the cell for visualisation.
  //if(mesh->labelCenterVis().empty() or (mesh->useParents() and mesh->parentCenterVis().empty()))
  mesh->updateCentersNormals();

  // cell centers
  typedef std::pair<int,Point3f> IntPoint3fP;
  IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());
  IntPoint3fAttr normals = (mesh->useParents()? mesh->parentNormal():mesh->labelNormal());


  // merge neighbor graph
  std::map<int, std::set<int> > cellNeighborMapParents;
  if(mesh->useParents()){
    std::set<int> meshLabels = findAllLabelsSet(S);
    // go through labels
    forall(int l, meshLabels){
      forall(int neighb, cellNeighborMap[l]){
        if(mesh->parents()[l] != mesh->parents()[neighb])
          cellNeighborMapParents[mesh->parents()[l]].insert(mesh->parents()[neighb]);
      }
    }
    cellNeighborMap = cellNeighborMapParents;
  }


    std::vector<double> used_areas;
    double area_normalization = 1;
    if(avgType == "area weighted geometric mean"){
        //Compute geo-mean of all areas
        forall(IntPoint3fP p, centers){
            if(cellAxis.find(p.first) == cellAxis.end())
                continue;
            used_areas.push_back(areaHeat[p.first]);
        }
        area_normalization = geoMeanVectorValue(used_areas);
    }



//Anything required here to work over parents, aot daughter cells?
 forall(IntPoint3fP p, centers){
    if(cellAxis.find(p.first) == cellAxis.end())
        continue;

    Matrix3d avg_pdg;
    double weight = 1;

    int cell = p.first; // cell label
    Point3d cellNormal(normals[cell]);
    //if(mesh->useParents()) cell = mesh->parents()[cell];

    //WEIGHT BY AREA AS WELL
    //if(avgType=="n")
    //    weight++;


    if(avgType=="area" || avgType=="area weighted geometric mean")
        weight=areaHeat[p.first]/area_normalization;
    double weight_sum = weight;
    forall(int neighb, cellNeighborMap[p.first]){
            double cur_weight=0;
            if(cellAxis.find(neighb) == cellAxis.end())
                continue;
            if(avgType == "geometric mean" || avgType == "n")
                cur_weight = 1;
            else if(avgType != "none")
                cur_weight=areaHeat[neighb]/area_normalization;

            weight_sum+=cur_weight;
    }



    if(avgType!="area weighted geometric mean" && avgType!="geometric mean")
       avg_pdg = weight/weight_sum*cellAxis[p.first].toMatrix();
    if(avgType == "area weighted geometric mean" || avgType == "geometric mean"){
          SymmetricTensor curAxis = cellAxis[p.first];
          Point3f e_vals = curAxis.evals();
//          std::cout<<"taking matrix power"<<std::endl;
//          std::cout<<e_vals<<std::endl;
          for(int i=0;i<3;i++)
                curAxis.evals()[i] = pow(fabs(e_vals[i]),weight/weight_sum);
//          std::cout<<weight<<" "<<area_normalization<<std::endl;
//          std::cout<<curAxis.evals()<<std::endl;
          avg_pdg = curAxis.toMatrix();
    }


/*
    SymmetricTensor test = cellAxis[p.first];
    std::cout<<"Testing to from code"<<std::endl;
    std::cout<<"Original"<<std::endl;
    std::cout<<test.ev1()<<std::endl<<test.ev2()<<std::endl<<test.ev3()<<std::endl;
    std::cout<<test.evals()<<std::endl;
    test.fromMatrix(test.toMatrix());
    std::cout<<"To/From"<<std::endl;
    std::cout<<test.ev1()<<std::endl<<test.ev2()<<std::endl<<test.ev3()<<std::endl;
    std::cout<<test.evals()<<std::endl;
*/


//    if(avgType == "area weighted geometric mean" || avgType == "geometric mean"){ //Compute sum of all weights, so we can take the root as we go

//    }

    forall(int neighb, cellNeighborMap[p.first]){
        if(cellAxis.find(neighb) == cellAxis.end())
            continue;
        double cur_weight = 1;
        if(avgType=="area" || avgType == "area weighted geometric mean")
            cur_weight=areaHeat[neighb]/area_normalization;

        cur_weight/=weight_sum;
        SymmetricTensor curAxis = cellAxis[neighb];
        //Assume all negative eigens are numerical errors, and thus just take abs-value
        if(avgType == "area weighted geometric mean" || avgType == "geometric mean"){
            Point3f e_vals = curAxis.evals();
            for(int i=0;i<3;i++)
                curAxis.evals()[i] = pow(fabs(e_vals[i]),cur_weight);
        }

        Point3d neighNormal(normals[neighb]);
        Matrix3d rotMat = calcRotMatrix(cellNormal, neighNormal);
        if(avgType!="area weighted geometric mean" && avgType!="geometric mean")
            avg_pdg+= cur_weight*(transpose(rotMat)*(curAxis.toMatrix())*rotMat);
        if(avgType=="geometric mean" || avgType=="area weighted geometric mean")
            avg_pdg= (transpose(rotMat)*(curAxis.toMatrix())*rotMat)*avg_pdg;

//        std::cout<<cur_weight<<":"<<weight_sum<<std::endl;
        weight+=cur_weight;
    }
    std::cout<<weight<<":"<<weight_sum<<std::endl;
    //PROJECT ON THE SURFACE
    if(projSurf){

    }
//    if(avgType=="n" || avgType=="area")
//        avg_pdg = (1.0f/weight)*avg_pdg;
    cellAxis[p.first].fromMatrix(avg_pdg);

//    if(avgType == "geometric mean" || avgType == "area weighted geometric mean"){
//            Point3f e_vals = cellAxis[p.first].evals();
//            for(int i=0;i<3;i++)
//                cellAxis[p.first].evals()[i] = pow(fabs(e_vals[i]),1.0f/weight);
//    }
    //HANDLE THE AVERAGING
  }
    //REDISPLAY GROWTH

    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(AveragePDGs);

}
