#include "DivisionAnalysis2D.hpp"

using namespace std;

namespace mgx { //namespace process



 bool FlatDivisionPlane2D::run(Mesh *m, Mesh *m2, bool display, double planeSize, bool resetMesh, int label)
  {

    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();
    m->updateCentersNormals();

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

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

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

    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]);
      }
      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();

    Point3d cellNormal(0,0,0);

    std::vector<Point3d> pointsSharedWall, pointsSharedWallExt;

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

    if(label > 0){
      vtxVec.clear();
      forall(const vertex& v, S){
        if(!m->useParents() and v->label == label)
          vtxVec2.push_back(v);
        if(m->useParents() and m->parents()[v->label] == label)
          vtxVec2.push_back(v);
      }
      forall(const vertex& v, vtxVec2){
        vtxVec.push_back(v);
        forall(vertex n, S.neighbors(v)){
          if(n->label == -1) vtxVec.push_back(n);
        }
      }
    }

    std::cout << "flat plane approx " << selectedParent << "/" << selParents.size() << "/" << vtxVec.size() << "/" << m->useParents() << "/" << label << std::endl;

    if(vtxVec.size() < 2) {
      return setErrorMessage("Please select exactly two neighboring cells with the same parent label");
    }

    forall(const vertex& v, vtxVec){
      
      int vLabel = m->parents()[v->label];
      if(v->label == -1) vLabel = -1;

      if(vLabel > -1 or v->margin) continue;

      bool vertexOk = true;

      forall(vertex n, S.neighbors(v)){
        int nLabel = m->parents()[n->label];
        if(n->label == -1) nLabel = -1;

        if(m->useParents() and (nLabel != selectedParent and nLabel != -1)) vertexOk=false;
        //if(!m->useParents() and selLabels.find(nLabel) == selLabels.end() and nLabel != -1) vertexOk=false;
      }
      if(vertexOk){
        pointsSharedWall.push_back(v->pos);
        cellNormal += v->nrml;
      } 
    }

    cellNormal /= pointsSharedWall.size();
    cellNormal /= norm(cellNormal);

    forall(Point3d p, pointsSharedWall){
      pointsSharedWallExt.push_back(p);
      pointsSharedWallExt.push_back(p+cellNormal);
      pointsSharedWallExt.push_back(p-cellNormal);
    }

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

    // enforce orthogonality to the surface
    Point3d orth = planeNrml ^ cellNormal;
    planeNrml = orth ^ cellNormal;

    CellDivisionAttr cda;
    cda.planeNrml = planeNrml;
    cda.planeDisplacement = 0; // TODO
    cda.planePos = planePos;
    cda.div2d = 1;
    cda.cellNormal = cellNormal;
    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, 0.9);
      m2->updateAll();
    }

    return true;

    m->updateAll();
  }
  REGISTER_PROCESS(FlatDivisionPlane2D);


  CellDivisionAttr convertCellDivAttr(CellDivision& cd)
  {
      CellDivisionAttr cda;
      cda.planeNrml = cd.divPlane.nrml;
      cda.planeDisplacement = cd.divPlane.displ;
      cda.daughterRatio = cd.volumeRatio;
      cda.planeArea = cd.divPlane.area;
      cda.planePos = cd.divPlane.pos;

      return cda;
  }

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

    Point3d cellNrml;
    forall(vertex v, verticesCell){
      cellNrml += v->nrml;
    }
    cellNrml /= verticesCell.size();
    cellNrml /= norm(cellNrml);

    if(norm(divPointNrml) < 0.99){
      divPointNrml = cellNrml;
    }

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

    std::set<Point3d> emptyMap;

    double sampling_res = 180.0;
    sampling_res=sampling_res/(double)divPlanes*(2.*centerDisplSteps+1.);
    generateSamplingDirections(sampling_res, centerDisplSteps, centerDisplStepSize, divPointNrml, emptyMap, false, samplingDirections);

    std::cout << "resolution/planes to be tested: " << sampling_res << "/" << samplingDirections.size() << std::endl;
    std::vector<CellDivision> possibleDivisions(samplingDirections.size());

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

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

    std::vector<vertex> borderVerticesOrdered;
    if(m->useParents()) findBorderVerticesOfCellsCombined(S, selLabels, verticesCell, borderVerticesOrdered);
    else findBorderVerticesOfCell(S, selectedLabel, verticesCell, borderVerticesOrdered);

    //m->updateAll();

    //std::cout << "testDiv2D bef " << borderVerticesOrdered.size() << "/" << samplingDirections.size() << std::endl;

    Point3d shortestPlaneNrml;
    double shortesPlaneDispl;
    //simulateDivision2D(S, vc, samplingDirections, shortestPlaneNrml, shortesPlaneDispl);
    simulateDivision2D(borderVerticesOrdered, divPointPos, divPointNrml, samplingDirections, possibleDivisions, shortestPlaneNrml, shortesPlaneDispl);

    std::vector<CellDivision> possibleDivisionsAct(1);
    if(testActual){
      std::vector<P3dDouble> samplingDirectionsAct;
      samplingDirectionsAct.push_back(std::make_pair(actual.planeNrml, 0.));
      simulateDivision2D(borderVerticesOrdered, actual.planePos, divPointNrml, samplingDirectionsAct, possibleDivisionsAct, shortestPlaneNrml, shortesPlaneDispl);
    }

    //std::cout << "testDiv2D " << possibleDivisions.size() << std::endl;

    std::vector<CellDivisionAttr> possibleDivisionsCDA;

    //int counter = 1;
    forall(auto cd, possibleDivisions){
      CellDivisionAttr cda = convertCellDivAttr(cd);
      //cda.planePos = divPointPos;
      cda.div2d = 1;
      cda.cellNormal = divPointNrml;
      cda.cellLabel = selectedLabel;

      if(testActual){
        CellDivision actualDiv = possibleDivisionsAct[0];
        if(actualDiv.divPlane.area != 0){
          cda.planeAreaRelActual = cda.planeArea / actualDiv.divPlane.area;
        } else {
          cda.planeAreaRelActual = -1;
        }
        
        actual.planeArea = actualDiv.divPlane.area;
        actual.daughterRatio = actualDiv.volumeRatio;
        cda.angleActual = angleVectors(cda.planeNrml, actual.planeNrml, false);
        actual.angleActual = 0;
        actual.planeAreaRelActual = 1;
      } else {
        cda.planeAreaRelActual = -1;
        cda.angleActual = -1;
      }

      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;

      if(testActual){
        actual.planeAreaRelShortest = actual.planeArea / minArea;
      }
    }

    return possibleDivisionsCDA;
  }

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

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

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

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

    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;

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

    selLabels.erase(-1);

    //if(selLabels.size() != 1 and !m->useParents()) return setErrorMessage("Please select exactly one cell");

    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;


    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);
        }
      }
    }



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

    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 " << cellCentroid << "/nrml " << divPointNrml << std::endl;
      divPointPos = cellCentroid;
    }



    FlatDivisionPlane2D *pFlat2D;
    if(!getProcess("Mesh/Division Analysis/Analysis 2D/Division Plane Approximation", pFlat2D))
      throw(QString("Division Analysis:: Unable to make Approximation process"));
    //if(displayActual) pFlat2D->parm("Display Plane") = "Yes";
    pFlat2D->parm("Display Plane") = "No";
    pFlat2D->parm("Reset Mesh") = "No";

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

    if(!actualOK) std::cout << "Couldn't compute actual plane for parent label: " << selectedParent << 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!");
    }
    bool testActual = actualOK;

    std::vector<CellDivisionAttr> possibleDivisions;

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

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


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

    return true;
  }
  REGISTER_PROCESS(TestDivisionPlanes2D);

 std::set<vertex> getVerticesOfCell(Mesh *m, int label, bool parent)
 {
   std::set<vertex> result, resultExtended;

   vvGraph& S = m->graph();

   forall(vertex v, S){
    if(v->label == label and !parent){
      result.insert(v);
    } else if(m->parents()[v->label] == label and parent){
      result.insert(v);
    } 
   }

  // add the border
   forall(vertex v, result){
     resultExtended.insert(v);
     forall(vertex n, S.neighbors(v)){
       resultExtended.insert(n);
     }
   }

   return resultExtended;
 }

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

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

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

    //AttrMap<int, CellDivisionAttr>& planeDataShortest = m->attributes().attrMap<int, CellDivisionAttr>("Cell Division Shortest Planes");

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

    m->updateCentersNormals();

    // 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 << "multi2d " << 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;

        TestDivisionPlanes2D *pDiv2D;
        if(!getProcess("Mesh/Division Analysis/Analysis 2D/Division Analysis", pDiv2D))
          throw(QString("Division Analysis:: Unable to make Approximation process"));

        bool res = pDiv2D->run(m, m2, pDiv2D->parm("Nr of Planes").toInt(), pDiv2D->parm("Steps Center Displacement").toInt(),
                    pDiv2D->parm("Stepsize Center Displacement").toDouble(), pDiv2D->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;

        TestDivisionPlanes2D *pDiv2D;
        if(!getProcess("Mesh/Division Analysis/Analysis 2D/Division Analysis", pDiv2D))
          throw(QString("Division Analysis:: Unable to make Approximation process"));

        bool res = pDiv2D->run(m, m2, pDiv2D->parm("Nr of Planes").toInt(), pDiv2D->parm("Steps Center Displacement").toInt(),
                    pDiv2D->parm("Stepsize Center Displacement").toDouble(), pDiv2D->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(DivisionAnalysisMulti2D);



}