//
// 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 TRIANGULATE_HPP
#define TRIANGULATE_HPP

#include <Misc.hpp>
#include <Subdivide.hpp>
#include <Attributes.hpp>

#include <ImageData.hpp>

namespace mgx
{

  /**
   * triangleOrientationCheck
   **********************
   * check the orientation of all triangles and flip them if necessary, hereby the first one is considered as properly oriented
   * NOTE: this function is not 100% bug free!
   * std::vector<Point3i>& triList : vector of tris
   */
  mgx_EXPORT bool triangleOrientationCheck(std::vector<Point3i>& triList, int idxRefTri = -1);

  /**
   * orderPolygonSegs
   **********************
   * put segments (pairs of Point3d) into a vector in the correct order (starting with the first seg)
   * example: input ((2,3),(5,1),(3,5),(1,2)) output (2,3,5,1)
   * std::vector<std::pair<Point3d,Point3d> > &polygonSegs : vector of segments
   */
  mgx_EXPORT std::vector<Point3d> orderPolygonSegs(std::vector<std::pair<Point3d,Point3d> > &polygonSegs, bool keepPolySegs = true);
  mgx_EXPORT std::vector<std::vector<Point3d> > orderPolygonSegsMulti(std::vector<std::pair<Point3d,Point3d> > &polygonSegs);

  /**
   * getPolygonNrml
   **********************
   * returns the normal of a polygon defined by its points by finding three points of the polygon
   * that are not in a line and calculating their normal
   * std::vector<Point3d> &pointsPolygon : the input point vector of polygon points
   * int& idxP3 : by default the first two points are used, for the third point it is checked whether 
   * its on the line defined by the first two, this is the index of the first third point not on the line
   */
  mgx_EXPORT Point3d getPolygonNrml(std::vector<Point3d>& pointsPolygon);
  mgx_EXPORT Point3d getPolygonNrml(std::vector<Point3d>& pointsPolygon, int& idxP3);

  /**
   * isDataPlanar
   **********************
   * returns true if inputVertices are in one plane or false if they are true 3D points
   * std::vector<Point3d> &inputVertices : the input point vector of the test points
   */
  mgx_EXPORT bool isDataPlanar(std::vector<Point3d>& inputVertices);

  /**
   * findPolygonPlane
   **********************
   * gives a fitted plane of a vector of points (the corner points of a polygon), calculates the fit using a PCA
   * std::vector<Point3d> &pointsPolygon : the input point vector
   * Point3d planePos : position of the plane
   * Point3d planeNrml : normal of the plane
   */
  mgx_EXPORT bool findPolygonPlane(std::vector<Point3d> &pointsPolygon, Point3d &planePos, Point3d &planeNrml);

  /**
   * inplaceMatrixMultiply
   * P := P*M
   */
  mgx_EXPORT void inplaceMatrixMultiply(std::vector<Point3d>& P, const Matrix3d& M);

  /**
   * projectPointsOnPlane
   **********************
   * project a vector of points onto a plane
   * std::vector<Point3d> &points : the point vector, input AND output
   * Point3d planePos : position of the plane
   * Point3d planeNrml : normal of the plane
   */
  mgx_EXPORT bool projectPointsOnPlane(std::vector<Point3d> &points, Point3d planePos, Point3d planeNrml);

  /**
   * calcRotMatrix
   **********************
   * calculates the rotation matrix given an initial plane normal and a final plane normal
   * Point3d nrmlStart : initial plane normal
   * Point3d nrmlEnd : final plane normal
   */
  mgx_EXPORT Matrix3d calcRotMatrix(Point3d nrmlStart, Point3d nrmlEnd, Point3d rot_axis = Point3d(0,0,0));

  /**
   * rotatePointsIntoXY
   **********************
   * rotate (planar) (polygon) points into the XY plane
   **********************
   * NOTE: the new position will be parallel to the XY plane, but not necessarily in z=0 position.
   * points have to have their original z value, if it is required to rotate them back into the original position
   * **********************
   * Point3d polygonNormal : normal of the polygon plane
   * std::vector<Point3d>& pointsPolygon : the point vector, input AND output
   * Matrix3d& inv_rot_polygon_plane : inverse of the resulting rotation matrix (in case it has to be rotated back)
   */
  mgx_EXPORT void rotatePointsIntoXY(Point3d polygonNormal, std::vector<Point3d>& pointsPolygon, Matrix3d& inv_rot_polygon_plane);

  /**
   * triangulatePolygon3D
   **********************
   * takes a vector of polygon corner points in a plane and triangulates the polygon
   * the input polygon can either be arbitrary with its points in order or convex with unordered points
   * parms
   * IMPORTANT NOTE: this function is not necessarily bug-free and could cause infinite loops
   * float maxArea : maximum triangle area of resulting triangulation
   * Point3d polygonNormal : normal of the input polygon
   * std::vector<Point3d> &pointsPolygon : points of the input polygon
   * std::vector<Point3i> &triList : output triangle list
   * std::vector<Point3d> &ptList : output point list
   * bool boundary : if true also allow to subdivide the border segments that are longer than sqrt(maxArea*2)
   * bool inOrder : if true the input points are in order; if false points will be ordered automatically assuming convexity
   */
  mgx_EXPORT bool triangulatePolygon3D(float maxArea, Point3d polygonNormal, std::vector<Point3d> &pointsPolygon, 
    std::vector<Point3i> &triList, std::vector<Point3d> &ptList, bool boundary, bool inOrder, bool pointsCenter);

  /**
   * orthogonalVector
   **********************
   * returns a vector orthogonal to the input vector
   * Point3d vec : the input vector
   */
  mgx_EXPORT Point3d orthogonalVector(Point3d vec);

  /**
   * calculates the nearest point of a tragetPoint on a bezier grid (bezier)
   *
   * \param targetPoint Point3d input: the target point
   * \param bezier std::vector<std::vector<Point3d> > input: the Bezier Grid array
   * \param bezIdx int output: index of nearest point on the line
   *
   * \ingroup CellAtlasUtils
   */
   mgx_EXPORT Point3d calcNearestPointOnBezierGrid(const Point3d& targetPoint, std::vector<std::vector<Point3d> >& bezier, Point2i& bezIdx);

  /**
   * projectPointOnPlane
   **********************
   *
   */
 // mgx_EXPORT Point3d projectPointOnPlane(const Point3d& P, const Point3d& Q, const Point3d& N);



  mgx_EXPORT Matrix3d calcProjectedGrowthTensor(Matrix3d tensor, Point3d planeNrml);

  mgx_EXPORT bool pointInPolygon(Point3d pointToTest, std::vector<Point3d>& polygon, bool debug = false);

/**
   * k means function for 2D data (lloyds algorithm)
   *
   * \param clusters std::vector<Point2d> input clusters (could be random if none available)
   * \param data std::vector<Point2d> data points
   * \param maxStep int max number of update steps, if negative it runs until less than 1% of points change
   *
   * \ingroup CellAtlasUtils
   */
  mgx_EXPORT std::vector<Point2d> lloydKMeans2D(std::vector<Point2d>& clusters, std::vector<Point2d>& data, int maxStep = -1);

  /**
   * k means function for n-D data (lloyds algorithm)
   * NOTE: dimensions have to match, otherwise function returns false (size of each inputClusters vector == dataPoints.size())
   *
   * \param inputClusters std::vector<std::vector<double> > vector of input clusters (could be random if none available)
   * \param dataPoints std::vector<std::map<int, double> > vector of label -> data maps
   * \param maxStep int max number of update steps, if negative it runs until less than 1% of points change
   *
   * \ingroup CellAtlasUtils
   */
  mgx_EXPORT bool lloydKMeans(std::vector<std::vector<double> >& inputClusters,
    std::vector<AttrMap<int, double> >& dataPoints, int maxStep = -1, double convergenceThreshold = 0.01);

  // interpolates the appropriate array index for the landscape function which has a 
  // resolution of newRange+1 
  mgx_EXPORT double interpolateArrayIndex(double value, double min, double max, double newRange);
  mgx_EXPORT void drawCross(QImage& image, int x, int y, int size, int crossSize, QRgb color);
  mgx_EXPORT void drawCircle(QImage& image, int x, int y, int size, QRgb color);

  /**
   * returns the value of a 2D gaussian function
   *
   * \param x double input value x dimension
   * \param y double input value y dimension
   * \param muX double mean value x dim
   * \param muY double mean value y dim
   * \param sigmaX double standard deviation value x dim
   * \param sigmaY double standard deviation value y dim
   *
   * \ingroup CellAtlasUtils
   */
  mgx_EXPORT double gauss2D(double x, double y, double muX, double muY, double sigmaX, double sigmaY);


  /**
   * calculates the nearest point of a tragetPoint on a bezier line (bMap)
   * as the bezier is discretized, it is approximated linearly between the nearest point and the closer neighbor point
   *
   * \param targetPoint Point3d input: the target point
   * \param bMap std::map<int, Point3d> input: the Bezier line
   * \param bezIdx int output: index of nearest point on the line
   * \param bezWeight double output: weight of the point
   *
   * \ingroup CellAtlasUtils
   */
  mgx_EXPORT Point3d calcNearestPointOnBezierLine(const Point3d& targetPoint, std::vector<Point3d>& bezierVec, int& bezIdx, double& bezWeight);
  mgx_EXPORT Point3d calcNearestPointOnBezierLine(const Point3d& targetPoint, std::map<int, Point3d>& bMap, int& bezIdx, double& bezWeight);
  mgx_EXPORT Point3d calcNearestPointOnBezierLine(const Point3d& targetPoint, std::map<int, Point3d>& bMap);

  mgx_EXPORT QRgb calcRGB(double maxValue, double currentValue);

  mgx_EXPORT QRgb getColorFromLabel(int label);

  mgx_EXPORT double angleVectors(Point3d v1, Point3d v2, bool directed = false);

}
#endif

