//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
//

#include <BezierProcess.hpp>
#include <Dir.hpp>
#include <QFileDialog>

#include <Triangulate.hpp>

namespace mgx
{
  bool LoadBezier::initialize(QWidget* parent)
  {
    QString filename = parm("Input File");
    if(filename.isEmpty())
      filename = QFileDialog::getOpenFileName(parent, "Choose Bezier to load", QDir::currentPath(),
                                              "Bezier file (*.bez);;All files (*.*)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".bez", Qt::CaseInsensitive))
      filename += ".bez";
    setParm("Input File",stripCurrentDir(filename));
    return true;
  }

  bool LoadBezier::run(const QString& filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }

    CuttingSurface* cs = cuttingSurface();
    cs->setMode(CuttingSurface::BEZIER);
    cs->bezier().loadBezier(filename);
    cs->showGrid();

    SETSTATUS("Loaded Bezier file: " << filename);
    return true;
  }
  REGISTER_PROCESS(LoadBezier);

  bool SaveBezier::initialize(QWidget* parent)
  {
    QString filename = parm("Output File");
    if(filename.isEmpty())
      filename = QFileDialog::getSaveFileName(parent, "Choose Bezier file to save", filename, "Bezier files (*.bez)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".bez", Qt::CaseInsensitive))
      filename += ".bez";
    setParm("Output File", stripCurrentDir(filename));
    return true;
  }

  bool SaveBezier::run(const QString& filename)
  {
    QFile file(filename);
    if(filename.isEmpty()) {
      setErrorMessage(QString("Unable to write to output file"));
      return false;
    }
    file.setFileName(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("Unable to write to output file"));
      return false;
    }

    CuttingSurface* cs = cuttingSurface();
    cs->bezier().saveBezier(filename);

    SETSTATUS("Saved Bezier file: " << filename);
    return true;
  }
  REGISTER_PROCESS(SaveBezier);

  bool NewBezier::run(uint bezPointsX, uint bezPointsY, float sizeX, float sizeY, uint linesX, uint linesY)
  {
    CuttingSurface* cs = cuttingSurface();// NB: this does not call the default constructor, just returns the pointer to existing cutting surface.
    cs->setMode(CuttingSurface::BEZIER);
    cs->bezier().setBezPoints(Point2u(bezPointsX, bezPointsY));
    cs->bezier().setBezSize(Point2d(sizeX, sizeY));
    cs->bezier().setBezLines(Point2u(linesX, linesY));
    cs->bezier().initBez();
    cs->showGrid();

    return true;
  }
  REGISTER_PROCESS(NewBezier);


  bool MeasureBezier::run(Mesh* m, int discretization)
  {

  typedef std::map<Point2i, Point3d> P2iP3dMap;

    // get the bezier grid
    CuttingSurface* cutSurf = cuttingSurface();
    P2iP3dMap bezGridMap, diffBezGridMapX, diffBezGridMapY;
    double bezStart = 0.0;
    double bezEnd = 1.0;
    int dataPointsBezier = discretization; // number of intermediate points
    double stepSize = (double)(bezEnd-bezStart)/(double)(dataPointsBezier-1);
    // correct directions, take into account the rotations/translations
    Matrix4d rotMatrixS1, rotMatrixCS;
    m->stack()->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    std::vector<Point3d> bMap;

    // create bezier map
    for(int i = 0; i < dataPointsBezier; i++){
      double u = bezStart + i*stepSize;

      Point3d pb(cutSurf->evalCoord(u,0)); // get bez world coordinate
      pb = multMatrix4Point3(mGLTot,pb); // correct rotations

      bMap.push_back(pb);

      for(int j = 0; j < dataPointsBezier; j++){
        double v = bezStart+j*stepSize;
        Point3d p(cutSurf->evalCoord(v,u)); // get bez world coordinate
        p = multMatrix4Point3(mGLTot,p); // correct rotations
        Point2i a(j,i);
        bezGridMap[a] = p;
      }
    }

    double area = 0;

    // create derivatives in X and Y direction of the bezier
    for(int i=0; i<dataPointsBezier-1; i++){
      for(int j=0; j<dataPointsBezier-1; j++){

        int ip = i;
        int jp = j;

        if(i == 0) ip++;
        if(j == 0) jp++;

        Point2i a (jp,ip);
        Point2i a1 (jp+1,ip);
        Point2i a2 (jp,ip+1);
        Point2i a3 (jp+1,ip+1);

        area += triangleArea(bezGridMap[a], bezGridMap[a1], bezGridMap[a3]) +
                triangleArea(bezGridMap[a], bezGridMap[a2], bezGridMap[a3]);
      }
    }

    double lineLength = 0;

    Point3d startP = bMap[0];

    // calc length of line
    forall(auto p, bMap){
      lineLength += norm(p - startP);
      startP = p;
    }

    std::cout << "Bezier Grid Area: " << area << " / Bezier Line Length: " << lineLength  << "/ Start-End Distance: " << norm(startP-bMap[0]) << std::endl;


    uint bezPoints = cutSurf->bezier().bezPoints().x();
    std::vector<Point3d>& bezVec = cutSurf->bezier().bezierV();

    if(bezPoints < 2) return true;

    uint counter = 0;
    double dis = 0;
    Point3d pLast = bezVec[0];
    for(uint i=1; i<bezPoints; i++){
      dis += norm(pLast - bezVec[i]);
      pLast = bezVec[i];
    }
    double dis2 = dis + norm(pLast - bezVec[0]);

    std::cout << "Distance through points: " << dis << "/ Distance points as circle:" << dis2 << std::endl;

    return true;
  }
  REGISTER_PROCESS(MeasureBezier);

  // file dialogue for save data
  bool MeasureBezierExport::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getSaveFileName(0, "Choose 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 MeasureBezierExport::run(Mesh* m, int discretization, QString filename)
  {

  typedef std::map<Point2i, Point3d> P2iP3dMap;

    // get the bezier grid
    CuttingSurface* cutSurf = cuttingSurface();
    P2iP3dMap bezGridMap, diffBezGridMapX, diffBezGridMapY;
    double bezStart = 0.0;
    double bezEnd = 1.0;
    int dataPointsBezier = discretization; // number of intermediate points
    double stepSize = (double)(bezEnd-bezStart)/(double)(dataPointsBezier-1);
    // correct directions, take into account the rotations/translations
    Matrix4d rotMatrixS1, rotMatrixCS;
    m->stack()->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    std::vector<Point3d> bMap;

    // create bezier map
    for(int i = 0; i < dataPointsBezier; i++){
      double u = bezStart + i*stepSize;

      Point3d pb(cutSurf->evalCoord(u,0)); // get bez world coordinate
      pb = multMatrix4Point3(mGLTot,pb); // correct rotations

      bMap.push_back(pb);

      for(int j = 0; j < dataPointsBezier; j++){
        double v = bezStart+j*stepSize;
        Point3d p(cutSurf->evalCoord(v,u)); // get bez world coordinate
        p = multMatrix4Point3(mGLTot,p); // correct rotations
        Point2i a(j,i);
        bezGridMap[a] = p;
      }
    }


    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 << "Bez Point,X,Y,Z,Distance Points" << endl;


    uint bezPoints = cutSurf->bezier().bezPoints().x();
    std::vector<Point3d>& bezVec = cutSurf->bezier().bezierV();

    if(bezPoints < 2) return true;

    uint counter = 1;
    double dis = 0;
    Point3d pLast = bezVec[0];
    out << counter << "," << pLast.x() << "," << pLast.y() << "," << pLast.z() << "," << endl;
    for(uint i=1; i<bezPoints; i++){
      dis = norm(pLast - bezVec[i]);
      pLast = bezVec[i];
      out << ++counter << "," << pLast.x() << "," << pLast.y() << "," << pLast.z() << "," << dis << endl;
    }

    return true;
  }
  REGISTER_PROCESS(MeasureBezierExport);


  // Collapse Bezier points into a line
  bool CollapseBezier::run()
  {
    // get current Bezier points
    CuttingSurface* cutSurf = cuttingSurface();
    uint bezPoints = cutSurf->bezier().bezPoints().x();
    std::vector<Point3d>& bezVec = cutSurf->bezier().bezierV();

    // collapse them into a line
    uint counter = 0;
    forall(Point3d &bezPoint, bezVec){
      bezPoint = bezVec[counter];
      counter++;
      if(counter >= bezPoints)
        counter = 0;
    }

    return true;
  }

  REGISTER_PROCESS(CollapseBezier);

  bool NewBezierRing::run(uint bezPointsX, float sizeX, float sizeY)
  {
    if(bezPointsX < 4) return setErrorMessage("Need at least 4 points.");

    CuttingSurface* cs = cuttingSurface();// NB: this does not call the default constructor, just returns the pointer to existing cutting surface.
    cs->setMode(CuttingSurface::BEZIER);
    int lengthY = 2;
    cs->bezier().setBezPoints(Point2u(bezPointsX, lengthY));
    cs->bezier().setBezSize(Point2d(sizeX, sizeY));
    cs->bezier().setBezLines(Point2u(15, 15));
    cs->bezier().initBez();

    CollapseBezier* cb;
    if(!getProcess("Misc/Bezier/Collapse Bezier Points", cb))
      throw(QString("BezierProcess:: Unable to make Collapse Bezier Points process"));
    cb->run();

    std::vector<Point3d>& bezVec = cs->bezier().bezierV();

    // collapse them into a line
    double segm = 2.*M_PI/(double)(bezPointsX-1);
    int c = 0;
    for(int i = 0; i < bezVec.size(); i++){
      bezVec[i].x() = sizeX * cos(c*segm);
      bezVec[i].y() = sizeY * sin(c*segm);
      c++;
      if(c>=bezPointsX) c=0;
    }

    cs->showGrid();

    return true;
  }
  REGISTER_PROCESS(NewBezierRing);

  bool NewBezierRingFromSelected::run(Mesh* m, uint bezPointsX, float sizeX, float sizeY, QString mode)
  {
    if(bezPointsX < 4) return setErrorMessage("Need at least 4 points.");

    // find selected cells
    std::vector<Point3d> selectedCellCenters;

    vvGraph& S = m->graph();
    std::set<int> selectedCells = findAllSelectedLabels(S);

    if(mode == "2D"){

      m->updateCentersNormals();
      auto &centers = m->labelCenterVis();

      forall(int l, selectedCells){
        if(centers.find(l) == centers.end()){
          std::cout << "Can't find center of cell with label " << l << std::endl;
          continue;
        }
        selectedCellCenters.push_back(Point3d(centers[l]));
      }

    } else { // mode 3D

      AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
      if(centroids.empty()) return setErrorMessage("No cell centroids found. Run Cell Analysis 3D process first.");


      forall(int l, selectedCells){
        if(centroids.find(l) == centroids.end()){
          std::cout << "Can't find centroid of cell with label " << l << std::endl;
          continue;
        }
        selectedCellCenters.push_back(centroids[l]);
      }

    }

    if(selectedCellCenters.empty()) return setErrorMessage("No cell centers found.");

    // do a pca for plane fitting
    Point3d planePos, planeNrml;
    findPolygonPlane(selectedCellCenters, planePos, planeNrml);

    // rotate cutting surface into the plane
    // Rotate polygon plane to align it with (x,y) plane
    Point3d z_axis(0, 0, 1);

    // get the rotation matrix
    Matrix3d inv_rot_polygon_plane = calcRotMatrix(z_axis, -planeNrml);

    Point3d avgPos(0,0,0);

    forall(Point3d p, selectedCellCenters){
      avgPos += p;
    }
    avgPos /= selectedCellCenters.size();



    CuttingSurface* cs = cuttingSurface();// NB: this does not call the default constructor, just returns the pointer to existing cutting surface.
    cs->setMode(CuttingSurface::BEZIER);
    int lengthY = 2;
    cs->bezier().setBezPoints(Point2u(bezPointsX, lengthY));
    cs->bezier().setBezSize(Point2d(sizeX, sizeY));
    cs->bezier().setBezLines(Point2u(15, 15));
    cs->bezier().initBez();

    CollapseBezier* cb;
    if(!getProcess("Misc/Bezier/Collapse Bezier Points", cb))
      throw(QString("BezierProcess:: Unable to make Collapse Bezier Points process"));
    cb->run();

    std::vector<Point3d>& bezVec = cs->bezier().bezierV();

    // collapse them into a line
    double segm = 2.*M_PI/(double)(bezPointsX-1);
    int c = 0;
    for(int i = 0; i < bezVec.size(); i++){
      bezVec[i].x() = sizeX * cos(c*segm);
      bezVec[i].y() = sizeY * sin(c*segm);

      bezVec[i] = bezVec[i] * inverse(inv_rot_polygon_plane);
      bezVec[i] += planePos;

      c++;
      if(c>=bezPointsX) c=0;
    }

    cs->showGrid();

    return true;
  }
  REGISTER_PROCESS(NewBezierRingFromSelected);


    bool FlipBezierLine::run(Mesh* m)
    {
      // get current Bezier points
      CuttingSurface* cutSurf = cuttingSurface();
      std::vector<Point3d>& bezVec = cutSurf->bezier().bezierV();
      std::vector<Point3d> bezVecNew(bezVec.size());

      for(int i=0; i<bezVec.size(); i++){
        bezVecNew[bezVec.size()-i-1] = bezVec[i];
      }
      bezVec = bezVecNew;


      return true;
    }
    REGISTER_PROCESS(FlipBezierLine);

}
