#include "DivisionAnalysisProcesses.hpp"

using namespace std;

namespace mgx { //namespace process

   double maxPlaneSimilarity(AttrMap<int, CellDivisionAttr>& cellDivs, CellDivisionAttr cda)
 {

  double minAngle = 90;// = 0;
  //std::cout << "size " << cellDivs.size() << std::endl;
  forall(auto p, cellDivs){
    double sim = angleVectors(p.second.planeNrml, cda.planeNrml);//180/3.1415*acos((p.second.planeNrml * cda.planeNrml)/ norm(p.second.planeNrml) / norm(cda.planeNrml));
    //std::cout << "sim " << sim << "/" << p.second.planeNrml << "/" << cda.planeNrml << std::endl;
    if(sim < minAngle) minAngle = sim;
  }
  return minAngle;
 }


 bool FilterPlanes::run(Mesh *m, Mesh *m2, QString planes, QString filter, int maxPlanes, double maxSimilarity, bool drawPlanes, bool resetMesh2,
  double planeSize, bool actual, bool selectActual, QString heatmap)
  {

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

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

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

    if(planes == "All Planes") planeDataFiltered.clear();

    AttrMap<int, CellDivisionAttr>& searchData = planeData;
    if(planes == "Filtered Planes") searchData = planeDataFiltered;

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


    if(searchData.empty()) return setErrorMessage("No plane data found.");

    AttrMap<int, CellDivisionAttr> newFiltered;

    std::set<int> labels;
    std::map<int, std::vector<CellDivisionAttr> > labelAttrMap;
    forall(auto p, searchData){
      int label = p.second.cellLabel;
      labels.insert(label);
      labelAttrMap[label].push_back(p.second);
    }


    forall(int l, labels){
      std::cout << "labels " << l << "/" << labelAttrMap[l].size() << std::endl;
    }


    forall(auto lAttr, labelAttrMap){
      std::map<double, CellDivisionAttr> sortedMap;
      AttrMap<int, CellDivisionAttr> newFilteredLabel;

      forall(CellDivisionAttr cda, lAttr.second){

        double planeValue = std::rand(); // random by default

        if(filter == "Min Area" or filter == "Max Area"){
          planeValue = cda.planeArea;
        } else if(filter == "Min Angle to Actual" or filter == "Max Angle to Actual"){
          planeValue = cda.angleActual;
        } else if(filter == "Min Volume Ratio" or filter == "Max Volume Ratio"){
          planeValue = cda.daughterRatio;
        }

        sortedMap[planeValue] = cda;

      }

      int counter = 0;
      auto itStart = sortedMap.begin();
      auto itEnd = sortedMap.end();
      if(filter == "Min Area" or filter == "Min Angle to Actual" or filter == "Min Volume Ratio"){

        for(; itStart!=itEnd; itStart++){
          if(newFilteredLabel.size() >= (uint)maxPlanes) continue;
          CellDivisionAttr newCDA = itStart->second;
          double minAngle = maxPlaneSimilarity(newFilteredLabel, newCDA);
          if(minAngle > maxSimilarity){
            //std::cout << "added plane " << minAngle << "/" << maxSimilarity << std::endl;
            newFilteredLabel[counter++] = newCDA;

          }

        }

      } else {
        itEnd--;
        for(; itStart!=itEnd; itEnd--){
          if(newFilteredLabel.size() >= (uint)maxPlanes) continue;
          CellDivisionAttr newCDA = itEnd->second;
          double minAngle = maxPlaneSimilarity(newFilteredLabel, newCDA);
          if(minAngle > maxSimilarity){
            //std::cout << "added plane " << minAngle << "/" << maxSimilarity << std::endl;
            newFilteredLabel[counter++] = newCDA;
          }
        }
      }

      forall(auto p, newFilteredLabel){
        std:cout << "Filtered " << p.first << "/" << p.second.planeArea << "/" << p.second.angleActual << "/" << p.second.daughterRatio << "/" << p.second.planeNrml << "/" << p.second.planePos << std::endl;
      }

      if(drawPlanes){

        int maxLabel = 0;
        forall(vertex v, m->graph()){
          if(maxLabel < v->label) maxLabel = v->label;
        }
        forall(vertex v, m2->graph()){
          if(maxLabel < v->label) maxLabel = v->label;
        }

        forall(auto p, newFilteredLabel){
          maxLabel++;
          DivisionPlane dp = DivisionPlane(planeSize, p.second.planeNrml, p.second.planePos, maxLabel);

          double planeValue = 0;

          if(heatmap == "Plane Area Abs"){
            planeValue = p.second.planeArea;
          } else if(heatmap == "Plane Area Rel Shortest"){
            planeValue = p.second.planeAreaRelShortest;
          } else if(heatmap == "Plane Area Rel Actual"){
            planeValue = p.second.planeAreaRelActual;
          } else if(heatmap == "Angle to Actual"){
            planeValue = p.second.angleActual;
          } else if(heatmap == "Daughter Ratio"){
            planeValue = p.second.daughterRatio;
          } else if(heatmap == "Distance Centroid"){
            planeValue = p.second.distanceCentroid;
          } else if(heatmap == "Distance Actual"){
            planeValue = p.second.distanceActual;
          } else if(heatmap == "Normalized Area"){
            planeValue = p.second.normalizedArea;
          }
          m2->labelHeat()[maxLabel] = planeValue;

          if(p.second.div2d != 1) drawPlane(*this, m2, dp);
          else drawPlane(*this, m2, dp, p.second.cellNormal, 0.95); // add some scaling for 2D planes
        }

        m2->heatMapBounds() = m2->calcHeatMapBounds();
        m2->setShowLabel("Label Heat");
        m2->updateAll();

      }

      counter = planeDataFiltered.size();
      forall(auto p, newFilteredLabel){
        planeDataFiltered[counter] = p.second;
        counter++;
      }
    }

    if(actual){
      forall(auto p, actualPlane){
        DivisionPlane dp = DivisionPlane(planeSize, p.second.planeNrml, p.second.planePos, p.first);
        if(p.second.div2d != 1) drawPlane(*this, m2, dp);
        else drawPlane(*this, m2, dp, p.second.cellNormal, parm("Plane Scaling 2D").toDouble()); // add some scaling for 2D planes
      
        double planeValue = 0;

        if(heatmap == "Plane Area Abs"){
          planeValue = p.second.planeArea;
        } else if(heatmap == "Plane Area Rel Shortest"){
          planeValue = p.second.planeAreaRelShortest;
        } else if(heatmap == "Plane Area Rel Actual"){
          planeValue = p.second.planeAreaRelActual;
        } else if(heatmap == "Angle to Actual"){
          planeValue = p.second.angleActual;
        } else if(heatmap == "Daughter Ratio"){
          planeValue = p.second.daughterRatio;
        } else if(heatmap == "Distance Centroid"){
          planeValue = p.second.distanceCentroid;
        } else if(heatmap == "Distance Actual"){
          planeValue = p.second.distanceActual;
        } else if(heatmap == "Normalized Area"){
          planeValue = p.second.normalizedArea;
        }
        m2->labelHeat()[p.first] = planeValue;

      }

    }
    if(selectActual){
      std::set<int> allParents;
      forall(auto p, m->parents()){
        allParents.insert(p.second);
      }
      forall(vertex v, S2){
        v->selected = false;
        if(allParents.find(v->label) == allParents.end()) continue;
        v->selected = true;
      }
      forall(vertex v, S2){
        if(!v->selected) continue;
        forall(vertex n, S2.neighbors(v)){
          n->selected = true;
        }
      }
    }
    m2->updateAll();

    return true;
  }
  REGISTER_PROCESS(FilterPlanes);

  REGISTER_PROCESS(ResetDivPlaneData);


 bool ComputeDivPlaneAngles::run(Mesh *m, Mesh *m2, QString plane1, QString plane2)
  {

    if(plane1 == plane2) return setErrorMessage("Plane1 and Plane2 must be different!");

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

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

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

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

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

    std::map<int, Point3d> mapPlane1;

    std::map<int, std::vector<Point3d> > mapVecPlane1;

    if(plane1 == "Actual"){
      // forall(auto p, normalsActual){
      //   mapPlane1[p.first] = p.second;
      // }
      forall(auto p, actualPlane){
        mapPlane1[p.first] = p.second.planeNrml;
      }
    } else if(plane1 == "Shortest"){
      FilterPlanes *pFilter;
      if(!getProcess("Mesh/Division Analysis/Display and Filter Planes", pFilter))
        throw(QString("Division Analysis:: Unable to make Filter Planes process"));

      pFilter->run(m, m2, "All Planes", "Min Area", 1, 1., false, false, 1.,false, false, "None");

      forall(auto p, planeDataFiltered){
        int label = p.second.cellLabel;
        mapPlane1[label] = p.second.planeNrml;
      }
    } else if(plane1 == "Filtered Best" or plane1 == "Filtered Worst"){
      forall(auto p, planeDataFiltered){
        mapPlane1[p.second.cellLabel] = Point3d(0,0,0);
        mapVecPlane1[p.second.cellLabel].push_back(p.second.planeNrml);
      }
    }
    std::cout << "s " << mapPlane1.size() << std::endl;

    forall(auto p, mapPlane1){
      int label = p.first;
      Point3d nrml1 = p.second;
      Point3d nrml2;
      std::cout << "l " << p.first << "/" << p.second << "/" << (customDirs.find(label) != customDirs.end()) << std::endl;
      std::cout << "c " << customDirs[label].ev1() << std::endl;
      if(plane2 == "Custom Dir X" and customDirs.find(label) != customDirs.end()){
        nrml2 = Point3d(customDirs[label].ev1());
      } else if(plane2 == "Custom Dir Y" and customDirs.find(label) != customDirs.end()){
        nrml2 = Point3d(customDirs[label].ev2());
      } else if(plane2 == "Custom Dir Z" and customDirs.find(label) != customDirs.end()){
        nrml2 = Point3d(customDirs[label].evals());
      } else if(plane2 == "Actual" and actualPlane.find(label) != actualPlane.end()){
        nrml2 = Point3d(actualPlane[label].planeNrml);
      } else {
        continue;
      }

      if(plane1 == "Shortest" or plane1 == "Actual"){
        angle[label] = angleVectors(nrml1, nrml2);
      } else { // check filtered
        double minAngle = 91;
        double maxAngle = -1;
        forall(Point3d nrmlComp, mapVecPlane1[label]){
          double angleCurrent = angleVectors(nrmlComp, nrml2);
          if(angleCurrent < minAngle) minAngle = angleCurrent;
          if(angleCurrent > maxAngle) maxAngle = angleCurrent;
        }

        if(plane1 == "Filtered Best") angle[label] = minAngle;
        if(plane1 == "Filtered Worst") angle[label] = maxAngle;
      }


      std::cout << "o " << p.first << "/" << p.second << "/" << nrml2 << "/" << angle[label] << std::endl;
    }

    m->labelHeat().clear();
    forall(auto p, angle){
      m->labelHeat()[p.first] = p.second;
    }

    m->setShowLabel("Label Heat");
    m->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(ComputeDivPlaneAngles);


}
