#ifndef CellDivision_H
#define CellDivision_H


#include <GraphAlgorithms.hpp>
#include <DivisionAnalysisUtils.hpp>
#include <Triangulate.hpp>
#include <MeshBuilder.hpp>


using namespace std;

namespace mgx
{


  // division plane properties
  struct DivPlane{
    Point3d nrml, pos;
    double displ;

    double area;

    // for drawing
    int label;
    double size;
    std::vector<vertex> vtxs;
    std::vector<Point3i> tris;

    DivPlane() {}
    DivPlane(Point3d nrml, Point3d pos, double displ): nrml(nrml), pos(pos), displ(displ) {}

  };


  class CellDivisionAttr{

    public:

    DivPlane divPlane;

    Point3d cellNormal;
    double volumeLeft, volumeRight, volumeRatio;
    double betw, rwBetw;
    double cutLowDegreeNeighbor;
    Point3d cellCenter;
    int neighbAll, neighbLeft, neighbRight;

    // new stuff
    Point3d planeNrml;
    Point3d planePos;
    double planeDisplacement;

    //Point3d cellNormal;

    int cellLabel;
    bool isParent;

    int parentLabel;
    int planeLabel;

    //double volumeRatio;
    double planeArea;
    double planeAreaRelShortest;
    double planeAreaRelActual;
    double daughterRatio;

    double distanceActual, distanceCentroid;

    double sharedTrisArea;

    //double rwBetw;
    //double centrBetw;

    double angleActual;

    //bool actualDiv; // delete
    double actualArea;

    double normalizedArea;

    Point3d surfNrml;


    int div2d;

    CellDivisionAttr() {}

      bool operator==(const CellDivisionAttr &other) const
      {
        if(planeNrml == other.planeNrml and planeDisplacement == other.planeDisplacement
            and cellLabel == other.cellLabel and parentLabel == other.parentLabel
            and volumeRatio == other.volumeRatio and planeArea == other.planeArea)
          return true;
        return false;
      }

  double 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;

  }


  };

    // Read/write attr map
    bool inline readAttr(CellDivisionAttr &m, const QByteArray &ba, size_t &pos)
    {
      return readChar((char *)&m, sizeof(CellDivisionAttr), ba, pos);
    }
    bool inline writeAttr(const CellDivisionAttr &m, QByteArray &ba)
    {
      return writeChar((char *)&m, sizeof(CellDivisionAttr), ba);
    }


  // new class for attr maps (single cell)
  class DivisionDataNew
  {
  public:

    std::vector<CellDivisionAttr> simulationData;

    DivisionDataNew() {}

      bool operator==(const DivisionDataNew &other) const
      {
        if(simulationData == other.simulationData)
          return true;
        return false;
      }


  };

    // Read/write attr map
    bool inline readAttr(DivisionDataNew &dd, const QByteArray &ba, size_t &pos)
    {

      uint sz;
      readAttr(sz, ba, pos);
      for(uint i = 0; i<sz; i++){
        CellDivisionAttr cda;
        readAttr(cda, ba, pos);
        dd.simulationData.push_back(cda);
      }

      return true;//readChar((char *)&m, sizeof(DivisionDataNew), ba, pos);
    }
    bool inline writeAttr(const DivisionDataNew &dd, QByteArray &ba)
    {

      uint sz = dd.simulationData.size();
      writeChar((char *)&sz, sizeof(uint), ba);
      for(uint i = 0; i < sz; i++){
        writeAttr(dd.simulationData[i], ba);
      }

      return true;//writeChar((char *)&m, sizeof(DivisionDataNew), ba);
    }



  class CellDivision{
    public:

    // the cells: c = cell to be divided, dL, dR: the left and right daughter
    cell c, dL, dR;

    DivPlane divPlane;

    Point3d cellNormal;

    std::vector<Triangle> trisLeft, trisRight, trisAll;
    std::set<Triangle> trisSplit;
    std::map<Triangle, std::pair<double,double> > splitTrisAreas;
    std::map<IntInt, double> nbhdMap;


    cellGraph C;

    //std::map<cell, std::vector<tri> > neighbTris;
    //std::map<cell, double> neighbWallArea;
    int labelLeft, labelRight;

    // map from a pair of vertices (edge in original vvgraph) to a new vertex that has to be introduced because the edge is cut by the division plane
    std::map<VtxVtx, vertex> segmentsToBeSplit;

    std::map<vertex,int> vtxDaughterMap; // 1 for left, -1 for right

    std::map<IntInt, double> neighborMapOld;
    std::map<Triangle, int> trisNeighbMap;

    // measures relevant for the division
    int neighbAll, neighbLeft, neighbRight;
    double volumeLeft, volumeRight, volumeRatio;
    double betw, rwBetw;
    std::map<IntInt, double> weights;
    int idxSamples;

    double cutLowDegreeNeighbor;

    // for the actual division
    Subdivide *subdiv;
    mgx::CellTissue::DivAlg divAlg;

    // new stuff to actually divide the cell
    vector<vector<Point3d> > pointsPolygonMulti;
    Point3d cellCenter;


    CellDivision() {}
    CellDivision(const cell& cellToDivide, DivPlane d, const cellGraph& C): c(cellToDivide), divPlane(d), C(C)
    {
      cellCenter = cellToDivide->pos;
    }
    CellDivision(P3dDouble plane, CellTissue& T, cell cellToDivide, int maxLabel)
    {
        T.updGeometry(c, 4);

        cellCenter = cellToDivide->pos;

        Point3d divPlanePoint = cellCenter + plane.second*plane.first;
        double distanceToCenter = plane.second;

        DivPlane p(plane.first, divPlanePoint, distanceToCenter);

        CellTissue Sim;
        copyCellGraph(T, Sim);

        c = cellToDivide;
        divPlane = p;
        C = Sim.C;
        labelLeft = maxLabel+2;
        labelRight = maxLabel+1;
    }

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

    // sorts tris of cell into 3 bins, all vtxs have one of two labels (left or right daughter cell)
    void findCellTris();
    void sortCellTris();
    // write data to vec for saving as csv
    void writeDataToVec(std::vector<double>& d);

    // update the neighborhoodmap (input for the graph measures)
    void createUpdatedNeighborhoodMap();

    // update the neighborhoodmap (input for the graph measures)
    void createUpdatedNeighborhoodMap(std::map<IntInt, double>& neighborMap, std::map<Triangle, int>& trisNeighbMap);

    vector<pair<Point3d, Point3d> > polygonCuttingSegs();

    // go through all trisToBeSplit, create a vector of all segments that are on the boundary of the new cell wall
    //vector<pair<Point3d, Point3d> > polygonCuttingSegsSimple();
    vector<vector<Point3d> > getCuttingPolygons();

    // split the segments of trisToBeSplit that are cut by the division plane
    void splitCuttingSegs();

    void divideCellInCellGraph();

    // simulate a division: calculate divPlaneArea & volumes, update cell graph (but not mesh)
    bool simulateDivision3D(std::map<IntInt, double> neighborMap, std::map<Triangle, int> trisNeighbMap);

    // do the proper division
    bool divideCell3D(CellTissue &T);

    double getValue(QString criteria);



    private:

    // initialization of variables
    void initialize();

    // returns true if the triangle vnm is cut by the divsion plane
    bool triToBeSplit(const vertex& v, const vertex& m, const vertex& n);

    // creates the meshes of the daughter cells and updates the meshes of the neighbor cells
    void createDaughterMeshes();


  };


std::vector<Triangle> findCellTris(vvGraph& S, int cellLabel);

   class CellDivisionMGXM{
     public:

    std::vector<Triangle> cellTrisAll, cellTrisShared, cellTrisOutside;

    std::set<Triangle> trisSplit;
    std::vector<Triangle> trisLeft, trisRight;

    std::set<vertex> cellVtxsOutside;
    std::set<vertex> vtxsLeft, vtxsRight;

    std::map<vertex, int> vtxDaughterMap;

    std::map<Triangle, std::pair<double,double> > splitTrisAreas;

    vector<vector<Point3d> > pointsPolygonMulti;

    int labelLeft;
    int labelRight;

    Point3d centroid;



    // map from a pair of vertices (edge in original vvgraph) to a new vertex that has to be introduced because the edge is cut by the division plane
    std::map<VtxVtx, vertex> segmentsToBeSplit;

     CellDivisionMGXM() {}

    void initialize(Mesh *m, int parentLabel)
    {
      // requirement: cellTrisAll is set!

      // find tris of shared wall and remove it
      cellTrisOutside = getOrRemoveCommonWall(m, cellTrisAll, true, true);
      //cellTrisShared = getOrRemoveCommonWall(m, cellTrisAll, false, true);

      double volumeTotal = 0;
      Point3d centroidP(0,0,0);

      forall(Triangle t, cellTrisOutside){
        vertex v = t.v[0];
        vertex n = t.v[1];
        vertex m = t.v[2];
        cellVtxsOutside.insert(v);
        cellVtxsOutside.insert(n);
        cellVtxsOutside.insert(m);


        double volume = signedTetraVolume(t.v[0]->pos, t.v[1]->pos, t.v[2]->pos);
        centroidP += volume * (t.v[0]->pos + t.v[1]->pos + t.v[2]->pos) / 4.0;
        volumeTotal += volume;
      }
      centroidP /= volumeTotal;

      AttrMap<int, Point3d>& centroidsP = m->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");
      centroidsP[parentLabel]=centroidP;
      std::cout << "pc  " << centroidP << "/" << parentLabel << std::endl;
    }


    void clearSimulatedData()
    {
      vtxDaughterMap.clear();
      vtxsLeft.clear();
      vtxsRight.clear();
      trisSplit.clear();
      trisLeft.clear();
      trisRight.clear();
      splitTrisAreas.clear();
    }

    CellDivisionAttr simulateDivision3D(P3dDouble divPlane, Point3d divPos)
    {
      clearSimulatedData();
      // requirement: initialize was ran!
      labelLeft = 1;
      labelRight = 2;

      Point3d divPlaneNrml = divPlane.first;
      Point3d divPlanePos = divPos + divPlane.second * divPlane.first;

      // check on which side the vtxs lie
      forall(const vertex& v, cellVtxsOutside){
        double sP;
        Point3d intersectP;
        planeLineIntersect(divPlanePos, divPlaneNrml, v->pos, v->pos-divPlaneNrml, sP, intersectP);
        if(sP<0){
          vtxDaughterMap[v] = labelLeft;
          vtxsLeft.insert(v);
        } else {
          vtxDaughterMap[v] = labelRight;
          vtxsRight.insert(v);
        }
      }

      std::cout << "vtxs " << vtxsLeft.size() << "/" << vtxsRight.size() << std::endl;

      // find the tris that will be split by the plane
      forall(Triangle& t, cellTrisOutside){
        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));
          }
        }
      }

      // find the ploygon cutting segments
      vector<pair<Point3d, Point3d> > segs = polygonCuttingSegs(divPlaneNrml, divPlanePos);

      // order the cutting segments
      pointsPolygonMulti = orderPolygonSegsMulti(segs);

      // compute the plane area
      double planeArea = divisionPlaneAreaMulti(divPlanePos, pointsPolygonMulti);

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

      std::cout << "divTest " << divPlanePos << "/" << divPlaneNrml << "/" << planeArea << "/"<< volumeLeft << "/" << volumeRight << std::endl;

      CellDivisionAttr cda;
      cda.div2d = 0;
      cda.planeNrml = divPlaneNrml;
      cda.planeDisplacement = divPlane.second;
      cda.daughterRatio = volumeRatio;
      cda.planeArea = planeArea;
      cda.planePos = divPlanePos;

      return cda;
    }

    // returns true if the triangle vnm is cut by the divsion plane
    bool triToBeSplit(const vertex& v, const vertex& m, const vertex& n);

    vector<pair<Point3d, Point3d> > polygonCuttingSegs(Point3d planeNrml, Point3d planeP);


   };

}
  #endif
