//
// This file is part of 3DCellAtlas.
// Copyright (C) 2015-2016 George W. Bassel and collaborators.
//
// If you use 3DCellAtlas in your work, please cite:
//   http://dx.doi.org/10.1105/tpc.15.00175
//
// 3DCellAtlas is an AddOn for MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2012 Richard S. Smith and collaborators.
//
// 3DCellAtlas and MorphoGraphX are free software, and are licensed under under the terms of the
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
//
#include <MeshProcessCellAnalysis3D.hpp>

#include <GraphUtils.hpp>
#include <MeshProcessLineage.hpp> // SetParent process

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

using namespace std;

namespace mgx
{
  bool ExplodeMesh::run(Mesh* m, double factorX, double factorY, double factorZ, bool onlySelected, bool reverse, bool bundleParents)
  {
    AttrMap<int, Point3d>& explDispl = m->attributes().attrMap<int, Point3d>("Explosion Displacement");
    AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    bool useParents = m->useParents();

    //if(useParents and bundleParents) return setErrorMessage("Bundle Parents option not available when parent labels are active.");

    if(centroids.empty())
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    if(!explDispl.empty()) {
      RestoreOriginalMesh re(*this);
      re.run(m, false);
    }
    if(reverse) {
      explDispl.clear();
    }

    vvGraph& S = m->graph();

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

    // std::set<int> selectedParents;

    // forall(int l, selectedCells) {
    //   selectedParents.insert(m->parents()[l]);
    // }

    //if(useParents) selectedParents = selectedCells;

    std::map<int, Point3d> labelDisplacementMap;

    typedef std::pair<int, Point3d> IntPoint3dP;
    if(bundleParents) {
      std::map<int, Point3d> parentCentroidMap;
      std::map<int, std::set<int> > parentLabelMap;
      typedef std::pair<int, std::set<int> > IntIntSetP;
      forall(IntInt p, m->parents()) {
        parentLabelMap[p.second].insert(p.first);
      }

      forall(IntIntSetP p, parentLabelMap) {
        if(p.second.empty()) continue;
        Point3d c(0,0,0);
        forall(int l, p.second) {
          //c += ca.rcp.rootData.cellCentroids[l];
          c += centroids[l];
        }
        c /= p.second.size();
        parentCentroidMap[p.first] = c;
      }

      forall(IntInt p, m->parents()) {
        if(onlySelected and selectedCells.find(p.first) == selectedCells.end()) continue;
        labelDisplacementMap[p.first].x() = parentCentroidMap[p.second].x() * factorX;
        labelDisplacementMap[p.first].y() = parentCentroidMap[p.second].y() * factorY;
        labelDisplacementMap[p.first].z() = parentCentroidMap[p.second].z() * factorZ;

        explDispl[p.first] += labelDisplacementMap[p.first];
      }
    } else {
      forall(IntPoint3dP p, centroids) {
        if(onlySelected and selectedCells.find(p.first) == selectedCells.end()) continue;
        labelDisplacementMap[p.first].x() = p.second.x() * factorX;
        labelDisplacementMap[p.first].y() = p.second.y() * factorY;
        labelDisplacementMap[p.first].z() = p.second.z() * factorZ;
        explDispl[p.first] += labelDisplacementMap[p.first];
      }
    }


    // move vertices
    forall(const vertex& v, S) {
      int l = v->label;
      //if(useParents) l = m->parents()[v->label];
      if(l < 1) {
        forall(const vertex& n, S.neighbors(v)) {
          int nl = n->label;
          //if(useParents) nl = m->parents()[n->label];
          if(nl > 0) {
            l = nl;
            break;
          }
          forall(const vertex& mv, S.neighbors(n)) {
          int ml = mv->label;
          //if(useParents) ml = m->parents()[mv->label];
            if(ml > 0) {
              l = ml;
              break;
            }
          }
        }
      }
      v->pos += explDispl[l];
    }

    // also for cell axis
    IntPoint3fAttr& centers = (m->useParents() ? m->parentCenterVis() : m->labelCenterVis());

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

    //m->setAxisView("Cell Axis");
    //m->updateTriangles();
    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(ExplodeMesh);

   bool TranslateByParent::run(Mesh* m, double factorX, double factorY, double factorZ, bool onlySelected, bool reverse, bool rot, Matrix4d rotMatrix)
  {

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

    std::cout << "here1 " << std::endl;

    if(centroids.empty())
      return setErrorMessage("No Attribute Map found, please run 'Heat Map/Analysis/Cell Analysis 3D' first.");

    if(!explDispl.empty()) {
      RestoreOriginalMesh re(*this);
      re.run(m, false);
    }
    if(reverse) {
      explDispl.clear();
    }

    bool useParents = m->useParents();

    vvGraph& S = m->graph();

    std::set<int> selectedCells = findAllSelectedLabels(S);
    std::map<int, int> parentTranslationFactorMap;

    std::map<int, Point3d> labelDisplacementMap;

    std::set<int> parentLabels;

    int counter = 1;

    forall(IntInt p, m->parents()) {
      if(p.second < 1) continue;
      if(parentTranslationFactorMap.find(p.second) != parentTranslationFactorMap.end()) continue;
      parentTranslationFactorMap[p.second] = 0;
      parentLabels.insert(p.second);
    }

    forall(int p, parentLabels) {
      parentTranslationFactorMap[p] = counter;
      //std::cout << "l " << p << "/" << counter << std::endl;
      counter++;
    }
    //std::cout << "fff " << factorX << "/" << factorY << "/" << factorZ << std::endl;

    typedef std::pair<int, Point3d> IntPoint3dP;
    forall(IntPoint3dP p, centroids) {
      if(onlySelected and selectedCells.find(p.first) == selectedCells.end()) continue;
      labelDisplacementMap[p.first].x() = factorX * parentTranslationFactorMap[m->parents()[p.first]];
      labelDisplacementMap[p.first].y() = factorY * parentTranslationFactorMap[m->parents()[p.first]];
      labelDisplacementMap[p.first].z() = factorZ * parentTranslationFactorMap[m->parents()[p.first]];
      explDispl[p.first] += labelDisplacementMap[p.first];
      if(rot) explDispl[p.first] += multMatrix4Point3(rotMatrix,labelDisplacementMap[p.first]);
    }

// move vertices
    forall(const vertex& v, S) {
      int l = v->label;
      //if(useParents) l = m->parents()[v->label];
      if(l < 1) {
        forall(const vertex& n, S.neighbors(v)) {
          int nl = n->label;
          //if(useParents) nl = m->parents()[n->label];
          if(nl > 0) {
            l = nl;
            break;
          }
          forall(const vertex& mv, S.neighbors(n)) {
          int ml = mv->label;
          //if(useParents) ml = m->parents()[mv->label];
            if(ml > 0) {
              l = ml;
              break;
            }
          }
        }
      }
      v->pos += explDispl[l];
    }

    // also for cell axis
    IntPoint3fAttr& centers = (m->useParents() ? m->parentCenterVis() : m->labelCenterVis());

    forall(auto& p, centers) {
      int l = p.first;
      if(useParents) l = m->parents()[l];
      p.second += Point3f(explDispl[l]);
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(TranslateByParent);


   bool RestoreOriginalMesh::run(Mesh* m, bool vis)
  {

    AttrMap<int, Point3d>& explDispl = m->attributes().attrMap<int, Point3d>("Explosion Displacement");

    if(explDispl.empty())
      return setErrorMessage("No explosion Attr Map present. Can't do anything.");

    bool useParents = m->useParents();

    vvGraph& S = m->graph();

    // forall(const vertex& v, S) {
    //   int l = v->label;

    //   v->pos -= explDispl[l];
    // }

    // move vertices
    forall(const vertex& v, S) {
      int l = v->label;
      //if(useParents) l = m->parents()[v->label];
      if(l < 1) {
        forall(const vertex& n, S.neighbors(v)) {
          int nl = n->label;
          //if(useParents) nl = m->parents()[n->label];
          if(nl > 0) {
            l = nl;
            break;
          }
          forall(const vertex& mv, S.neighbors(n)) {
          int ml = mv->label;
          //if(useParents) ml = m->parents()[mv->label];
            if(ml > 0) {
              l = ml;
              break;
            }
          }
        }
      }
      v->pos -= explDispl[l];
    }

    // also for cell axis
    IntPoint3fAttr& centers = (m->useParents() ? m->parentCenterVis() : m->labelCenterVis());

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

    if(vis) {
      explDispl.clear();
      m->updateAll();
    }

    return true;
  }
  REGISTER_PROCESS(RestoreOriginalMesh);


   bool SetParentExploded::run(Mesh* m, int label)
  {

    RestoreOriginalMesh re(*this);
    re.run(m, false);

    vvGraph& S = m->graph();
    std::set<int> selectedLabels = findAllSelectedLabels(m->graph());

    // update the position of the selected vertices
    AttrMap<int, Point3d>& explDispl = m->attributes().attrMap<int, Point3d>("Explosion Displacement");

    std::map<int, Point3d> parentDispl;
    forall(auto p, explDispl) {
      int pLabel = m->parents()[p.first];
      if(parentDispl.find(pLabel) == parentDispl.end()) parentDispl[pLabel] = p.second;
      else {
        if(parentDispl[pLabel] != p.second) std::cout << "WARNING: Not a parent based explosion, locations could be wrong." << std::endl;
      }
    }

    // call Set Parent
    SetParent sp(*this);
    sp.run(m, label);


    if(explDispl.empty())
      return true;

    forall(int l, selectedLabels) {
      std::cout << "label " << l << "/ " << label << "/" << explDispl[l] << "/" << parentDispl[label] << std::endl;
      explDispl[l] = parentDispl[label];
    }

    // move vertices
    forall(const vertex& v, S) {
      int l = v->label;
      //if(useParents) l = m->parents()[v->label];
      if(l < 1) {
        forall(const vertex& n, S.neighbors(v)) {
          int nl = n->label;
          //if(useParents) nl = m->parents()[n->label];
          if(nl > 0) {
            l = nl;
            break;
          }
          forall(const vertex& mv, S.neighbors(n)) {
          int ml = mv->label;
          //if(useParents) ml = m->parents()[mv->label];
            if(ml > 0) {
              l = ml;
              break;
            }
          }
        }
      }
      v->pos += explDispl[l];
    }

    m->updateAll();
    return true;
  }
  REGISTER_PROCESS(SetParentExploded);


  bool LabelCellWalls::run(Mesh *mesh1, int extend, bool onlyAttr)
  {
    vvGraph& S = mesh1->graph();

    AttrMap<vertex, int>& vertexCellWallLabel = mesh1->attributes().attrMap<vertex, int>("Cell Wall Labels");
    vertexCellWallLabel.clear();

    AttrMap<vertex, int>& vertexCellLabel = mesh1->attributes().attrMap<vertex, int>("Cell Labels");
    vertexCellLabel.clear();

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

    //mesh1->parents().clear();

    // write labels to parents

    //typedef std::set<vertex> tri;

    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;

    // go through all tris, check what cells they share
    // vtx -> tri
    // e.g. tri t1 -> cell 2 and 3

    // build vertex IntInt map

    // label vtx according to their entry in the map

    AttrMap<vertex, int>& vtxCounter = mesh1->attributes().attrMap<vertex, int>("Vertex Cell Counter");
    vtxCounter.clear();

    AttrMap<int, IntInt>& wallLabelCellLabelsMap = mesh1->attributes().attrMap<int, IntInt>("Wall Label Cell Labels");
    wallLabelCellLabelsMap.clear();

    std::set<int> outerLabels;

    typedef std::pair<IntInt, int> IntIntPIntP;
    forall(IntIntPIntP p, info.allWalls) {
      wallLabelCellLabelsMap[p.second] = p.first;
    }

    AttrMap<int, int>& neighbor1Map = mesh1->attributes().attrMap<int, int>("Measure Label Int NeighborCell1");
    AttrMap<int, int>& neighbor2Map = mesh1->attributes().attrMap<int, int>("Measure Label Int NeighborCell2");
    neighbor1Map.clear();
    neighbor2Map.clear();

    forall(auto p, wallLabelCellLabelsMap) {
      neighbor1Map[p.first] = p.second.first;
      neighbor2Map[p.first] = p.second.second;
    }

    forall(const vertex& v, S) {
      int oldLabel = v->label;
      if(!onlyAttr) v->label = vtxWallMap[v];

      vertexCellWallLabel[v] = vtxWallMap[v];
      if(!onlyAttr) mesh1->parents()[v->label] = oldLabel;
      vtxCounter[v] = vtxCellCounterMap[v];
      if(vtxCellCounterMap[v] == 1 and v->label>0) outerLabels.insert(v->label);
      //if(v->label == 0) { cout << "zero label2 " << vtxWallMap[v] << endl;
      //v->selected = true;
      //}
    }

    for(int i =0; i< extend; i++) {
      std::map<vertex, int> newLabel;
      typedef std::pair<vertex, int> VtxIntPair;

      mesh1->updateAll();

      if(extend) {
        forall(const vertex& v, S) {
          if(v->label > -1) continue;
          std::set<int> neighborLabels;
          forall(const vertex& n, S.neighbors(v)) {
            if(outerLabels.find(n->label) == outerLabels.end() and n->label > 0) neighborLabels.insert(n->label);
          }
          if(neighborLabels.size() == 1) newLabel[v] = *(neighborLabels.begin());
        }
      }
      forall(VtxIntPair p, newLabel) {
          if(!onlyAttr)
            p.first->label = p.second;
          vertexCellWallLabel[p.first] = p.second;
      }
    }

    mesh1->parents().erase(-1);

    mesh1->updateAll();
    return true;
  }
  REGISTER_PROCESS(LabelCellWalls);


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

    AttrMap<vertex, int>& vertexCellLabel = mesh1->attributes().attrMap<vertex, int>("Cell Labels");

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

    mesh1->updateAll();
    return true;
  }
  REGISTER_PROCESS(RestoreCellLabels);


  bool FindCellWallCenters::run(Mesh *mesh1, Mesh *mesh2, QString mode, bool writeMesh)
  {
    vvGraph& S = mesh1->graph();
    vvGraph& S2 = mesh2->graph();

    double epsilon = 1E-8;

    AttrMap<vertex, int>& vertexCellWallLabel = mesh1->attributes().attrMap<vertex, int>("Cell Wall Labels");
    if(vertexCellWallLabel.empty()) {
      return setErrorMessage("Run the Label Cell Walls process first!");
    }

    AttrMap<IntInt, Point3d>& wallCenters = mesh1->attributes().attrMap<IntInt, Point3d>("Cell Wall Centers");
    wallCenters.clear();
    AttrMap<int, IntInt>& wallLabelCellLabelsMap = mesh1->attributes().attrMap<int, IntInt>("Wall Label Cell Labels");

    // create label vertex map
    std::map<int, std::set<vertex> > labelVertexMap;
    forall(const vertex& v, S) {
      int l = vertexCellWallLabel[v];
      if(l < 1) continue;
      labelVertexMap[l].insert(v);
    }

    if(mode == "Average Tri Weighted") {

      forall(auto p, labelVertexMap) {
        int label = p.first;
        std::set<vertex> vtxs = p.second;
        double totalArea = 0;
        Point3d pos(0,0,0);

        forall(vertex v, vtxs) {
          forall(vertex n, S.neighbors(v)) {
            vertex m = S.nextTo(v,n);
            if(!S.uniqueTri(v,n,m)) continue;
            double tArea = triangleArea(v->pos, n->pos, m->pos);
            totalArea += tArea;
            pos += tArea * (v->pos + n->pos + m->pos)/3.;
          }
        }
        if(totalArea > 0)
          pos /= totalArea;

        vertex v;
        v->pos = pos;

        if(writeMesh) S2.insert(v);
        wallCenters[wallLabelCellLabelsMap[label]] = pos;
      }

    } else {

      forall(auto p, labelVertexMap) {
        int label = p.first;
        std::set<vertex> vtxs = p.second;

        std::map<double, vertex> distanceVtxMap;
        std::unordered_map<vertex, double> vtxDistanceMap;
        std::set<vertex> finishedVtxs;

        forall(vertex v, vtxs) {
          double minDis = 1E20;
          forall(vertex n, S.neighbors(v)) {
            if(n->label == -1) {
              double dis = norm(v->pos - n->pos);
              if(dis<minDis) minDis = dis;
            }
          }

          if(minDis < 1E20) {
            while(distanceVtxMap.find(minDis) != distanceVtxMap.end()) {
              minDis += epsilon;
            }
            distanceVtxMap[minDis] = v;
            vtxDistanceMap[v] = minDis;
            finishedVtxs.insert(v);
          }
        }

        vertex vLast;
        double disLast;
        while(!distanceVtxMap.empty()) {

          // find the lowest value in map
          auto p = distanceVtxMap.begin();
          double disVtx = p->first;
          vertex v = p->second;

          // go to all neighbors of vertex and replace distance with lowest possible value
          forall(const vertex& n, S.neighbors(v)) {
            if(finishedVtxs.find(n) != finishedVtxs.end() or n->label != v->label) continue;
            double newDis = disVtx + norm(v->pos - n->pos);
            if(vtxDistanceMap.find(n) == vtxDistanceMap.end()) { // new vertex: create entry
              while(distanceVtxMap.find(newDis) != distanceVtxMap.end()) {
                newDis += epsilon;
              }
              distanceVtxMap[newDis] = n;
              vtxDistanceMap[n] = newDis;
            } else { // existing vertex: update distance
              if(newDis < vtxDistanceMap[n]) {
                distanceVtxMap.erase(vtxDistanceMap[n]);


                distanceVtxMap[newDis] = n;
                vtxDistanceMap[n] = newDis;
              }
            }

          }
          // erase the vertex
          vLast = v;
          disLast = disVtx;
          finishedVtxs.insert(v);
          distanceVtxMap.erase(disVtx);
          vtxDistanceMap.erase(v);
        }

        vertex v;
        v->pos = vLast->pos;

        if(writeMesh) S2.insert(v);
        wallCenters[wallLabelCellLabelsMap[label]] = vLast->pos;
      }

    }

    mesh2->updateAll();
    return true;
  }
  REGISTER_PROCESS(FindCellWallCenters);


  bool RotateCamera::run(Stack *s1, Mesh *m1, Stack *s2, Mesh *m2, double xRot, double yRot, double zRot, int steps, int breakStep, bool cameraStack2)
  {

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

    double xRotStep = xRot * M_PI / 180. / (double)steps;
    double yRotStep = yRot * M_PI / 180. / (double)steps;
    double zRotStep = zRot * M_PI / 180. / (double)steps;

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

    rotMatrixStep = Matrix4d::identity();

    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 = rotMatrixS1;
    updatedMat2 = rotMatrixS2;

    for(int i = 0; i < steps; i++) {
      if(!progressAdvance()) {
        s1->getFrame().setFromMatrix(rotMatrixS1.data());
        m1->updateAll();
        s2->getFrame().setFromMatrix(rotMatrixS2.data());
        m2->updateAll();
        updateState();
        updateViewer();
        return false;
      }
      auto startTime = std::chrono::system_clock::now();

      updatedMat = updatedMat * rotMatrixStep;
      if(cameraStack2) updatedMat2 = updatedMat2 * rotMatrixStep;

      s1->getFrame().setFromMatrix(updatedMat.data());
      s2->getFrame().setFromMatrix(updatedMat2.data());

      updateState();
      updateViewer();

      sleep_until(startTime + nanoseconds(breakStep*1000));

    }


    return true;
  }
  REGISTER_PROCESS(RotateCamera);


  bool LabelNuclei::run(Mesh* m, Mesh* m2)
  {
    AttrMap<int, Point3d>& centroidsCells = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    AttrMap<int, Point3d>& centroidsNuclei = m2->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    if(centroidsCells.empty()) return setErrorMessage("No Centroids in active mesh");
    if(centroidsNuclei.empty()) return setErrorMessage("No Centroids in other mesh");

    // get label trinagle map
    vvGraph& S1 = m->graph();
    vvGraph& S2 = m2->graph();

    std::unordered_map<int, std::vector<Triangle> > labelTriangleMap;
    forall(vertex v, S1) {
      forall(vertex n, S1.neighbors(v)) {
        vertex m = S1.nextTo(v,n);
        if(!S1.uniqueTri(v,n,m)) continue;
        Triangle t(v,n,m);
        labelTriangleMap[v->label].push_back(t);
      }
    }

    Point3d dirVec(1,0,0);
    Point3d origin(0,0,0);

    std::map<int, int> newNucleusLabels;

    std::vector<int> vecCentroidsNucleiLabel;
    std::vector<Point3d> vecCentroidsNucleiPos;

    forall(auto p, centroidsNuclei) {
      vecCentroidsNucleiLabel.push_back(p.first);
      vecCentroidsNucleiPos.push_back(p.second);
    }

    // shoot ray from nucleus centroid, check how many intersections
    #pragma omp parallel for
    for(uint i = 0; i < vecCentroidsNucleiLabel.size(); i++) {
    //forall(auto p, centroidsNuclei) {

      int label = vecCentroidsNucleiLabel[i]; //p.first;
      Point3d pos = vecCentroidsNucleiPos[i]; //p.second;

      std::map<int, int> labelIntersections;
      std::set<int> insideCell;

      //std::cout << "nucleus " << label << "/" << pos << std::endl;

      if(norm(pos - origin) < 1E-5) continue;


      forall(auto p2, centroidsCells) {
        if(norm(pos - origin) < 1E-5) continue;
        Point3d coord;
        int stopCount = 0;
        findIntersectPoint(pos, dirVec, labelTriangleMap[p2.first], stopCount, coord);

        if(stopCount % 2 == 1) {
          insideCell.insert(p2.first);
        }
      }

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

      if(insideCell.size() == 1) {
        newNucleusLabels[label] = *(insideCell.begin());
      } else if(insideCell.size() == 0) {
        newNucleusLabels[label] = 0;
      } else{
        std::cout << "more than a single cell found for label " << label << std::endl;
      }
    }

    m2->parents().clear();

    forall(vertex v, S2) {
      m2->parents()[v->label] = newNucleusLabels[v->label];
    }

    m2->updateAll();
    return true;

  }
  REGISTER_PROCESS(LabelNuclei);


  bool DistanceNuclei::run(Mesh* m, Mesh* m2)
  {
    AttrMap<int, Point3d>& centroidsCells = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    AttrMap<int, Point3d>& centroidsNuclei = m2->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    if(m2->useParents())
      centroidsNuclei = m2->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");

    AttrMap<int, double>& distanceNuclei = m->attributes().attrMap<int, double>("Measure Label Double DistanceNuclei");
    distanceNuclei.clear();

    AttrMap<int, Point3d>& distanceCustom = m->attributes().attrMap<int, Point3d>("Measure Label Vector DistanceCustom");
    distanceCustom.clear();

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

    if(centroidsCells.empty()) return setErrorMessage("No Centroids in active mesh");
    if(centroidsNuclei.empty()) return setErrorMessage("No Centroids in other mesh");


    std::map<int, Point3d> labelNucleiDirMap;

    forall(auto p, centroidsNuclei) {
      int label = p.first;//m2->parents()[p.first];
      if(label < 1) continue;
      if(centroidsCells.find(label) == centroidsCells.end()) continue;

      labelNucleiDirMap[label] = p.second - centroidsCells[label];

      m->labelHeat()[label] = norm(labelNucleiDirMap[label]);

      if(customDirs.find(label) == customDirs.end()) continue;

      Point3d customDistance(0,0,0);
      customDistance.x() = (labelNucleiDirMap[label] * Point3d(customDirs[label].ev1()));
      customDistance.x() /= norm(customDirs[label].ev1());
      customDistance.y() = (labelNucleiDirMap[label] * Point3d(customDirs[label].ev2()));
      customDistance.y() /= norm(customDirs[label].ev2());
      customDistance.z() = (labelNucleiDirMap[label] * Point3d(customDirs[label].evals()));
      customDistance.z() /= norm(customDirs[label].evals());

      std::cout << "dis " << label << "/" << norm(labelNucleiDirMap[label]) << "/" << customDistance << std::endl;

      distanceCustom[label] = customDistance;
      distanceNuclei[label] = norm(labelNucleiDirMap[label]);
    }

    m->updateAll();
    return true;
  }
  REGISTER_PROCESS(DistanceNuclei);


  void labelsInRange(double coordPos, double threshold, std::set<int>& labelSet, std::map<double, int> coordToLabel)
  {
    std::map<double, int>::iterator itlow, itup;

    itlow = coordToLabel.lower_bound(coordPos - threshold);
    itup = coordToLabel.upper_bound(coordPos + threshold);

    std::set<int> labelsInRange, newLabelsSet;

    for(std::map<double, int>::iterator it = itlow; it != itup; it++) {
      if(norm((*it).first - coordPos) < threshold)
        labelsInRange.insert((*it).second);
    }

    forall(int l, labelsInRange) {
      if(labelSet.find(l) == labelSet.end()) continue;
      newLabelsSet.insert(l);
    }
    labelSet = newLabelsSet;
  }


  bool SelectDuplicatedNuclei::run(Mesh* m, double disThreshold)
  {

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

    std::map<double, int> xToLabel, yToLabel, zToLabel;

    std::set<int> labelsToSelect;

    std::set<std::set<int> > clusters;

    // build position maps
    forall(auto p, centroidsNuclei) {
      int l = p.first;
      Point3d pos = p.second;

      double xPos = pos.x();
      double yPos = pos.y();
      double zPos = pos.z();

      while(xToLabel.find(xPos) != xToLabel.end()) {
        xPos += 1E-8;
      }
      while(yToLabel.find(yPos) != yToLabel.end()) {
        yPos += 1E-8;
      }
      while(zToLabel.find(zPos) != zToLabel.end()) {
        zPos += 1E-8;
      }
      xToLabel[xPos] = l;
      yToLabel[yPos] = l;
      zToLabel[zPos] = l;
    }

    std::cout << "maps built" << std::endl;

    std::vector<int> vecCentroidsNucleiLabel;
    std::vector<Point3d> vecCentroidsNucleiPos;
    std::set<int> labelsAll;

    forall(auto p, centroidsNuclei) {
      vecCentroidsNucleiLabel.push_back(p.first);
      vecCentroidsNucleiPos.push_back(p.second);
      //labelsAll.insert(p.first);
    }

    // check for close centroids
    //#pragma omp parallel for
    for(uint i = 0; i < vecCentroidsNucleiLabel.size(); i++) {

      int l = vecCentroidsNucleiLabel[i]; //p.first;
      Point3d pos = vecCentroidsNucleiPos[i]; //p.second;

      forall(auto p, centroidsNuclei) {
        labelsAll.insert(p.first);
      }
      labelsAll.erase(l);

      //std::cout << "loop1 " << l << "/" << labelsAll.size() <<std::endl;
      labelsInRange(pos.x(), disThreshold, labelsAll, xToLabel);
      //std::cout << "loop2 " << l << "/" << labelsAll.size() <<std::endl;
      labelsInRange(pos.y(), disThreshold, labelsAll, yToLabel);
      //std::cout << "loop3 " << l << "/" << labelsAll.size() <<std::endl;
      labelsInRange(pos.z(), disThreshold, labelsAll, zToLabel);


      labelsAll.insert(l);

      #pragma omp critical
      clusters.insert(labelsAll);
      forall(int l, labelsAll) {
        labelsToSelect.insert(l);
      }

    }

    std::cout << "labels found" << std::endl;

    labelsToSelect.clear();

    forall(std::set<int> c, clusters) {
    //  std::cout << "c " << c.size() << std::endl;
      for(std::set<int>::iterator it = c.begin(); it != c.end(); it++) {
        if(it != c.begin())
          labelsToSelect.insert((*it));
      }
    }

    forall(vertex v, m->graph()) {
      v->selected = false;
      if(labelsToSelect.find(v->label) == labelsToSelect.end()) continue;
      v->selected = true;
    }

    m->updateAll();
    return true;
  }

  REGISTER_PROCESS(SelectDuplicatedNuclei);
}
