//
// 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 <Progress.hpp>
#include <LobynessPlugin.hpp>
#include <MeshProcessCellMesh.hpp>
#include <Triangulate.hpp>
#include <triangle.h>
#include <GraphUtils.hpp>
#include <limits>
#include <fstream>

namespace mgx 
{

  bool LobynessPlugin::run(Mesh *m,Mesh *m2, const QString &measure, const QString &iteration_string,const QString &fix_mesh_string, const QString &save_LEC_string,const QString &LEC_CSV_string)
  {  
    //epsilon = epsilon_string.toDouble();
    smooth_its = iteration_string.toInt();
    fix_mesh = stringToBool(fix_mesh_string);
    save_LEC = stringToBool(save_LEC_string);
    LEC_CSV = LEC_CSV_string;
  
#ifdef DEBUG_LOBYNESS_PLUGIN
    debug_mesh = m2;
    //MakeMeshCircle ext(*this);
    //ext.run(debug_mesh, 0.01, 0.01, 100.0, 20, false, 0, true); //A hack to clear the debug_mesh...for now
    //debug_mesh->clear();
#endif
   if(fix_mesh && m->tissue().meshType() == "MGXM"){
	FixBorderTriangles fixB(*this);
	fixB.run(m);
  	FixMeshCorners fixC(*this);
	fixC.run(m,true,false,5);  //run fix corners with default parameters (minus selecting bad vertices)
   }

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


   if(m->tissue().meshType() == "MGXM")
   	mgxmToMgxc(S, T, 1);
   else
	  T = S;

   IntFloatAttr &heat_val = m->labelHeat();
   heat_val.clear();


   if(measure == "Convex(Perimeter)")
     std::cout<<"Calculating convexity (perimeter)"<<std::endl;
   if(measure == "Convex(Area)")
     std::cout<<"Calculating convexity (area)"<<std::endl;
   if(measure == "Visibility(Stomata)")
     std::cout<<"Calculating visibility (Stomata)"<<std::endl;
   if(measure == "Visibility(Pavement)")
     std::cout<<"Calculating visibility (Pavement)"<<std::endl;
   if(measure == "Circularity")
     std::cout<<"Calculating circularity"<<std::endl;
   if(measure == "LargestEmptySpace"){
     std::cout<<"Calculating largest empty space"<<std::endl;
     lec_data.clear();
   }
   if(measure == "Skeleticity")
   	 std::cout<<"Calculating skeleticity"<<std::endl;

  #pragma omp parallel for schedule(dynamic)
  for(uint i = 0; i < T.size(); i++){
    const vertex &v = T[i];
     if(v->label==-1)
       continue;
     std::vector<Point3d> npos;
     Point3d ppos;
     Point3d pnorm;
     Matrix3d prot;
     forall(const vertex &s, T.neighbors(v))
       npos.push_back(s->pos);

     findPolygonPlane(npos,ppos,pnorm);
     projectPointsOnPlane(npos,ppos,pnorm);
     rotatePointsIntoXY(pnorm,npos,prot);
#ifdef DEBUG_LOBYNESS_PLUGIN
     debug_rot = prot;
     debug_z = npos[0].z();
     debug_inv_x = false;
#endif 

     for(uint i=0;i<npos.size();i++)
       npos[i].z()=0;
     if(Area(npos)<0)
       for(uint i=0;i<npos.size();i++){
         npos[i].x()*=-1;
         debug_inv_x = true;
         //debug_rot*=-1;
       }
       //std::cout<<"Smoothing"<<std::endl;
       SmoothContour(npos,smooth_its);
       //std::cout<<"Calculating heat val";
       float cell_val=0;
       if(measure == "Convex(Perimeter)")
         cell_val = ConvexityPerimeter(npos);
       if(measure == "Convex(Area)")
         cell_val = ConvexityArea(npos);
       if(measure == "Visibility(Stomata)")
         cell_val = VisibilityIndex(ConnectivityMap(npos));
       if(measure == "Visibility(Pavement)")
         cell_val = 1.0-VisibilityIndex(ConnectivityMap(npos));
       if(measure == "LargestEmptySpace"){
	       if(save_LEC)
         	cell_val = LargestEmptySpace(npos,v->label);
	       else
         	cell_val = LargestEmptySpace(npos);			 
       }
       if(measure == "Circularity")
         cell_val = Circularity(npos);
       //std::cout<<"Storing heat val"<<std::endl;
       if(measure == "Skeleticity")
       	 cell_val = Skeleticity(npos);
        #pragma omp critical
       heat_val[v->label] = cell_val;
       //break;
     }


    m->heatMapBounds() = m->calcHeatMapBounds();
    m->updateTriangles();
    m->setShowLabel("Label Heat");

    if(save_LEC){
	std::fstream csvout(LEC_CSV.toLatin1().data(),std::fstream::out);
	csvout<<"Label"<<","<<"Radius"<<","<<"X position"<<","<<"Y position"<<","<<"Z position"<<std::endl;
	for(unsigned int i=0;i<lec_data.size();i++){
		LEC_output cur_cell = lec_data[i];
		csvout<<cur_cell.label<<","<<cur_cell.radius<<","<<cur_cell.centre.x()<<","<<cur_cell.centre.y()<<","<<cur_cell.centre.z()<<std::endl;

	}
	csvout.close();
    }
 
    return true;
}

void LobynessPlugin::MakeNeighMap(vvGraph &T, IntFloatAttr &heat_val, std::map<int,std::map<int, double> > &neigh_vals, Mesh *m)
{
  forall(const vertex &v, T){
    if(v->label==-1) //!m->uniqueTri(v, s, n))
      continue;
    neigh_vals[v->label]= std::map<int,double>();
    forall(const vertex &s, T.neighbors(v)) {
      const vertex &n = T.nextTo(s,T.nextTo(s,v)); 
      if(n->label==-1)
        continue;
      neigh_vals[v->label][n->label]=  heat_val[n->label];  
    }
  }
}

float LobynessPlugin::Circularity(std::vector<Point3d> contour){
  return pow(Perimeter(contour),2)/4.0/M_PI/(Area(contour));
}

float LobynessPlugin::ConvexityPerimeter(std::vector<Point3d> contour){
  return Perimeter(contour)/Perimeter(ConvexHull(contour));
}

float LobynessPlugin::ConvexityArea(std::vector<Point3d> contour){
  return Area(ConvexHull(contour))/Area(contour);
}

float LobynessPlugin::Skeleticity(std::vector<Point3d> contour){
	std::vector<Point3d> chull = ConvexHull(contour);
	float avgLen = Perimeter(contour)/(float)(contour.size());
	float avgCHLen = Perimeter(chull)/(float)(chull.size());
	std::cout<<avgCHLen<<std::endl;
	while(avgCHLen>avgLen){
		std::vector<Point3d> refined_cont;
		for(uint i=0;i<chull.size();i++){
			refined_cont.push_back(chull[i]);
			if(norm(chull[i]-chull[(i+1)%chull.size()])>avgLen){
				refined_cont.push_back( 0.5*(chull[i]+chull[(i+1)%chull.size()]) );
			}
		}
		chull = refined_cont;
		avgCHLen = Perimeter(chull)/(float)(chull.size());
		std::cout<<chull.size()<<":"<<avgLen<<":"<<avgCHLen<<std::endl;
	}

   return LargestEmptySpace(contour)/LargestEmptySpace(ConvexHull(contour));
}

float LobynessPlugin::SignedTriArea(Point3d p1,Point3d p2,Point3d p3)
{
  Point3d d21 = p2-p1;
  Point3d d32 = p3-p2;
  return (d21.x()*d32.y() - d21.y()*d32.x())/2.0;
}

float LobynessPlugin::Perimeter(std::vector<Point3d> contour)
{
  float length=0;
  for(uint i=1;i<contour.size();i++)
    length+=norm(contour[i]-contour[i-1]);
  if(contour.size())
    length+=norm(contour[0]-contour[contour.size()-1]);
  return length;
}


float LobynessPlugin::Area(std::vector<Point3d> contour)
{
  float area=0;
  for(uint i=1;i<contour.size();i++)
    area+=SignedTriArea(Point3d(0,0,0),contour[i-1],contour[i]);
  area+=SignedTriArea(Point3d(0,0,0),contour[contour.size()-1],contour[0]);
  return area;
}

std::vector<Point3d> LobynessPlugin::ConvexHull(std::vector<Point3d> contour)
{
  std::vector<Point3d> chull;
  bool convex = false;
  while(!convex && contour.size()>3){
  convex = true;
  for(uint i=0;i<contour.size();i++){
    int ind1 = (i+contour.size()-1)%contour.size();
    int ind2 = i;
    int ind3 = (i+contour.size()+1)%contour.size(); 
    if(SignedTriArea(contour[ind1],contour[ind2],contour[ind3])>0)
      chull.push_back(contour[i]);
    else
      convex = false;
  }
    contour = chull;
    chull.clear();
  }
  return contour;
}

float LobynessPlugin::IntPoint(Point3d p1,Point3d p2, Point3d q1,Point3d q2)
{
  Point3d d= p2-p1;
  Point3d n(-(q2.y()-q1.y()),q2.x()-q1.x(),0);  
  float t = -(p1 - q1)*n/(d*n); 
  return t;
}

bool LobynessPlugin::isVisible(int i,int j,const std::vector<Point3d> &contour)
{
  bool visible = true;
//  srand(10);
  for(uint k=0;k<contour.size();k++){
    int ind1 = k;
    int ind2 = (k+1)%contour.size();
    if(ind1 == i || ind2==i || ind1==j || ind2==j)
      continue;
    float t1 = IntPoint(contour[ind1],contour[ind2],contour[i],contour[j]);
    float t2 = IntPoint(contour[i],contour[j],contour[ind1],contour[ind2]);
    if(t1>0 && t1<1 && t2>0 && t2<1)
      visible = false;
    }

  //if(rand()%500)
  //  visible = false;

  return visible;

}

bool LobynessPlugin::ThroughInterior(int i,int j,const std::vector<Point3d> &contour)
{
  int il = (i-1+contour.size())%contour.size();
  int ir = (i+1)%contour.size();
  int jl = (j-1+contour.size())%contour.size();
  int jr = (j+1)%contour.size();
  if(SignedTriArea(contour[il],contour[i],contour[j])<0 && SignedTriArea(contour[j],contour[i],contour[ir])<0 )
    return false;
  if(SignedTriArea(contour[jl],contour[j],contour[i])<0 && SignedTriArea(contour[i],contour[j],contour[jr])<0 )
    return false;

  return true;
}

//Epsilon for the inclusion of short hops outside the contour
std::vector<std::vector<double> > LobynessPlugin::ConnectivityMap(std::vector<Point3d> contour)
{
  std::vector<std::vector<double> > map;
  map.resize(contour.size());
  for(uint i=0;i<contour.size();i++)
    map[i].resize(contour.size());

  for(uint i=0;i<contour.size();i++)
    for(uint j=0;j<contour.size();j++){
      if(ThroughInterior(i,j,contour) and isVisible(i,j,contour)){
        map[i][j] = 1;
      }
      else
        map[i][j] = 0;
    }
  return map;
}

void LobynessPlugin::SmoothContour(std::vector<Point3d> &contour,int its,float rate)
{
  std::vector<Point3d> ncontour;
  ncontour.resize(contour.size());
  for(;its>0;its--){
   for(uint i=0;i<contour.size();i++){
    int ind1 = (i+contour.size()-1)%contour.size();
    int ind2 = i;
    int ind3 = (i+1)%contour.size();
    ncontour[i]= rate*(contour[ind1]-2.0*contour[ind2]+contour[ind3])+contour[ind2];
  }
    contour = ncontour;
  }
}

float LobynessPlugin::VisibilityIndex(std::vector<std::vector<double> > map)
{
  float visible=0;
  float total=0;
  for(uint i=0;i<map.size();i++)
    for(uint j=0;j<map[i].size();j++){
      if(map[i][j]>0 || i==j)
        visible++;
      total++;
    }
  return visible/total; //there will be some small errors due to numerics and perhaps incorrect handling of neighbors and reflexive points
}

float LobynessPlugin::LargestEmptySpace(std::vector<Point3d> contour,int label)
{
  if(contour.size()<3)
    return 0;
  float emp_space = 0;
  std::vector<Point2f> pList;
  std::vector<Point2i> segList;
  std::vector<Point2f> ptList;
  std::vector<Point3i> triList;
  std::vector<Point2i> segtList;

  for(uint i=0;i<contour.size();i++){
    Point2f cpos(contour[i].x(),contour[i].y());
    Point2i seg_ind(i,(i+1)%contour.size());
    pList.push_back(cpos);
    segList.push_back(seg_ind);
  }
#ifdef DEBUG_LOBYNESS_PLUGIN
/*
//Maps 2D cells to their 3D position
  for(int i=0;i<contour.size();i++){
    vertex v;
    v->pos = contour[i];
    if(debug_inv_x)
           v->pos.x()*=-1;
    v->pos.z() = -debug_z;
    v->pos = inverse(debug_rot)*v->pos;
    v->pos.z()*=-1;
    debug_mesh->graph().insert(v);
  }
  forall(const vertex &v, debug_mesh->graph())
    if(!vtxExists[v])
      vtxExists[v] = true;
*/
#endif

  triangulatePolygon(10000, pList, segList, ptList, triList, segtList, false, false);
  Point3d debug_pos;
  for(uint i=0;i<triList.size();i++){
    Point3i tind = triList[i];
    Point2f p1 = ptList[tind.x()]; 
    Point2f p2 = ptList[tind.y()];
    Point2f p3 = ptList[tind.z()];
    Point2f a = p1-p2;
    Point2f b = p3-p2;
    Point2f c = p1-p3;
    float rad = norm(a)*norm(b)*norm(c);
    float denom = sqrt(fabs(a.x()*b.y() - a.y()*b.x()));
    denom*=denom*2.0;
    rad/=denom;
    //rad = sqrt(rad);
    if(rad>emp_space){
      denom = (a%b)*(a%b)*2.0;
      emp_space = rad;
      float a1 = (b*b)*(a*c)/denom;
      float a2 = (c*c)*(a*b)/denom;
      float a3 = -(a*a)*(c*b)/denom;

      debug_pos.x() = a1*p1.x()+a2*p2.x()+a3*p3.x();
      debug_pos.y() = a1*p1.y()+a2*p2.y()+a3*p3.y();
      debug_pos.z() = 0;

      //debug_pos.x() = 1.0/3.0*(p1.x()+p2.x()+p3.x());//a1*p1.x()+a2*p2.x()+a3*p3.x();
      //debug_pos.y() = 1.0/3.0*(p1.y()+p2.y()+p3.y());
      //debug_pos.z() = 0;


    }
  }

#ifdef DEBUG_LOBYNESS_PLUGIN
/*
   MakeMeshCircle ext(*this);
    //ext.run(debug_mesh, emp_space*0.01, emp_space*0.01, 100.0, 20, false, 0, false);
   ext.run(debug_mesh, emp_space, emp_space, 100.0, 20, false, 0, false);
      forall(const vertex &v, debug_mesh->graph()){
        if(!vtxExists[v]){
          vtxExists[v] = true;
          v->pos+=debug_pos;
          if(debug_inv_x)
      v->pos.x()*=-1;
      v->pos.z() = -debug_z;
      v->pos = inverse(debug_rot)*v->pos;
      v->pos.z()*=-1;
          //v->pos = inverse(debug_rot)*(v->pos + debug_pos);
          //v->pos.x()*=-1;
//          v->pos.z()*=-1;
          v->label = rand();
        }
      }
*/
      //emp_space = 100;
#endif

  if(save_LEC){
	LEC_output cur_cell;
	cur_cell.label = label;
	cur_cell.radius = emp_space;
	cur_cell.centre = debug_pos;
        if(debug_inv_x)
      		cur_cell.centre.x()*=-1;
      	cur_cell.centre.z() = -debug_z;
        cur_cell.centre = inverse(debug_rot)*cur_cell.centre;
        cur_cell.centre.z()*=-1;
	lec_data.push_back(cur_cell);
    vertex v;
    v->pos = cur_cell.centre;
    debug_mesh->graph().insert(v);

  }

  return emp_space;
}


/*float LobynessPlugin::ExteriorDistance(int i,int j,std::vector<Point3d> contour){
  float int_dist = 0;
  for(int k=0;k<contour.size();k++){
    int ind1 = k;
    int ind2 = (k+1)%contour.size();
    if(ind1 == i || ind2==i || ind1==j || ind2==j)
      continue;
    float t1 = IntPoint(contour[ind1],contour[ind2],contour[i],contour[j]);
    float t2 = IntPoint(contour[i],contour[j],contour[ind1],contour[ind2]);
    Point3d int_point = t1*(contour[ind2]-contour[ind1])+contour[ind1];
    if(t1>=0 && t1<1 && t2>=0 && t2<1){
      float sign = SignedTriArea(contour[ind1],contour[ind2],contour[i])>0 ? 1:-1;
      int_dist += norm(contour[i]-int_point)*sign;
    }
  }
  return int_dist;
}

//This function hasn't been fully tested
 void LobynessPlugin::NeighAverage(std::map<int, double> &avg_heat_val, std::map<int,std::map<int, double> > &neigh_vals, Mesh *m){
      forall(const vertex &v, T){
          if(v->label==-1)
                continue;    
          double cell_val=0;
          std::map<int,double> cvalues=neigh_vals[v->label];
          for(std::map<int,double>::iterator cval = cvalues.begin(); cval != cvalues.end();++cval)
            cell_val+=cval->second;
          cell_val/=(double)cvalues.size();
          avg_heat_val[v->label]=cell_val;
     }
 }

//This function hasn't been fully tested 
 void LobynessPlugin::NeighStdDev(std::map<int, double> &std_heat_val,std::map<int, double> &avg_heat_val, std::map<int,std::map<int, double> > &neigh_vals, Mesh *m){
     forall(const vertex &v, T){
          if(v->label==-1)
                continue;    
          double cell_val=0;
          std::map<int,double> cvalues=neigh_vals[v->label];
          double avg = avg_heat_val[v->label];
          for(std::map<int,double>::iterator cval = cvalues.begin(); cval != cvalues.end();++cval)
            cell_val+=(cval->second-avg)*(cval->second-avg);
          cell_val/=(double)cvalues.size();
       cell_val = sqrt(cell_val);
          std_heat_val[v->label]=cell_val;
      }

 }

*/
  REGISTER_PROCESS(LobynessPlugin);

  bool LobynessMeasureConvPeri::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";

    LobynessPlugin lp(*this);
    lp.run(m,m2,"Convex(Perimeter)",0,"No","No",fileName);

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

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureConvPeri);

  bool LobynessMeasureConvArea::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.run(m,m2,"Convex(Area)",0,"No","No",fileName);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Solidarity");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureConvArea);

  bool LobynessMeasureVisSto::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.run(m,m2,"Visibility(Stomata)",0,"No","No",fileName);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Visibility Stomata");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureVisSto);

  bool LobynessMeasureVisPav::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.run(m,m2,"Visibility(Pavement)",0,"No","No",fileName);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Visibility Pavement");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureVisPav);

  bool LobynessMeasureSkel::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.run(m,m2,"Skeleticity",0,"No","No",fileName);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Skeleticity");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
//  REGISTER_PROCESS(LobynessMeasureSkel);


  bool LobynessMeasureCirc::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.run(m,m2,"Circularity",0,"No","No",fileName);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Circularity");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureCirc);

  bool LobynessMeasureLES::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.run(m,m2,"LargestEmptySpace",0,"No","No",fileName);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Largest Empty Space");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureLES);


  bool LobynessMeasureLESOutput::run(Mesh *m, Mesh *m2)
  {
    if(m->graph().empty()) return false;
//    QString fileName = "";
    LobynessPlugin lp(*this);
    lp.parent = parent;
    lp.run(m,m2,"LargestEmptySpace",0,"No","Yes",LEC_CSV);

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double Lobeyness/Output Largest Empty Space");
    attrData.clear();

    forall(IntFloatPair p, m->labelHeat())
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(LobynessMeasureLESOutput);


  // Select the vertices belonging to the LES, based on center position
	// and radius given in the csv output file. 
  bool LobynessMeasureLESSelect::run(Mesh *m)
  {
    // Open csv file
    QFile file(LEC_CSV);
    if(!file.open(QIODevice::ReadOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(LEC_CSV));
      return false;
    }
    QTextStream ts(&file);
    QString line = ts.readLine();
    QStringList fields = line.split(",");
		if(fields.size() < 5) {
      setErrorMessage(QString("File '%1' should be a CSV file with at least 5 columns").arg(LEC_CSV));
      return false;
    }
		// Read the headers
		QHash<QString, int> fieldMap;
    int field_radius, field_x, field_y, field_z;
    int field_label;

    field_label = fields.indexOf("Label");
    if(field_label < 0)
      return setErrorMessage("Error, no field 'Label' in the CSV file");
    field_radius = fields.indexOf("Radius");
    if(field_radius < 0)
      return setErrorMessage("Error, no field 'Radius' in the CSV file");
    field_x = fields.indexOf("X position");
    if(field_x < 0)
      return setErrorMessage("Error, no field 'X position' in the CSV file");
    field_y = fields.indexOf("Y position");
    if(field_y < 0)
      return setErrorMessage("Error, no field 'Y position' in the CSV file");
    field_z = fields.indexOf("Z position");
    if(field_z < 0)
      return setErrorMessage("Error, no field 'Z position' in the CSV file");
    
		// Read the values of csv file, put the center position and radius in maps
		IntPoint3fMap oldCenter, newCenter; 
		IntFloatMap radius; 
		
		int nb_fields = fields.size();
    int linenum = 1;
    while(ts.status() == QTextStream::Ok) {
      ++linenum;
      fields = ts.readLine().split(",");
      if(fields.empty() or fields.size() == 1)
        break;
      if(fields.size() < nb_fields)
        return setErrorMessage(
          QString("Error on line %1: not enough fields (Expected: %2, Read: %3)").arg(linenum).arg(nb_fields).arg(
            fields.size()));

      bool ok;
      int label = fields[field_label].toInt(&ok);
      if(!ok) {
        setErrorMessage(
          QString("Error line %1: invalid label '%2'").arg(linenum).arg(fields[field_label].trimmed()));
        return false;
      }

			double radiusLES = fields[field_radius].toDouble(&ok);
      if(!ok) {
        setErrorMessage(
          QString("Error line %1: invalid label '%2'").arg(linenum).arg(fields[field_radius].trimmed()));
        return false;
      }
  
      Point3f centerLES;
      centerLES.x() = fields[field_x].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid x '%2'").arg(linenum).arg(fields[field_x]));
      centerLES.y() = fields[field_y].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid y '%2'").arg(linenum).arg(fields[field_y]));
      centerLES.z() = fields[field_z].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid z '%2'").arg(linenum).arg(fields[field_z]));
    
      radius[label] = radiusLES;
      oldCenter[label] = centerLES;
			newCenter[label] = centerLES;
    }

    // Find the center of the LES on the cell surface
		// NB the position can be a bit different if we don't use exactly the same mesh 
		// as the one used to compute the LES

		IntFloatMap dist;
		// If we have a normal mesh, label is stored in vertices
		if(m->meshType() == "MGXM"){
	    const vvGraph& S = m->graph();
		  // find LES center position on the mesh
	  	forall(const vertex &v, S){
        int label = v->label;
	  		if(label <= 0 or oldCenter.find(label) == oldCenter.end())
	  		  continue; 
	  		float d = norm(Point3f(v->pos) - oldCenter[label]);
	  		if(dist.find(label) != dist.end()) {// if not first run
	  		  if(d < dist[label])
	  		    newCenter[label] = Point3f(v->pos);
          dist[label] = d;
        }
	  	}
	  	// select the vertices
	  	forall(const vertex &v, S){
        int label = v->label;
	  		if(label <= 0 or oldCenter.find(label) == oldCenter.end())
	  		  continue; 
	  		float d = norm(Point3f(v->pos) - newCenter[label]);
	  		if(d < radius[label])
	  		  v->selected = true;
	  		else
	  		  v->selected = false;
	  	}	  
		} 
		// If we have a tissue mesh, label stored in cell centers
		else if(m->meshType() == "MGX3D"){
		  const cellGraph &Cells = m->cells();
		  // find LES center position on the mesh
      forall(const cell &c, Cells){
			  int label = c->label;
				if(oldCenter.find(label) == oldCenter.end())
				  continue; 
        forall(const vertex &v, c->S){
	  	  	float d = norm(Point3f(v->pos) - oldCenter[label]);
	  	  	if(dist.find(label) != dist.end()) {// if not first run
	  	  	  if(d < dist[label])
	  	  	    newCenter[label] = Point3f(v->pos);
            dist[label] = d;
          }
				}
	  	  // select the vertices
	  	  forall(const vertex &v, c->S){
	  	  	float d = norm(Point3f(v->pos) - newCenter[label]);
	  	  	if(d < radius[label])
	  	  	  v->selected = true;
	  	  	else
	  	  	  v->selected = false;
	  	  }	
			}
		}

		m->setShowMesh();
		m->setShowMeshLines(true);
		m->updateSelection();
		return true; 
  }
  REGISTER_PROCESS(LobynessMeasureLESSelect);

  typedef std::vector<Point3d> Point3dVec;
  typedef BoundingBox<3, double> BoundingBox3d;

  // RSS I think this one is already in Geometry, we should move all of them there.
  double signedTriArea(const Point3d &p1, const Point3d &p2, const Point3d &p3)
  {
    Point3d d21 = p2-p1;
    Point3d d32 = p3-p2;
    return (d21.x()*d32.y() - d21.y()*d32.x())/2.0;
  }

  Point3dVec convexHull(const Point3dVec &contour)
  {
    Point3dVec chull, cont = contour;
    bool convex = false;
    while(!convex && cont.size() > 3) {
      convex = true;
      for(uint i = 0; i < cont.size(); i++) {
        int ind1 = (i+cont.size()-1)%cont.size();
        int ind2 = i;
        int ind3 = (i+cont.size()+1)%cont.size(); 
        if(signedTriArea(cont[ind1], cont[ind2], cont[ind3]) > 0)
          chull.push_back(cont[i]);
        else
          convex = false;
      }
      cont = chull;
      chull.clear();
    }
    return cont;
  }

  // Rotate contour by a given theta
  Point3dVec rotateCont(Point3dVec contour, double theta)
  {
    forall(auto &p, contour)
      p = Point3d(cos(theta)*p.x()-sin(theta)*p.y(), sin(theta)*p.x()+cos(theta)*p.y(), 0.0);
  
    return contour;
  }

  // Compute the bounding box for a contour
  BoundingBox3d calcBBox(const Point3dVec &contour)
  {
    BoundingBox3d bBox;
    forall(const auto &p, contour)
      bBox |= p;
  
    return bBox;
  }
  
  // `Compute avg orientation based on vector from centroid to each point. Converts vectors to orientations and averages as "circular data"
//  Point3dVec rotateMinAxisApprox(Point3dVec contour)
//  {
//    Point3d centroid;
//    for(const auto &p : contour)
//      centroid += p;
//    centroid = centroid*(1.0/(double)contour.size());
//  
//    for(auto &p : contour)
//      p -= centroid;
//  
//    Point3d avgDir;
//    for(const auto &p : contour) {
//      double theta = atan2(p.y(), p.x()) + 2.0 * M_PI;
//      avgDir += Point3d(cos(2.0*theta), sin(2.0*theta), 0.0);
//    }
//    double avgTheta = -(atan2(avgDir.y(),avgDir.x())+2.0*M_PI)/2.0+M_PI/2.0;
//    return rotateCont(contour, avgTheta);
//  }
  Point3d findCentroid(const Point3dVec &contour) 
  {
    Point3d centroid;
    forall(const auto &p, contour)
      centroid += p;
    centroid = centroid*(1.0/(double)contour.size());

    return centroid;
  }
  
  Point3dVec translateContour(Point3dVec contour, const Point3d translate)
  {
    forall(auto &p, contour)
      p += translate;

    return contour;
  }

  Point3dVec centerContour(const Point3dVec &contour)
  {
    Point3d centroid = findCentroid(contour);
    return translateContour(contour, -centroid);
  }

  // Exhaustively try all rotations
  double findMinAxis(Point3dVec contour)
  {
    // Min axis of convex hull and contour must agree, so just work on chull
    Point3dVec chull = convexHull(centerContour(contour)); 

    double minTheta = 0;
    double minAxis = 1e20;
    Point3dVec rotatedHull = chull;
    for(double theta=0; theta<2*M_PI+0.1; theta+=0.001) { // Larger theta steps could probably be used
      rotatedHull = rotateCont(chull, theta);
      auto bBox = calcBBox(rotatedHull);
      double axis = bBox[1].x() - bBox[0].x();
      if(axis < minAxis){
        minAxis = axis;
        minTheta = theta;
      }
    }
    return minTheta;
  }

  Point3dVec rotateMinAxis(Point3dVec contour)
  {
    double theta = findMinAxis(contour);
    return rotateCont(contour, theta);
  }

  // Exhaustively try all rotations
  double findMinRectangle(Point3dVec contour)
  {
    // Min axis of convex hull and contour must agree, so just work on chull
    Point3dVec chull = convexHull(centerContour(contour)); 

    double minTheta = 0;
    double minArea= 1e20;
    Point3dVec rotatedHull = chull;
    for(double theta=0; theta<2*M_PI+0.1; theta+=0.001) { // Larger theta steps could probably be used
      rotatedHull = rotateCont(chull, theta);
      auto bBox = calcBBox(rotatedHull);
      double area = (bBox[1].x() - bBox[0].x()) * (bBox[1].y() - bBox[0].y());
      if(area < minArea){
        minArea = area;
        minTheta = theta;
      }
    }
    return minTheta;
  }

  Point3dVec rotateMinRectangle(Point3dVec contour)
  {
    double theta = findMinRectangle(contour);
    return rotateCont(contour, theta);
  }

  double polyArea(const Point3dVec &contour)
  {
    double area = 0;
    for(uint i = 1; i < contour.size(); i++)
      area += signedTriArea(Point3d(0,0,0), contour[i-1], contour[i]);
    area += signedTriArea(Point3d(0,0,0), contour[contour.size()-1], contour[0]);
    return area;
  }
  

  bool MeasureRectangularity::run(Mesh *m, IntFloatAttr& heatMap)
  {
    heatMap.clear();

    if(m->meshType() != "MGXC" and m->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    if(m->graph().empty()) return false;

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


    if(m->tissue().meshType() == "MGXM")
      mgxmToMgxc(S, T, 1);
    else
      T = S;

     IntFloatAttr &heat_val = m->labelHeat();
     heat_val.clear();


    #pragma omp parallel for
    for(uint i = 0; i < T.size(); i++){
      const vertex &v = T[i];
       if(v->label==-1)
         continue;
       std::vector<Point3d> contour;
       Point3d ppos;
       Point3d pnorm;
       Matrix3d prot;
       forall(const vertex &s, T.neighbors(v))
         contour.push_back(s->pos);

       findPolygonPlane(contour,ppos,pnorm);
       projectPointsOnPlane(contour,ppos,pnorm);
       rotatePointsIntoXY(pnorm,contour,prot);

      // Flip if inverted
      if(polyArea(contour) < 0)
        for(uint j = 0; j < contour.size(); j++)
          contour[j].x() *= -1;

       auto bBox = calcBBox(rotateMinRectangle(contour));
       heatMap[v->label] = polyArea(contour)/((bBox[1].x() - bBox[0].x()) * (bBox[1].y() - bBox[0].y()));
     }

    m->labelHeat() = heatMap;
    m->setShowLabel("Label Heat");
    m->heatMapUnit() = "";
    m->heatMapBounds() = m->calcHeatMapBounds();
    m->updateTriangles();

    AttrMap<int, double>& attrData = m->attributes().attrMap<int, double>("Measure Label Double /Lobeyness/Rectangularity");
    attrData.clear();

    forall(IntFloatPair p, heatMap)
      attrData[p.first] = p.second;

    return true;
  }
  REGISTER_PROCESS(MeasureRectangularity);

}

