#include "CellDivision.hpp"


using namespace std;

namespace mgx
{

  // return if triangle has to be split (upon cell division)
  bool CellDivision::triToBeSplit(const vertex& v, const vertex& m, const vertex& n)
  {
    if(vtxDaughterMap[v] == vtxDaughterMap[m] and vtxDaughterMap[v] == vtxDaughterMap[n]) return false;
    return true;
  }

  // label existing vertices according to their position regarding the division plane to find which daughter cell a vertex will belong to
  void CellDivision::newLabelDaughters()
  {

    forall(const vertex& v, c->S){
      double sP;
      Point3d intersectP;
      planeLineIntersect(divPlane.pos, divPlane.nrml, v->pos, v->pos-divPlane.nrml, sP, intersectP);
      if(sP<0){
        vtxDaughterMap[v] = dL->label;
        dL->S.insert(v);
      } else {
        vtxDaughterMap[v] = dR->label;
        dR->S.insert(v);
      }
    }
  }

  // sorts tris of cell into 3 bins, all vtxs have one of two labels (left or right daughter cell)
  void CellDivision::findCellTris()
  {
    std::vector<Triangle> tAll;

    forall(const vertex& v, c->S){
      forall(const vertex& n, c->S.neighbors(v)){
        vertex m = c->S.nextTo(v,n);
        if(c->S.uniqueTri(v,n,m)){
          Triangle t(v,n,m);
          // tri t;
          // t[0] = v;
          // t[1] = n;
          // t[2] = m;
          tAll.push_back(t);
        }
      }
    }
    trisAll = tAll;

  }

  // sorts tris of cell into 3 bins, all vtxs have one of two labels (left or right daughter cell)
  void CellDivision::sortCellTris()
  {

    forall(Triangle& t, trisAll){
      if(triToBeSplit(t.v[0],t.v[1],t.v[2])){
        trisSplit.insert(t);
      } else {
        if(vtxDaughterMap[t.v[0]] == labelLeft){
          trisLeft.push_back(t);
          splitTrisAreas[t] = make_pair(triangleArea(t.v[0]->pos, t.v[1]->pos, t.v[2]->pos),0);
        } else {
          trisRight.push_back(t);
          splitTrisAreas[t] = make_pair(0,triangleArea(t.v[0]->pos, t.v[1]->pos, t.v[2]->pos));
        }
      }
    }
  }

  // calculates the segment where the division plane cuts through the tri
  //res.first ist the point segment ab is cut, res.second: a\c
  pair<Point3d, Point3d> cuttingSegmentOfTri(Point3d oddTriPoint, Point3d triP2, Point3d triP3, Point3d planeNrml, Point3d planeP)
  {
    pair<Point3d, Point3d> res;

    double sP;
    Point3d intersectP;
    planeLineIntersect(planeP, planeNrml, oddTriPoint, triP2, sP, intersectP);

    double sP2;
    Point3d intersectP2;
    planeLineIntersect(planeP, planeNrml, oddTriPoint, triP3, sP2, intersectP2);

    // just for debug
    if(sP > 1 or sP < 0 or sP2 > 1 or sP2 < 0) cout << "sth wrong in cutting seg " << sP << "/" << sP2 << "/" << intersectP << "/" << intersectP2 << endl;

    res.first = intersectP;
    res.second = intersectP2;
    return res;
  }


  // trisToBeSplit will always have 2 points with a common label and the third point with a different label
  vector<pair<Point3d, Point3d> > CellDivision::polygonCuttingSegs()
  {
    Point3d planeNrml = divPlane.nrml;
    Point3d planeP = divPlane.pos;

    vector<pair<Point3d, Point3d> > res;

    std::map<VtxVtx, vertex> segmentsToBeSplitNew;

    Point3d avgP (0,0,0);

    vertex m; // midpoint of div plane
    m->pos = planeP;


    // find the single vtx
    forall(const Triangle& t, trisSplit){
      pair<Point3d, Point3d> seg;
      // order tri corner vertices
      vertex a,b,c;
      if(vtxDaughterMap[t.v[0]] != vtxDaughterMap[t.v[1]] and vtxDaughterMap[t.v[1]] == vtxDaughterMap[t.v[2]]){ // t0 is the odd one
        a = t.v[0]; b = t.v[1]; c = t.v[2];
      } else if(vtxDaughterMap[t.v[1]] != vtxDaughterMap[t.v[0]] and vtxDaughterMap[t.v[0]] == vtxDaughterMap[t.v[2]]){ // t1
        a = t.v[1]; b = t.v[2]; c = t.v[0];
      } else if(vtxDaughterMap[t.v[2]] != vtxDaughterMap[t.v[0]] and vtxDaughterMap[t.v[0]] == vtxDaughterMap[t.v[1]]){ // t2
        a = t.v[2]; b = t.v[0]; c = t.v[1];
      } else {
        cout << "sth wrong in polygonCuttingSegs" << endl;
      }
      // calc the cutting segment of the tri
      seg = cuttingSegmentOfTri(a->pos, b->pos, c->pos, planeNrml, planeP);

      vertex d,e;
      d->pos = seg.first;
      e->pos = seg.second;

      Triangle t1(a,d,e);
      Triangle t2(d,b,e);
      Triangle t3(e,b,c);

      Triangle t4(e,d,m);
      Triangle t5(d,e,m);

      double a1 = triangleArea(a->pos, d->pos, e->pos);
      double a2 = triangleArea(d->pos, b->pos, e->pos);
      double a3 = triangleArea(e->pos, b->pos, c->pos);

      if(vtxDaughterMap[a] == labelRight){
        trisRight.push_back(t1);
        trisRight.push_back(t4);
        splitTrisAreas[t] = make_pair(a1, a2+a3);

        trisLeft.push_back(t2);
        trisLeft.push_back(t3);
        trisLeft.push_back(t5);

      } else {
        trisLeft.push_back(t1);
        trisLeft.push_back(t4);
        splitTrisAreas[t] = make_pair(a2+a3, a1);
        trisRight.push_back(t2);
        trisRight.push_back(t3);
        trisRight.push_back(t5);
      }
      // legacy from old procedure
      vertex v1;
      v1->pos = seg.first;
      vertex v2;
      v2->pos = seg.second;
      VtxVtx ab = make_pair(a,b);
      if(b<a) ab = make_pair(b,a);
      VtxVtx ac = make_pair(a,c);
      if(c<a) ac = make_pair(c,a);
      segmentsToBeSplitNew[ab] = v1;
      segmentsToBeSplitNew[ac] = v2;

      res.push_back(seg);
    }
    segmentsToBeSplit = segmentsToBeSplitNew;
    return res;
  }
/*
//THIS IS THE OLD PROCEDURE
  // trisToBeSplit will always have 2 points with a common label and the third point with a different label
  // go through all trisToBeSplit, create a vector of all segments that are on the boundary of the new cell wall
  vector<pair<Point3d, Point3d> > CellDivision::polygonCuttingSegsSimple()
  {
    vector<pair<Point3d, Point3d> > res;

    std::map<VtxVtx, vertex> segmentsToBeSplitNew;

    Point3d planeNrml = divPlane.nrml;
    Point3d planeP = divPlane.pos;

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

    // find the single vtx
    forall(const tri& t, trisSplit){
      pair<Point3d, Point3d> seg;
      vertex a,b,c;
      if(vtxDaughterMap[t[0]] != vtxDaughterMap[t[1]] and vtxDaughterMap[t[1]] == vtxDaughterMap[t[2]]){ // t0 is the odd one
        a = t[0]; b = t[1]; c = t[2];
      } else if(vtxDaughterMap[t[1]] != vtxDaughterMap[t[0]] and vtxDaughterMap[t[0]] == vtxDaughterMap[t[2]]){ // t1
        a = t[1]; b = t[2]; c = t[0];
      } else if(vtxDaughterMap[t[2]] != vtxDaughterMap[t[0]] and vtxDaughterMap[t[0]] == vtxDaughterMap[t[1]]){ // t2
        a = t[2]; b = t[0]; c = t[1];
      } else {
        cout << "sth wrong in polygonCuttingSegs" << endl;
      }

      seg = cuttingSegmentOfTri(a->pos, b->pos, c->pos, planeNrml, planeP);
      res.push_back(seg);
      vertex v1;
      v1->pos = seg.first;
      vertex v2;
      v2->pos = seg.second;
      VtxVtx ab = make_pair(a,b);
      if(b<a) ab = make_pair(b,a);
      VtxVtx ac = make_pair(a,c);
      if(c<a) ac = make_pair(c,a);
      segmentsToBeSplitNew[ab] = v1;
      segmentsToBeSplitNew[ac] = v2;
    }

    segmentsToBeSplit = segmentsToBeSplitNew;
    return res;
  }
*/


  void addTrisOnLongSide2(MeshBuilder& mb, const vertex &v, const vertex &n, const vertex &m, std::map<VtxVtx, vertex> &segmentsToBeSplit,
    Point3d cellCenter, bool splitN)
  {

    Point3d avgDir = (v->pos + n->pos + m->pos) /3.0 - cellCenter;
    VtxVtx vn, vm, nm;
    if(v<n) vn = std::make_pair(v,n);
    else vn = std::make_pair(n,v);
    if(v<m) vm = std::make_pair(v,m);
    else vm = std::make_pair(m,v);
    if(n<m) nm = std::make_pair(n,m);
    else nm = std::make_pair(m,n);

    vertex& vnV = segmentsToBeSplit[vn];
    vertex& vmV = segmentsToBeSplit[vm];
    //vertex nmV = segmentsToBeSplit[nm];
    if(vnV->pos.x() == 0 or vmV->pos.x() == 0)
      std::cout << "bingo1 " << std::endl;
    else {
      if(splitN){
        mb.addTriCheckOrient(vnV, n, vmV, avgDir, true);
        mb.addTriCheckOrient(n, m, vmV, avgDir, true);
      } else {
        mb.addTriCheckOrient(vnV, n, m, avgDir, true);
        mb.addTriCheckOrient(vnV, m, vmV, avgDir, true);
      }
    }
  }

  void addTrisOnShortSide2(MeshBuilder& mb, const vertex &v, const vertex &n, const vertex &m, std::map<VtxVtx, vertex> &segmentsToBeSplit,
    Point3d cellCenter)
  {

    Point3d avgDir = (v->pos + n->pos + m->pos) /3.0 - cellCenter;
    VtxVtx vn, vm, nm;
    if(v<n) vn = std::make_pair(v,n);
    else vn = std::make_pair(n,v);
    if(v<m) vm = std::make_pair(v,m);
    else vm = std::make_pair(m,v);
    if(n<m) nm = std::make_pair(n,m);
    else nm = std::make_pair(m,n);

    vertex& vnV = segmentsToBeSplit[vn];
    vertex& vmV = segmentsToBeSplit[vm];
    //vertex nmV = segmentsToBeSplit[nm];
    if(vnV->pos.x() == 0 or vmV->pos.x() == 0)
      std::cout << "bingo2 " << std::endl;
    else {
      mb.addTriCheckOrient(v, vnV, vmV, avgDir, true);
    }
  }

  void CellDivision::createDaughterMeshes()
  {

    dL->label = labelLeft;
    dR->label = labelRight;

    MeshBuilder mbLeft;
    MeshBuilder mbRight;

    forall(std::vector<Point3d> pointsPolygon, pointsPolygonMulti){
      //pointsPolygon = p;
      std::vector<Point3i> triList;
      std::vector<Point3d> ptList;
      bool boundary = false;
      double maxArea = 1.0;

      //std::cout << "before " << pointsPolygon.size() << "//" << polygonSegsNew2.size() << std::endl;

      triangulatePolygon3D(maxArea, divPlane.nrml, pointsPolygon, triList, ptList, boundary, true, true);

      //std::cout << "after " << ptList.size() << "/" << triList.size() << std::endl;

      double eps = 0.0001;

    // add tris of triangulated plane to the list of vtxs and tris of the new cell
      for(uint i = 0; i < triList.size(); ++i) {
        Point3d a(ptList[triList[i].x()].x(), ptList[triList[i].x()].y(), ptList[triList[i].x()].z());
        Point3d b(ptList[triList[i].z()].x(), ptList[triList[i].z()].y(), ptList[triList[i].z()].z());
        Point3d c(ptList[triList[i].y()].x(), ptList[triList[i].y()].y(), ptList[triList[i].y()].z());

        vertex v1;
        v1->pos = a;
        vertex v2;
        v2->pos = b;
        vertex v3;
        v3->pos = c;

        // if very close to the div plane polygon border move the vertex a small amount
        forall(const Vtxvtxp& p, segmentsToBeSplit){
          if(norm(a-p.second->pos) < eps){
            v1 = p.second;
          }
          if(norm(b-p.second->pos) < eps){
            v2 = p.second;
          }
          if(norm(c-p.second->pos) < eps){
            v3 = p.second;
          }
        }

        mbLeft.addTri(v1, v2, v3);
        mbRight.addTri(v2, v1, v3);
      }
    }

    // determine vertices in left and right daughter cell (and relabel them to remember for later)
    //std::map<vertex,int> vtxDaughterMap; // 1 for left, -1 for right
    double epsi = -0.0001;
    forall(const vertex& v, c->S){
      double sP;
      Point3d intersectP;
      //planeLineIntersect(cellCenter, divPlane.nrml, v->pos, v->pos-divPlane.nrml, sP, intersectP);
      planeLineIntersect(divPlane.pos, divPlane.nrml, v->pos, v->pos-divPlane.nrml, sP, intersectP);

      if(norm(intersectP-v->pos)<epsi){ // vtxs too close to boundary all go to the right cell
        std::cout << "intersect " << intersectP << "/" << v->pos << std::endl;
        v->pos = intersectP;
        vtxDaughterMap[v] = -1;
        v->label = labelRight;
      } else if(sP<0){
        v->label = labelLeft;
        vtxDaughterMap[v] = 1;
      } else {
        v->label = labelRight;
        vtxDaughterMap[v] = -1;
      }
    }

    // check all neighbors of c, whether they share a vertex of

    int idxRefTriLeft = -1;
    int idxRefTriRight = -1;

    // get triangles of open halves
    // connect division plane and "open" halves
    forall(const vertex& v, c->S){
      forall(const vertex& n, c->S.neighbors(v)){
        vertex m = c->S.nextTo(v,n);
        if(c->S.uniqueTri(v,n,m)){
          if(vtxDaughterMap[v] == vtxDaughterMap[m] and vtxDaughterMap[v] == vtxDaughterMap[n]){ // first the easy triangles (that will stay the same)
            if(vtxDaughterMap[v]==1){
              // addtri left
              mbLeft.addTri(v, n, m, true);
              if(idxRefTriLeft == -1) idxRefTriLeft = mbLeft.triVec.size();
            } else if(vtxDaughterMap[v]==-1){//} else if(v->label == labelRight) {
              // addtri right
              mbRight.addTri(v, n, m, true);
              if(idxRefTriRight == -1) idxRefTriRight = mbLeft.triVec.size();
            }
          } else { // now the triangles that have to be divided
            if(vtxDaughterMap[v]==1 and vtxDaughterMap[n]==-1 and vtxDaughterMap[m]==-1){ // LRR
              addTrisOnShortSide2(mbLeft, v, n, m, segmentsToBeSplit, divPlane.pos);
              addTrisOnLongSide2(mbRight, v, n, m, segmentsToBeSplit, divPlane.pos, true);
            } else if(vtxDaughterMap[v]==-1 and vtxDaughterMap[n]==1 and vtxDaughterMap[m]==1){ // RLL
              addTrisOnLongSide2(mbLeft, v, n, m, segmentsToBeSplit, divPlane.pos, true);
              addTrisOnShortSide2(mbRight, v, n, m, segmentsToBeSplit, divPlane.pos);
            } else if(vtxDaughterMap[v]==1 and vtxDaughterMap[n]==1 and vtxDaughterMap[m]==-1){ // LLR
              addTrisOnLongSide2(mbLeft, m, v, n, segmentsToBeSplit, divPlane.pos, true);
              addTrisOnShortSide2(mbRight, m, v, n, segmentsToBeSplit, divPlane.pos);
            } else if(vtxDaughterMap[v]==-1 and vtxDaughterMap[n]==-1 and vtxDaughterMap[m]==1){ // RRL
              addTrisOnShortSide2(mbLeft, m, v, n, segmentsToBeSplit, divPlane.pos);
              addTrisOnLongSide2(mbRight, m, v, n, segmentsToBeSplit, divPlane.pos, true);
            } else if(vtxDaughterMap[v]==1 and vtxDaughterMap[n]==-1 and vtxDaughterMap[m]==1){ // LRL
              addTrisOnLongSide2(mbLeft, n, m, v, segmentsToBeSplit, divPlane.pos, true);
              addTrisOnShortSide2(mbRight, n, m, v, segmentsToBeSplit, divPlane.pos);
            } else if(vtxDaughterMap[v]==-1 and vtxDaughterMap[n]==1 and vtxDaughterMap[m]==-1){ // RLR
              addTrisOnShortSide2(mbLeft, n, m, v, segmentsToBeSplit, divPlane.pos);
              addTrisOnLongSide2(mbRight, n, m, v, segmentsToBeSplit, divPlane.pos, true);
            }
          }
            if(v->pos.x() == 0 or n->pos.x() == 0 or m->pos.x() == 0)
              std::cout << "tri " << v->pos << "/" << n->pos << "/" << m->pos << std::endl;
        }

      }

    }

    //std::cout << "final " << Vleft.size() << "/" << Tleft.size() << "/" << Vright.size() << "/" << Tright.size() << std::endl;

    // now take care of the neighbors
    forall(const cell& cn, C.neighbors(c)){
      cout << "cell " << cn->label << endl;
      MeshBuilder mbNeighbor;
      forall(const vertex& v, cn->S){
        forall(const vertex& n, cn->S.neighbors(v)){
          vertex m = cn->S.nextTo(v,n);
          if(!cn->S.uniqueTri(v,n,m)) continue;

          VtxVtx vn = make_pair(v,n);
          if(n<v) vn = make_pair(n,v);
          VtxVtx vm = make_pair(v,m);
          if(m<v) vm = make_pair(m,v);
          VtxVtx nm = make_pair(n,m);
          if(m<n) nm = make_pair(m,n);

          Point3d avgDir = (n->pos - v->pos) % (m->pos - v->pos);//(v->pos + n->pos + m->pos) /3.0 - cn->pos;

          // case 1: no edge of tri is influenced by division -> just take the tri as it is
          if(segmentsToBeSplit.find(vn) == segmentsToBeSplit.end() and segmentsToBeSplit.find(vm) == segmentsToBeSplit.end()
            and segmentsToBeSplit.find(nm) == segmentsToBeSplit.end()){
            mbNeighbor.addTri(v, n, m);
          } else {

            // case 2: two edges of tri is influenced by division (these are all shared tris of the common wall that are cut by the plane)
            if(segmentsToBeSplit.find(vn) != segmentsToBeSplit.end() and segmentsToBeSplit.find(vm) != segmentsToBeSplit.end()){
//cout << "1" << endl;
            addTrisOnShortSide2(mbNeighbor,v, n, m, segmentsToBeSplit, cn->pos);
            addTrisOnLongSide2(mbNeighbor,v, n, m, segmentsToBeSplit, cn->pos, false);
            } else if(segmentsToBeSplit.find(vn) != segmentsToBeSplit.end() and segmentsToBeSplit.find(nm) != segmentsToBeSplit.end()){
//cout << "2" << endl;
            addTrisOnShortSide2(mbNeighbor,n, m, v, segmentsToBeSplit, cn->pos);
            addTrisOnLongSide2(mbNeighbor,n, m, v, segmentsToBeSplit, cn->pos, false);
            } else if(segmentsToBeSplit.find(vm) != segmentsToBeSplit.end() and segmentsToBeSplit.find(nm) != segmentsToBeSplit.end()){
//cout << "3" << endl;
            addTrisOnShortSide2(mbNeighbor,m, v, n, segmentsToBeSplit, cn->pos);
            addTrisOnLongSide2(mbNeighbor,m, v, n, segmentsToBeSplit, cn->pos, false);
            // case 3: one edges of tri is influenced by division (these are tris that only share an edge with a tri that is cut by the plane)
            } else if(segmentsToBeSplit.find(vn) != segmentsToBeSplit.end()){
              //cout << "4" << endl;
              mbNeighbor.addTriCheckOrient(v, segmentsToBeSplit[vn], m, avgDir);
              mbNeighbor.addTriCheckOrient(segmentsToBeSplit[vn], n, m, avgDir);
            } else if(segmentsToBeSplit.find(vm) != segmentsToBeSplit.end()){
              //cout << "5" << endl;
              mbNeighbor.addTriCheckOrient(n, segmentsToBeSplit[vm], v, avgDir);
              mbNeighbor.addTriCheckOrient(segmentsToBeSplit[vm], n, m, avgDir);
            } else if(segmentsToBeSplit.find(nm) != segmentsToBeSplit.end()){
              //cout << "6" << endl;
              mbNeighbor.addTriCheckOrient(n, segmentsToBeSplit[nm], v, avgDir);
              mbNeighbor.addTriCheckOrient(segmentsToBeSplit[nm], m, v, avgDir);
            } else {
              //cout << "777777777777777" << endl;
            }
          }

        }
      }
      vvGraph Snew;
      //triangleOrientationCheck(mbNeighbor.triVec);
      if(!Mesh::meshFromTriangles(Snew, mbNeighbor.vtxVec, mbNeighbor.triVec)){
        std::cout << "Error, cannot add all the triangles (Neighbors)" << std::endl;
      }
      cn->S = Snew;
    }

/*
    // correct orientation
    forall(vertex& v, mbLeft.vtxVec){

    }
    Point3d avgDir = (mbLeft.triVec[0]->pos + mbLeft.triVec[1]->pos + mbLeft.triVec[2]->pos) /3.0 - cellCenter;
      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);
      }
*/
    // check triVec for connected components, only keep largest

    //std::cout << "size " << mbLeft.triVec.size() << "/" << mbRight.triVec.size() << std::endl;

    mbLeft.eraseSmallComponents();
    mbRight.eraseSmallComponents();

    //std::cout << "size2 " << mbLeft.triVec.size() << "/" << mbRight.triVec.size() << std::endl;

    triangleOrientationCheck(mbLeft.triVec, idxRefTriLeft);
    triangleOrientationCheck(mbRight.triVec, idxRefTriRight);


    vvGraph test;

    if(!Mesh::meshFromTriangles(dL->S, mbLeft.vtxVec, mbLeft.triVec)){
      std::cout << "Error, cannot add all the triangles (Left)";
    }


    if(!Mesh::meshFromTriangles(dR->S, mbRight.vtxVec, mbRight.triVec)){
      std::cout << "Error, cannot add all the triangles (Right)";
    }

    // erase vtxs that dont have edges
    std::set<vertex> toErase;
    forall(const vertex& v, dL->S){
      int n = 0;
      forall(const vertex& m, dL->S.neighbors(v)){
        n++;
      }
      if(n==0) toErase.insert(v);
    }
    forall(const vertex& v, toErase){
      dL->S.erase(v);
    }

    toErase.clear();
    forall(const vertex& v, dR->S){
      int n = 0;
      forall(const vertex& m, dR->S.neighbors(v)){
        n++;
      }
      if(n==0) toErase.insert(v);
    }
    forall(const vertex& v, toErase){
      dR->S.erase(v);
    }

  }

  // update the cell graph: insert both daughter cells will all their neighbors and delete father cell
  void CellDivision::divideCellInCellGraph()
  {

    std::map<int, cell> labelCellMap;

    neighbAll = 0;
    neighbLeft = 0;
    neighbRight = 0;

    forall(const cell& cn, C.neighbors(c)){
      labelCellMap[cn->label] = cn;
      neighbAll++;
    //  cout << "l" << cn->label << endl;
    }

    forall(const IntIntDouble& p, nbhdMap){
      if(p.first.first == dL->label and p.second > 0.00001 and p.first.second != dR->label){
        C.insertEdge(labelCellMap[p.first.second], dL);
        C.insertEdge(dL, labelCellMap[p.first.second]);
      //  cout << "c" << p.first.second << "/" << p.second << endl;
      }
      if(p.first.first == dR->label and p.second > 0.00001 and p.first.second != dL->label){
        C.insertEdge(labelCellMap[p.first.second], dR);
        C.insertEdge(dR, labelCellMap[p.first.second]);
      //  cout << "c" << p.first.second << "/" << p.second << endl;
      }
    }

    forall(const cell& cn, C.neighbors(dL)){
      neighbLeft++;
    }
    forall(const cell& cn, C.neighbors(dR)){
      neighbRight++;
    }

    std::map<Point3d,int> posLabelMap;
    forall(const vertex& v, dL->S){
      posLabelMap[v->pos] = dL->label;
    }

    forall(const vertex& v, dR->S){
      if(posLabelMap.find(v->pos) != posLabelMap.end()){
        posLabelMap[v->pos] = -2;
      } else {
        posLabelMap[v->pos] = dR->label;
      }
    }

    std::set<cell> neighborLeft, neighborRight;
    // now the neighbor cells
    forall(const cell& cn, C.neighbors(c)){

      forall(const vertex& v, cn->S){
        if(posLabelMap.find(v->pos) != posLabelMap.end()){
          if(posLabelMap[v->pos] == -2 or posLabelMap[v->pos] == dL->label){
            neighborLeft.insert(cn);
          }
          if(posLabelMap[v->pos] == -2 or posLabelMap[v->pos] == dR->label){
            neighborRight.insert(cn);
          }
        }
      }
    }

    forall(const cell& c, neighborLeft){
      C.insertEdge(c, dL);
      C.insertEdge(dL, c);
    }

    forall(const cell& c, neighborRight){
      C.insertEdge(c, dR);
      C.insertEdge(dR, c);
    }
  }

  vector<vector<Point3d> > CellDivision::getCuttingPolygons(){
    //vector<vector<Point3d> > pointsPolygonMulti;

    sortCellTris();

    vector<pair<Point3d, Point3d> > segs = polygonCuttingSegs();
    //pointsPolygonMulti = orderPolygonSegsMulti(segs);
    return orderPolygonSegsMulti(segs);
  }

  // simulate a division:
  // check on which side of the div plane are the vtxs and update cell graph
  bool CellDivision::simulateDivision3D(std::map<IntInt, double> neighborMap, std::map<Triangle, int> trisNeighbMap)
  {
    // create daughter cells (without inserting them into the tissue)

    initialize();
    // determine vertices in left and right daughter cell

    forall(const TriInt& p, trisNeighbMap){
      trisAll.push_back(p.first);
    }

    newLabelDaughters();

    //sortCellTris();
    //estimateVolumeDaughterCells(trisLeft, trisRight, trisSplit);

    //vector<pair<Point3d, Point3d> > segs = polygonCuttingSegs();
    //vector<Point3d> pointsPolygon = orderPolygonSegs(segs);
    pointsPolygonMulti = getCuttingPolygons();//= orderPolygonSegsMulti(segs);

    divPlane.area = divisionPlaneAreaMulti(divPlane.pos, pointsPolygonMulti);
    //forall(std::vector<Point3d> p, pointsPolygonMulti){
    //  divPlane.area += divisionPlaneArea(divPlane.pos, p);
    //}
    //divPlane.area = divisionPlaneArea(divPlane.pos, pointsPolygon);

    createUpdatedNeighborhoodMap(neighborMap, trisNeighbMap);

    // update the cell graph
    // insert new cells as vtxs
    C.insert(dL);
    C.insert(dR);

    // add edges in between them
    C.insertEdge(dL, dR);
    C.insertEdge(dR, dL);

    // add edges to neighbors
    divideCellInCellGraph();

    volumeLeft = calcVolume(trisLeft);
    volumeRight = calcVolume(trisRight);
    volumeRatio = calcRatio(volumeLeft, volumeRight);

    return true;
  }

  void CellDivision::initialize(){

    // get the labels for the daughter cells (find highest label and add 1 and 2)
    int maxLabel = 0;
    forall(const cell& c, C){
      if(maxLabel < c->label) maxLabel = c->label;
    }
    labelLeft = maxLabel+1;
    labelRight = maxLabel+2;

    // to avoid complications in where vertices touch but are not cenntected
    separateIdenticalVtxs(c);

    // create daughter cells (without inserting them into the tissue)
    cell daughterLeft, daughterRight;
    dL = daughterLeft;
    dR = daughterRight;

    dL->label = labelLeft;
    dR->label = labelRight;

    // get the center of the cell to divide
    cellCenter = c->pos;

  }

  // Divide a cell
  bool CellDivision::divideCell3D(CellTissue &T)
  {

    std::cout << "divide1" << std::endl;
    // update cell properties one last time
    T.updGeometry(c, 4);

    // set labels, create cells
    initialize();

    // get all tris of cell
    findCellTris();

    // determine vertices in left and right daughter cell
    newLabelDaughters();

    // determine triangles that are cut by the division plane and the open polygons
    pointsPolygonMulti = getCuttingPolygons();

    // calculate the division plane area
    divPlane.area = divisionPlaneAreaMulti(divPlane.pos, pointsPolygonMulti);

    // create the actual meshes of the daughter cells
    createDaughterMeshes();

    // insert the daughters
    C.insert(dL);
    C.insert(dR);

    // update the cell graph
    C.insertEdge(dL, dR);
    C.insertEdge(dR, dL);
    divideCellInCellGraph();

    // get rid of the father cell
    C.erase(c);

    // set the new cell graph
    T.C = C;

    // update the cell tissue
    T.updateAllGraph();
    T.updGeometry();

    return true;
  }

  // write data to vec for saving as csv
  void CellDivision::writeDataToVec(std::vector<double>& d)
  {
    d.push_back(divPlane.nrml.x());
    d.push_back(divPlane.nrml.y());
    d.push_back(divPlane.nrml.z());
    d.push_back(divPlane.displ);
    d.push_back(divPlane.area);
    d.push_back(betw);
    d.push_back(rwBetw);
    d.push_back(calcRatio(volumeLeft, volumeRight));
    d.push_back(volumeLeft);
    d.push_back(volumeRight);
    d.push_back(cutLowDegreeNeighbor);
  }

  // update the neighborhoodmap after a division (or simulated division) (needed input for the graph measures)
  void CellDivision::createUpdatedNeighborhoodMap()
  {
    createUpdatedNeighborhoodMap(neighborMapOld,trisNeighbMap);
  }

  // update the neighborhoodmap after a division (or simulated division) (needed input for the graph measures)
  void CellDivision::createUpdatedNeighborhoodMap(std::map<IntInt, double>& neighborMap, std::map<Triangle, int>& trisNeighbMap)
  {
    std::map<IntInt, double> newNeighMap;

    double totalSum = 0;
    double totalSumF = 0;

    // copy old map (unchanged, new cells have new labels)
    forall(const IntIntDouble& p, neighborMap){
      newNeighMap[p.first] = p.second;
    }

    forall(const TriInt& p, trisNeighbMap){
      if(p.second != -1){ // tri in trisNeighbMap with a valid neighbor cell
        totalSum+=triangleArea(p.first.v[0]->pos, p.first.v[1]->pos, p.first.v[2]->pos);
        if(splitTrisAreas[p.first].first > 0){ // tri area that belongs to left cell
          IntInt w = make_pair(p.second, dL->label);
          IntInt w2 = make_pair(dL->label, p.second);
          newNeighMap[w]+=splitTrisAreas[p.first].first;
          newNeighMap[w2]+=splitTrisAreas[p.first].first;
          totalSumF+=splitTrisAreas[p.first].first;
        }
        if(splitTrisAreas[p.first].second > 0){ // tri area that belongs to right cell
          IntInt w = make_pair(p.second, dR->label);
          IntInt w2 = make_pair(dR->label, p.second);
          newNeighMap[w]+=splitTrisAreas[p.first].second;
          newNeighMap[w2]+=splitTrisAreas[p.first].second;
          totalSumF+=splitTrisAreas[p.first].second;
        }
      }
    }

    if(totalSum - totalSumF > 0.01) cout << "sums " << totalSum << "/" << totalSumF << endl;

    //cout << "sums " << totalSum << "/" << totalSumF << endl;

    // the area between the two daughter cells
    IntInt w = make_pair(dR->label, dL->label);
    IntInt w2 = make_pair(dL->label, dR->label);

    newNeighMap[w]=divPlane.area;
    newNeighMap[w2]=divPlane.area;

    nbhdMap = newNeighMap;
  }

  double CellDivision::getValue(QString criteria)
  {
    if(criteria == "Wall Area")
      return divPlane.area;
    else if(criteria == "Volume Ratio")
      return volumeRatio;
    else if(criteria == "Betweenness Centrality" or criteria == "Betweenness Centrality (weighted)")
      return betw;
    else if(criteria == "Random Walk Betweenness" or criteria == "Random Walk Betweenness (weighted)")
      return rwBetw;
    else if(criteria == "Neighbor Count")
      return neighbLeft+neighbRight;
    else
      return 0;

  }


  // return if triangle has to be split (upon cell division)
  bool CellDivisionMGXM::triToBeSplit(const vertex& v, const vertex& m, const vertex& n)
  {
    if(vtxDaughterMap[v] == vtxDaughterMap[m] and vtxDaughterMap[v] == vtxDaughterMap[n]) return false;
    return true;
  }


    // trisToBeSplit will always have 2 points with a common label and the third point with a different label
    vector<pair<Point3d, Point3d> > CellDivisionMGXM::polygonCuttingSegs(Point3d planeNrml, Point3d planeP)
    {

      vector<pair<Point3d, Point3d> > res;

      std::map<VtxVtx, vertex> segmentsToBeSplitNew;

      Point3d avgP (0,0,0);

      vertex m; // midpoint of div plane
      m->pos = planeP;


      // find the single vtx
      forall(const Triangle& t, trisSplit){
        pair<Point3d, Point3d> seg;
        // order tri corner vertices
        vertex a,b,c;
        if(vtxDaughterMap[t.v[0]] != vtxDaughterMap[t.v[1]] and vtxDaughterMap[t.v[1]] == vtxDaughterMap[t.v[2]]){ // t0 is the odd one
          a = t.v[0]; b = t.v[1]; c = t.v[2];
        } else if(vtxDaughterMap[t.v[1]] != vtxDaughterMap[t.v[0]] and vtxDaughterMap[t.v[0]] == vtxDaughterMap[t.v[2]]){ // t1
          a = t.v[1]; b = t.v[2]; c = t.v[0];
        } else if(vtxDaughterMap[t.v[2]] != vtxDaughterMap[t.v[0]] and vtxDaughterMap[t.v[0]] == vtxDaughterMap[t.v[1]]){ // t2
          a = t.v[2]; b = t.v[0]; c = t.v[1];
        } else {
          cout << "sth wrong in polygonCuttingSegs" << endl;
        }
        // calc the cutting segment of the tri
        seg = cuttingSegmentOfTri(a->pos, b->pos, c->pos, planeNrml, planeP);

        vertex d,e;
        d->pos = seg.first;
        e->pos = seg.second;

        Triangle t1(a,d,e);
        Triangle t2(d,b,e);
        Triangle t3(e,b,c);

        Triangle t4(e,d,m);
        Triangle t5(d,e,m);

        double a1 = triangleArea(a->pos, d->pos, e->pos);
        double a2 = triangleArea(d->pos, b->pos, e->pos);
        double a3 = triangleArea(e->pos, b->pos, c->pos);

        if(vtxDaughterMap[a] == labelRight){
          trisRight.push_back(t1);
          trisRight.push_back(t4);
          splitTrisAreas[t] = make_pair(a1, a2+a3);

          trisLeft.push_back(t2);
          trisLeft.push_back(t3);
          trisLeft.push_back(t5);

        } else {
          trisLeft.push_back(t1);
          trisLeft.push_back(t4);
          splitTrisAreas[t] = make_pair(a2+a3, a1);
          trisRight.push_back(t2);
          trisRight.push_back(t3);
          trisRight.push_back(t5);
        }
        // legacy from old procedure
        vertex v1;
        v1->pos = seg.first;
        vertex v2;
        v2->pos = seg.second;
        VtxVtx ab = make_pair(a,b);
        if(b<a) ab = make_pair(b,a);
        VtxVtx ac = make_pair(a,c);
        if(c<a) ac = make_pair(c,a);
        segmentsToBeSplitNew[ab] = v1;
        segmentsToBeSplitNew[ac] = v2;

        res.push_back(seg);
      }
      segmentsToBeSplit = segmentsToBeSplitNew;
      return res;
    }

  // sorts tris of cell into 3 bins, all vtxs have one of two labels (left or right daughter cell)
  std::vector<Triangle> findCellTris(vvGraph& S, int cellLabel)
  {
    std::vector<Triangle> trisCell;

    std::set<vertex> cellVtxs;

    forall(const vertex& v, S){
      if(v->label == cellLabel) cellVtxs.insert(v);
    }

    forall(const vertex& v, cellVtxs){
      forall(const vertex& n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
        if(S.uniqueTri(v,n,m)){
          Triangle t(v,n,m);
          trisCell.push_back(t);
        }
      }
    }
    return trisCell;
  }

}
