#include "DivisionAnalysis.hpp"
#include "Information.hpp"
#include "Progress.hpp"

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

#include <CuttingSurface.hpp>

#include <QFile>
#include <QTextStream>
#include <QFileDialog>

#include <Triangulate.hpp>
#include <triangle.h>

#include <MeshBuilder.hpp>

#include <DivisionAnalysisUtils.hpp>

#include <GraphAlgorithms.hpp>
#include <GraphUtils.hpp>

#include <KrylovMethods.hpp>
#include <DistMatrix.hpp>
#include <DistObject.hpp>

#include <MeshProcessStructure.hpp>
#include <MeshProcessMulti.hpp>

#include <CellTissue.hpp>

#include <ctime>

using namespace std;

namespace mgx { //namespace process
//{

  CellDivisionAttr convertToCDA(CellDivision& cd){
    CellDivisionAttr cda;

    cda.divPlane.nrml = cd.divPlane.nrml;
    cda.divPlane.pos = cd.divPlane.pos;
    cda.divPlane.displ = cd.divPlane.displ;
    cda.volumeRatio = cd.volumeRatio;
    cda.divPlane.area = cd.divPlane.area;

    return cda;
  }

 double planeSimiliarity(double planePos1, Point3d planeNrml1, double planePos2, Point3d planeNrml2)
 {
  double scalarNrmls = planeNrml1 * planeNrml2;
  if((scalarNrmls>0 and (planePos1 == planePos2)) or (scalarNrmls<=0 and (planePos1 == -planePos2))){
    return norm(scalarNrmls);
  } else {
    return norm(scalarNrmls);
  }
 }


// // stuff from cell tissue
//   void extractTriangles(vvGraph& S, std::vector<CDiv::Triangle>& triangles)
//   {

//     //const auto& c_faces = C.bounds(c);
//     triangles.clear();

//     forall(const vertex& v, S){
//       forall(const vertex& n, S.neighbors(v)){
//         vertex m = S.nextTo(v,n);
//         if(S.uniqueTri(v,n,m)){
//           CDiv::Point p1 (v->pos.x(),v->pos.y(),v->pos.z());
//           CDiv::Point p2 (n->pos.x(),n->pos.y(),n->pos.z());
//           CDiv::Point p3 (m->pos.x(),m->pos.y(),m->pos.z());
//           CDiv::Triangle tr(p1, p2, p3);
//           triangles.push_back(tr);
//         }
//       }
//     }

//   }
/*
  void extractTriangles(vvGraph& S, std::vector<tri>& trianglesMGX)
  {

    //const auto& c_faces = C.bounds(c);
    trianglesMGX.clear();

    forall(const vertex& v, S){
      forall(const vertex& n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
        if(S.uniqueTri(v,n,m)){
          tri t;
          t[0] = v;
          t[1] = n;
          t[2] = m;
          trianglesMGX.push_back(t);
        }
      }
    }

  }*/
/*
  void copyTriangles(std::vector<tri>& trianglesMGX, std::vector<CDiv::Triangle>& triangles)
  {

    forall(const tri& t, trianglesMGX){
      CDiv::Point p1 (t[0]->pos.x(),t[0]->pos.y(),t[0]->pos.z());
      CDiv::Point p2 (t[1]->pos.x(),t[1]->pos.y(),t[1]->pos.z());
      CDiv::Point p3 (t[2]->pos.x(),t[2]->pos.y(),t[2]->pos.z());
      CDiv::Triangle tr(p1, p2, p3);
      triangles.push_back(tr);
    }

  }
*/
/*    void addTriCheckOrient(std::map<Point3d, int>& vMap, const vertex &a, const vertex &b, const vertex &c,
                                           std::vector<vertex> &vtxs, std::vector<Point3i> &tris, Point3d &nrml)
    {
      Point3d triNrml = (b->pos-a->pos) % (c->pos-a->pos);
      if(triNrml * nrml > 0){
        addTri(vMap,a,b,c,vtxs,tris);
      } else {
        addTri(vMap,b,a,c,vtxs,tris);
      }

    }*/
/*
  void trianglesFromLabelledVertices(vvGraph& S, std::vector<vertex>& vertices, int label, std::vector<Point3i>& tris, std::vector<Point3d> vtxsPos)
  {

    // Map for caching vertex positions.
    std::map<Point3d, int> pMap;

    // get all unique tris
    forall(const vertex &v, vertices){
      forall(const vertex &n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
        if(S.uniqueTri(v,n,m) and (v->label == label and n->label == label and m->label == label)){

          Point3i currentTri(0,0,0);

          if(pMap[v->pos] <= 0){
            pMap[v->pos] = tris.size()+1;
          }
          currentTri.x() = pMap[v->pos];


        }
      }
    }

  }
*/
  //   std::vector<CDiv::Object_and_primitive_id> objPrimCGAL(const CDiv::Tree& tree, Point3d point, Point3d normal)
  //   {

  //   CDiv::Point cp(point.x(), point.y(), point.z()); // cell center
  //   CDiv::Vector d(normal.x(), normal.y(), normal.z()); // cutting plane

  //   // Construct the planar intersection query
  //   CDiv::Plane plane_query(cp,d);

  //   std::vector<CDiv::Object_and_primitive_id> intersections;
  //   tree.all_intersections(plane_query, std::back_inserter(intersections));
  //   return intersections;
  //   }

  //   bool triPlaneIntersect()
  //   {

  //   }


  //   std::vector<std::pair<Point3d,Point3d> > findDivisionPlane(std::vector<Triangle> tris, P3dDouble plane,
  //                                                              std::map<VtxVtx, vertex>& segmentsToBeSplit)
  // {
  //   std::vector<std::pair<Point3d,Point3d> > res;



  //   return res;
  // }


  //   std::vector<std::pair<Point3d,Point3d> > findDivisionPlane(const vvGraph& cS, std::vector<CDiv::Object_and_primitive_id> intersections,
  //                                                              std::map<VtxVtx, vertex>& segmentsToBeSplit, bool calcSegs)
  // {

  //   std::map<VtxVtx, vertex> temp;
  //   std::vector<std::pair<Point3d,Point3d> > result;

  //   //CDiv::Point cp(cellCenter.x(), cellCenter.y(), cellCenter.z()); // cell center
  //   //CDiv::Vector d(dir.x(), dir.y(), dir.z()); // cutting plane

  //   // Construct the planar intersection query
  //   //CDiv::Plane plane_query(cp,d);

  //   //#pragma omp critical
  //   //tree.all_intersections(plane_query, std::back_inserter(intersections));

  //   CDiv::Segment segment;
  //   CDiv::Point point;
  //   //DivisionLines div;
  //   forall(const CDiv::Object_and_primitive_id& obj, intersections)
  //   {
  //     //auto edges = C.bounds(f);
  //     if(CGAL::assign(segment, obj.first)){

  //       Point3d p1(segment.source());
  //       Point3d p2(segment.target());

  //       if(calcSegs){
  //       forall(const vertex& v, cS){
  //         forall(const vertex& n, cS.neighbors(v)){
  //           if(distLinePoint(v->pos, n->pos, p1, true) < 0.0001){
  //             VtxVtx vv;
  //             if(v<n) vv = std::make_pair(v,n);
  //             else vv = std::make_pair(n,v);
  //             vertex vp1;
  //             vp1->pos = p1;
  //             temp[vv] = vp1;
  //           }
  //           if(distLinePoint(v->pos, n->pos, p2, true) < 0.0001){
  //             VtxVtx vv;
  //             if(v<n) vv = std::make_pair(v,n);
  //             else vv = std::make_pair(n,v);
  //             vertex vp2;
  //             vp2->pos = p2;
  //             temp[vv] = vp2;
  //           }
  //           //CDiv::Point sp1(p1.x(), p1.y(), p1.z());
  //           //CDiv::Point sp2(p2.x(), p2.y(), p2.z());
  //           //CDiv::Segment seg(sp1, sp2);
  //           //if(CGAL::do_intersect(seg, plane_query)){

  //           //}
  //         }
  //       }
  //       }
  //       // go through all edges, find

  //       //cout << "p segment " << p1 << "/" << p2 << std::endl;
  //       std::pair<Point3d,Point3d> p = std::make_pair(p1,p2);
  //       result.push_back(p);
  //       //result.push_back(p2);

  //     } else if(CGAL::assign(point, obj.first)){
  //       Point3d p(point);
  //       //cout << "p single  " << p << endl;
  //       //result.push_back(p);
  //       //cout << "############################################################" << endl;
  //       //cout << "WARNING ####@@@@@@@@@@#### -- intersection is a single point" << endl;
  //       //cout << "############################################################" << endl;
  //     }
  //   }
  //   segmentsToBeSplit = temp;
  //   return result;
  // }

  // // finds the smallest division plane (custom point)
  // P3dDouble findSmallestDivisionPlane(const cell &c, std::vector<P3dDouble> &samplingDirections, std::map<P3dDouble, double>& resultMap){

  //   //std::cout << "find " << std::endl;
  //   // cell center
  //   Point3d cellCenter = c->pos;

  //   // find the division plane
  //   Point3d divPlaneNrml (0,1,0);

  //   // prepare CGAL for getting the line segments
  //   std::vector<CDiv::Triangle> triangles;
  //   extractTriangles(c->S, triangles);
  //   CDiv::Tree tree(triangles.begin(), triangles.end());

  //   std::vector<std::pair<Point3d,Point3d> > polygonSegs;
  //   std::vector<Point3d> pointsPolygon;
  //   std::map<VtxVtx, vertex> segmentsToBeSplit;

  //   P3dDouble minSampleDir;
  //   std::map<P3dDouble, double> planeAreaMap;

  //   std::map<P3dDouble, std::vector<CDiv::Object_and_primitive_id> > objPrimCGALMap;

  //   //progressStart("Test Division Planes - Find Smallest/Shortest Division Plane (CGAL object)", 0);
  //   forall(P3dDouble p, samplingDirections){
  //     objPrimCGALMap[p] = objPrimCGAL(tree, cellCenter+p.second*p.first, p.first);
  //   }

  //   double minArea = HUGE_VAL;
  //   //progressStart("Test Division Planes - Find Smallest/Shortest Division Plane (calc areas)", 0);
  //   #pragma omp parallel for
  //   for(size_t i=0; i<samplingDirections.size(); i++){//const P3dDouble& p, samplingDirections){

  //     const P3dDouble& p = samplingDirections[i];
  //     std::vector<std::pair<Point3d,Point3d> > polygonSegsTemp = findDivisionPlane(c->S, objPrimCGALMap[p], segmentsToBeSplit,false);
  //     if(polygonSegsTemp.size() > 0){
  //       std::vector<Point3d> pointsPolygonTemp = orderPolygonSegs(polygonSegsTemp);
  //       planeAreaMap[p] = divisionPlaneArea(cellCenter+p.second*p.first, pointsPolygonTemp);
  //     }
  //     if(planeAreaMap[p]<minArea and planeAreaMap[p]>0){
  //      minArea = planeAreaMap[p];
  //       minSampleDir = p;
  //     }
  //   }

  //   resultMap = planeAreaMap;
  //   return minSampleDir;
  // }

// end stuff from cell tissue
  /*
  void extractTriangles(vvGraph& S, std::vector<vertex>& cellVertices, std::vector<CDiv::Triangle>& triangles)
  {

    //const auto& c_faces = C.bounds(c);
    triangles.clear();

    forall(const vertex& v, cellVertices){
      forall(const vertex& n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
        if(S.uniqueTri(v,n,m)){
          CDiv::Point p1 (v->pos.x(),v->pos.y(),v->pos.z());// = CDiv::point(v->pos);
          CDiv::Point p2 (n->pos.x(),n->pos.y(),n->pos.z());// = CDiv::point(n->pos);
          CDiv::Point p3 (m->pos.x(),m->pos.y(),m->pos.z());// = CDiv::point(m->pos);
          CDiv::Triangle tr(p1, p2, p3);
          triangles.push_back(tr);
        }
      }
    }
  }
*/
/*
 // uses CGAL to find the (unordered) border vertices of the division plane
 // planeP, dir: point and normal of division plane
 std::vector<Point3d> findDivisionPlane(Point3d planeP, const Point3d& dir, const CDiv::Tree& tree)
  {

    std::vector<Point3d> result;

    CDiv::Point cp(planeP.x(), planeP.y(), planeP.z()); // plane point (cell center or a bit displaced)
    CDiv::Vector d(dir.x(), dir.y(), dir.z()); // cutting plane

    // Construct the planar intersection query
    CDiv::Plane plane_query(cp,d);

    std::vector<CDiv::Object_and_primitive_id> intersections;
    tree.all_intersections(plane_query, std::back_inserter(intersections));

    CDiv::Segment segment;
    CDiv::Point point;
    //DivisionLines div;
    forall(const CDiv::Object_and_primitive_id& obj, intersections)
    {
      //auto edges = C.bounds(f);
      if(CGAL::assign(segment, obj.first)){

        Point3d p1(segment.source());
        Point3d p2(segment.target());

        CDiv::Point sp1(p1.x(), p1.y(), p1.z());
        CDiv::Point sp2(p2.x(), p2.y(), p2.z());
        CDiv::Segment seg(sp1, sp2);
        if(CGAL::do_intersect(seg, plane_query)){
        }

        // go through all edges, find

        //cout << "p segment " << p1 << "/" << p2 << endl;
        result.push_back(p1);
        result.push_back(p2);

      } else if(CGAL::assign(point, obj.first)){
        Point3d p(point);
        //cout << "p single  " << p << endl;
        //result.push_back(p);
        //cout << "############################################################" << endl;
        //cout << "WARNING ####@@@@@@@@@@#### -- intersection is a single point" << endl;
        //cout << "############################################################" << endl;
      }
    }
    return result;
  }

bool DivideCell3D::run(Mesh *m)
  {
    //if(_meshType != "MGX3D") return setErrorMessage("error: wrong mesh type");

    vvGraph& S = m->graph();
    std::map<int, std::vector<vertex> > labelVerticesMap;

    forall(const vertex& v, S){
      labelVerticesMap[v->label].push_back(v);
      v->selected = false;
    }

    std::vector<Point3d> samplingDirections;

    int label = 1;

    // get cell center
    Point3d cellCenter (0,0,0);

    forall(const vertex& v, labelVerticesMap[label]){
      cellCenter += v->pos;
    }
    cellCenter/=labelVerticesMap[label].size();

    forall(const vertex& v, labelVerticesMap[label]){
      samplingDirections.push_back((v->pos - cellCenter)/norm(v->pos - cellCenter));
    }

    cout << "cellCenter  " << cellCenter << endl;

    // find division plane

    // correct division plane

    // divide the cell
    Point3d divPlanePos (0,0,0);
    Point3d divPlaneNrml (0,1,0);


    //d = randomDirection();

    std::vector<CDiv::Triangle> triangles;
    //std::vector<
    //std::vector<face> faces;
    //extractTriangles(c, triangles, faces);
    extractTriangles(S, labelVerticesMap[label], triangles);
    CDiv::Tree tree(triangles.begin(), triangles.end());

    //return findDivisionPlane(c, d, 0, triangles, faces, tree);
    std::vector<Point3d> d = findDivisionPlane(cellCenter,divPlaneNrml, tree);

    double minArea = HUGE_VAL;
    Point3d minSampleDir;
    forall(const Point3d& p, samplingDirections){
      d = findDivisionPlane(cellCenter, p, tree);
      double area = divisionPlaneArea(cellCenter, d);
      cout << "area " << area << "/" << p << endl;
      if(area<minArea){
        minArea = area;
        minSampleDir = p;
      }
    }
    cout << "min  " << minArea << "/" << minSampleDir << endl;
    divPlaneNrml = minSampleDir;

    d = findDivisionPlane(cellCenter,divPlaneNrml, tree);

    std::set<vertex> daughterLeft, daughterRight;
    // which vertex is above and which below?
    forall(const vertex& v, labelVerticesMap[label]){
      double sP;
      Point3d intersectP;
      planeLineIntersect(cellCenter, divPlaneNrml, v->pos, v->pos-divPlaneNrml, sP, intersectP);
      if(abs(sP)<0.01){
        daughterLeft.insert(v);
        daughterRight.insert(v);
      } else if(sP<0){
        daughterLeft.insert(v);
      } else {
        daughterRight.insert(v);
      }
    }

    forall(const vertex& v, daughterLeft){
      v->selected = true;
    }

    m->updateAll();
    updateViewer();
    updateState();

    return true;
  }
  REGISTER_PROCESS(DivideCell3D);
  */

  bool DivideCell3DMGX3D::run(Mesh *mesh, int label, QString division, double x, double y, double z, double displacement/*bool newM*//*, QString fragmentHandling*/)
  {

    std::vector<P3dDouble> samplingDirections;
    //Subdivide* s = NULL;
    //mgx::CellTissue::DivAlg divAlg;

    cell cellToDivide;
    bool found = false;

    CellTissue &T = mesh->tissue();
    forall(const cell& c, T.C){
      if(c->label == label){
        cellToDivide = c;
        found = true;
        cout << "divide cell " << c->label << "/center " << c->pos << "/" << /*samplingDirections[0].first*/ endl;
      }
    }
    if(!found)
      return setErrorMessage("Cell doesn't exist!");

    if(division == "shortest wall"){ // shortest Wall, generate sampling directions and test them
      double sampling_res = 10;
      generateSamplingDirections(sampling_res, 1, displacement, samplingDirections);
      //divAlg = mgx::CellTissue::stringToDivAlg("SHORTEST WALL");

      std::map<P3dDouble, double> divPlaneSizeMap;

      P3dDouble divPlane = make_pair(Point3d(std::rand(), std::rand(), std::rand()), displacement);//TODO find proper smallest plane! findSmallestDivisionPlane(cellToDivide, samplingDirections, divPlaneSizeMap);
      samplingDirections.push_back(divPlane);
    } else { // random or predefined direction
      if(division == "random direction"){
        P3dDouble p = make_pair(Point3d(std::rand(), std::rand(), std::rand()), displacement);
        samplingDirections.push_back(p);
      } else if(division == "predefined direction"){
        if(x==y and x==z and x==0){ // cant divide with (0,0,0)
          z=1;
        }
        samplingDirections.push_back(make_pair(Point3d(x,y,z),displacement));
      }
        samplingDirections[0].first/=norm(samplingDirections[0].first); // normalize
        //divAlg = mgx::CellTissue::stringToDivAlg("PREDEFINED DIRECTION");
    }


    cout << "divide cell " << cellToDivide->label << "/center " << cellToDivide->pos << "/plane " << samplingDirections[0].first << endl;
    DivPlane p(samplingDirections[0].first, cellToDivide->pos + displacement*samplingDirections[0].first, displacement);
    CellDivision cd (cellToDivide, p, T.C);
    cd.divideCell3D(T);
    // this was the old method
    //if(!divideCell3D(T, cellToDivide, s, divAlg, samplingDirections)) return setErrorMessage("wrong mesh type, MGX3D required");

    cout << "Process ended, now update" << endl;
    T.updateAllGraph();
    mesh->updateAll();

    return true;
  }
  REGISTER_PROCESS(DivideCell3DMGX3D);

   int countNeighbours(cellGraph &C, const cell &c1)
   {
    int result = 0;

     forall(const cell &cn, C.neighbors(c1)){
       result++;
     }
     return result;
   }

   int countNeighbours(cellGraph &C, cell &c1, cell &c2)
   {
    int result = 0;

     forall(const cell &cn, C.neighbors(c1)){
       result++;
     }
     forall(const cell &cn, C.neighbors(c2)){
       result++;
     }
     return result;
   }


// check if the value is in the reach to be considered optimal
bool checkForExtremaQuick(const CellDivision& cd, QString direction, double value, double tolerance,
  std::vector<CellDivision>& optPlanes, double& optValue)
{
  if(direction == "Min"){
    if(optValue > value){ // new optimum
      optValue = value;
      return true;
    }
  } else if(direction == "Max"){ // in range
    if(optValue < value){ // new optimum
      optValue = value;
      return true;
    }
  }
  return false;
}


// check if the value is in the reach to be considered optimal
bool checkForExtrema(const CellDivision& cd, QString direction, double value, double tolerance,
  std::vector<CellDivision>& optPlanes, double& optValue)
{

  if(direction == "Min"){
    if(value > optValue and norm(optValue/value-1)<tolerance){ // in range
      optPlanes.push_back(cd);
    } else if(optValue > value){ // new optimum
      //std::vector<CellDivision> optPlanesTemp;
      //optPlanesTemp.push_back(cd);
      optValue = value;
      return true;
      //optPlanes = optPlanesTemp;
    }
  } else if(direction == "Max"){ // in range
    if(value < optValue and norm(value/optValue-1)<tolerance){
      optPlanes.push_back(cd);
    } else if(optValue < value){ // new optimum
      //std::vector<CellDivision> optPlanesTemp;
      //optPlanesTemp.push_back(cd);
      optValue = value;
      //optPlanes = optPlanesTemp;
      return true;
    }
  } else { // random
      //int randPlaneNr = std::rand() % optPlanes.size();
      std::vector<CellDivision> optPlanesTemp;
      optPlanesTemp.push_back(cd);
  }
  return false;
}

  void filterPlanes(std::vector<CellDivision>& inputDivisions, std::vector<CellDivision>& outputDivisions,
  double& optValue, double tolerance, P3dDouble& bestDivPlane, QString crit, QString direction, double planeSim = 1.0)
  {


    std::vector<CellDivision> primaryPlanes;
    std::vector<double> primaryPlanesValues;
    optValue = HUGE_VAL;
    if(direction == "Max") optValue = -HUGE_VAL;

    std::map<double, CellDivision> inputDivisionsSorted;

    // first go through all divisions and find the ones with the best betweness measure
    forall(CellDivision& cd, inputDivisions){
      double value = cd.getValue(crit);
      checkForExtremaQuick(cd, direction, value, tolerance, primaryPlanes, optValue);
      inputDivisionsSorted[value] = cd;
    }

    if(direction == "Min"){
      for(auto it = inputDivisionsSorted.begin(); it != inputDivisionsSorted.end(); it++){
        CellDivision cd = it->second;
        //double value = cd.getValue(primary);
        // in range of tolerance?
        if(norm(optValue/it->first-1)>tolerance) continue;
        // check similarity with existing planes
        bool tooSimilar = false;
        forall(CellDivision& cd2, primaryPlanes){
          if(planeSimiliarity(cd.divPlane.displ, cd.divPlane.nrml, cd2.divPlane.displ, cd2.divPlane.nrml) > planeSim){
            tooSimilar = true;
          }
        }
        if(!tooSimilar){
          primaryPlanes.push_back(cd);
        }
      }
    } else if(direction == "Max"){
      for(auto it = inputDivisionsSorted.end(); it != inputDivisionsSorted.begin();){
        it--;
        CellDivision cd = it->second;
        //double value = cd.getValue(primary);
        // in range of tolerance?
        if(norm(it->first/optValue-1)>tolerance) continue;
        // check similarity with existing planes
        bool tooSimilar = false;
        forall(CellDivision& cd2, primaryPlanes){
          if(planeSimiliarity(cd.divPlane.displ, cd.divPlane.nrml, cd2.divPlane.displ, cd2.divPlane.nrml) > planeSim){
            tooSimilar = true;
          }
        }
        if(!tooSimilar){
          primaryPlanes.push_back(cd);
        }
      }
    } else if(direction == "Random"){
      int idx = std::rand() % inputDivisions.size();
      if(idx < 0) idx = 0;
      if((uint)idx > inputDivisions.size() - 1) idx = inputDivisions.size();
      std::cout << "Random! " << idx << "/" << inputDivisions.size() << std::endl;
      primaryPlanes.push_back(inputDivisions[idx]);
    }

    bestDivPlane = make_pair(primaryPlanes[0].divPlane.nrml, primaryPlanes[0].divPlane.displ);
    outputDivisions = primaryPlanes;
  }


 // new routine for display planes
 // finds the lowest/highest value of certain measures in CellDivisions
 void findBestPlane(std::vector<CellDivision>& possibleDivisions, std::vector<CellDivision>& optPlanes,
  double& optValue, double tolerance, P3dDouble& bestDivPlane, QString primary, QString secondary, QString primaryDir,
  QString secondaryDir, double planeSim = 1.0)
 {
  std::vector<CellDivision> filteredPlanes;
  filterPlanes(possibleDivisions, filteredPlanes, optValue, tolerance, bestDivPlane, primary, primaryDir, planeSim);

  if(secondary != "None"){
    filterPlanes(filteredPlanes, optPlanes, optValue, tolerance, bestDivPlane, secondary, secondaryDir, planeSim);
  } else {
    optPlanes = filteredPlanes;
  }

 }

  void calcBetweennessMeasures(std::vector<CellDivision>& possibleDivisions, std::vector<cellGraph>& cellGraphsToTest, std::vector<cCPair>& cellPairsToTest,
   std::map<P3dDouble, int>& divPlaneCellGraphIdxMap, std::vector<P3dDouble>& samplingDirections, bool constWeight, bool dijk, int percGPU)
  {

    unsigned concurentThreadsSupported = omp_get_max_threads();
    //int nthreads;
    int tid;
    //cout << "bla " << percGPU << "/" << concurentThreadsSupported << endl;

    // now with openmp, calculate betweenness
    //progressStart("Test Division Planes - Calculate Betweenness", 0);
    AttrMap<int, double> betweenessMap, rwBetweenessMap;

    if(!dijk){
      cout << cellGraphsToTest.size() << " different cell graphs will be considered" << endl;
      for(size_t i = 0; i<cellGraphsToTest.size(); i++){
        betweenessMap[i] = 0;
      }
      std::map<IntInt, double> equalWeights;

      #pragma omp parallel for schedule(guided)
      for(size_t i = 0; i<cellGraphsToTest.size(); i++){
        tid = omp_get_thread_num();
        std::map<cell, double> cellBetwMap, cellRWBetwMap;
        cCPair currentPair = cellPairsToTest[i];
        double num = (double)(tid+1)/(double)concurentThreadsSupported*100;
        bool useGPU = false;
        if(num <= percGPU or percGPU>=100) useGPU = true;
        calcBetweenness(cellGraphsToTest[i], cellBetwMap);
        calcRandWalkBetweenness(cellGraphsToTest[i], equalWeights, cellRWBetwMap, useGPU, false);
        betweenessMap[i] = cellBetwMap[currentPair.first] + cellBetwMap[currentPair.second];
        rwBetweenessMap[i] = cellRWBetwMap[currentPair.first] + cellRWBetwMap[currentPair.second];
      }

      for(size_t i = 0; i<possibleDivisions.size(); i++){
        CellDivision& cd = possibleDivisions[i];
        cd.betw = betweenessMap[divPlaneCellGraphIdxMap[samplingDirections[cd.idxSamples]]];
        cd.rwBetw = rwBetweenessMap[divPlaneCellGraphIdxMap[samplingDirections[cd.idxSamples]]];
      }

    } else { // with weights (dijkstra)

      // weights
      #pragma omp parallel for
      for(size_t i = 0; i<possibleDivisions.size(); i++){
        CellDivision& cd = possibleDivisions[i];
        cd.weights = generateWeights(cd.nbhdMap, constWeight);
      }
      // betw centr
      #pragma omp parallel for schedule(dynamic)
      for(size_t i = 0; i<possibleDivisions.size(); i++){
        std::map<cell, double> cellBetwMap;
        CellDivision& cd = possibleDivisions[i];
        if(!calcBetweenness(cd.C, cd.weights, cellBetwMap)) cout << cd.divPlane.nrml << "/" << cd.divPlane.displ << endl;
        cd.betw = cellBetwMap[cd.dL] + cellBetwMap[cd.dR];
      }
      progressStart("Test Division Planes - Calculate Random Walk Betweenness", 0);
      // rw betw
      #pragma omp parallel for schedule(dynamic) private(tid)
      for(size_t i = 0; i<possibleDivisions.size(); i++){
        tid = omp_get_thread_num();
        //nthreads = omp_get_num_threads();
        double num = (double)(tid+1)/(double)concurentThreadsSupported*100;
        bool useGPU = false;
        if(num <= percGPU or percGPU>=100) useGPU = true;
        //cout << "bla " << tid << "/" << useGPU << "/" << num << endl;
        std::map<cell, double> cellRWBetwMap;
        CellDivision& cd = possibleDivisions[i];
        calcRandWalkBetweenness(cd.C, cd.weights, cellRWBetwMap, useGPU, false);
        cd.rwBetw = cellRWBetwMap[cd.dL] + cellRWBetwMap[cd.dR];
      }

    }
  }

   bool createSamplingDirs(AttrMap<int, Point3d>& surfaceDirs, CellTissue& T, cell& cellToDivide, DivisionParms& p,
    std::vector<P3dDouble>& samplingDirections)
  {

    std::cout << "perp " << p.perpendicular << "/" << p.perpendicularNeighbor << "/" << p.perpendicularConstrained << std::endl;

    //std::cout << "datp " << dat.cellNormalMap.size() << "/" << cellToDivide->label << std::endl;

    std::set<Point3d> perpNeighborPlanes;
    std::unordered_map<Point3d, std::vector<Point3d> > perpNeighborWalls;
    std::unordered_map<Point3d, Point3d> perpNeighborPoints;
    //if(perpendicularNeighbor or cuttingNeighborWall){
      // go through all neighbor cells, count their number of neighbors, choose the lowest ones and approximate the plane

      std::set<cell> minCells;

      int minNum = 1E6;

      forall(const cell& c, T.C.neighbors(cellToDivide)){
        int count = 0;
        forall(const cell& n, T.C.neighbors(c)){
          count++;
        }
        //std::cout << "neigh " << c->label << "/" << count << std::endl;
        if(count < minNum){
          std::set<cell> minCellsNew;
          minCellsNew.insert(c);
          minCells = minCellsNew;
          minNum = count;
        } else if(count == minNum){
          minCells.insert(c);
        }
      }

      std::cout << "neighbor cells with lowest degree: " << std::endl;
      forall(const cell& c, minCells){
        std::cout << c->label << " ";
        std::vector<Point3d> pointsPolygon = getSharedVtxs(cellToDivide,c);
        Point3d planePos, planeNrml;
        findPolygonPlane(pointsPolygon, planePos, planeNrml);
        perpNeighborPlanes.insert(planeNrml);
        perpNeighborWalls[planeNrml] = pointsPolygon;

      }
      std::cout << std::endl;

      typedef std::pair<Point3d, std::vector<Point3d> > P3dVecP3dPair;
      forall(P3dVecP3dPair p, perpNeighborWalls){
        if(p.second.empty()) continue;
        Point3d avg(0.,0.,0.);
        int pointCount = 0;
        forall(Point3d p3d, p.second){
          avg+=p3d;
          pointCount++;
        }
        avg/=pointCount;
        perpNeighborPoints[p.first] = avg;
        //std::cout << "wallCenter " << avg << "/" << p.first << std::endl;
      }


    // generate the division planes to test

    //std::cout << "perps " << p.perpendicular << "/" << p.perpendicularNeighbor << "/" << p.perpendicularConstrained << "/" << p.cuttingNeighborWall << std::endl;

    double sampling_res = 180.0;
    int numberOfPlanes = p.planeNr;

    do{


      if(p.perpendicular or p.perpendicularNeighbor){
        Point3d surfaceNormal = Point3d(0,0,0);
        if(p.perpendicular){
          surfaceNormal = surfaceDirs[cellToDivide->label];
          if(surfaceNormal == Point3d(0,0,0)){
            std::cout << "Error: No nearest Surface point found. Analyze the Surface first or set Perpendicular to Surface to No" << std::endl;
            return false;
          }
        }
        sampling_res=sampling_res/(double)numberOfPlanes*(2*p.stepsCenter+1);
        std::set<Point3d> perpNeighborPlanesToBeUsed;
        if(p.perpendicularNeighbor) perpNeighborPlanesToBeUsed = perpNeighborPlanes;



        generateSamplingDirections(sampling_res, p.stepsCenter, p.stepSizeCenter,
          surfaceNormal, perpNeighborPlanesToBeUsed, p.perpendicularConstrained, samplingDirections);

        if(p.perpendicular and p.perpendicularConstrained){
          addCenterDisplacementTestCut(p.stepsCenter, p.stepSizeCenter, cellToDivide->pos, surfaceNormal, perpNeighborPlanes,
                                      perpNeighborPoints, perpNeighborWalls, samplingDirections);
        }
      } else { // no restrictions
        sampling_res=sampling_res/std::pow((double)numberOfPlanes/(double)(2*p.stepsCenter+1), 1/3.);
        generateSamplingDirections(sampling_res, p.stepsCenter, p.stepSizeCenter, samplingDirections);
      }

      //std::cout << "Debug " << perpendicular << "/" << perpendicularNeighbor << "/" << perpendicularConstrained << "/" << cuttingNeighborWall << "/" << samplingDirections.size() << std::endl;

      if(p.cuttingNeighborWall){

        forall(P3dDouble p, samplingDirections){
          //std::cout << "planes " << p.first << "/" << p.second << std::endl;
        }

        std::vector<P3dDouble> samplingDirectionsNew;
        //std::cout << "cut " << perpNeighborPlanes.size() << std::endl;
        for(size_t i = 0; i< samplingDirections.size(); i++){
          P3dDouble plane = samplingDirections[i];
          bool found = false;
          forall(Point3d p, perpNeighborPlanes){
            //std::cout << "plane " << plane.first << "/" << plane.second << "//" << perpNeighborWalls[p].size() << std::endl;
            if(!found and planePolygonIntersect(plane.first, cellToDivide->pos+plane.first*plane.second, perpNeighborWalls[p])){
              samplingDirectionsNew.push_back(plane);
              found = true;
            }
          }
        }
        samplingDirections = samplingDirectionsNew;
                forall(P3dDouble p, samplingDirections){
          //std::cout << "planesKeep " << p.first << "/" << p.second << std::endl;
        }
      }

      if(numberOfPlanes > 50000){
        std::cout << "Either too many planes to be tested or unable to find suitable planes (due to rectrictions?)" << std::endl;
        return false;
      }
      if(samplingDirections.size() == 0) numberOfPlanes *= 10;
    } while(samplingDirections.size() == 0);

    p.perpNeighborWalls = perpNeighborWalls;
    p.perpNeighborPlanes = perpNeighborPlanes;

    cout << "Test division of cell with label " << cellToDivide->label << endl;
    cout << samplingDirections.size() << " possible division planes will be tested (delta angel " << sampling_res << ")" << endl;

    return true;

  }

// - check if cell with label exists and return it
// - calculate the maxLabel of the mesh
// - create a neigbormap
// - check if cells are isolated
// - update geometry for accurate volume calculations
bool checkUpdateMesh(CellTissue& T, cell& cellToDivide, int label, int& maxLabel, std::map<IntInt, double>& neighborMap)
{

    //cell cellToDivide;

    //std::map<IntInt, double> neighborMap;
    T.getNeighborhoodMap3D(neighborMap);

    //cout << "neighborhood " << double(clock()-begin)/CLOCKS_PER_SEC << endl;

    maxLabel = 0;

    if(!findCell(T.C, cellToDivide, label, maxLabel)){
      std::cout << "Cell with label " << label << " doesn't exist" << std::endl;
      return false;//setErrorMessage(QString("Cell with label %1 does not exist!").arg(label));
    }

    // look for isolated cells with no neighbors
    // this can cause trouble for the graph betweeness measures (matrix inversion!)
    std::set<cell> isolatedCells;
    forall(const cell& c, T.C){
      int count = 0;
      forall(const cell& n, T.C.neighbors(c)){
        count++;
      }
      if(count == 0 and c->label != label) isolatedCells.insert(c);
    }

    if(isolatedCells.size() > 0){
      cout << "There are isolated cells: ";
      forall(const cell& c, isolatedCells){
        cout << " " << c->label << " /";
      }
      cout << endl;
      //std::cout << "Error: There are isolated cells" << std::endl;
      return false;
    }

    // update the cell
    T.updGeometry(cellToDivide, 4);

  return true;
}

   //bool TestDivisionPlanes::testDivisionPlanesFunc(Mesh *mesh, Mesh *mesh2, int label, int divPlanes, bool perpendicular, bool perpendicularNeighbor,
   //  int cellCenterVarSteps, double cellCenterVarStepSize, bool displayResult, QString filename, bool constWeight, bool dijk,
   //  int percGPU, bool calcBetw, bool overwrite, bool cellCenterSelVertex, bool perpendicularConstrained, bool cuttingNeighborWall, Point3d divisionPoint)
   bool TestDivisionPlanes::testDivisionPlanesFunc(Mesh *mesh, Mesh *mesh2, int label, bool displayResult, QString filename, bool constWeight, bool dijk,
     int percGPU, bool calcBetw, bool overwrite, QString cellCenterSelVertex, DivisionParms& par, Point3d divisionPoint)

  {

    AttrMap<int, Point3d>& surfaceDirections = mesh->attributes().attrMap<int, Point3d>("Surface Directions");
    if(surfaceDirections.empty()){
      std::cout << "No Surface Normals found." << std::endl;
    }

    AttrMap<int, CellDivisionAttr>& attrDivData = mesh->attributes().attrMap<int, CellDivisionAttr>("Cell Division Simulated Planes");
    //attrDivData.clear();

    //AttrMap<int, CellDivisionAttr>& actualPlane = mesh->attributes().attrMap<int, CellDivisionAttr>("Cell Division Actual Plane");

    std::cout << "dat " << dat.cellNormalMap.size() << std::endl;

    DivisionAnalysisData datNew;
    datNew.cellNormalMap = dat.cellNormalMap;
    datNew.cellBetwMap = dat.cellBetwMap;
    datNew.cellCurrentFlowMap = dat.cellCurrentFlowMap;

    dat = datNew;

    int divPlanes = par.planeNr;

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

    //std::cout << "center " << cellCenterSelVertex << "/" << perpendicularConstrained << "/" << perpendicularNeighbor << std::endl;

    progressStart("Test Division Planes - Initialize", 0);
    // parameter check
    if(divPlanes<=0) divPlanes = 1;

    CellTissue &T = mesh->tissue();

    //clock_t begin = clock();
    //cout << "start Div Plane Testing " << double(clock()-begin)/CLOCKS_PER_SEC << endl;

    cell cellToDivide;
    int maxLabel;

    std::map<IntInt, double> neighborMap;
    if(!checkUpdateMesh(T, cellToDivide, label, maxLabel, neighborMap)){
      if(cellToDivide->label > 0){
        std::cout << "Error: There are isolated cells (Can't calculate betweeness measures)" << std::endl;
        if(calcBetw) return false;
      } else {
        return false;
      }
    }
    double sizeCell = avgCellSizeMGX3D(cellToDivide);

    std::vector<P3dDouble> samplingDirections;
    createSamplingDirs(surfaceDirections, T, cellToDivide, par, samplingDirections);

    if(samplingDirections.empty()){
      std::cout << "Error: No possible division planes generated." << std::endl;
      return false;
    }

    std::map<CellSetPair, P3dDouble> testedGraphs;
    std::map<P3dDouble, int> divPlaneCellGraphIdxMap;
    std::vector<cellGraph> cellGraphsToTest;
    std::vector<cCPair> cellPairsToTest;

    //std::map<cell, std::vector<tri> > neighbTris;
    std::map<Triangle, int> trisNeighbMap; // contains all tris of the cellToDivide

    cellTriNeighbors(T.C, cellToDivide, trisNeighbMap);

    // forall division planes:
    progressStart("Test Division Planes - Simulate Divisions", 0);

    std::vector<CellDivision> possibleDivisions(samplingDirections.size());

    //std::cout << "cell center " << divisionPoint << std::endl;
    if(cellCenterSelVertex == "Manual"){
      //divisionPoint = cellToDivide->pos;
      VtxVec selectedV = mesh->selectedVertices();
      std::cout << "selected: " << selectedV.size() << std::endl;
      forall(const vertex& v, selectedV){
      //forall(const vertex& v, T.S){
        if(!v->selected) continue;
        divisionPoint = v->pos;
        std::cout << "selected vertex found!" << std::endl;
      }

      forall(const cell& c, T.C){
        forall(const vertex& v, c->S){
          if(!v->selected) continue;
          divisionPoint = v->pos;
          std::cout << "selected vertex found here!" << std::endl;
        }
      }

      for(size_t i = 0; i< samplingDirections.size(); i++){
        std::cout << "before " << samplingDirections[i].first << "/" << samplingDirections[i].second << std::endl;
        double sN;
        Point3d intersectN;
        P3dDouble plane = samplingDirections[i];
        Point3d dir = divisionPoint-plane.first;
        planeLineIntersect(cellToDivide->pos+plane.first*plane.second, plane.first, divisionPoint, dir, sN, intersectN);
        std::cout << "p1 " << divisionPoint << "/" << cellToDivide->pos << std::endl;
        std::cout << "p2 " << sN << "/" << intersectN << std::endl;
        if(dir * plane.first > 0)
          samplingDirections[i].second = norm(intersectN - divisionPoint);
        else
          samplingDirections[i].second = -norm(intersectN - divisionPoint);

        std::cout << "after " << samplingDirections[i].first << "/" << samplingDirections[i].second << std::endl;
      }
    }
    std::cout << "cell div point" << divisionPoint << std::endl;

    //cout << "simulate divisions " << double(clock()-begin)/CLOCKS_PER_SEC << endl;
    #pragma omp parallel for
    for(size_t i = 0; i< samplingDirections.size(); i++){
      if(!progressAdvance()) continue;

      P3dDouble plane = samplingDirections[i];
      CellDivision cd (plane, T, cellToDivide, maxLabel);

      cd.cutLowDegreeNeighbor = false;

      // set cutLowDegreeNeighbor true if any of the neighbor planes is cut
      forall(Point3d p, par.perpNeighborPlanes){
        if(planePolygonIntersect(plane.first, cellToDivide->pos+plane.first*plane.second, par.perpNeighborWalls[p])) cd.cutLowDegreeNeighbor = true;
      }

      // // simulate division
      cd.simulateDivision3D(neighborMap, trisNeighbMap);

      if(calcBetw){
        if(!dijk){ // combine identical cell graphs

          std::set<cell> neighborDl, neighborDr;

          // update cell graph
          forall(const cell& c, cd.C.neighbors(cd.dL)){
            neighborDl.insert(c);
          }
          forall(const cell& c, cd.C.neighbors(cd.dR)){
            neighborDr.insert(c);
          }
          CellSetPair currentPair = std::make_pair(neighborDl,neighborDr);
          CellSetPair currentPairRev = std::make_pair(neighborDr,neighborDl);

          // write necessary cellGraph to map
          #pragma omp critical
          if(testedGraphs.find(currentPair) == testedGraphs.end()){
            testedGraphs[currentPair] = plane;
            divPlaneCellGraphIdxMap[plane] = cellGraphsToTest.size();
            cellGraphsToTest.push_back(cd.C);
            cellPairsToTest.push_back(std::make_pair(cd.dL,cd.dR));

          } else {
            divPlaneCellGraphIdxMap[plane] = divPlaneCellGraphIdxMap[testedGraphs[currentPair]];
          }

        }
      }
        possibleDivisions[i] = cd;
//cout << "done" << double(clock()-begin)/CLOCKS_PER_SEC << endl;
    }

    // find the shortest wall
    P3dDouble shortestDivPlane;
    double shortest = HUGE_VAL;
    forall(const CellDivision& cd, possibleDivisions){
      //cout << "cd " << cd.divPlane.area << "/" << cd.divPlane.nrml << "/" << cd.divPlane.displ << endl;
      if(cd.divPlane.area < shortest){
        P3dDouble p = make_pair(cd.divPlane.nrml, cd.divPlane.displ);
        shortest = cd.divPlane.area;
        shortestDivPlane = p;
      }
    }

    cout << "total planes:" << samplingDirections.size() << endl;

    cout << "pos div " << possibleDivisions.size() << endl;

    //std::map<P3dDouble, double> divPlaneSizeMap = divPlaneSizeMap2;
    //P3dDouble shortestDivPlane = findSmallestDivisionPlane(cellToDivide, samplingDirections, divPlaneSizeMap);

    cout << "smallest all  :" << shortestDivPlane.first << "/" << shortestDivPlane.second << endl;

    dat.divPlaneShortest = DivisionPlane(sizeCell, shortestDivPlane.first, cellToDivide->pos + shortestDivPlane.first*shortestDivPlane.second, 2);

    //cout << "betweenness " << double(clock()-begin)/CLOCKS_PER_SEC << endl;
    if(calcBetw){

      // calculate the betweeness measures for all possible divisions
      progressStart("Test Division Planes - Calc Betweenness", 0);
      calcBetweennessMeasures(possibleDivisions, cellGraphsToTest, cellPairsToTest, divPlaneCellGraphIdxMap, samplingDirections, constWeight, dijk, percGPU);

      progressStart("Test Division Planes - Finalize", 0);

      double optValue, optValueRW;
      std::vector<CellDivision> optPlanes, optPlanesRW;
      P3dDouble bestDivPlane, bestDivPlaneRW;
      double toleranceBetw = 0.0001;

      // find the best (smallest) planes of the optimal graph measures planes
      //findBestPlane(possibleDivisions, optPlanes, optValue, toleranceBetw, bestDivPlane, false);
      //findBestPlane(possibleDivisions, optPlanesRW, optValueRW, toleranceBetw, bestDivPlaneRW, true);

      QString primary = "Betweenness Centrality";
      if(!constWeight) primary = "Betweenness Centrality (weighted)";
      QString primaryRW = "Random Walk Betweenness";
      if(!constWeight) primaryRW = "Random Walk Betweenness (weighted)";
      QString secondary = "Wall Area";
      QString primaryDir = "Min";
      QString secondaryDir = "Min";

      findBestPlane(possibleDivisions, optPlanes, optValue, toleranceBetw, bestDivPlane, primary, secondary, primaryDir, secondaryDir);
      findBestPlane(possibleDivisions, optPlanesRW, optValueRW, toleranceBetw, bestDivPlaneRW, primaryRW, secondary, primaryDir, secondaryDir);

      cout << "opt planes (betw centr): " << optPlanes.size() << endl;
      cout << "betw centr  :" << optValue << endl;
      cout << "opt planes (rand walk betw): " << optPlanesRW.size() << endl;
      cout << "rand walk betw  :" << optValueRW << endl;
      cout << "smallest best (betw centr):" << bestDivPlane.first << "/" << bestDivPlane.second << endl;
      cout << "smallest best (rand walk betw):" << bestDivPlaneRW.first << "/" << bestDivPlaneRW.second << endl;

      dat.divPlaneShortestBest = DivisionPlane(sizeCell, bestDivPlane.first, cellToDivide->pos + bestDivPlane.first*bestDivPlane.second, 1);
      dat.divPlaneShortestBestRW = DivisionPlane(sizeCell, bestDivPlaneRW.first, cellToDivide->pos + bestDivPlaneRW.first*bestDivPlaneRW.second, 3);

    }

    // now generate and display the division planes in the non-active mesh
    if(displayResult){
      cout << "display result" << endl;
      if(overwrite){
        mesh->reset();
        cout << "display result " << __LINE__ << endl;
        //drawPlane(*this, mesh, dat.divPlaneFlat);
        cout << "display result " << __LINE__ << endl;
        drawPlane(*this, mesh, dat.divPlaneShortest);
        cout << "display result " << __LINE__ << endl;
        if(calcBetw) drawPlane(*this, mesh, dat.divPlaneShortestBest);
        if(calcBetw) drawPlane(*this, mesh, dat.divPlaneShortestBestRW);
        cout << "display result " << __LINE__ << endl;
        mesh->updateAll();
        cout << "display result " << __LINE__ << endl;
      } else {
        mesh2->reset();
        cout << "display result " << __LINE__ << endl;
        //drawPlane(*this, mesh2, dat.divPlaneFlat);
        cout << "display result " << __LINE__ << endl;
        drawPlane(*this, mesh2, dat.divPlaneShortest);
        cout << "display result " << __LINE__ << endl;
        if(calcBetw) drawPlane(*this, mesh2, dat.divPlaneShortestBest);
        if(calcBetw) drawPlane(*this, mesh2, dat.divPlaneShortestBestRW);
        cout << "display result " << __LINE__ << endl;
        mesh2->updateAll();
        cout << "display result " << __LINE__ << endl;
      }
    }
cout << "done  " << possibleDivisions.size() << endl;

    int lCounter = 1;
    for(uint i=0; i<possibleDivisions.size(); i++){
      possibleDivisions[i].divPlane.label = lCounter++;
      //cout << "planes " << possibleDivisions[i].divPlane.nrml << endl;
    }
    dat.divData = possibleDivisions;
    if(filename != "" and possibleDivisions.size()>0){
      cout << "save file" << endl;
      QString header = "divPlaneNrml X,divPlaneNrml Y,divPlaneNrml Z,divPlaneNrml D,area,betw centr,rand walk betw,vol ratio,vol left,vol right";
      vector< vector<double> > data;

      forall(CellDivision& cd, possibleDivisions){
        std::vector<double> d;
        cd.writeDataToVec(d);
        data.push_back(d);
      }

      if(!writeCsvFile(filename, header, data))
        std::cout << "Error writing file: " << std::endl;
        //return setErrorMessage(QString("Error writing file: ").arg(filename));
    }
    // cout << "attr map  " << possibleDivisions.size() << endl;
    // int counter = 0;
    // actualPlane.clear();
    // forall(auto cd, possibleDivisions){
    //   CellDivisionAttr cda;
    //   cda.planeNrml = cd.divPlane.nrml;
    //   cda.planePos = cd.divPlane.pos;
    //   cda.planeDisplacement = cd.divPlane.displ;
    //   cda.div2d = 1;
    //   cda.cellNormal = Point3d(0,0,0); // TODO
    //   //if(!m->useParents()) cda.cellLabel = selectedLabel;
    //   //else cda.parentLabel = selectedParent;

    //   cda.volumeRatio = cd.volumeRatio;
    //   cda.planeArea = cd.divPlane.area;
    //   cout << "cd  " << cda.planeArea << "/" << cda.planeDisplacement << "/" << actualPlane.size() << endl;
    //   if(!actualPlane.empty()){
    //     cda.angleActual = 180/3.1415*acos(actualPlane[0].planeNrml * cda.planeNrml / norm(actualPlane[0].planeNrml) / norm(cda.planeNrml));
    //     if(cda.angleActual > 90) cda.angleActual = 180 - cda.angleActual;
    //   }

    //   cout << "cd  2" << endl;

    //   attrDivData[counter] = cda;
    //   counter++;
    // }

cout << "done  " << endl;

    datTest = dat;

    updateState();
    updateViewer();
    //cout << "done " << double(clock()-begin)/CLOCKS_PER_SEC << endl;
    return true;
  }


   //bool TestDivisionPlanes::run(Mesh *m, Mesh *m2, int label1, int label2, int divPlanes, bool perpendicular, bool perpendicularNeighbor,
   //  int cellCenterVarSteps, double cellCenterVarStepSize, bool displayResult, QString filename, bool constWeight,
   //  bool dijk, int percGPU, bool calcBetw, bool cellCenterSelVertex, bool perpendicularConstrained, bool cuttingNeighborWall)
   bool TestDivisionPlanes::run(Mesh *m, Mesh *m2, DivisionParms& par, int label1, int label2, bool displayResult, QString filename, bool constWeight,
     bool dijk, int percGPU, bool calcBetw, QString cellCenterSelVertex)

  {

    AttrMap<int, vertex>& divPoint = m->attributes().attrMap<int, vertex>("Division Analysis Custom Point");

    Point3d divisionPoint = divPoint[0]->pos;

    if(label2 > 0){
       // copy mesh,
      CopySwapMeshes copy(*this);
      if(!copy.run("1 -> 2")) return setErrorMessage(QString("Error while copying mesh."));

      if(cellCenterSelVertex == "Center of Actual"){
        CellTissue &T2 = m2->tissue();
        cell c1, c2;
        bool found1=false, found2=false;
        forall(const cell& c, T2.C){
          if(c->label == label1){
            c1 = c;
            found1 = true;
          }
          if(c->label == label2){
            c2 = c;
            found2 = true;
          }
        }
        std::vector<Point3d> sharedVtx = getSharedVtxs(c1, c2);

        if(!found1 or !found2 or sharedVtx.empty()){
          divisionPoint = Point3d(0,0,0);
        } else {
          divisionPoint = Point3d(0,0,0);
          forall(Point3d p, sharedVtx){
            divisionPoint += p;
          }
          divisionPoint/=sharedVtx.size();
        }

      }



      // join cells on copy
      JoinCellsMGX3D join(*this);
      if(!join.run(m2,label1, label2)) return setErrorMessage(QString("Error while joining cells. Make sure they exist and are neighbors."));

      // now run the standard process on the mesh with the joined cells
      //TestDivisionPlanes div(*this);



      if(cellCenterSelVertex != "Centroid"){
        if(!testDivisionPlanesFunc(m2, m, label1, displayResult, filename,
          constWeight, dijk, percGPU, calcBetw, true, cellCenterSelVertex, par, divisionPoint))
          return setErrorMessage(QString("Error doing the test divisions."));
      } else {
        if(!testDivisionPlanesFunc(m2, m, label1, displayResult, filename,
          constWeight, dijk, percGPU, calcBetw, true, cellCenterSelVertex, par))
          return setErrorMessage(QString("Error doing the test divisions."));
      }
    } else {
      if(cellCenterSelVertex != "Centroid"){
        if(!testDivisionPlanesFunc(m, m2, label1, displayResult, filename,
          constWeight, dijk, percGPU, calcBetw, false, cellCenterSelVertex, par, divisionPoint))
          return setErrorMessage(QString("Error doing the test divisions."));
        } else {
        if(!testDivisionPlanesFunc(m, m2, label1, displayResult, filename,
          constWeight, dijk, percGPU, calcBetw, false, cellCenterSelVertex, par))
          return setErrorMessage(QString("Error doing the test divisions."));
      }
    }




    return true;
  }
  REGISTER_PROCESS(TestDivisionPlanes);

  // simulate a division along a specified division plane
  CellDivision simDiv(P3dDouble plane, CellTissue& T, cell cellToDivide, int maxLabel)
  {
    //Point3d cellCenter = cellToDivide->pos;

    std::map<IntInt, double> neighborMap;
    T.getNeighborhoodMap3D(neighborMap);

    std::map<Triangle, int > trisNeighbMap; // contains all tris of the cellToDivide
    cellTriNeighbors(T.C, cellToDivide, trisNeighbMap);

    // create CellDivision object and simulate division
    CellDivision cd(plane, T, cellToDivide, maxLabel);
    cd.simulateDivision3D(neighborMap, trisNeighbMap);

    return cd;
  }

  // post division (with two cells)
   bool TestDivisionPlanesCellAtlas::run(Mesh *m, Mesh *m2, int label1, int label2, bool displayResult, RootDivData& result, bool copyToM2)
  {

    //CellTissue &T = m->tissue();
    CellTissue &T2 = m2->tissue();

    AttrMap<int, Point3d>& dirRadAttr = m->attributes().attrMap<int, Point3d>("Measure Label Vector DirRadial");
    AttrMap<int, Point3d>& dirLongAttr = m->attributes().attrMap<int, Point3d>("Measure Label Vector DirLongitudinal");
    AttrMap<int, Point3d>& dirCircAttr = m->attributes().attrMap<int, Point3d>("Measure Label Vector DirCircumferential");

    AttrMap<int, double>& coordRadAttr = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Radial");
    AttrMap<int, double>& coordLongAttr = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Longitudinal");
    AttrMap<int, double>& coordCircAttr = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Circumferential");

    if(dirRadAttr.empty() or dirLongAttr.empty() or dirCircAttr.empty() or coordRadAttr.empty() or coordLongAttr.empty() or coordCircAttr.empty())
      return setErrorMessage("Can't find Cell Atlas Attributes. Run Cell Atlas first!");

     // copy mesh
    if(copyToM2){
      m2->reset();
      CopySwapMeshes copy(*this);
      if(!copy.run("1 -> 2")) return setErrorMessage(QString("Error while copying mesh."));
    }

    // convert to mgx3d if required
    //Information::out << "start" << endl;

    const QString meshType = m2->tissue().meshType();
    //Information::out << "meshtype in " << meshType.toLocal8Bit().constData() << endl;
    if(meshType != "MGX3D"){
    //Information::out << "convert" << endl;
      if(!T2.toMgx3d())
        throw(QString("Unexpected error converting to normal (MGX3D) mesh"));
    }

    // get the actual plane
    CellDivisionAttr cda;
    FlatDivisionPlane approxRealPlane(*this);
    approxRealPlane.run(m2, m2, label1, label2, false, dat, cda);

    // join cells on copy
    JoinCellsMGX3D join(*this);
    if(!join.run(m2,label1, label2)){
      Information::out << "Problem with joining cells" << label1 << " and " << label2 << ". Make sure they exist and are neighbors." << endl;
      return setErrorMessage(QString("Error while joining cells. Make sure they exist and are neighbors."));
    }

    Point3d dirRad, dirLong, dirCirc;
    double coordLong = 0., coordRad = 0., coordCirc = 0.;
    int count = 0;

    forall(auto p, dirRadAttr){
      int currentLabel = p.first;//(*cellAtlasAttr)[i].cellLabel;
      if(currentLabel != label1 and currentLabel != label2) continue;
      count++;

      dirRad += dirRadAttr[currentLabel] * 0.5;
      dirLong += dirLongAttr[currentLabel] * 0.5;
      dirCirc += dirCircAttr[currentLabel] * 0.5;

      coordLong += coordLongAttr[currentLabel] * 0.5;
      coordRad += coordRadAttr[currentLabel] * 0.5;
      coordCirc += coordCircAttr[currentLabel] * 0.5;

    }
    if(count != 2){
      Information::out << "not exactly one cell with each of the labels" << endl;
      return false;
    }

    P3dDouble planeRad = make_pair(dirRad,0);
    P3dDouble planeLong = make_pair(dirLong,0);
    P3dDouble planeCirc = make_pair(dirCirc,0);

    P3dDouble planeReal = make_pair(dat.divPlaneFlat.planeNrml, 0);

    Point3d cellCenter;
    int maxLabel;
    cell cellToDivide;
    findCell(T2.C, cellToDivide, label1, maxLabel);
//Information::out << "found cell " << cellToDivide->label << endl;
    CellDivision cdRad = simDiv(planeRad, T2, cellToDivide, maxLabel);
    CellDivision cdLong = simDiv(planeLong, T2, cellToDivide, maxLabel);
    CellDivision cdCirc = simDiv(planeCirc, T2, cellToDivide, maxLabel);

    CellDivision cdReal = simDiv(planeReal, T2, cellToDivide, maxLabel);
//Information::out << "created divisions" << endl;
    if(displayResult){
      dat.customPlanes.push_back(DivisionPlane(10.0, planeRad.first, cellToDivide->pos + planeRad.first*planeRad.second, 1));
      dat.customPlanes.push_back(DivisionPlane(10.0, planeLong.first, cellToDivide->pos + planeLong.first*planeLong.second, 2));
      dat.customPlanes.push_back(DivisionPlane(10.0, planeCirc.first, cellToDivide->pos + planeCirc.first*planeCirc.second, 3));

      dat.customPlanes.push_back(DivisionPlane(10.0, planeReal.first, cellToDivide->pos + planeReal.first*planeReal.second, 4));

      m2->reset();
      DisplayOptimalPlanes disp(*this);
      disp.run(m,m2,true,"","","","",dat, 0.01, 100, false,1,10.0);

    }

    //Information::out << "dirs " << dirRad << "/" << dirLong << "/" << dirCirc;

    result.dirRad = cdRad;
    result.dirLong = cdLong;
    result.dirCirc = cdCirc;
    result.dirReal = cdReal;

    result.coordLong = coordLong;
    result.coordRad = coordRad;
    result.coordCirc = coordCirc;

    return true;
  }
  REGISTER_PROCESS(TestDivisionPlanesCellAtlas);

  bool saveRootDivFile(QString filename, std::vector<RootDivData> result, Mesh *m)
  {
    // file for root cell data
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream out(&file);

    out << "Cell Label 1, Cell Label 2,Parent label,Longitudinal Coord, Radial Coord, Circumferential Coord";
    out << ",Actual Plane X,Actual Plane Y,Actual Plane Z,Actual Plane Area,Actual Vol Left,Actual Vol Right,Actual Vol Ratio";
    out << ",Long X,Long Y,Long Z,Long Plane Area,Long Vol Left,Long Vol Right,Long Vol Ratio,Long Angle";
    out << ",Rad X,Rad Y,Rad Z,Rad Plane Area,Rad Vol Left,Rad Vol Right,Rad Vol Ratio,Rad Angle";
    out << ",Circ X,Circ Y,Circ Z,Circ Plane Area,Circ Vol Left,Circ Vol Right,Circ Vol Ratio,Circ Angle";
    out << endl;

    int bestFitLong = 0;
    int bestFitRad = 0;
    int bestFitCirc = 0;

    AttrMap<int, double>& angleLongAttr = m->attributes().attrMap<int, double>("Measure Label Double AngleLong");
    AttrMap<int, double>& angleRadAttr = m->attributes().attrMap<int, double>("Measure Label Double AngleRad");
    AttrMap<int, double>& angleCircAttr = m->attributes().attrMap<int, double>("Measure Label Double AngleCirc");
    AttrMap<int, Point3d>& planeActualAttr = m->attributes().attrMap<int, Point3d>("Measure Label Vector ActualPlane");
    AttrMap<int, double>& planeActualAreaAttr = m->attributes().attrMap<int, double>("Measure Label Double ActualPlaneArea");
    AttrMap<int, double>& planeActualVolRatioAttr = m->attributes().attrMap<int, double>("Measure Label Double ActualVolumeRatio");

    angleLongAttr.clear();
    angleRadAttr.clear();
    angleCircAttr.clear();
    planeActualAttr.clear();
    planeActualAreaAttr.clear();
    planeActualVolRatioAttr.clear();

    forall(RootDivData rd, result){

      int parentLabel = m->parents()[rd.label1];

      double scalarLong = angleVectors(rd.dirLong.divPlane.nrml, rd.dirReal.divPlane.nrml);//abs(rd.dirLong.divPlane.nrml * rd.dirReal.divPlane.nrml);
      double scalarRad = angleVectors(rd.dirRad.divPlane.nrml, rd.dirReal.divPlane.nrml);//abs(rd.dirRad.divPlane.nrml * rd.dirReal.divPlane.nrml);
      double scalarCirc = angleVectors(rd.dirCirc.divPlane.nrml, rd.dirReal.divPlane.nrml);//abs(rd.dirCirc.divPlane.nrml * rd.dirReal.divPlane.nrml);

      out << rd.label1 << "," << rd.label2 << "," << m->parents()[rd.label1] << ",";
      out << rd.coordLong << "," << rd.coordRad << "," << rd.coordCirc << ",";
      out << rd.dirReal.divPlane.nrml.x() << ","  << rd.dirReal.divPlane.nrml.y() << ","  << rd.dirReal.divPlane.nrml.z() << "," << rd.dirReal.divPlane.area << ",";
      out << rd.dirReal.volumeLeft << ","  << rd.dirReal.volumeRight << "," << rd.dirReal.volumeRatio << ",";

      out << rd.dirLong.divPlane.nrml.x() << ","  << rd.dirLong.divPlane.nrml.y() << ","  << rd.dirLong.divPlane.nrml.z() << "," << rd.dirLong.divPlane.area << ",";
      out << rd.dirLong.volumeLeft << ","  << rd.dirLong.volumeRight << "," << rd.dirLong.volumeRatio << "," << scalarLong << ",";

      out << rd.dirRad.divPlane.nrml.x() << ","  << rd.dirRad.divPlane.nrml.y() << ","  << rd.dirRad.divPlane.nrml.z() << "," << rd.dirRad.divPlane.area << ",";
      out << rd.dirRad.volumeLeft << ","  << rd.dirRad.volumeRight << "," << rd.dirRad.volumeRatio << "," << scalarRad << ",";

      out << rd.dirCirc.divPlane.nrml.x() << ","  << rd.dirCirc.divPlane.nrml.y() << ","  << rd.dirCirc.divPlane.nrml.z() << "," << rd.dirCirc.divPlane.area << ",";
      out << rd.dirCirc.volumeLeft << ","  << rd.dirCirc.volumeRight << "," << rd.dirCirc.volumeRatio << "," << scalarCirc << ",";

      out << endl;

      if(scalarLong < scalarRad and scalarLong < scalarCirc) bestFitLong++;
      else if(scalarRad < scalarLong and scalarRad < scalarCirc) bestFitRad++;
      else bestFitCirc++;

      std::cout << "f: " << rd.label1 << "/" << rd.label2 << "//" << m->parents()[rd.label1] << "/" << m->parents()[rd.label2] << "//" << scalarLong << std::endl;

      angleLongAttr[m->parents()[rd.label1]] = scalarLong;
      angleRadAttr[m->parents()[rd.label1]] = scalarRad;
      angleCircAttr[m->parents()[rd.label1]] = scalarCirc;
      planeActualAttr[m->parents()[rd.label1]] = rd.dirReal.divPlane.nrml;
      planeActualAreaAttr[m->parents()[rd.label1]] = rd.dirReal.divPlane.area;
      planeActualVolRatioAttr[m->parents()[rd.label1]] = rd.dirReal.volumeRatio;
    }

    Information::out << result.size() << " cell pairs were tested." << endl;
    Information::out << "Best fitting Long: " << bestFitLong << "; Rad: " << bestFitRad << "; Circ: " << bestFitCirc << endl;

    return true;
  }

 bool RootDivisionAnalysis::run(Mesh *m, Mesh *m2, bool displayResult, QString filename)
  {

    vvGraph& S = m->graph();
    std::map<int, std::set<int> > parentCellLabelMap;
    typedef std::pair<int, std::set<int> > IntIntSetPair;

    // find all pairs
    forall(const vertex& v, S){
      if(m->parents()[v->label] > 0){
        parentCellLabelMap[m->parents()[v->label]].insert(v->label);
      }
    }


    CopySwapMeshes copy(*this);
    if(!copy.run("1 -> 2")) return setErrorMessage(QString("Error while copying mesh."));

    CellTissue& T2 = m2->tissue();

    const QString meshType = m2->tissue().meshType();
    if(meshType != "MGX3D"){
      //std::cout << "meshtype " << meshType.toLocal8Bit().constData() << std::endl;
      if(!T2.toMgx3d())
        throw(QString("Unexpected error converting to MGX3D mesh"));
    }

    std::vector<RootDivData> result;

    std::cout << "size2 " << parentCellLabelMap.size() << std::endl;

    forall(const IntIntSetPair& p, parentCellLabelMap){
      // run the div analysis for all pairs of same parent labels
      std::cout << "sizeP " << p.first << "/" << p.second.size() << std::endl;
      if(p.second.size() == 2){

        RootDivData d;

        d.label1 = *p.second.begin();
        d.label2 = *(++p.second.begin());

        std::cout << "Processing cells " << d.label1 << " and " << d.label2 << std::endl;
        TestDivisionPlanesCellAtlas divAnalysis(*this);
        if(divAnalysis.run(m, m2, d.label1, d.label2, false, d, false)){
          result.push_back(d);
          std::cout << "Done!" << std::endl;
        } else {
          std::cout << "Division analysis failed!" << std::endl;
        }

      } else {
        std::cout << "The following cells are ignored as they have a parent label that does not have 2 occurances: ";
        forall(int i, p.second) cout << i << " ";
        std::cout << std::endl;
      }
    }

      //int count = 0;
      forall(RootDivData rd, result){
        int color = 0;
        dat.customPlanes.push_back(DivisionPlane(10.0, rd.dirLong.divPlane.nrml, rd.dirLong.divPlane.pos, ++color));
        dat.customPlanes.push_back(DivisionPlane(10.0, rd.dirRad.divPlane.nrml, rd.dirRad.divPlane.pos, ++color));
        dat.customPlanes.push_back(DivisionPlane(10.0, rd.dirCirc.divPlane.nrml, rd.dirCirc.divPlane.pos, ++color));
        dat.customPlanes.push_back(DivisionPlane(15.0, rd.dirReal.divPlane.nrml, rd.dirReal.divPlane.pos, ++color));
       // cout++;
      }

    if(displayResult){
      m2->reset();
      DisplayOptimalPlanes disp(*this);
      disp.run(m,m2,true,"","","","",dat, 0.01, 100, false,1.0,10.0);
    }

    if(filename != "")
      saveRootDivFile(filename, result, m);

    return true;
  }
  REGISTER_PROCESS(RootDivisionAnalysis);


 bool DivisionAnalysisMulti3DOld::run(Stack* s1, Store* s, Mesh *m, Mesh *m2, bool singleCells, bool doubleCells, bool useStack,
  bool displayShortest, bool displayActual)
  {

    

    AttrMap<int, CellDivisionAttr>& attrData = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Actual Plane");
    attrData.clear();
    AttrMap<int, DivisionDataNew>& divDataNew = m->attributes().attrMap<int, DivisionDataNew>("Cell Division Data");
    divDataNew.clear();

    AttrMap<int, double>& vols = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Volume");

    vvGraph& S = m->graph();
    std::map<int, std::set<int> > parentCellLabelMap;
    typedef std::pair<int, std::set<int> > IntIntSetPair;

    // create parents -> cells map
    forall(const vertex& v, S){
      // std::cout << "ll " << v->label << "/" << m->parents()[v->label] << std::endl;
      if(m->parents()[v->label] > 0){
        parentCellLabelMap[m->parents()[v->label]].insert(v->label);
      }
    }

    std::cout << "multi3d " << parentCellLabelMap.size() << std::endl;

    std::vector<RootDivData> result;
    forall(const IntIntSetPair& p, parentCellLabelMap){
      // run the div analysis for all pairs of same parent labels
      if(p.second.size() == 2 and doubleCells){

        // copy mesh to mesh 2 to back it up
        CopySwapMeshes copy(*this);
        if(!copy.run("1 -> 2")) return setErrorMessage(QString("Error while copying mesh."));

        CellTissue& T2 = m2->tissue();

        const QString meshType = m2->tissue().meshType();
        if(meshType != "MGX3D"){
          //std::cout << "meshtype " << meshType.toLocal8Bit().constData() << std::endl;
          if(!T2.toMgx3d())
            throw(QString("Unexpected error converting to MGX3D mesh"));
        }

        RootDivData d;

        d.label1 = *p.second.begin();
        d.label2 = *(++p.second.begin());
        std::cout << "Processing cells " << d.label1 << " and " << d.label2 << "/Parent "<< p.first <<  std::endl;

        CellDivisionAttr actual;

        if(useStack){
          Stack* s1 = currentStack();
          Store* store1 = s1->currentStore();

          ActualPlaneStack aps(*this);
          aps.run(s1, store1, m2, m2, d.label1, d.label2, false, dat, actual);
        } else {
          FlatDivisionPlane fdp(*this);
          fdp.run(m2, m2, d.label1, d.label2, false, dat, actual);
        }

        actual.volumeRatio = calcRatio(vols[d.label1],vols[d.label2]);

        std::cout << "actual " << "/" << actual.divPlane.pos << "/" << "/" << actual.divPlane.area << std::endl;

        std::cout << "Div Planes " << attrData.size() << "/vols " << actual.volumeRatio << std::endl;

        // join cells on copy
        JoinCellsMGX3D join(*this);
        if(!join.run(m2,d.label1, d.label2)){
          Information::out << "Problem with joining cells" << d.label1 << " and " << d.label2 << ". Make sure they exist and are neighbors." << endl;
          return setErrorMessage(QString("Error while joining cells. Make sure they exist and are neighbors."));
        }

        // find centroid of joined cell

        DivisionParms dp;

        cell cellToDivide;
        int maxLabel;
        findCell(T2.C, cellToDivide, d.label1, maxLabel);
        Point3d centerPos = cellToDivide->pos;
        Point3d projected = projectPointOnPlane(centerPos, actual.divPlane.pos, actual.divPlane.nrml);
        double disPlane = norm(projected - centerPos);
        Point3d vectDispl = projected - centerPos;
        if(vectDispl * actual.divPlane.nrml > 0) disPlane = -disPlane;
        P3dDouble planeActual = make_pair(actual.divPlane.nrml,disPlane);
        CellDivision cdActual = simDiv(planeActual, T2, cellToDivide, maxLabel);

        actual.divPlane.area = cdActual.divPlane.area;
        actual.divPlane.displ = disPlane;
        attrData[p.first] = actual;

        std::cout << "actualstuff " << projected << "/" << centerPos << "/" << actual.divPlane.pos << "/" << disPlane << "/" << actual.divPlane.area << std::endl;

        dp.perpendicularNeighbor = false;
        dp.perpendicularConstrained = false;
        dp.cuttingNeighborWall = false;
        dp.perpendicular = false;

        dp.stepsCenter = 0;
        dp.stepSizeCenter = 0;

        dp.planeNr = 2000;

        TestDivisionPlanes t(*this);
        t.testDivisionPlanesFunc(m2, m2, d.label1, false, "", false, false, 0, false, false, "Centroid", dp);
        //t.testDivisionPlanesFunc(m, m2, cellToDivide->label, false, "", constWeight, dijk, 0, calcBetw, false, "Centroid", p);
        //Mesh *mesh, Mesh *mesh2, int label, bool displayResult, QString filename, bool constWeight, bool dijk,
        //  int percGPU, bool calcBetw, bool overwrite, QString cellCenterSelVertex, DivisionParms& par, Point3d divisionPoint = Point3d(0,0,0)
        //DivisionAnalysisData dat;
        dat = t.datTest;

        DivisionDataNew divDataAll;

        // convert cell division data to store in new attr map
        forall(CellDivision cd, dat.divData){
          CellDivisionAttr cda = convertToCDA(cd);
          //std::cout << "converted? " << cda.divPlane.nrml << std::endl;
          divDataAll.simulationData.push_back(cda);
        }

        divDataNew[p.first] = divDataAll;

        std::cout << "Planes " << dat.divData.size() << std::endl;

      } else if(p.second.size() == 1 and singleCells){
        RootDivData d;
        d.label1 = *p.second.begin();
      } else {
        std::cout << "The following cells are ignored (invalid number of parent occurances: ";
        forall(int i, p.second) cout << i << " ";
        std::cout << std::endl;
      }
    }


    //if(displayResult){
    m2->reset();
    m2->updateAll();
    //   DisplayOptimalPlanes disp(*this);
    //   disp.run(m,m2,true,"","","","",dat, 0.01, 100, false,1.0,10.0);
    //}

    return true;
  }
  REGISTER_PROCESS(DivisionAnalysisMulti3DOld);

 // bool DivisionAnalysisMulti::run(Mesh *m, Mesh *m2, bool runPre, bool runPost, bool displayResult, QString filename)
 //  {

 //    vvGraph& S = m->graph();
 //    std::map<int, std::set<int> > parentCellLabelMap;
 //    typedef std::pair<int, std::set<int> > IntIntSetPair;

 //    // find all pairs and single
 //    forall(const vertex& v, S){
 //      if(m->parents()[v->label] > 0){
 //        parentCellLabelMap[m->parents()[v->label]].insert(v->label);
 //      }
 //    }

 //    //CopySwapMeshes copy(*this);
 //    //if(!copy.run("1 -> 2")) return setErrorMessage(QString("Error while copying mesh."));

 //    //CellTissue& T2 = m2->tissue();
 //    CellTissue& T = m->tissue();

 //    const QString meshType = m2->tissue().meshType();
 //    if(meshType != "MGX3D"){
 //      //std::cout << "meshtype " << meshType.toLocal8Bit().constData() << std::endl;
 //      if(!T.toMgx3d())
 //        throw(QString("Unexpected error converting to MGX3D mesh"));
 //    }

 //    std::vector<RootDivData> result;
 //    forall(const IntIntSetPair& p, parentCellLabelMap){
 //      // run the div analysis for all pairs of same parent labels
 //      if((p.second.size() == 2 and runPost) or (p.second.size() == 1 and runPre)){

 //        RootDivData d;

 //        d.label1 = *p.second.begin();

 //        if(p.second.size() == 1) d.label2 = 0;
 //        else d.label2 = *(++p.second.begin());

 //        std::cout << "Processing cells " << d.label1 << " and " << d.label2 << std::endl;

 //        TestDivisionPlanes *divTest;
 //        if(!getProcess("Mesh/Division Analysis/Analysis 3D/B Cell Division Analysis MGX3D", divTest))
 //         throw(QString("DivisionAnalysisMulti:: Unable to make TestDivisionPlanes process"));

 //        divTest->setParm("Label Cell 1", QString::number(d.label1));
 //        divTest->setParm("Label Cell 2 (optional)", QString::number(d.label2));
 //        divTest->run();
 //      } else {
 //        std::cout << "The following cells are ignored as they have a parent label that does not have 1/2 occurances: ";
 //        forall(int i, p.second) cout << i << " ";
 //        std::cout << std::endl;
 //      }
 //    }

 //    return true;
 //  }
 //  REGISTER_PROCESS(DivisionAnalysisMulti);

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

 //  bool SaveCellGraph::run(Mesh *m, QString filename, double minArea)
 //  {
 //    // first get the cell graph
 //    std::map<int, vertex> labelCellVtxMap;

 //    std::map<IntInt, double> neighborMap;

 //    CellTissue &T = m->tissue();

 //    forall(const cell& c, T.C){
 //      Point3d cellCenter (0,0,0);
 //      int count = 0;
 //      forall(const vertex& v, c->S){
 //        cellCenter+=v->pos;
 //        count++;
 //      }
 //      cellCenter/=count;
 //      vertex v;
 //      v->pos = cellCenter;
 //      labelCellVtxMap[c->label] = v;
 //    }

 //    // now save the cell graph

 //    QFile file(filename);
 //    if(!file.open(QIODevice::WriteOnly))
 //    {
 //      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
 //      return false;
 //    }
 //    QTextStream out(&file);

 //    out << "Cell label,Neighbour Label,Wall Area" << endl;

 //    std::set<Tri> setOfTris;

 //    std::map<vertex, bool> vtxShared;
 //    std::set<vertex> setOfVtx;

 //    T.getNeighborhoodMap3D(neighborMap);

 //    forall(const IntIntDouble &p, neighborMap){
 //      if(p.second > minArea){
 //        out << p.first.first << "," << p.first.second << "," << p.second << endl;
 //        cout << p.first.first << "," << p.first.second << "/" << p.second << endl;
 //      }
 //    }

 //    return true;
 //  }
 //  REGISTER_PROCESS(SaveCellGraph);

  bool DisplayCellGraph::run(Mesh *m, Mesh *m2, double minArea)
  {

    std::map<IntInt, double> neighborMap;

    std::map<int, vertex> labelCellVtxMap;

    CellTissue &T = m->tissue();
    vvGraph &S = m2->graph();
    m2->reset();
    forall(const cell& c, T.C){
      Point3d cellCenter (0,0,0);
      int count = 0;
      forall(const vertex& v, c->S){
        cellCenter+=v->pos;
        count++;
      }
      cellCenter/=count;
      vertex v;
      v->pos = cellCenter;
      labelCellVtxMap[c->label] = v;
      S.insert(v);
    }

    std::map<vertex, bool> vtxShared;
    std::set<vertex> setOfVtx;

    T.getNeighborhoodMap3D(neighborMap);

    forall(const IntIntDouble &p, neighborMap)
      if(p.second > minArea)
        S.insertEdge(labelCellVtxMap[p.first.first], labelCellVtxMap[p.first.second]);

    m2->updateAll();
    m2->setShowSurface(false);
    m2->setMeshView("All");

    return true;
  }
  REGISTER_PROCESS(DisplayCellGraph);


  bool DisplayCellGraph2D::run(Mesh *m, Mesh *m2, double minArea)
  {

    std::map<IntInt, double> neighMap;
    vvGraph &S1 = m->graph();
    m->updateCentersNormals();

    IntPoint3fAttr& labelCenterVis = m->labelCenterVis();
    IntPoint3fAttr& parentCenterVis = m->parentCenterVis();

    if(m->meshType() == "MGXM" || m->meshType() == "MGXC"){
      neighborhood2D(S1, neighMap);
    }

    std::map<int, vertex> labelVtxMap;
    vvGraph &S = m2->graph();
    m2->reset();

    if(m->useParents()){
      forall(auto p, parentCenterVis){
        vertex v;
        v->pos = Point3d(p.second);
        v->label = m->parents()[p.first];
        S.insert(v);
        labelVtxMap[p.first] = v;
      }
    } else {
      forall(auto p, labelCenterVis){
        vertex v;
        v->pos = Point3d(p.second);
        v->label = p.first;
        S.insert(v);
        labelVtxMap[p.first] = v;
      }
    }


    forall(const IntIntDouble &p, neighMap){
      int label1 = p.first.first;
      int label2 = p.first.second;
      // if(label1 < label2) continue;
      if(m->useParents()){
        label1 = m->parents()[label1];
        label2 = m->parents()[label2];
      }
      if(p.second > minArea)
        S.insertEdge(labelVtxMap[label1], labelVtxMap[label2]);
    }

    m2->updateAll();
    m2->setShowSurface(false);
    m2->setMeshView("All");

    return true;
  }
  REGISTER_PROCESS(DisplayCellGraph2D);


  bool DisplayCellGraph3D::run(Mesh *m, Mesh *m2, double minArea)
  {

    std::map<IntInt, double> neighMap;
    vvGraph &S1 = m->graph();

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

    if(m->meshType() == "MGXM"){
      neighborhood2D(S1, neighMap);
      double tolerance = 0.001;
      NhbdGraphInfo info;
      neighborhoodGraph(S1, tolerance, info, false);
      neighMap = info.sharedWallArea;
    }

    std::map<int, vertex> labelVtxMap;
    vvGraph &S = m2->graph();
    m2->reset();

    // if(m->useParents()){
    //   forall(auto p, centroidsP){
    //     vertex v;
    //     v->pos = Point3d(p.second);
    //     v->label = m->parents()[p.first];
    //     S.insert(v);
    //     labelVtxMap[p.first] = v;
    //   }
    // } else {
      forall(auto p, centroids){
        vertex v;
        v->pos = Point3d(p.second);
        v->label = p.first;
        S.insert(v);
        labelVtxMap[p.first] = v;
      }
    // }


    forall(const IntIntDouble &p, neighMap){
      int label1 = p.first.first;
      int label2 = p.first.second;
      // if(label1 < label2) continue;
      // if(m->useParents()){
      //   label1 = m->parents()[label1];
      //   label2 = m->parents()[label2];
      // }
      if(p.second > minArea)
        S.insertEdge(labelVtxMap[label1], labelVtxMap[label2]);
    }

    m2->updateAll();
    m2->setShowSurface(false);
    m2->setMeshView("All");

    return true;
  }
  REGISTER_PROCESS(DisplayCellGraph3D);



 bool JoinCellsMGX3D::run(Mesh *m, int label1, int label2)
  {
    CellTissue &T = m->tissue();

    cell c1, c2;

    // check if both cells exist and neighbor each other
    bool valid = false;
    forall(const cell &c, T.C){
      if(c->label == label1){
        c1 = c;
        forall(const cell &d, T.C.neighbors(c)){
          if(d->label == label2){
            valid = true;
            c2 = d;
          }
        }
      }
    }

    if(!valid) return setErrorMessage("The specified cells cannot be joined! Check if they really exist and whether they share a common border.");

    //joinCellsMGX3D(c1, c2, T);

    vvGraph Temp;

    bool simple = false;

    // join them
    calcJoinedCell(c1, c2, Temp, simple);

    // add new cell, fix neighborhood and delete old cells
    //cell joinedCell = c1;
    //joinedCell->label = c1->label;

    calcCellGraphJoinedCells(T.C, c1, c2);

    // take cell 1 as joined cell, delete cell 2
    c1->S = Temp;
    //T.C.erase(c1);
    T.C.erase(c2);
    //T.C.insert(joinedCell);

    separateIdenticalVtxs(c1);

    T.updateAllGraph();
    T.updGeometry();


    m->updateAll();
    updateViewer();
    updateState();

    //const QString meshType = m->tissue().meshType();
    //Information::out << "meshtype after join " << meshType.toLocal8Bit().constData() << endl;

    return true;
  }
  REGISTER_PROCESS(JoinCellsMGX3D);


 bool MergeCellsMGXM3D::run(Mesh *m, int label1, int label2, bool forceMerge)
  {

    vvGraph& S = m->graph();

    int cell1 = -1;
    int cell2 = -1;
    // find the two cells (parameter or selection)
    if(label1 > 0 and label2 > 0){
      cell1 = label1;
      cell2 = label2;
    } else {
      // get selected cells
      std::set<int> selCells = findAllSelectedLabels(S);
      if(selCells.size() != 2) return setErrorMessage("Please select exactly two cells.");
      cell1 = *selCells.begin();
      cell2 = *(++selCells.begin());
    }

    std::cout << "Merge Cells " << cell1 << " and " << cell2 << std::endl;

    // get neighborhoodgraph

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

    std::map<int, std::vector<vertex> > cellVtxs = generateLabelVertexMap(S);

    // create a map from vtx to vtx

    //trivector tris1 = cellTriangles[cell1];
    //triVector tris2 = cellTriangles[cell2];

    double threshold = 0.0001;

    std::unordered_map<vertex, vertex> vtxVtxMap, vtxVtxMap2;

    // make a map of the duplicated vertices
    forall(vertex v, cellVtxs[cell1]){
      //find nearest in other cell, check threshold
      vertex minV;
      double minDis = HUGE_VAL;
      forall(vertex v2, cellVtxs[cell2]){
        double dis = norm(v->pos - v2->pos);
        if(dis < minDis){
          minDis = dis;
          minV = v2;
        }
      }

      if(minDis < threshold){
        vtxVtxMap[v] = minV;
        vtxVtxMap2[minV] = v;
      }
    }

    triVector newTris;
    triVector trisShared;
    std::set<vertex> sharedVtxs;



    // get shared triangles between the two cells
    forall(Triangle t, cellTriangles[cell1]){

      vertex v1 = t.v[0];
      vertex v2 = t.v[1];
      vertex v3 = t.v[2];

      if(vtxVtxMap.find(v1) == vtxVtxMap.end() or vtxVtxMap.find(v2) == vtxVtxMap.end() or vtxVtxMap.find(v3) == vtxVtxMap.end()){
        // tri to be kept
        newTris.push_back(t);
      } else {
        // tri shared, to be deleted
        trisShared.push_back(t);
        sharedVtxs.insert(v1);
        sharedVtxs.insert(v2);
        sharedVtxs.insert(v3);
      }

    }




    // get shared triangles between the two cells
    forall(Triangle t, cellTriangles[cell2]){

      vertex v1 = t.v[0];
      vertex v2 = t.v[1];
      vertex v3 = t.v[2];

      if(vtxVtxMap2.find(v1) == vtxVtxMap2.end() or vtxVtxMap2.find(v2) == vtxVtxMap2.end() or vtxVtxMap2.find(v3) == vtxVtxMap2.end()){
        // tri to be kept
        newTris.push_back(t);
      } else {
        // tri shared, to be deleted
        trisShared.push_back(t);
        sharedVtxs.insert(vtxVtxMap2[v1]);
        sharedVtxs.insert(vtxVtxMap2[v2]);
        sharedVtxs.insert(vtxVtxMap2[v3]);
      }

    }

    if(trisShared.empty()) return setErrorMessage("The two cells don't share triangles!");

    std::vector<Point3i> triVec;
    std::vector<vertex> vtxVec;

    std::map<vertex, int> vtxPosMap;

    forall(Triangle t, newTris){
      vertex v1 = t.v[0];
      vertex v2 = t.v[1];
      vertex v3 = t.v[2];

      if(vtxVtxMap.find(v1) != vtxVtxMap.end() and sharedVtxs.find(v1) != sharedVtxs.end()) v1 = vtxVtxMap[v1];
      if(vtxVtxMap.find(v2) != vtxVtxMap.end() and sharedVtxs.find(v2) != sharedVtxs.end()) v2 = vtxVtxMap[v2];
      if(vtxVtxMap.find(v3) != vtxVtxMap.end() and sharedVtxs.find(v3) != sharedVtxs.end()) v3 = vtxVtxMap[v3];

      if(vtxPosMap.find(v1) == vtxPosMap.end()){
        vtxPosMap[v1] = vtxVec.size();
        vtxVec.push_back(v1);
      }
      if(vtxPosMap.find(v2) == vtxPosMap.end()){
        vtxPosMap[v2] = vtxVec.size();
        vtxVec.push_back(v2);
      }
      if(vtxPosMap.find(v3) == vtxPosMap.end()){
        vtxPosMap[v3] = vtxVec.size();
        vtxVec.push_back(v3);
      }

      v1->label = cell1;
      v2->label = cell1;
      v3->label = cell1;


      triVec.push_back(Point3i(vtxPosMap[v1],vtxPosMap[v2],vtxPosMap[v3]));

    }


    vvGraph Temp = S;

    // erase all vertices of the cells
    forall(vertex v, cellVtxs[cell1]){
      S.erase(v);
    }
    // erase all vertices of the cells
    forall(vertex v, cellVtxs[cell2]){
      S.erase(v);
    }

    //assembleCellFromWalls(S, vtxVec, triVec, sharedVtxs, threshold, 0, false);



    // add only the non-shared triangles back to the mesh and label the merged cell
    if(!Mesh::meshFromTriangles(S, vtxVec, triVec)){
      if(!forceMerge) S = Temp;
      m->updateAll();
      return setErrorMessage("Merge Error occured!");
    }

    forall(Triangle t, newTris){

      vertex v1 = t.v[0];
      vertex v2 = t.v[1];
      vertex v3 = t.v[2];

      v1->label = cell1;
      v2->label = cell1;
      v3->label = cell1;
    }

    // update the mesh
    m->updateAll();


    return true;
  }
  REGISTER_PROCESS(MergeCellsMGXM3D);


 // bool JoinCellsMGX2D::run(Mesh *m, int label1, int label2)
 //  {

 //    vvGraph& S = m->graph();
 //    CellTissue& T = m->tissue();


 //    if(!T.toMgxm())
 //      throw(QString("Unexpected error converting to normal (MGXM) mesh"));

 //    vertex v1, v2;
 //    bool found1 = false;
 //    bool found2 = false;

 //    std::set<vertex> border;
 //    std::set<vertex> commonVert, endOfBorder;

 //    forall(const vertex& v, S){
 //      if(v->label == label1){ v1 = v; found1 = true; }
 //      if(v->label == label2){ v2 = v; found2 = true; }
 //    }

 //    if(!found1 or !found2) throw(QString("Couldnt find the cells"));

 //    forall(const vertex& n, S.neighbors(v1)){
 //      forall(const vertex& n2, S.neighbors(v2)){
 //        if(n==n2) border.insert(n);
 //      }
 //    }

 //    forall(const vertex& n, border){
 //      bool found = false;
 //      forall(const vertex& m, S.neighbors(n)){
 //        if((m->label > -1 and m->label != label1 and m->label != label2) or n->margin){
 //          endOfBorder.insert(n);
 //          found = true;
 //        }
 //      }
 //      if(!found) commonVert.insert(n);
 //    }

 //    if(endOfBorder.size() == 0) return setErrorMessage("The specified cells cannot be joined! Check if they really exist and whether they share a common border.");

 //    if(commonVert.size() > 0) {
 //      // just merge the vtx
 //      forall(const vertex& v, commonVert){
 //        v->selected = true;
 //      }
 //    } else {
 //      // split the border and merge
 //      Subdivide *sDiv = 0;
 //      vertex b1 = *(endOfBorder.begin());
 //      vertex b2 = *(++endOfBorder.begin());

 //      b1->selected = true;
 //      b2->selected = true;
 //      bool res;
 //      if(S.uniqueLine(b1, b2)) res = splitEdge(S, b1, b2, sDiv);
 //      else res = splitEdge(S, b2, b1, sDiv);
 //      if(!res) std::cout << "sth wrong" << std::endl;

 //      b1->selected = false;
 //      b2->selected = false;
 //    }

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

 //    MergeVertices mVtx(*this);
 //    mVtx.run(m);

 //    m->updateAll();
 //    updateViewer();
 //    updateState();

 //    if(!T.toMgx2d(-1))
 //      throw(QString("Unexpected error converting to MGX2D mesh"));

 //    // consistently label the resulting cell to label1
 //    cell newCell;
 //    int maxL;
 //    if(findCell(T.C, newCell, label2, maxL)) newCell->label = label1;

 //    m->updateAll();
 //    updateViewer();
 //    updateState();

 //    return true;
 //  }
 //  REGISTER_PROCESS(JoinCellsMGX2D);

 bool FlatDivisionPlane::run(Mesh *m, Mesh *m2, int label1, int label2, bool display, DivisionAnalysisData &dat, CellDivisionAttr& cdaExp)
  {

    AttrMap<int, CellDivisionAttr>& attrData = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Actual Plane");

    int parentLabel1 = m->parents()[label1];
    int parentLabel2 = m->parents()[label2];

    if(parentLabel1 != parentLabel2) return setErrorMessage("Cells don't have the same parent labels.");

    CellTissue &T = m->tissue();

    cell c1, c2;

    // check if both cells exist and neighbor each other
    bool valid = false;
    forall(const cell &c, T.C){
      if(c->label == label1){
        c1 = c;
        forall(const cell &d, T.C.neighbors(c)){
          if(d->label == label2){
            valid = true;
            c2 = d;
          }
        }
      }
    }

    if(!valid) return setErrorMessage("cells not valid!");

    std::vector<Point3d> pointsPolygon = getSharedVtxs(c1,c2);
    Point3d planePos, planeNrml;
    findPolygonPlane(pointsPolygon, planePos, planeNrml);

    dat.divPlaneFlat = DivisionPlane(15, planeNrml, planePos, 9);

    std::set<Triangle> ts = getSharedTris(c1,c2);
    double actualArea = calcTriArea(ts);

    std::cout << "Flat " << ts.size() << "/" << actualArea << std::endl;

    CellDivisionAttr cda;
    cda.planeNrml = planeNrml;
    cda.planeDisplacement = 0; // TODO
    cda.planePos = planePos;
    cda.sharedTrisArea = actualArea;
    attrData[parentLabel1] = cda;
    cdaExp = cda;

    if(display){
      drawPlane(*this, m2, dat.divPlaneFlat);
      m2->updateAll();
    }

    return true;
  }
  REGISTER_PROCESS(FlatDivisionPlane);



 bool ActualPlaneStack::run(Stack* s1, Store* store1, Mesh *m, Mesh *m2, int label1, int label2, bool display, DivisionAnalysisData &dat, CellDivisionAttr &cdaExp)
  {

    AttrMap<int, CellDivisionAttr>& attrData = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Actual Plane");
    //attrData.clear();

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

    // go through the stack and find voxels of the two cells
    int parentLabel1 = m->parents()[label1];
    int parentLabel2 = m->parents()[label2];

    if(parentLabel1 != parentLabel2) return setErrorMessage("Cells don't have the same parent labels.");

    CellTissue &T = m->tissue();

    cell c1, c2;

    // check if both cells exist and neighbor each other
    bool valid = false;
    forall(const cell &c, T.C){
      if(c->label == label1){
        c1 = c;
        forall(const cell &d, T.C.neighbors(c)){
          if(d->label == label2){
            valid = true;
            c2 = d;
          }
        }
      }
    }

    if(!valid) return setErrorMessage("cells not valid!");


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

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

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

    int xDim = imgSize.x();
    int yDim = imgSize.y();
    int zDim = imgSize.z();

    Point3f voxelSizes = s1->step();

    double voxelVolume = voxelSizes.x() * voxelSizes.y() * voxelSizes.z();

    //Point3i lower(0,0,0);
    //Point3i upper(xDim-1, yDim-1, zDim-1);

    //double totalSignal = 0;

    int neighborhoodSize = 1;

    std::vector<Point3d> borderPoints;

    //#pragma omp parallel for schedule(guided)
    for(int z=neighborhoodSize; z<zDim-neighborhoodSize; z++){
      for(int x=neighborhoodSize; x<xDim-neighborhoodSize; x++){
        for(int y=neighborhoodSize; y<yDim-neighborhoodSize; y++){
          Point3i ps(x,y,z);
          int label = data[s1->offset(ps)];
          if(label != label1 and label != label2) continue;

          bool border = false;

          for(int zn=-neighborhoodSize; zn<neighborhoodSize; zn++){
            for(int xn=-neighborhoodSize; xn<neighborhoodSize; xn++){
              for(int yn=-neighborhoodSize; yn<neighborhoodSize; yn++){
                Point3i psNew(x+xn,y+yn,z+zn);
                int labelNew = data[s1->offset(psNew)];

                if(labelNew != label1 and labelNew != label2) continue;
                if(labelNew == label) continue;

                border = true;


              }
            }
          }
          if(border)
            borderPoints.push_back(Point3d(s1->imageToWorld(ps)));
        }
      }
    }

    // forall(Point3d p, borderPoints){
    //   vertex v;
    //   v->pos = p;
    //   S2.insert(v);
    // }

    // m2->updateAll();


    if(borderPoints.empty()) return setErrorMessage("No shared wall found");

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

    dat.divPlaneFlat = DivisionPlane(15, planeNrml, planePos, 9);

    // std::set<Triangle> ts = getSharedTris(c1,c2);
    // double actualArea = calcTriArea(ts);

    std::cout << "Flat " << planeNrml << "/" << planePos << std::endl;

    std::set<Triangle> ts = getSharedTris(c1,c2);
    double actualArea = calcTriArea(ts);

    CellDivisionAttr cda;
    cda.planeNrml = planeNrml;
    cda.planeDisplacement = 0; // TODO
    cda.planePos = planePos;
    cda.divPlane.pos = planePos;
    cda.divPlane.nrml = planeNrml;
    cda.sharedTrisArea = actualArea;
    attrData[parentLabel1] = cda;
    cdaExp = cda;

    if(display){
      drawPlane(*this, m2, dat.divPlaneFlat);
      m2->updateAll();
    }

    return true;
  }
  REGISTER_PROCESS(ActualPlaneStack);


 bool AnalyzeSurface::run(Mesh *m, Mesh *m2, DivisionAnalysisData &dat)
  {
    CellTissue &T = m->tissue();
    const vvGraph& S2 = m2->graph();

    AttrMap<int, Point3d>& surfaceDirections = m->attributes().attrMap<int, Point3d>("Surface Directions");
    surfaceDirections.clear();

    std::map<int, Point3d> newCellNormalMap;

    // forall cells: find the nearest surface vertex
    forall(const cell& c, T.C){
      vertex minV;
      Point3d perpDir (0,0,0);
      T.updGeometry(c, 4);
      double minDis = 1E20;
      forall(const vertex &v, S2){
        double dis = norm(v->pos - c->pos);
        if(dis < minDis){
          minDis = dis;
          minV = v;
        }
      }
      perpDir = minV->nrml;
      newCellNormalMap[c->label] = perpDir/norm(perpDir);
      surfaceDirections[c->label] = newCellNormalMap[c->label];
      minV->selected = true;
    }

    dat.cellNormalMap = newCellNormalMap;

/*
    // calculate betweenness centrality
    if(betw){
      if(weights){
        std::map<IntInt, double> neighborMap;
        T.getNeighborhoodMap3D(neighborMap);
        neighborMap = generateWeights(neighborMap, false);
        calcBetweenness(T.C, neighborMap, dat.cellBetwMap);
      } else {
        calcBetweenness(T.C, dat.cellBetwMap);
      }
    }
    // calculate current flow betweenness
    if(currentFlow){
      if(weightsCur){
        std::map<IntInt, double> neighborMap;
        T.getNeighborhoodMap3D(neighborMap);
        neighborMap = generateWeights(neighborMap, !weightsCur);
        calcRandWalkBetweenness(T.C, neighborMap, dat.cellCurrentFlowMap, useGPU, false);
      } else {
        std::map<IntInt, double> equalWeights;
        calcRandWalkBetweenness(T.C, equalWeights, dat.cellCurrentFlowMap, useGPU, false);
      }
    }
*/
    return true;
  }
  REGISTER_PROCESS(AnalyzeSurface);




 bool SetDivisionPoint::run(Mesh *m)
  {

    AttrMap<int, vertex>& divPoint = m->attributes().attrMap<int, vertex>("Division Analysis Custom Point");

    std::vector<vertex> selVtxs = m->selectedVertices();

    if(selVtxs.size() > 1) return setErrorMessage("Please select exactly one vertex!");

    divPoint[0] = selVtxs[0];

    return true;
  }
  REGISTER_PROCESS(SetDivisionPoint);


 bool DisplayOptimalPlanes::run(Mesh *m, Mesh *m2, bool custom, QString primary, QString secondary,
  QString primaryDir, QString secondaryDir, DivisionAnalysisData &dat, double tolerance, int maxPlanes, bool resetM2, double planeSim, double planeSize)
  {

    //cout << "display" << endl;
    if(resetM2) m2->reset();

    int nrOfPlanes = 0;

    if(custom){
      forall(DivisionPlane d, dat.customPlanes){
        drawPlane(*this, m2, d);
        nrOfPlanes++;
        if(nrOfPlanes>=maxPlanes) return true;
      }
    } else {

      std::vector<CellDivision> optPlanes;
      double optValue;
      P3dDouble bestDivPlane;

      // first sort for the primary criterium
      // then take best results and sort for secondary criterium
      findBestPlane(dat.divData, optPlanes, optValue, tolerance, bestDivPlane, primary, secondary, primaryDir, secondaryDir, planeSim);
      std::cout << "best " << bestDivPlane.first << "/" << bestDivPlane.second << "//" << optValue << endl;

      std::map<double, CellDivision> optPlanesSorted;
      std::pair<double, CellDivision> DoubleCDP;
      forall(CellDivision cd, optPlanes){
        double value = cd.getValue(primary);
        optPlanesSorted[value] = cd;
      }


      for(auto it = optPlanesSorted.begin();it != optPlanesSorted.end(); it++){
        if(nrOfPlanes>=maxPlanes) continue;
        CellDivision cd = it->second;
        nrOfPlanes++;
        DivisionPlane dp = DivisionPlane(planeSize, cd.divPlane.nrml, cd.divPlane.pos, cd.divPlane.label);
        m2->labelHeat()[cd.divPlane.label] = it->first;

        if(norm(cd.cellNormal) < 1E-5) drawPlane(*this, m2, dp);
        else drawPlane(*this, m2, dp, cd.cellNormal, 0.92); // add some scaling for 2D planes
        //drawPlane(*this, m2, dp);
        m2->heatMapBounds() = m2->calcHeatMapBounds();
        m2->setShowLabel("Label Heat");
        m2->updateAll();
        std::cout << "plane #" << cd.divPlane.label << ": " << cd.divPlane.nrml << "/" << cd.divPlane.pos << "/ values /" << cd.getValue(primary) << "/" << cd.getValue(secondary) << endl;
      }
/*
      forall(CellDivision cd, optPlanes){
        DivisionPlane dp = DivisionPlane(planeSize, cd.divPlane.nrml, cd.divPlane.pos, 1);
        drawPlane(*this, m2, dp);
        nrOfPlanes++;
        std::cout << "plane #" << nrOfPlanes << ": " << cd.divPlane.nrml << "/" << cd.divPlane.pos << endl;
        if(nrOfPlanes>=maxPlanes) return true;
      }
      */
    }

    return true;
  }
  REGISTER_PROCESS(DisplayOptimalPlanes);


 bool DisplayOptimalPlanesSimple::run(Mesh *m, Mesh *m2, QString divRule, double maxDis, int maxPlanes, double planeSim, double planeSize, bool resetMesh)
  {
    DisplayOptimalPlanes dop(*this);
    if(divRule == "Shortest Wall"){
      return dop.run(m,m2,false,"Wall Area","None","Min","Min",dat,maxDis,maxPlanes,resetMesh,planeSim,planeSize);
    } else if(divRule == "Best Volume Ratio"){
      return dop.run(m,m2,false,"Volume Ratio","None","Min","Min",dat,maxDis,maxPlanes,resetMesh,planeSim,planeSize);
    } else if(divRule == "Longest Wall"){
      return dop.run(m,m2,false,"Wall Area","None","Max","Min",dat,maxDis,maxPlanes,resetMesh,planeSim,planeSize);
    } else {

    }

    return true;
  }
  REGISTER_PROCESS(DisplayOptimalPlanesSimple);


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

    m2->reset();

    forall(CellDivision cd, dat.divData){
      DivisionPlane dp = DivisionPlane(10., cd.divPlane.nrml, cd.divPlane.pos, 1);
      drawPlane(*this, m2, dp);
    }

    return true;
  }
  REGISTER_PROCESS(DisplayAllPlanes);

 bool AnalysisDivPlanes::run(Mesh *m, Mesh *m2, bool flat, bool shortest, bool shortestBest, DivisionAnalysisData &dat)
  {

    m2->reset();

    if(flat){
      drawPlane(*this, m2, dat.divPlaneFlat);
    }
    if(shortest){
      drawPlane(*this, m2, dat.divPlaneShortest);
    }
    if(shortestBest){
      drawPlane(*this, m2, dat.divPlaneShortestBest);
    }

    return true;
  }
  REGISTER_PROCESS(AnalysisDivPlanes);

 // bool DisplayHeatMaps::run(QString choice)
 //  {

 //    Mesh *m1 = currentMesh();// mesh(0);

 //    std::map<cell, double> heatMap;

 //    if(choice == "Betweenness Centrality"){
 //      heatMap = dat.cellBetwMap;
 //    } else if(choice == "Current Flow Betweenness"){
 //      heatMap = dat.cellCurrentFlowMap;
 //    } else {
 //      //option = -1;
 //      //warning << QString("Invalid type '%1'").arg(outputType);
 //    }

 //    double maxValue = -1E20, minValue = 1E20;

 //    m1->wallHeat().clear();
 //    m1->labelHeat().clear();
 //    forall(const cellDPair& p, heatMap){
 //      m1->labelHeat()[p.first->label] = p.second;
 //      if(maxValue < p.second) maxValue = p.second;
 //      if(minValue > p.second) minValue = p.second;
 //    }

 //    m1->heatMapBounds() = Point2f(minValue, maxValue);
 //    m1->updateTriangles();
 //    m1->showHeat();
 //   // SETSTATUS(QString("Loaded heat map with values ranging from %1 to %2").arg(minValue).arg(maxValue));

 //    return true;
 //  }
 //  REGISTER_PROCESS(DisplayHeatMaps);



    typedef std::pair<cell, cell> cellPair;
    typedef std::pair<cell, cellPair> cellTriplet;

    cellTriplet createCellTriplet(const cell& c1, const cell& c2, const cell& c3)
    {
      cellTriplet ct;

      if(c1 < c2 and c1 < c3){
        ct.first = c1;
        ct.second = std::make_pair(c2,c3);
        if(c3<c2) ct.second = std::make_pair(c3,c2);
      } else if(c2 < c1 and c2 < c3){
        ct.first = c2;
        ct.second = std::make_pair(c1,c3);
        if(c3<c1) ct.second = std::make_pair(c3,c1);
      } else if(c3 < c2 and c3 < c1){
        ct.first = c3;
        ct.second = std::make_pair(c1,c2);
        if(c2<c1) ct.second = std::make_pair(c2,c1);
      }

      return ct;
    }

    cellPair createCellPair(const cell& c1, const cell& c2)
    {
      cellPair cp;

      if(c1 < c2){
        cp = std::make_pair(c1,c2);
      } else {
        cp = std::make_pair(c2,c1);
      }
      return cp;
    }

 // create a map of label,label and their shared wall Area
 void getNeighborhoodMapOutside3D(CellTissue& T, std::map<IntInt, double>& neighborMap, std::map<int, double>& outsideCellWall)
 {

  typedef std::pair<vertex, vertex> VtxVtx;
  typedef std::pair<VtxVtx, vertex> Tri;

  std::map<IntInt, double> result;
  std::map<int, double> resultOut;

  std::map<Tri, std::vector<int> > triLabelsMap;
  typedef std::pair<Tri, std::vector<int> > TriLabelPair;

  // go through all cells
  forall(const cell& c, T.C){
    // write triangles of cell to a map
    forall(const vertex &v, c->S){
      forall(const vertex &n, c->S.neighbors(v)){
        vertex m = c->S.nextTo(v,n);
        if(T.S.uniqueTri(v,n,m)){
          VtxVtx vv = std::make_pair(v,n);
          Tri tt = std::make_pair(vv,m);
          triLabelsMap[tt].push_back(c->label);
        }
        vertex m2 = c->S.prevTo(v,n); // needed as triangles are oriented the other way in neighboring cells
        if(T.S.uniqueTri(v,n,m2)){
          VtxVtx vv = std::make_pair(v,n);
          Tri tt = std::make_pair(vv,m2);
          triLabelsMap[tt].push_back(c->label);
        }
      }
    }
  }

  forall(const TriLabelPair &p, triLabelsMap){
    //if(p.second.size()>2) std::cout << "sth weird happened" << std::endl;
    if(p.second.size()==2){
      std::pair<int, int> neighborPair1 = std::make_pair(p.second[0], p.second[1]);
      std::pair<int, int> neighborPair2 = std::make_pair(p.second[1], p.second[0]);
      double triArea = triangleArea(p.first.first.first->pos, p.first.first.second->pos, p.first.second->pos)/2.; // divide by 2 as everything was counted twice
      result[neighborPair1] += triArea;
      result[neighborPair2] += triArea;
    } else if(p.second.size()==1){
      double triArea = triangleArea(p.first.first.first->pos, p.first.first.second->pos, p.first.second->pos)/2.; // divide by 2 as everything was counted twice
      resultOut[p.second[0]] += triArea;
    }
  }

  neighborMap = result;
  outsideCellWall = resultOut;

  }

 bool OutsideCellWall::run(Mesh *m1)
 {

    CellTissue& T = m1->tissue();

    std::map<IntInt, double> neighbMap;
    std::map<int, double> outsMap;

    getNeighborhoodMapOutside3D(T, neighbMap, outsMap);

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

    std::map<int, double> cellWallAreaMap;

    forall(const cell& c, T.C){
      // write triangles of cell to a map
      forall(const vertex &v, c->S){
        forall(const vertex &n, c->S.neighbors(v)){
          vertex m = c->S.nextTo(v,n);
          if(T.S.uniqueTri(v,n,m)){
            double triArea = triangleArea(v->pos, n->pos, m->pos); // divide by 2 as everything was counted twice
            cellWallAreaMap[c->label]+=triArea;
          }
        }
      }
    }

    forall(IntDouble p, outsMap){
      m1->labelHeat()[p.first] = p.second/cellWallAreaMap[p.first];
    }

    m1->heatMapBounds() = m1->calcHeatMapBounds();
    m1->showHeat();
    m1->updateTriangles();

    return true;
 }
   REGISTER_PROCESS(OutsideCellWall);

  bool orderpolygonVtx(std::set<vertex>& vtxs)
  {
    std::set<vertex> vtxsLeftToDo;

    forall(const vertex& v, vtxs){
      vtxsLeftToDo.insert(v);
    }

    while(!vtxsLeftToDo.empty()){

      const vertex& v = *(vtxsLeftToDo.begin());

      vtxsLeftToDo.erase(v);

    }

    return true;
  }


 // file dialogue for save data
 bool CreateJunctionGraph::initialize(QWidget* parent)
  {
    QString filename = parm("Save to file");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Save to file", filename);
    return true;
  }

    typedef std::pair<vertex, int> VtxIntP;


  struct JunctionType{

  std::set<int> neighborCells;
  int count;
  Point3d pos;
  vertex v;

};


  typedef std::pair<std::set<int>, JunctionType> JunctionInfo;
  typedef std::pair<int,int> IntInt;
    typedef std::pair<IntInt, std::set<vertex> > IntIntVtxSetP;
    typedef std::pair<Point3i, std::set<vertex> > P3iVtxSetP;


 bool writeJunctionGraphToFiles(QString filename, std::map<vertex, int>& vtxIdxMap, std::vector<std::pair<vertex,vertex> >& edgeList,
  std::map<vertex, JunctionType>& vtxsJunctions)
 {

     std::cout << "Writing Files.....";

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

    out << "Junction Vertex ID,Pos X,Pos Y,Pos Z,Adjacent Cells"<< endl;


    forall(VtxIntP p, vtxIdxMap){
      out << p.second << "," << p.first->pos.x() << "," << p.first->pos.y() << ","  << p.first->pos.z() << ",";
      forall(int c, vtxsJunctions[p.first].neighborCells)
        out << c << ",";
      out << endl;
    }

    filename.remove(".csv", Qt::CaseInsensitive);
    filename += "_edges";
    filename += ".csv";

    QFile file2(filename);
    if(!file2.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream out2(&file2);

    out2 << "Wall Edge ID,Source Junction Vertex ID,Target Junction Vertex ID,Length"<< endl;

    typedef std::pair<vertex, vertex> VtxVtxP;
    int counter = 1;
    forall(VtxVtxP p, edgeList){
      out2 << counter << "," << vtxIdxMap[p.first] << "," << vtxIdxMap[p.second] << "," << norm(p.first->pos-p.second->pos) << endl;
      counter++;
    }
  std::cout << "....Done!" << std::endl;
  return true;
 }


  std::set<std::pair<int,int> > findAllCellPairs(const std::set<int>& junctionCells)
  {
    std::set<std::pair<int,int> > cellPairsVec;

    forall(int c1, junctionCells){
      forall(int c2, junctionCells){
        if(c1 == c2) continue;
        IntInt p = std::make_pair(c1, c2);
        if(c2<c1) p = std::make_pair(c2, c1);
        cellPairsVec.insert(p);
      }
    }

    return cellPairsVec;
  }

  std::set<Point3i> findAllCellTriplets(const std::set<int>& junctionCells)
  {
    std::set<Point3i> cellTripletVec;

    forall(int c1, junctionCells){
      forall(int c2, junctionCells){
        if(c1 == c2) continue;
        forall(int c3, junctionCells){
          if(c1 == c3 or c2 == c3) continue;

          Point3i triplet(c1,c2,c3); // 123
          if(c2 < c1){ // 21x, x21 or 2x1
            if(c1 < c3) triplet = Point3i(c2,c1,c3);
            else if(c3 < c2) triplet = Point3i(c3,c2,c1);
            else triplet = Point3i(c2,c3,c1);
          } else if(c3 < c1){
            triplet = Point3i(c3,c1,c2);
          } else {
            if(c3 < c2) triplet = Point3i(c1,c3,c2);
          }
          cellTripletVec.insert(triplet);
        }
      }
    }

    return cellTripletVec;
  }

  void getRedundantTuple5(const std::set<int>& tuple6, std::vector<std::set<int> >& tuples5)
  {

    std::set<int>::iterator it = tuple6.begin();

    int i1,i2,i3,i4,i5,i6;
    i1 = *(it);
    it++;
    i2 = *(it);
    it++;
    i3 = *(it);
    it++;
    i4 = *(it);
    it++;
    i5 = *(it);
    it++;
    i6 = *(it);
    it++;

    std::set<int> no1,no2,no3,no4,no5,no6;
    no1.insert(i2);no1.insert(i3);no1.insert(i4);no1.insert(i5);no1.insert(i6);
    no2.insert(i1);no2.insert(i3);no2.insert(i4);no2.insert(i5);no2.insert(i6);
    no3.insert(i2);no3.insert(i1);no3.insert(i4);no3.insert(i5);no3.insert(i6);
    no4.insert(i2);no4.insert(i3);no4.insert(i1);no4.insert(i5);no4.insert(i6);
    no5.insert(i2);no5.insert(i3);no5.insert(i4);no5.insert(i1);no5.insert(i6);
    no6.insert(i1);no6.insert(i2);no6.insert(i3);no6.insert(i4);no6.insert(i5);

    tuples5.push_back(no1);
    tuples5.push_back(no2);
    tuples5.push_back(no3);
    tuples5.push_back(no4);
    tuples5.push_back(no5);
    tuples5.push_back(no6);

  }

  void getRedundantTuple4(const std::set<int>& tuple5, std::vector<std::set<int> >& tuples4)
  {

    std::set<int>::iterator it = tuple5.begin();

    int i1,i2,i3,i4,i5;
    i1 = *(it);
    it++;
    i2 = *(it);
    it++;
    i3 = *(it);
    it++;
    i4 = *(it);
    it++;
    i5 = *(it);
    it++;

    std::set<int> no1,no2,no3,no4,no5;
    no1.insert(i2);no1.insert(i3);no1.insert(i4);no1.insert(i5);
    no2.insert(i1);no2.insert(i3);no2.insert(i4);no2.insert(i5);
    no3.insert(i2);no3.insert(i1);no3.insert(i4);no3.insert(i5);
    no4.insert(i2);no4.insert(i3);no4.insert(i1);no4.insert(i5);
    no5.insert(i2);no5.insert(i3);no5.insert(i4);no5.insert(i1);

    tuples4.push_back(no1);
    tuples4.push_back(no2);
    tuples4.push_back(no3);
    tuples4.push_back(no4);
    tuples4.push_back(no5);

  }

  void deleteRedundantJunctions45(std::map<std::set<int>, JunctionType>& junctionsMap)
  {

    std::set<std::set<int> > toBeDeleted6;
    forall(const JunctionInfo& j, junctionsMap){
      if(j.first.size() != 6) continue;

      std::vector<std::set<int > > tuples5;
      getRedundantTuple5(j.first, tuples5);
      for(size_t i = 0; i<tuples5.size(); i++){
        toBeDeleted6.insert(tuples5[i]);
        std::vector<std::set<int > > tuples4;
        getRedundantTuple4(j.first, tuples4);
        for(size_t j = 0; j<tuples4.size(); j++){
          toBeDeleted6.insert(tuples4[j]);
        }

      }
    }

    forall(std::set<int> s, toBeDeleted6){
      if(junctionsMap.find(s) != junctionsMap.end()){
        junctionsMap.erase(s);
      }
    }

    // find 5-way junctions and kick out redundant 4-way junctions
    std::set<std::set<int> > toBeDeleted;

    forall(const JunctionInfo& j, junctionsMap){
      if(j.first.size() != 5) continue;

      std::vector<std::set<int > > tuples4;
      getRedundantTuple4(j.first, tuples4);
      for(size_t i = 0; i<tuples4.size(); i++)
        toBeDeleted.insert(tuples4[i]);

    }

    forall(std::set<int> s, toBeDeleted){
      if(junctionsMap.find(s) != junctionsMap.end()){
        junctionsMap.erase(s);
      }
    }

  }

 void buildJunctionGraph(vvGraph& S,
  std::map<std::set<int>, JunctionType>& junctionsMap,
  std::map<vertex, JunctionType>& vtxsJunctions,
  std::vector<std::pair<vertex,vertex> >& edgeList,
  std::map<vertex, int >& vtxIdxMap,
  std::map<IntInt, std::set<vertex> >& cellPairJunctionMap)
  {

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

   // std::vector<int> counter(30);

    // calculated junction pos and add them to vvGraph
    forall(const JunctionInfo& j, junctionsMap){
      //vertex v;
      //v->pos = j.second.pos;
      S.insert(j.second.v);
      vtxIdxMap[j.second.v] = vtxIdxMap.size();
      vtxsJunctions[j.second.v] = j.second;
    }
    //std::cout << "mesh " << S.size() << std::endl;
    //std::map<IntInt, std::set<vertex> > cellPairJunctionMap;

    std::map<Point3i, std::set<vertex> > cellTripletJunctionMap;

    // go through junctions, create map cellPair->junction, connect them
    forall(const JunctionInfo& j, junctionsMap){
      std::set<IntInt> cellPairs = findAllCellPairs(j.first);
      std::set<Point3i> cellTriplet = findAllCellTriplets(j.first);
      forall(const IntInt& p, cellPairs){
        cellPairJunctionMap[p].insert(junctionsMap[j.first].v);
      }
      forall(const Point3i& p, cellTriplet){
        cellTripletJunctionMap[p].insert(junctionsMap[j.first].v);
      }
    }

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

    //forall(const IntIntVtxSetP& p, cellPairJunctionMap){
    forall(const P3iVtxSetP& p, cellTripletJunctionMap){
      //if(p.first.x() != 1457 and p.first.y()!= 1457 and p.first.z()!= 1457) continue;
      forall(const vertex& v, p.second){
        forall(const vertex& n, p.second){
          if(v == n) continue;
          S.insertEdge(v,n);
          std::pair<vertex,vertex> vv = std::make_pair(v,n);
          edgeList.push_back(vv);
        }
      }
    }
    //std::cout << "edges " << edgeList.size() << std::endl;
  }

  void buildMeshFromJunctions(vvGraph& S, vvGraph& S2, std::set<int>& allLabels,
    std::map<IntInt, std::set<vertex> >& cellPairJunctionMap,
    std::map<int, Point3d>& labelPosMap)
  {

    std::unordered_map<int, MeshBuilder> labelMeshBMap;

    forall(int l, allLabels){
      MeshBuilder mb;
      labelMeshBMap[l] = mb;
    }

    forall(const IntIntVtxSetP& p, cellPairJunctionMap){

      int c1 = p.first.first;
      int c2 = p.first.second;

      //if(c2 !=1457 and c2 != 1457) continue;

      // order the vertices
      std::set<vertex> vtxs = p.second;
      //std::cout << "vtx " << c1 << "/" << c2 << "//" << vtxs.size() << std::endl;

      // create segments
      std::vector<std::pair<Point3d,Point3d> > polygonSegs;
      std::map<Point3d, vertex> P3dVertexMap;
      std::set<std::pair<vertex,vertex> > segsDone;
      forall(const vertex& v, vtxs){
        forall(const vertex& n, S.neighbors(v)){
          if(vtxs.find(n) == vtxs.end()) continue;
          std::pair<vertex,vertex> vv = std::make_pair(v,n);
          if(n<v) vv = std::make_pair(n,v);
          if(segsDone.find(vv) != segsDone.end()) continue;

          segsDone.insert(vv);
          polygonSegs.push_back(std::make_pair(v->pos,n->pos));
          P3dVertexMap[v->pos] = v;
          P3dVertexMap[n->pos] = n;

        }
      }

      // order the segments
      //std::cout << "segs " << c1 << "/" << c2 << "//" << polygonSegs.size() << std::endl;
      if(polygonSegs.size() < 3) continue;
      std::vector<std::vector<Point3d> > orderedSegs = orderPolygonSegsMulti(polygonSegs);

      forall(std::vector<Point3d>& wall, orderedSegs){

        // add a vertex to the center
        Point3d centerPos(0,0,0);
        forall(Point3d& p, wall)
          centerPos+=p;
        centerPos/=wall.size();

        vertex center;
        center->pos = centerPos;

        Point3d prev = wall[wall.size()-1];
        forall(Point3d& p, wall){
          // add the triangles with right orientation
          Point3d dirC1 = labelPosMap[c1]-centerPos;
          double s1 = ((p - prev) ^ (centerPos - prev))*(dirC1);
          vertex w1,w2,w3;
          w1->pos = prev;
          w2->pos = p;
          w3->pos = centerPos;

          vertex w4,w5,w6;
          w4->pos = prev;
          w5->pos = p;
          w6->pos = centerPos;

          if(s1 > 0){
            labelMeshBMap[c1].addTri(w2, w1, w3);
            labelMeshBMap[c2].addTri(w4, w5, w6);
          } else {
            labelMeshBMap[c1].addTri(w1, w2, w3);
            labelMeshBMap[c2].addTri(w5, w4, w6);
          }
          prev = p;
        }

      }

    }
      //m2->reset();
      //S2.clear();

    forall(int c, allLabels){
      //if(c != 1457) continue;
      if(labelMeshBMap[c].triVec.size() == 0) continue;

      forall(const vertex& v, labelMeshBMap[c].vtxVec){
        v->label = c;
      }
      //std::cout << "orient" << std::endl;
      triangleOrientationCheck(labelMeshBMap[c].triVec);
      //std::cout << "orient out" << std::endl;
      if(!Mesh::meshFromTriangles(S2, labelMeshBMap[c].vtxVec, labelMeshBMap[c].triVec)){
        std::cout << "Error, cannot add all the triangles" << std::endl;
      }
    }

  }

 bool CreateJunctionGraph::run(Mesh *m1, Mesh *m2, QString filename, bool createGraph, bool createMesh)
  {
    if(m1->tissue().meshType() != "MGX3D") return setErrorMessage(QString("Mesh Type must be MGX3D"));

    m2->reset();
    CellTissue& T = m1->tissue();
    vvGraph S2;// = m2->graph();

    std::map<std::set<int>, JunctionType> junctionsMap;
    std::map<vertex, JunctionType> vtxsJunctions;

    std::set<int> allLabels;

    std::map<int, Point3d> labelPosMap;
    std::map<int, int> labelVoxelCounter;

    // remember network for writing to the files at the end
    std::map<vertex, int> vtxIdxMap;
    std::vector<std::pair<vertex,vertex> > edgeList;
    //std::map<vertex, std::set<int> > junctionNeighborCellsMap;


    std::map<vertex, std::set<int> > vtxCellsMap;
    //typedef std::pair<vertex, std::set<int> > VtxIntSetP;

    std::map<int, cell> labelCellMap;

    // create a map of vtxs to their set of cells
    forall(const cell& c, T.C){
      forall(const vertex& v, c->S){
        vtxCellsMap[v].insert(c->label);
      }
      labelCellMap[c->label] = c;
      labelPosMap[c->label] = c->pos;
      allLabels.insert(c->label);
    }

    /*forall(VtxIntSetP p, vtxCellsMap){
      if(p.second.size()>=3){
        vertex v;
        v->pos = p.first->pos;
      }
    }*/

    std::map<cell, std::set<int> > cellNeighborsMap;

    // forall cells find neighbors
    forall(const cell& c, T.C){
      forall(const cell& n, T.C.neighbors(c)){
        cellNeighborsMap[c].insert(n->label);
      }
    }

    std::map<vertex, std::set<cell> > junctionNeighborMap;
    std::set<std::set<cell> > finishedTuples;

    std::set<std::set<cell> > finishedTris;

    std::map<std::set<cell>, vertex> neighborJunctionMap;

    std::set<vertex> innerJ, outerJ;

    std::set<cellTriplet> allTris, allTrisInTetras, allOutsideTris;

    std::cout << "find junctions" << std::endl;

    // forall cells
    forall(const cell& c, T.C){


      // forall neighbors create (c1,c2) pairs
      forall(const cell& n, T.C.neighbors(c)){
        //cellPair cp = std::make_pair(c,n);
        //if(n<c) cp = std::make_pair(n,c);

        // forall common neighbors create (c1,c2,c3) triplets
        forall(const cell& m, T.C.neighbors(n)){
          if(m == c) continue;
          if(cellNeighborsMap[c].find(m->label) == cellNeighborsMap[c].end()) continue; // m no neighbor of c

          cellTriplet ct = createCellTriplet(c,n,m);
          allTris.insert(ct);

          // m has to be a neighbor of both c and n

          forall(const cell& b, T.C.neighbors(m)){
            // forall common neighbors add point in m2
            if(b == c) continue;
            if(b == n) continue;
            if(cellNeighborsMap[c].find(b->label) == cellNeighborsMap[c].end()) continue;
            if(cellNeighborsMap[n].find(b->label) == cellNeighborsMap[n].end()) continue;

            // b has to be a neighbor of c, n and m

            std::set<cell> currentTuple;
            Point3d pos = c->pos + n->pos + m->pos + b->pos;
            double count = 4;
            currentTuple.insert(c);
            currentTuple.insert(n);
            currentTuple.insert(m);
            currentTuple.insert(b);

            // check for 5 way junction
            forall(const cell& a, T.C.neighbors(b)){
              if(a == c or a == n or a == m) continue;
              if(cellNeighborsMap[c].find(a->label) == cellNeighborsMap[c].end()) continue;
              if(cellNeighborsMap[n].find(a->label) == cellNeighborsMap[n].end()) continue;
              if(cellNeighborsMap[m].find(a->label) == cellNeighborsMap[m].end()) continue;

              currentTuple.insert(a);
              pos+=a->pos;
              count++;


              //std::cout << " 5-way junction!" << c->label << "/" << n->label <<"/"<< m->label << "/" << b->label << "/" << a->label << std::endl;
              forall(const cell& z, T.C.neighbors(a)){
                if(z == c or z == n or z == m or z == b) continue;
                if(cellNeighborsMap[c].find(z->label) == cellNeighborsMap[c].end()) continue;
                if(cellNeighborsMap[n].find(z->label) == cellNeighborsMap[n].end()) continue;
                if(cellNeighborsMap[m].find(z->label) == cellNeighborsMap[m].end()) continue;
                if(cellNeighborsMap[b].find(z->label) == cellNeighborsMap[b].end()) continue;
                //std::cout << " 6-way junction!" << c->label << "/" << n->label <<"/"<< m->label << "/" << b->label << "/" << a->label << "/" << z->label << std::endl;
                currentTuple.insert(z);
                pos+=z->pos;
                count++;



                forall(const cell& y, T.C.neighbors(z)){
                  if(y == c or y == n or y == m or y == b or y == a) continue;
                  if(cellNeighborsMap[c].find(y->label) == cellNeighborsMap[c].end()) continue;
                  if(cellNeighborsMap[n].find(y->label) == cellNeighborsMap[n].end()) continue;
                  if(cellNeighborsMap[m].find(y->label) == cellNeighborsMap[m].end()) continue;
                  if(cellNeighborsMap[b].find(y->label) == cellNeighborsMap[b].end()) continue;
                  if(cellNeighborsMap[a].find(y->label) == cellNeighborsMap[a].end()) continue;
                  std::cout << " 7-way junction!" << c->label << "/" << n->label <<"/"<< m->label << "/" << b->label << "/" << a->label << "/" << z->label << std::endl;
                }
              }

            }



            if(finishedTuples.find(currentTuple) != finishedTuples.end()) continue;

            finishedTuples.insert(currentTuple);

            //std::cout << "neighbor" << c->label << "/" << n->label <<"/"<< m->label << "/" << b->label << std::endl;

            vertex v;
            v->pos = pos/count;
            S2.insert(v);
            int currentIdx = vtxIdxMap.size()+1;
            vtxIdxMap[v] = currentIdx;
            junctionNeighborMap[v] = currentTuple;
            neighborJunctionMap[currentTuple] = v;

          }
        }
      }
    }


    std::cout << "find cell walls" << std::endl;
    std::map<cellTriplet, std::set<vertex> > cellWalls;
    std::map<cellTriplet, vertex> centerOfTetraMap;

    std::map<IntInt, std::set<vertex> > cellWallsPairwise;


    forall(const vertex& v, S2){
      cell c1,c2,c3,c4;
      std::set<cell>::iterator it = junctionNeighborMap[v].begin();
     // std::cout << "vv " << v->pos << std::endl;
      if(it == junctionNeighborMap[v].end()) continue;
      //std::cout << "take " << std::endl;
      c1 = *(it);
      it++;
      c2 = *(it);
      it++;
      c3 = *(it);
      it++;
      c4 = *(it);
//std::cout << "done " << std::endl;
      cellTriplet ct1,ct2,ct3,ct4;
      ct1 = createCellTriplet(c1,c2,c3);
      ct2 = createCellTriplet(c1,c2,c4);
      ct3 = createCellTriplet(c1,c3,c4);
      ct4 = createCellTriplet(c2,c3,c4);

      cellWalls[ct1].insert(v);
      cellWalls[ct2].insert(v);
      cellWalls[ct3].insert(v);
      cellWalls[ct4].insert(v);

      allTrisInTetras.insert(ct1);
      allTrisInTetras.insert(ct2);
      allTrisInTetras.insert(ct3);
      allTrisInTetras.insert(ct4);

      centerOfTetraMap[ct1] = v;
      centerOfTetraMap[ct2] = v;
      centerOfTetraMap[ct3] = v;
      centerOfTetraMap[ct4] = v;

      cellPair cp1,cp2,cp3,cp4,cp5,cp6;
      cp1 = createCellPair(c1,c2);
      cp2 = createCellPair(c1,c3);
      cp3 = createCellPair(c1,c4);
      cp4 = createCellPair(c2,c3);
      cp5 = createCellPair(c2,c4);
      cp6 = createCellPair(c3,c4);


      IntInt cp1i = std::make_pair(c1->label, c2->label);
      if(c1->label > c2->label) cp1i = std::make_pair(c2->label, c1->label);

      IntInt cp2i = std::make_pair(c1->label, c3->label);
      if(c1->label > c3->label) cp2i = std::make_pair(c3->label, c1->label);

      IntInt cp3i = std::make_pair(c1->label, c4->label);
      if(c1->label > c4->label) cp3i = std::make_pair(c4->label, c1->label);

      IntInt cp4i = std::make_pair(c2->label, c3->label);
      if(c2->label > c3->label) cp4i = std::make_pair(c3->label, c2->label);

      IntInt cp5i = std::make_pair(c2->label, c4->label);
      if(c2->label > c4->label) cp5i = std::make_pair(c4->label, c2->label);

      IntInt cp6i = std::make_pair(c3->label, c4->label);
      if(c3->label > c4->label) cp6i = std::make_pair(c4->label, c3->label);

      cellWallsPairwise[cp1i].insert(v);
      cellWallsPairwise[cp2i].insert(v);
      cellWallsPairwise[cp3i].insert(v);
      cellWallsPairwise[cp4i].insert(v);
      cellWallsPairwise[cp5i].insert(v);
      cellWallsPairwise[cp6i].insert(v);


//std::cout << "done2 " << std::endl;
    }

    typedef std::pair<cellTriplet, std::set<vertex> > cellTripletVtxSetP;

    std::cout << "connect junctions" << std::endl;

    std::map<cellTriplet, vertex> triVtxMap;

    forall(cellTripletVtxSetP p, cellWalls){
      // find outside tris of tetras
      if(p.second.size() == 1){
        vertex v;
        v->pos = (p.first.first->pos + p.first.second.first->pos + p.first.second.second->pos)/3.;
        S2.insert(v);
        junctionNeighborMap[v].insert(p.first.first);
        junctionNeighborMap[v].insert(p.first.second.first);
        junctionNeighborMap[v].insert(p.first.second.second);
        int currentIdx = vtxIdxMap.size()+1;
        vtxIdxMap[v] = currentIdx;
        S2.insertEdge(v,centerOfTetraMap[p.first]);
        S2.insertEdge(centerOfTetraMap[p.first],v);
        edgeList.push_back(std::make_pair(v,centerOfTetraMap[p.first]));
        edgeList.push_back(std::make_pair(centerOfTetraMap[p.first],v));
        allOutsideTris.insert(p.first);
        triVtxMap[p.first] = v;
      } else { // "inside" triangles, connect them
        forall(const vertex& v, p.second){
          forall(const vertex& n, p.second){
            if(v==n) continue;
            S2.insertEdge(v,n);
            edgeList.push_back(std::make_pair(v,n));
          }
        }
      }
    }



    //std::map<cell, std::vector<std::set<vertex> > > cellWallsMap;

    // outside junctions that are not part of a tetra
    forall(cellTriplet ct, allTris){
      if(allTrisInTetras.find(ct) == allTrisInTetras.end()){
        vertex v;
        v->pos = (ct.first->pos + ct.second.first->pos + ct.second.second->pos)/3.;
        S2.insert(v);
        junctionNeighborMap[v].insert(ct.first);
        junctionNeighborMap[v].insert(ct.second.first);
        junctionNeighborMap[v].insert(ct.second.second);
        int currentIdx = vtxIdxMap.size()+1;
        vtxIdxMap[v] = currentIdx;
        allOutsideTris.insert(ct);
        triVtxMap[ct] = v;
      }
    }

    std::map<cellPair, std::set<vertex> > cellWallsOutside;

    // go through all outside tris make map of edges
    forall(cellTriplet ct, allOutsideTris){
      cellPair cp1,cp2,cp3;
      cp1 = std::make_pair(ct.first, ct.second.first);
      cp2 = std::make_pair(ct.first, ct.second.second);
      cp3 = std::make_pair(ct.second.first, ct.second.second);
      cellWallsOutside[cp1].insert(triVtxMap[ct]);
      cellWallsOutside[cp2].insert(triVtxMap[ct]);
      cellWallsOutside[cp3].insert(triVtxMap[ct]);
    }

    typedef std::pair<cellPair, std::set<vertex> > cellPairVtxSetP;

    forall(cellPairVtxSetP p, cellWallsOutside){
      // find outside tris
      if(p.second.size() == 1){
      } else {
        forall(const vertex& v, p.second){
          //v->selected = true;
          forall(const vertex& n, p.second){
            if(v==n) continue;
            S2.insertEdge(v,n);
            edgeList.push_back(std::make_pair(v,n));
          }
        }
      }
    }

    std::cout << "all finished" << std::endl;
    if(createGraph){
      vvGraph& m2Graph = m2->graph();
      m2Graph = S2;
      m2->updateAll();
    }

    if(createMesh){
      vvGraph S;
      buildMeshFromJunctions(S2, S, allLabels, cellWallsPairwise, labelPosMap);
      m1->reset();
      m1->graph() = S;
      m1->updateAll();
      /*std::map<cell, MeshBuilder> cellMeshBuilderMap;

      forall(const cell& c, T.C){
        MeshBuilder mb(1E-5);
        cellMeshBuilderMap[c] = mb;
      }


      forall(cellPairVtxSetP p, cellWallsPairwise){

        const cell& c1 = p.first.first;
        const cell& c2 = p.first.second;

        // order the vertices
        std::set<vertex> vtxs = p.second;

        // create segments
        std::vector<std::pair<Point3d,Point3d> > polygonSegs;
        std::map<Point3d, vertex> P3dVertexMap;
        std::set<std::pair<vertex,vertex> > segsDone;
        forall(const vertex& v, vtxs){
          forall(const vertex& n, S2.neighbors(v)){
            if(vtxs.find(n) == vtxs.end()) continue;
            std::pair<vertex,vertex> vv = std::make_pair(v,n);
            if(n<v) vv = std::make_pair(v,n);
            if(segsDone.find(vv) != segsDone.end()) continue;

            segsDone.insert(vv);
            polygonSegs.push_back(std::make_pair(v->pos,n->pos));
            P3dVertexMap[v->pos] = v;
            P3dVertexMap[n->pos] = n;

          }
        }

        // order the segments
        std::vector<std::vector<Point3d> > orderedSegs = orderPolygonSegsMulti(polygonSegs);


        forall(std::vector<Point3d>& wall, orderedSegs){

          // add a vertex to the center
          Point3d centerPos(0,0,0);
          forall(Point3d& p, wall)
            centerPos+=p;
          centerPos/=wall.size();

          vertex center;
          center->pos = centerPos;

          Point3d prev = wall[wall.size()-1];
          forall(Point3d& p, wall){
            // add the triangles with right orientation
            Point3d dirC1 = c1->pos-centerPos;
            double s1 = ((p - prev) ^ (centerPos - prev))*(dirC1);
            vertex w1,w2,w3;
            w1->pos = prev;
            w2->pos = p;
            w3->pos = centerPos;
            if(s1 > 0) cellMeshBuilderMap[c1].addTri(w2, w1, w3);
            else cellMeshBuilderMap[c1].addTri(w1, w2, w3);

            vertex w4,w5,w6;
            w4->pos = prev;
            w5->pos = p;
            w6->pos = centerPos;

            Point3d dirC2 = c2->pos-centerPos;
            double s2 = ((p - prev) ^ (centerPos - prev))*(dirC2);
            if(s2 > 0) cellMeshBuilderMap[c2].addTri(w5, w4, w6);
            else cellMeshBuilderMap[c2].addTri(w4, w5, w6);

            prev = p;
          }

        }

      }
      vvGraph S;
      forall(const cell& c, T.C){

        forall(const vertex& v, cellMeshBuilderMap[c].vtxVec){
          v->label = c->label;
        }

        if(!Mesh::meshFromTriangles(S, cellMeshBuilderMap[c].vtxVec, cellMeshBuilderMap[c].triVec)){
          std::cout << "Error, cannot add all the triangles" << std::endl;
        }
      }
      m1->graph() = S;*/
    }

    forall(VtxIntP p, vtxIdxMap){
      forall(const cell& c, junctionNeighborMap[p.first]){
        vtxsJunctions[p.first].neighborCells.insert(c->label);
      }
    }

    writeJunctionGraphToFiles(filename, vtxIdxMap, edgeList, vtxsJunctions);

    return true;
  }
  REGISTER_PROCESS(CreateJunctionGraph);

 // file dialogue for save data
 bool CreateJunctionGraphStack::initialize( QWidget* parent)
  {
    QString filename = parm("Save to file");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Save to file", filename);
    return true;
  }


 bool findJunctionInStack(const Stack* s, const HVecUS& data, Point3i pos, Point3i size, std::set<int>& neighbors)
 {
  bool res = false;
  std::set<int> neighborsNew;
  int l = data[s->offset(pos)];
  neighborsNew.insert(l);
/* // go through 3x3 neighborhood -> too many junctions found
  for(int i = -1; i<=1; i++){
    for(int j =-1; j<=1; j++){
      for(int k =-1; k<=1; k++){
        Point3i vx = Point3i(i+pos.x(),j+pos.y(),k+pos.z());
        int label = data[s->offset(vx)];
        //std::cout << "in " << l << "/" << label << "/" << pos << "//" << vx << std::endl;
        neighborsNew.insert(label);
      }
    }
  }
  */
  // go through 2x2 neighborhood
  for(int i = 0; i<=1; i++){
    for(int j =0; j<=1; j++){
      for(int k =0; k<=1; k++){
        Point3i vx = Point3i(i+pos.x(),j+pos.y(),k+pos.z());
        int label = data[s->offset(vx)];
        neighborsNew.insert(label);
      }
    }
  }

  /* // go through direct neighborhood -> not enough junctions found
  for(int i = -1; i<=1; i=i+2){
    Point3i vx = Point3i(i+pos.x(),pos.y(),pos.z());
    Point3i vy = Point3i(pos.x(),i+pos.y(),pos.z());
    Point3i vz = Point3i(pos.x(),pos.y(),i+pos.z());
    neighborsNew.insert(data[s->offset(vx)]);
    neighborsNew.insert(data[s->offset(vy)]);
    neighborsNew.insert(data[s->offset(vz)]);
  }
  */
   if(neighborsNew.size() > 3) res = true;
   //if(neighborsNew.size() > 5) res = true;
   neighbors = neighborsNew;

   return res;
 }


 bool CreateJunctionGraphStack::run(const Stack* s, Mesh *m1, Mesh *m2, QString filename, bool createGraph, bool createMesh)
  {

    Point3i imgSize(s->size());
    const Store* sw = s->main();
    const HVecUS& data = sw->data();

    int x = imgSize.x();
    int y = imgSize.y();
    int z = imgSize.z();

    int start = 1;
    int stop = 1;

    std::map<std::set<int>, JunctionType> junctionsMap;
    std::map<vertex, JunctionType> vtxsJunctions;

    std::set<int> allLabels;

    std::map<int, Point3d> labelPosMap;
    std::map<int, int> labelVoxelCounter;

    // remember network for writing to the files at the end
    std::map<vertex, int> vtxIdxMap;
    std::vector<std::pair<vertex,vertex> > edgeList;
    //std::map<vertex, std::set<int> > junctionNeighborCellsMap;

    //#pragma omp parallel for
    for(int i = start; i<x-stop; i++){
      for(int j = start; j<y-stop; j++){
        for(int k = start; k<z-stop; k++){
          Point3f p = s->imageToWorld(Point3f(i,j,k));
          Point3i vx = Point3i(i,j,k);
          if(data[s->offset(vx)]==0) continue;
          allLabels.insert(data[s->offset(vx)]);
          //Point3d pOutWorld(s2->imageToWorld(pOut));
          //std::cout << "voxel " << data[s->offset(vx)] << std::endl;
          std::set<int> neighbors;
          if(findJunctionInStack(s,data,vx,imgSize,neighbors)){
            if(junctionsMap.find(neighbors) == junctionsMap.end()){
              JunctionType j;
              junctionsMap[neighbors] = j;
              junctionsMap[neighbors].count = 0;
              junctionsMap[neighbors].pos = Point3d(0,0,0);
              vertex v;
              junctionsMap[neighbors].v = v;
            }
            junctionsMap[neighbors].neighborCells = neighbors;
            junctionsMap[neighbors].count++;
            junctionsMap[neighbors].pos+=Point3d(vx.x()+0.5,vx.y()+0.5,vx.z()+0.5);
            //junctions[neighbors].first++;
            //junctions[neighbors].second+=Point3d(vx.x(),vx.y(),vx.z());
            //std::cout << "j " << "/" << neighbors.size() << std::endl;
            //forall(int zzz, neighbors) std::cout << " a " << zzz;
            //std::cout << "j " << "/" << neighbors.size() << std::endl;
          }

          Point3d pd = Point3d(p.x(), p.y(), p.z());
          labelPosMap[data[s->offset(vx)]]+=pd;
          labelVoxelCounter[s->offset(vx)]++;

        }
      }
    }
    typedef std::pair<int,Point3d> IntP3dP;

    forall(IntP3dP p, labelPosMap){
      p.second /= labelVoxelCounter[p.first];
    }

    // check if 6 and 5 way junctions exist and delete redundant 5 and 6 way junctions
    deleteRedundantJunctions45(junctionsMap);

    std::cout << "total junctions " << junctionsMap.size() << std::endl;

    std::vector<int> counter(8);

    // calculate the junction positions by averaging across all found points of identical junctions
    for(std::map<std::set<int>, JunctionType>::iterator it = junctionsMap.begin(); it!=junctionsMap.end(); it++){
      //std::cout << "p1" << it->second.pos << std::endl;
      Point3f imgP(it->second.pos / (double)it->second.count);
      imgP = s->imageToWorld(imgP);
      it->second.pos = Point3d(imgP.x(), imgP.y(), imgP.z());
      it->second.v->pos = it->second.pos;
      counter[it->first.size()]++;
      //std::cout << "p2" << it->second.pos << std::endl;
    }
    int k = 0;
    forall(int i, counter){
      std::cout << "junctions of size " << k++ << ": " << i << std::endl;
    }

    // create the junction graph
    vvGraph S;
    std::map<IntInt, std::set<vertex> > cellPairJunctionMap;
    buildJunctionGraph(S, junctionsMap, vtxsJunctions, edgeList, vtxIdxMap, cellPairJunctionMap);

    if(createGraph){
      //std::cout << "create graph " << S.size() << std::endl;
      m1->reset();
      m1->graph() = S;
      //std::cout << "create graph " << m1->graph().size() << std::endl;
      m1->updateAll();
    }

    // create the mesh
    if(createMesh){
      m2->reset();
      vvGraph& S2 = m2->graph();
      buildMeshFromJunctions(S, S2, allLabels, cellPairJunctionMap, labelPosMap);
      m2->updateAll();
    }

    writeJunctionGraphToFiles(filename, vtxIdxMap, edgeList, vtxsJunctions);

    return true;
  }
    REGISTER_PROCESS(CreateJunctionGraphStack);

 // Post-Division - 2 cells as input, they will be joined in a temporary mesh and
 //bool OrganDivisionSimulation::run(Mesh *m, Mesh *m2, QString order, QString primary, QString secondary, QString primaryDir, QString secondaryDir,
 //       double tolerance, int planeNr, bool perpendicular, bool perpendicularNeighbor, int stepsCenter, double stepSizeCenter,
 //       bool perpendicularConstrained, bool cuttingNeighborWall)
 bool OrganDivisionSimulation::run(Mesh *m, Mesh *m2, QString order, DivisionParms& p, double tolerance)
  {

    std::cout << "par " << p.perpendicularNeighbor << "/" << p.perpendicularConstrained << std::endl;

    // find parent labels
    //vvGraph& S = m->graph();
    CellTissue& T = m->tissue();

    std::set<cell> cellsToDivide;

    // find the cells with parent > 0 that have to be divided
    forall(const cell& c, T.C){
      if(m->parents()[c->label] > 0){
        cellsToDivide.insert(c);
        //m->parents()[c->label] = 0;
      }
    }

    if(order != "Parent Label Min" and order != "Neighbor Count Max")
      return setErrorMessage("Wrong Division Order");

    bool calcBetw = false;
    if(p.primary == "Betweenness Centrality" or p.primary == "Random Walk Betweenness"
      or p.primary == "Betweenness Centrality (weighted)" or p.primary == "Random Walk Betweenness (weighted)"
      or p.secondary == "Betweenness Centrality" or p.secondary == "Random Walk Betweenness"
      or p.secondary == "Betweenness Centrality (weighted)" or p.secondary == "Random Walk Betweenness (weighted)") calcBetw = true;

    bool constWeight = true;
    if(calcBetw and (p.primary == "Betweenness Centrality (weighted)" or p.primary == "Random Walk Betweenness (weighted)"
      or p.secondary == "Betweenness Centrality (weighted)" or p.secondary == "Random Walk Betweenness (weighted)")) constWeight = false;
    bool dijk = true;


    // loop through all cellsToDivide, current one deleted at the end of the loop
    while(cellsToDivide.size()>0){

      // find the cell with most neighbors or lowest parent label
      cell cellToDivide;
      int maxN = -1;
      int minL = 1E9;

      forall(const cell& c, cellsToDivide){
        if(order == "Neighbor Count Max"){
          int neighbNr = countNeighbours(T.C, c);
          if(neighbNr > maxN){
            maxN = neighbNr;
            cellToDivide = c;
          }
        } else if(order == "Parent Label Min"){
          if(m->parents()[c->label] < minL and m->parents()[c->label] > 0){
            minL = m->parents()[c->label];
            cellToDivide = c;
          }
        }
      }

      // do divison plane analysis
      TestDivisionPlanes t(*this);
      t.testDivisionPlanesFunc(m, m2, cellToDivide->label, false, "", constWeight, dijk, 0, calcBetw, false, "Centroid", p);

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

      // find the optimal plane
      std::vector<CellDivision> optPlanes;
      double optValue;
      P3dDouble bestDivPlane;

      // first sort for the primary criterium
      // then take best results and sort for secondary criterium
      findBestPlane(dat.divData, optPlanes, optValue, tolerance, bestDivPlane, p.primary, p.secondary, p.primaryDir, p.secondaryDir);
      std::cout << "best plane " << bestDivPlane.first << "/" << bestDivPlane.second << "/ of total planes: " << optPlanes.size() << endl;

      // divide it
      std::vector<P3dDouble> samplingDirections;
      //Subdivide* s;
      //mgx::CellTissue::DivAlg divAlg;
      //divAlg = mgx::CellTissue::stringToDivAlg("PREDEFINED DIRECTION");
      samplingDirections.push_back(bestDivPlane);

      //if(newMethod){
        double displacement = bestDivPlane.second;
        std::cout << "pl  " << cellToDivide->pos << "/" << displacement << std::endl;

        DivPlane p(bestDivPlane.first, cellToDivide->pos + displacement*bestDivPlane.first, displacement);
        CellDivision cd (cellToDivide, p, T.C);
        cd.divideCell3D(T);

      //} else {
        //if(!divideCell3D(T, cellToDivide, s, divAlg, samplingDirections)) return setErrorMessage("wrong mesh type, MGX3D required");
      //}
      // erase it from cellsToDivide
      cellsToDivide.erase(cellToDivide);
      T.updateAllGraph();
    }

    forall(const cell& c, T.C){
      m->parents()[c->label] = 0;
    }
    m->updateAll();
    return true;
  }
  REGISTER_PROCESS(OrganDivisionSimulation);






 bool DivisionPlaneHeat::run(Mesh *m, QString option)
  {

    // get the division planes


    // collor them accoring to the option

    if(option == "Plane Area"){

    } else if(option == "Angle from Actual"){

    } else if(option == "Angle from Selected"){

    } else if(option == "Volume Ratio"){

    }


    return true;
  }
  REGISTER_PROCESS(DivisionPlaneHeat);

 bool MeasureBetweenessCentrality::run(Mesh *m, IntFloatAttr& heatMap, QString weightType)
  {

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

    vvGraph& S = m->graph();
    std::map<IntInt, double> neighWallLengthMap;
    QString neighborhoodWeightType = "";
    if(weightType == "Euclidean" or weightType == "1 / Euclidean") neighborhoodWeightType = "Euclidean";

    // calc neighborhhod graph
    neighborhood2D(S, neighWallLengthMap, neighborhoodWeightType);

    std::map<int, double> cellBetweenessMap;


    std::map<IntInt, double> neighWeightMap;

    forall(auto p, neighWallLengthMap){
      if(p.second < 1E-5) continue;
      if(weightType == "None") neighWeightMap[p.first] = 1;
      else if(weightType == "Wall Area") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Wall Area") neighWeightMap[p.first] = 1./p.second;
      else if(weightType == "Euclidean") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Euclidean") neighWeightMap[p.first] = 1./p.second;
    }

    // calc betweenness
    calcBetweenness(neighWeightMap, cellBetweenessMap);


    heatMap.clear();
    forall(auto p, cellBetweenessMap){
      heatMap[p.first] = p.second;
      attrData[p.first] = p.second;
    }

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

    return true;
  }
  REGISTER_PROCESS(MeasureBetweenessCentrality);


 bool MeasureBetweenessCentrality3D::run(Mesh *m, IntFloatAttr& heatMap, QString weightType)
  {

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

    std::map<IntInt, double> neighMap;
    vvGraph& S = m->graph();

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

    neighMap = info.sharedWallArea;

    std::map<IntInt, double> neighWeightMap;

    forall(auto p, neighMap){
      if(p.second < 1E-5) continue;
      if(weightType == "None") neighWeightMap[p.first] = 1;
      else if(weightType == "Wall Area") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Wall Area") neighWeightMap[p.first] = 1./p.second;
      else if(weightType == "Euclidean") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Euclidean") neighWeightMap[p.first] = 1./p.second;
    }


    std::map<int, double> cellBetweenessMap;

    // calc betweenness
    calcBetweenness(neighWeightMap, cellBetweenessMap);

    heatMap.clear();
    forall(auto p, cellBetweenessMap){
      heatMap[p.first] = p.second;
    }

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

    return true;
  }
  REGISTER_PROCESS(MeasureBetweenessCentrality3D);

 bool MeasureBetweenessCurrentFlow::run(Mesh *m, IntFloatAttr& heatMap, QString weightType)
  {

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Network/Betweenness Current Flow");
    attrData.clear();

    vvGraph& S = m->graph();
    std::map<IntInt, double> neighMap;
    //QString weightType;

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

    // calc neighborhhod graph
    neighborhood2D(S, neighMap, weightType);

    std::map<IntInt, double> neighWeightMap;

    forall(auto p, neighMap){
      if(p.second < 1E-5) continue;
      if(weightType == "None") neighWeightMap[p.first] = 1;
      else if(weightType == "Wall Area") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Wall Area") neighWeightMap[p.first] = 1./p.second;
      else if(weightType == "Euclidean") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Euclidean") neighWeightMap[p.first] = 1./p.second;
    }

    std::map<int, double> currentFlowMap;

    // calc betweenness
    calcRandWalkBetweenness(neighWeightMap, currentFlowMap, false);

    heatMap.clear();
    forall(auto p, currentFlowMap){
      heatMap[p.first] = p.second;
      attrData[p.first] = p.second;
    }

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


    return true;
  }
  REGISTER_PROCESS(MeasureBetweenessCurrentFlow);


 bool MeasureBetweenessCurrentFlow3D::run(Mesh *m, IntFloatAttr& heatMap, QString weightType)
  {

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Network/Betweenness Current Flow");
    attrData.clear();


    vvGraph& S = m->graph();
    std::map<IntInt, double> neighMap;

    // calc neighborhhod graph
    NhbdGraphInfo info;
    double tolerance = 0.0001;
    neighborhoodGraph(S, tolerance, info); // now in graphutils

    neighMap = info.sharedWallArea;

    std::map<IntInt, double> neighWeightMap;

    forall(auto p, neighMap){
      if(p.second < 1E-5) continue;
      if(weightType == "None") neighWeightMap[p.first] = 1;
      else if(weightType == "Wall Area") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Wall Area") neighWeightMap[p.first] = 1./p.second;
      else if(weightType == "Euclidean") neighWeightMap[p.first] = p.second;
      else if(weightType == "1 / Euclidean") neighWeightMap[p.first] = 1./p.second;
    }


    std::map<int, double> currentFlowMap;

    // calc betweenness
    calcRandWalkBetweenness(neighWeightMap, currentFlowMap, false);

    heatMap.clear();
    forall(auto p, currentFlowMap){
      heatMap[p.first] = p.second;
      attrData[p.first] = p.second;
    }

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


    return true;
  }
  REGISTER_PROCESS(MeasureBetweenessCurrentFlow3D);


 bool CreatePlaneFromSelectedVtxs::run(Mesh *m, Mesh *m2, double size, bool overwrite)
  {

    AttrMap<int, Point3d>& planePositions = m2->attributes().attrMap<int, Point3d>("Circular Planes Pos");
    AttrMap<int, Point3d>& planeNormals = m2->attributes().attrMap<int, Point3d>("Circular Planes Nrml");


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

    std::vector<Point3d> vtxPos;

    if(overwrite){
      m2->reset();
      planePositions.clear();
      planeNormals.clear();
    }

    forall(const vertex& v, S){
      if(!v->selected) continue;
      vtxPos.push_back(v->pos);
    }

    // compute PCA of plane
    Point3d planePos, planeNrml;
    findPolygonPlane(vtxPos, planePos, planeNrml);

    int maxLabel = 1;
    forall(const vertex& v, S2){
      if(v->label >= maxLabel) maxLabel = v->label + 1;
    }

    DivisionPlane d = DivisionPlane(size, planeNrml, planePos, maxLabel);

    planePositions[maxLabel] = planePos;
    planeNormals[maxLabel] = planeNrml;

    drawPlane(*this, m2, d);
    m2->updateAll();


    return true;
  }
  REGISTER_PROCESS(CreatePlaneFromSelectedVtxs);


 bool AngleTwoPlanes::run(Mesh *m)
  {

    AttrMap<int, Point3d>& planeNormals = m->attributes().attrMap<int, Point3d>("Circular Planes Nrml");

    // get selected vtxs
    vvGraph& S = m->graph();

    std::set<int> selectedLabels;

    forall(const vertex& v, S){
      if(!v->selected) continue;
      if(v->label < 1) continue;
      selectedLabels.insert(v->label);
    }

    if(selectedLabels.size() != 2) return setErrorMessage("Select exactly two planes.");

    int label1 = *selectedLabels.begin();
    int label2 = *(++selectedLabels.begin());

    double angle = angleVectors(planeNormals[label1], planeNormals[label2]);

    std::cout << "Angle between Planes " << label1 << " and " << label2 << ": " << angle << std::endl;


    return true;
  }
  REGISTER_PROCESS(AngleTwoPlanes);

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


 bool ExportPlaneData::run(Mesh *m, QString filename, QString mode)
  {

    QString attrName = "Cell Division Simulated Planes Filtered";

    if(mode == "Actual Planes")
      attrName = "Cell Division Actual Plane";
    else if(mode == "Simulated Planes")
      attrName = "Cell Division Simulated Planes";

    AttrMap<int, CellDivisionAttr>& attrMap = m->attributes().attrMap<int, CellDivisionAttr>(attrName);

    if(attrMap.empty()) return setErrorMessage("Attr Map empty");

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

    out << "Cell Parent Label,Plane Area,Plane Area Relative Shortest,Plane Area Relative Actual,Normalized Area,Angle Actual,Daughter Ratio,Plane Normal X,Plane Normal Y,Plane Normal Z";
    out << endl;

    forall(auto p, attrMap){
      CellDivisionAttr cda = p.second;
      out << cda.cellLabel << "," << cda.planeArea << "," << cda.planeAreaRelShortest << "," << cda.planeAreaRelActual
      << "," << cda.normalizedArea
      << "," << cda.angleActual << "," << cda.daughterRatio << ",";
      out << cda.planeNrml.x() << "," << cda.planeNrml.y() << "," << cda.planeNrml.z();
      out << endl;
    }




    return true;
  }
  REGISTER_PROCESS(ExportPlaneData);




 bool FlatDivisionPlane3D::run(Mesh *m, Mesh *m2, bool display, double planeSize, bool resetMesh, int label)
  {
    std::cout << "FlatDivPlane: " << display << "/" << planeSize << "/" << resetMesh << "/" << label << std::endl;
    if(planeSize < 0.0001 and display) return setErrorMessage("Plane Size must be > 0");

    AttrMap<int, CellDivisionAttr>& attrData = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Actual Plane");

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

    vvGraph& S = m->graph();

    if(!m->useParents()) return setErrorMessage("This process only works with activated parent labels!");

    std::set<int> selLabels;

    if(label < 1)
      selLabels = findAllSelectedLabels(S);
    else {
      forall(auto p, m->parents()){
        if(p.second == label) selLabels.insert(p.first);
      }
    }

    selLabels.erase(-1);

    forall(int l, selLabels){
      std::cout << "labels " << l << std::endl;
    }

    //std::cout << "FlatDivPlane: sel1 " << selLabels.size() << std::endl;

    if(selLabels.size() != 2) return setErrorMessage("Please select exactly two neighboring cells with the same parent label");

    std::set<int> selParents;
    if(label < 1){ // find the parent label

      forall(int l, selLabels){
        selParents.insert(m->parents()[l]);
      }
      //std::cout << "FlatDivPlane: sel2 " << selParents.size() << std::endl;
      if(selParents.size() != 1) return setErrorMessage("Please select exactly one parent label");
    }
    if(label > 0){ // parent label is set
      selParents.insert(label);
    }

    int selectedParent = *selParents.begin();
    //std::cout << "FlatDivPlane: sel3 " << selectedParent << std::endl;

    Point3d cellNormal(0,0,0);

    std::vector<vertex> cellVtxs = getVtxsOfCells(m, selectedParent, true);

    std::vector<Point3d> pointsSharedWall;
    std::vector<vertex> vtxsSharedWall = getOrRemoveCommonWall(m, cellVtxs, false, true);

    forall(vertex v, vtxsSharedWall){
      pointsSharedWall.push_back(Point3d(v->pos));
    }

    std::cout << "size " << pointsSharedWall.size() << std::endl;

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

    CellDivisionMGXM cd;
    cd.cellTrisAll = getTrisOfCells(m, selectedParent, true);
    //std::cout << "cdddddddd " << selectedParent << "/"  << cd.cellTrisAll.size() << std::endl;
    cd.initialize(m,selectedParent);

    P3dDouble divPlane = std::make_pair(planeNrml, 0);
    CellDivisionAttr cda = cd.simulateDivision3D(divPlane, planePos);

    cda.planeNrml = planeNrml;
    cda.planeDisplacement = 0; // TODO
    cda.planePos = planePos;
    cda.div2d = 1;
    cda.cellLabel = selectedParent;
    attrData[selectedParent] = cda;
    normals[selectedParent] = planeNrml;

    if(resetMesh){
      m2->reset();
    }

    if(display){
      DivisionPlane dp = DivisionPlane(planeSize, cda.planeNrml, cda.planePos, selectedParent);
      drawPlane(*this, m2, dp, cellNormal);
      m2->updateAll();
    }

    m->updateAll();

    return true;

  }
  REGISTER_PROCESS(FlatDivisionPlane3D);


     bool createSamplingDirsNew(int planeNr, bool perpSurf, Point3d surfNrml, int stepsCentr, double stepSizeCentr, std::vector<P3dDouble>& samplingDirections)
    {

      // generate the division planes to test

        if(perpSurf){ // perpendicular to surfNrml direction
          if(surfNrml == Point3d(0,0,0)){
          std::cout << "Error: No nearest Surface point found. Analyze the Surface first or set Perpendicular to Surface to No" << std::endl;
          return false;
        }
          // sampling_res=sampling_res/(double)numberOfPlanes*(2*stepsCentr+1);
          //
          //
          // generateSamplingDirections(sampling_res, stepsCentr, stepSizeCentr,
          //   surfNrml, perpNeighborPlanesToBeUsed, false, samplingDirections);

          // if(p.perpendicular and p.perpendicularConstrained){
          //   addCenterDisplacementTestCut(p.stepsCenter, p.stepSizeCenter, cellToDivide->pos, surfaceNormal, perpNeighborPlanes,
          //                               perpNeighborPoints, perpNeighborWalls, samplingDirections);
          // }
        } else { // no restrictions
          generateSamplingDirectionsNew(planeNr, stepsCentr, stepSizeCentr, samplingDirections);
        }

      cout << samplingDirections.size() << " possible division planes will be tested" << endl;

      return true;

    }

  std::vector<CellDivisionAttr> testDivision3D(Mesh *m, std::set<vertex>& verticesCell, int divPlanes, int centerDisplSteps, double centerDisplStepSize,
    Point3d divPointPos, Point3d surfaceNrml, std::set<int>& selLabels, int selectedLabel, bool testActual, CellDivisionAttr& actual)
  {
    vvGraph& S = m->graph();

    // generate the division planes to test
    std::vector<P3dDouble> samplingDirections;

    createSamplingDirsNew(divPlanes, false, Point3d(0,0,0), centerDisplSteps, centerDisplStepSize, samplingDirections);
    std::vector<CellDivisionAttr> possibleDivisionsCDA;

    //std::cout << "testDiv2D bef2 " << verticesCell.size() << "/" << selLabels.size() << std::endl;

forall(int l, selLabels){
  std::cout << "label " << l << std::endl;
}

    // get the cell tris with out the common wall
    CellDivisionMGXM cd;
    cd.cellTrisAll = getTrisOfCells(m, selectedLabel, true);
    std::cout << "cdddddddd " << selectedLabel << "/"  << cd.cellTrisAll.size() << std::endl;
    cd.initialize(m, selectedLabel);

    forall(P3dDouble divPlane, samplingDirections){
      CellDivisionAttr cda = cd.simulateDivision3D(divPlane, divPointPos);
      possibleDivisionsCDA.push_back(cda);
    }

    double minArea = HUGE_VAL;
    double maxArea = -HUGE_VAL;
    forall(CellDivisionAttr d, possibleDivisionsCDA){
      // find min and max of values
      if(d.planeArea < minArea) minArea = d.planeArea;
      if(d.planeArea > maxArea) maxArea = d.planeArea;
    }

    //std::cout << "minAA " << minArea << "/" << maxArea << "/" << possibleDivisionsCDA.size() << std::endl;

    forall(CellDivisionAttr& d, possibleDivisionsCDA){
      d.planeAreaRelShortest = d.planeArea / minArea;
      d.cellLabel = selectedLabel;

      if(testActual){
        actual.planeAreaRelShortest = actual.planeArea / minArea;
        actual.angleActual = 0;
        d.planeAreaRelActual = d.planeArea / actual.planeArea;
        d.angleActual = angleVectors(d.planeNrml, actual.planeNrml, false);
      } else {
        d.planeAreaRelActual = -1;
        d.angleActual = -1;
      }
    }


    return possibleDivisionsCDA;
  }


 bool TestDivisionPlanes3D::run(Mesh *m, Mesh *m2, int divPlanes, int centerDisplSteps, double centerDisplStepSize, QString divCenter, int parentLabel)
  {

    std::cout << "Div analysis s2" << std::endl;

    if(!m->useParents()) return setErrorMessage("This process only works with activated parent labels!");
    std::cout << "Div analysis s21 " << parentLabel << std::endl;
    AttrMap<int, CellDivisionAttr>& attrData = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Simulated Planes");
    //attrData.clear();
    std::cout << "Div analysis s22 " << parentLabel << std::endl;
    AttrMap<int, CellDivisionAttr>& actualPlane = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Actual Plane");
    std::cout << "Div analysis s23 " << parentLabel << std::endl;
    //AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    AttrMap<int, Point3d>& centroidsP = m->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");

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

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

    vertex divPoint;
    //Point3d divPointNrml(0,0,0);
    Point3d divPointPos;

    // find the label if not specfified
    std::set<int> selLabels;

    std::cout << "Div analysis s24 " << parentLabel << std::endl;
    if(parentLabel > 0){
      forall(auto p, m->parents()){
        if(p.second == parentLabel){
          selLabels.insert(p.first);
        }
      }
    } else {
      selLabels = findAllSelectedLabels(S);
    }
    std::cout << "Div analysis s21" << std::endl;
    selLabels.erase(-1);

    std::set<int> selParents;
    //if(m->useParents()){

    forall(int l, selLabels){
      selParents.insert(m->parents()[l]);
    }
    if(selParents.size() != 1) return setErrorMessage("Please select exactly one parent cell");
    if(selLabels.size() != 1 and selLabels.size() != 2) return setErrorMessage("Please select exactly one parent cell with 1 or 2 cell labels");
    //}

    int selectedLabel = *selLabels.begin();
    int selectedParent = *selParents.begin();
    if(parentLabel > 0) selectedParent = parentLabel;


    std::cout << "Div analysis s3" << std::endl;
    if(divCenter == "Custom Vertex"){
      AttrMap<int, vertex>& divPointMap = m->attributes().attrMap<int, vertex>("Division Analysis Custom Point");

      if(divPointMap.empty()) return setErrorMessage("Set a custom division point first!");

      divPoint = divPointMap[0];
      //divPointNrml = divPoint->nrml;
      divPointPos = divPoint->pos;

      std::cout << "custom division point " << divPointPos << /*"/nrml " << divPointNrml <<*/ std::endl;

    }

    std::set<vertex> verticesCell;

    if(parentLabel < 1){
      forall(const vertex& v, m->selectedVertices()){
        verticesCell.insert(v);
      }
    }
    else{
      forall(const vertex& v, S){
        forall(int l, selLabels){
          if(v->label == l) verticesCell.insert(v);
        }
      }
      forall(const vertex& v, verticesCell){
        forall(const vertex& n, S.neighbors(v)){
          if(n->label > 0) continue;
          verticesCell.insert(n);
        }
      }
    }



    std::cout << "Div analysis s4" << std::endl;
    // get all cell vertices from the selection
    //std::set<vertex> verticesCell;

    // forall(const vertex& v, verticesCell){
    //   // verticesCell.insert(v);
    //   if(divCenter != "Custom Vertex"){
    //     divPointNrml += v->nrml;
    //   }
    //   //v->selected = true;
    // }

    Point3d cellCentroid = centroidsP[selectedParent];

    if(divCenter != "Custom Vertex"){
      //divPointNrml /= verticesCell.size();
      //if(m->useParents()) cellCentroid = Point3d(m->parentCenterVis()[selectedParent]);
      //else cellCentroid = Point3d(m->labelCenterVis()[selectedLabel]);

      //std::cout << "cell center division point " << centroids[selectedParent] << /*"/nrml " << divPointNrml <<*/ std::endl;
      //divPointPos = centroids[selectedParent];
      std::cout << "not implemented yet" << std::endl;
    }

    std::cout << "find actual plane" << std::endl;
    FlatDivisionPlane3D *pFlat3D;
    if(!getProcess("Mesh/Division Analysis/Analysis 3D/Division Plane Approximation", pFlat3D))
      throw(QString("Division Analysis:: Unable to make Approximation process"));
    pFlat3D->parm("Display Plane") = "No";
    pFlat3D->parm("Reset Mesh") = "No";

    std::cout << "run1 " << actualPlane.size() << std::endl;
    bool actualOK = pFlat3D->run(m, m2, stringToBool(pFlat3D->parm("Display Plane")), 0.1, stringToBool(pFlat3D->parm("Reset Mesh")), selectedParent);

    if(!actualOK) std::cout << "Couldn't compute actual plane for parent label: " << selectedParent << std::endl;

    // get vtxs of cell
    std::vector<vertex> cellVtxs = getVtxsOfCells(m, selectedParent, true);

    std::vector<Point3d> pointsSharedWall;
    std::vector<vertex> vtxsCellNoSharedWall = getOrRemoveCommonWall(m, cellVtxs, true, false);
    forall(vertex v, S){
      v->selected = false;
    }
    forall(vertex v, vtxsCellNoSharedWall){
      v->selected = true;
    }
    std::cout << "blas " << "/" <<cellVtxs.size() << "/" << vtxsCellNoSharedWall.size() << std::endl;


    if(actualOK and divCenter == "Center Actual"){
      CellDivisionAttr& actualP = actualPlane[selectedParent];
      divPointPos = actualP.planePos;
      std::cout << "actual center: " << divPointPos << std::endl;
    } else if(!actualOK and divCenter == "Center Actual"){
      return setErrorMessage("Couldn't compute actual plane and can't use it's center for Div Analysis!");
    }

    std::vector<CellDivisionAttr> possibleDivisions;

    Point3d divPointNrml;

    // if(actualOK){
    //   CellDivisionAttr& t = actualPlane[selectedParent];
    //   possibleDivisions = testDivision3D(m, verticesCell, divPlanes, centerDisplSteps, centerDisplStepSize, divPointPos, divPointNrml, selLabels, selectedParent, testActual, t);
    // } else {
    //   CellDivisionAttr t;
    //   possibleDivisions = testDivision3D(m, verticesCell, divPlanes, centerDisplSteps, centerDisplStepSize, divPointPos, divPointNrml, selLabels, selectedParent, testActual, t);
    // }

    std::vector<P3dDouble> samplingDirections;
    createSamplingDirsNew(divPlanes, false, Point3d(0,0,0), centerDisplSteps, centerDisplStepSize, samplingDirections);

    std::vector<CellDivisionAttr> possibleDivisionsCDA;

    //std::cout << "testDiv2D bef2 " << verticesCell.size() << "/" << selLabels.size() << std::endl;

    forall(int l, selLabels){
    std::cout << "label " << l << std::endl;
    }

    // get the cell tris with out the common wall
    CellDivisionMGXM cd;
    cd.cellTrisAll = getTrisOfCells(m, selectedParent, true);
    std::cout << "cdddddddd " << selectedParent << "/"  << cd.cellTrisAll.size() << std::endl;
    cd.initialize(m, selectedParent);

    if(divCenter == "Centroid"){
      divPointPos = centroidsP[selectedParent];
    }

    forall(P3dDouble divPlane, samplingDirections){
      CellDivisionAttr cda = cd.simulateDivision3D(divPlane, divPointPos);
      possibleDivisionsCDA.push_back(cda);
    }

    double minArea = HUGE_VAL;
    double maxArea = -HUGE_VAL;
    forall(CellDivisionAttr d, possibleDivisionsCDA){
      // find min and max of values
      if(d.planeArea < minArea) minArea = d.planeArea;
      if(d.planeArea > maxArea) maxArea = d.planeArea;
    }

    //std::cout << "minAA " << minArea << "/" << maxArea << "/" << possibleDivisionsCDA.size() << std::endl;

    forall(CellDivisionAttr& d, possibleDivisionsCDA){
      d.planeAreaRelShortest = d.planeArea / minArea;
      d.cellLabel = selectedLabel;
      if(maxArea > minArea) {
        d.normalizedArea = (d.planeArea - minArea) / (maxArea - minArea);
        if(d.normalizedArea < 0) d.normalizedArea = 0;
        if(d.normalizedArea > 1) d.normalizedArea = 1;
      } else {
        d.normalizedArea = -1;
      }

      if(actualOK){
        CellDivisionAttr& actual = actualPlane[selectedParent];
        actual.planeAreaRelShortest = actual.planeArea / minArea;
        actual.angleActual = 0;
        d.planeAreaRelActual = d.planeArea / actual.planeArea;
        d.angleActual = angleVectors(d.planeNrml, actual.planeNrml, false);
        if(maxArea > minArea) {
          actual.normalizedArea = (actual.planeArea - minArea) / (maxArea - minArea);
          if(d.normalizedArea < 0) d.normalizedArea = 0;
          if(d.normalizedArea > 1) d.normalizedArea = 1;
        } else {
          d.normalizedArea = -1;
        }
    
      } else {
        d.planeAreaRelActual = -1;
        d.angleActual = -1;
        d.normalizedArea = -1;
      }
    }
    possibleDivisions = possibleDivisionsCDA;


    int counter = attrData.size();
    forall(auto cda, possibleDivisions){
      Point3d pointPlane = projectPointOnPlane(actualPlane[selectedParent].planePos, cda.planePos, cda.planeNrml);
      cda.distanceActual = norm(pointPlane - actualPlane[selectedParent].planePos);

      pointPlane = projectPointOnPlane(cellCentroid, cda.planePos, cda.planeNrml);
      cda.distanceCentroid = norm(pointPlane - cellCentroid);

      cda.cellLabel = selectedParent;

      attrData[counter] = cda;
      counter++;
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(TestDivisionPlanes3D);


   bool DivisionAnalysisMulti3D::run(Mesh *m, Mesh *m2, bool singleCells, bool doubleCells)
    {

      std::cout << "Div analysis Multi" << std::endl;

      if(!m->useParents()) return setErrorMessage("This process only works with activated parent labels!");

      // find all cells
      vvGraph& S = m->graph();
      std::map<int, std::set<int> > parentCellLabelMap;
      typedef std::pair<int, std::set<int> > IntIntSetPair;

      // create parents -> cells map
      forall(const vertex& v, S){
        // std::cout << "ll " << v->label << "/" << m->parents()[v->label] << std::endl;
        if(m->parents()[v->label] > 0){
          parentCellLabelMap[m->parents()[v->label]].insert(v->label);
        }
      }

      std::cout << "multi3d " << parentCellLabelMap.size() << std::endl;

      forall(const IntIntSetPair& p, parentCellLabelMap){
        int currentParent = p.first;
        // run the div analysis for all pairs of same parent labels
        if(p.second.size() == 2 and doubleCells){

          std::cout << "Test Cell Pair with parent label " << currentParent << std::endl;

          TestDivisionPlanes3D *pDiv3D;
          if(!getProcess("Mesh/Division Analysis/Analysis 3D/Division Analysis", pDiv3D))
            throw(QString("Division Analysis:: Unable to make 3D Division Analysis process"));

          bool res = pDiv3D->run(m, m2, pDiv3D->parm("Nr of Planes").toInt(), pDiv3D->parm("Steps Center Displacement").toInt(),
                      pDiv3D->parm("Stepsize Center Displacement").toDouble(), pDiv3D->parm("Division Point"), currentParent);

          if(!res) std::cout << "Problems when analyzing parent labelled cell pair " << currentParent <<". Cells Ignored!" << std::endl;
          else std::cout << "Cell pair " << currentParent << " done!" << std::endl;

        // run div analysis for single cells
        } else if(p.second.size() == 1 and singleCells){

          std::cout << "Test Cell with parent label " << currentParent << std::endl;

          TestDivisionPlanes3D *pDiv3D;
          if(!getProcess("Mesh/Division Analysis/Analysis 3D/Division Analysis", pDiv3D))
            throw(QString("Division Analysis:: Unable to make 3D Division Analysis process"));

          bool res = pDiv3D->run(m, m2, pDiv3D->parm("Nr of Planes").toInt(), pDiv3D->parm("Steps Center Displacement").toInt(),
                      pDiv3D->parm("Stepsize Center Displacement").toDouble(), pDiv3D->parm("Division Point"), currentParent);

          if(!res) std::cout << "Problems when analyzing parent labelled cell " << currentParent <<". Cell Ignored!" << std::endl;
          else std::cout << "Cell " << currentParent << " done!" << std::endl;

        }
      }

      return true;
    }
    REGISTER_PROCESS(DivisionAnalysisMulti3D);



}
