//
// 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.
//
#ifndef MESH_PROCESSS_PDG_HPP
#define MESH_PROCESSS_PDG_HPP

#include <Process.hpp>

#include <MeshProcessCellMesh.hpp>

namespace mgx
{
  mgxBase_EXPORT SymmetricTensor customDirectionTensor(SymmetricTensor origTensor, Point3f customX, Point3f customY, bool shearB);
  mgxBase_EXPORT double lengthCustomDirTensor3D(SymmetricTensor origTensor, Point3f customDir);

  ///\addtogroup MeshProcess
  ///@{
  /**
   * \class CorrespondenceJunctions ProcessPDG.hpp <MeshProcessPDG.hpp>
   */
  class mgxBase_EXPORT CorrespondenceJunctions : public Process
  {
  public:
    CorrespondenceJunctions(const Process& process) : Process(process)
    {
	  setName("Mesh/Cell Axis/PDG/Check Correspondence");
	  setDesc("Find matching cell junctions between 2 meshes based on parent labeling.\n"
	  "Both meshes are simplified with Make Cells to keep only the cell junctions and centers.\n"
	  "SAVE your meshes before running!");
	  setIcon(QIcon(":/images/CorrespondenceJunctions.png"));

	  addParm("Show correspondence","Draw connecting lines between corresponding junctions.","No",booleanChoice());
	  addParm("Convert to Cell Mesh","Convert the meshes into cell meshes","No",booleanChoice());
	  addParm("New","New","No",booleanChoice());
	}

    bool run()
    {
      if(!checkState().mesh(MESH_NON_EMPTY, 0).mesh(MESH_NON_EMPTY, 1))
        return false;
      Mesh* mesh1, *mesh2;
      if(currentMesh() == mesh(0)) {
        mesh1 = mesh(0);
        mesh2 = mesh(1);
      } else if(currentMesh() == mesh(1)) {
        mesh2 = mesh(0);
        mesh1 = mesh(1);
      } else
        return false;

      bool res = run(mesh1, mesh2, stringToBool(parm("Show correspondence")), stringToBool(parm("Convert to Cell Mesh")), stringToBool(parm("New")));
      return res;
    }

    bool run(Mesh* mesh1, Mesh* mesh2, bool ShowVVCorrespondence, bool convert, bool newMethod = false);  //Added convert parameter to enable choice of converting into a cell mesh
    bool checkNhbd(Mesh* mesh1, Mesh* mesh2, vvGraph &T1, vvGraph &T2, bool ShowVVCorrespondence);  //Original Correspondence Junctions

    bool checkNhbdNew(Mesh* mesh1, Mesh* mesh2, vvGraph &T1, vvGraph &T2, bool ShowVVCorrespondence);

  };

typedef std::unordered_map<vertex, std::set<int> > JunctionCellsMap;
typedef std::unordered_map<int, std::set<vertex> > CellJunctionsMap;
typedef std::unordered_map<int, std::vector<vertex> > CellJunctionsInOrderMap;

typedef std::pair<vertex, std::set<int> > VertLabelsP;
typedef std::pair<int, std::set<vertex> > CellJunctionsP;
typedef std::pair<int, std::vector<vertex> > CellJunctionsInOrderP;

  // mesh2 = later timepoint with parent labels
  // mesh1 = earlier timepoint without parents
  struct JunctionInformation{
    JunctionCellsMap junctionLabels1, junctionLabels2; // map from (junction) vertex to set of neighbor labels
    CellJunctionsMap cellJunctions1, cellJunctions2; // map from cell label to vertices of its junctions
    CellJunctionsInOrderMap cellJunctionInOrder1;
    std::unordered_map<vertex, vertex> J1J2, J2J1;

    std::map<int, Point3d> labelParentCenterMap;
    std::map<int, Point3d> labelDaughtersCenterMap;
  };

  mgxBase_EXPORT void findCellJunctionsInOrder(Mesh* mesh, vvGraph& S, JunctionInformation& jInfo);
  mgxBase_EXPORT void findCellJunctions(Mesh* mesh, vvGraph& S, JunctionInformation& jInfo);
  mgxBase_EXPORT void findCellJunctionsParents(Mesh* mesh, vvGraph& S, JunctionInformation& jInfo);

  /**
   * \class GrowthDirections ProcessPDG.hpp <MeshProcessPDG.hpp>
   */
  class mgxBase_EXPORT GrowthDirections : public Process
  {
  public:
    GrowthDirections(const Process& process) : Process(process)
    {
	  setName("Mesh/Cell Axis/PDG/Compute Growth Directions");
      setDesc("Compute PDGs based on correspondence between junctions.\n"
      "Adapted from Goodall and Green, 'Quantitative Analysis of Surface Growth.'Botanical Gazette (1986) \n"
      "and Dumais and Kwiatkowska, 'Analysis of surface growth in shoot apices.'Plant Journal (2002)");
	  setIcon(QIcon(":/images/PDG.png"));
		}
    bool run()
    {
      if(!checkState().mesh(MESH_NON_EMPTY, 0).mesh(MESH_NON_EMPTY, 1))
        return false;

      Mesh* mesh1, *mesh2;
      if(currentMesh() == mesh(0)) {
        mesh1 = mesh(0);
        mesh2 = mesh(1);
      } else if(currentMesh() == mesh(1)) {
        mesh2 = mesh(0);
        mesh1 = mesh(1);
      } else
        return false;

      bool res = run(mesh1, mesh2);
      return res;
    }

    bool run(Mesh* mesh1, Mesh* mesh2);
  };

  /**
   * \class DisplayPDGs ProcessPDG.hpp <MeshProcessPDG.hpp>
   */
  class mgxBase_EXPORT DisplayPDGs : public Process
  {
  public:
    enum ParmNames { pDisplayHeatMap,	pScaleHeatMap, pRangeHeatLow, pRangeHeatHigh, pDisplayAxis,
		  pColorPos, pColorNeg, pAxisWidth, pAxisScale, pAxisOffset, pAnisotropyThreshold, pCustomShear, pNumParms};

    DisplayPDGs(const Process& process) : Process(process)
    {
	  setName("Mesh/Cell Axis/PDG/Display Growth Directions");
      setDesc("Display the principle growth directions on a vertex level");
      setIcon(QIcon(":/images/PDG.png"));

	  addParm("Heatmap","Display stretch ratio values in max or min direction as a color map.\n"
	  "stretch ratio = (length after/length before). No deformation means stretch ratio = 1.","StretchMax", QStringList() << "None" << "StretchMax" << "StretchMin" << "Aniso=StretchMax/StretchMin" << "StretchMax-StretchMin" << "ProductStretches" << "StretchCustomX" << "StretchCustomY" << "StretchCustomX-StretchCustomY" << "CustomX Growth Proportion" << "CustomY Growth Proportion" << "AnisoCustom" << "Shear" << "AnisoRatio");	// 0
	  addParm("ScaleHeat","Scale heat map","Auto", QStringList() << "None" << "Auto" << "Manual");
	  addParm("Heat min","High bound heat map","1");
	  addParm("Heat max","Low bound heat map","3");
	  addParm("Show Axis","Draw directions as vectors, scaled by strain.\n"
	  "strain = (stretch ratio - 1). No deformation means strain = 0.","Both", QStringList() << "Both" << "StrainMax" << "StretchMax-StretchMin" << "StrainMin" << "BothCustom" << "StrainCustomX" << "StrainCustomY" << "None");	// 4
	  addParm("Color +","Color used for expansion (strain > 0)","white", QColor::colorNames());
	  addParm("Color -","Color used for shrinkage (strain < 0)","red", QColor::colorNames());
	  addParm("Line Width","Line Width","2.0");
	  addParm("Line Scale","Length of the vectors = Scale * Strain.","2.0");
	  addParm("Line Offset","Draw the vector ends a bit tilted up for proper display on surfaces.","0.1");
	  addParm("Threshold","Minimal value of anisotropy (= stretchMax/StretchMin) required for drawing PDGs. \n"
	  "Use a value above 1.","0.0");
	  addParm("Custom Dir Shear","Custom Direction Shear","No",booleanChoice());
	  addParm("Min Dis Vtx","Min Dis Vtx","1.0");
	}

	 bool initialize(QWidget* parent);

		// Process the parameters
    void processParms();

    bool run()
    {
      if(!checkState().mesh(MESH_NON_EMPTY))
        return false;
      return run(currentMesh());
    }

    bool run(Mesh* mesh)
    {
      processParms();
      return run(mesh, DisplayHeatMap, ScaleHeatMap, RangeHeat, DisplayAxis,
			  ColorPos, ColorNeg, AxisWidth, AxisScale, AxisOffset, AnisotropyThreshold, CustomShear);
    }


    bool run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
            const QString displayAxis, const QColor& colorPos, const QColor& colorNeg, float axisWidth, float axisScale,
            float axisOffset, float anisotropyThreshold, bool customShear);




  private:
		QString DisplayHeatMap;
		QString ScaleHeatMap;
		Point2d RangeHeat;
		QString DisplayAxis;
		QColor ColorPos;
		QColor ColorNeg;
		float AxisWidth;
		float AxisScale;
		float AxisOffset;
		float AnisotropyThreshold;
		bool CustomShear;
  };
  ///@}


  /**
   * \class DisplayPDGs ProcessPDG.hpp <MeshProcessPDG.hpp>
   */
  class mgxBase_EXPORT DisplayPDGsVertex : public Process
  {
  public:
    enum ParmNames { pDisplayHeatMap, pScaleHeatMap, pRangeHeatLow, pRangeHeatHigh, pDisplayAxis,
      pColorPos, pColorNeg, pAxisWidth, pAxisScale, pAxisOffset, pAnisotropyThreshold, pCustomShear, pVtxMinDis, pNumParms};

    DisplayPDGsVertex(const Process& process) : Process(process)
    {
	  setName("Mesh/Cell Axis/Deformation Gradient/Display Vertex Gradient");
      setDesc("Display the principle growth directions on a vertex level");
	  setIcon(QIcon(":/images/PDG.png"));

      addParm("Heatmap","Display stretch ratio values in max or min direction as a color map.\n"
      "stretch ratio = (length after/length before). No deformation means stretch ratio = 1.","StretchMax", QStringList() << "None" << "StretchMax" << "StretchMin" << "Aniso" << "StretchMax-StretchMin" << "ProductStretches" << "StretchCustomX" << "StretchCustomY" << "StretchCustomX-StretchCustomY" << "CustomX Growth Proportion" << "CustomY Growth Proportion" << "AnisoCustom" << "Shear" << "AnisoRatio");	// 0
	  addParm("ScaleHeat","Scale heat map","Auto", QStringList() << "None" << "Auto" << "Manual");
	  addParm("Heat min","High bound heat map","1");
	  addParm("Heat max","Low bound heat map","3");
	  addParm("Show Axis","Draw directions as vectors, scaled by strain.\n"
	  "strain = (stretch ratio - 1). No deformation means strain = 0.","Both", QStringList() << "Both" << "StrainMax" << "StretchMax-StretchMin" << "StrainMin" << "BothCustom" << "StrainCustomX" << "StrainCustomY" << "None");
	  addParm("Color +","Color used for expansion (strain > 0)","white", QColor::colorNames());
	  addParm("Color -","Color used for shrinkage (strain < 0)","red", QColor::colorNames());
	  addParm("Line Width","Line Width","2.0");	// 7
	  addParm("Line Scale","Length of the vectors = Scale * Strain.","2.0");
	  addParm("Line Offset","Draw the vector ends a bit tilted up for proper display on surfaces.","0.1");
	  addParm("Threshold","Minimal value of anisotropy (= stretchMax/StretchMin) required for drawing PDGs.\n"
	  "Use a value above 1.","0.0");
	  addParm("Custom Dir Shear","Custom Direction Shear","No",booleanChoice());
	  addParm("Min Dis Vtx","Min Dis Vtx","1.0");
	}

    bool initialize(QWidget* parent);

    // Process the parameters
    void processParms();

    bool run()
    {
      if(!checkState().mesh(MESH_NON_EMPTY))
        return false;
      return run(currentMesh());
    }

    bool run(Mesh* mesh)
    {
      processParms();
      return run(mesh, DisplayHeatMap, ScaleHeatMap, RangeHeat, DisplayAxis,
        ColorPos, ColorNeg, AxisWidth, AxisScale, AxisOffset, AnisotropyThreshold, CustomShear, vtxMinDis);
    }

    bool run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
      const QString displayPDGs, const QColor& colorPos, const QColor& colorNeg, float axisLineWidth,
      float scaleAxisLength, float axisOffset, float anisotropyThreshold, bool customShear, double vtxMinDis);

  private:
    QString DisplayHeatMap;
    QString ScaleHeatMap;
    Point2d RangeHeat;
    QString DisplayAxis;
    QColor ColorPos;
    QColor ColorNeg;
    float AxisWidth;
    float AxisScale;
    float AxisOffset;
    float AnisotropyThreshold;
    bool CustomShear;
    double vtxMinDis;
  };


  class mgxBase_EXPORT AveragePDGs : public Process
  {
  public:
    AveragePDGs(const Process& process) : Process(process)
    {
    setName("Mesh/Cell Axis/PDG/Average Growth Directions");
    setDesc("Averages growth directions of neighboring cells");
    setIcon(QIcon(":/images/PDG.png"));

    addParm("Project Directions on Surface","Project Directions on Surface","Yes",booleanChoice());
    addParm("Weighting","Weighting","none",QStringList() << "none" << "n"<<"area"<<"geometric mean"<<"area weighted geometric mean");

    }

    bool run()
    {
      if(!checkState().mesh(MESH_NON_EMPTY))
        return false;
      return run(currentMesh(),stringToBool(parm("Project Directions on Surface")),parm("Weighting"));
    }
    bool run(Mesh *m,bool projSurf,QString avgType);

  };

}
#endif
