lloyds algorithm in l1 metric

Mesh/meshGFaceLloyd.cpp

@@ -9,6 +9,18 @@
 #include "meshGFace.h"
 #include "BackgroundMesh.h"
 
+/*List of parameters to modify*/
+//Line 76 : the variable "number" contains the number of points to remove.
+//Line 165 : the variable "test" is printed by the script.
+//Line 190 : the variable "angle" specifies the orientation of the mesh.
+//Line 287 : the variable "num" contains the number of iteration of the bisection method
+//used to find the fixed point.
+//Line 447 : it is the function "optimize" that needs to be called by the script.
+//Line 455 : put "lloyd(face)" in commentary to avoid calling the l1 lloyds method.
+//Line 456 : put "recombineIntoQuads(face,1,1)" in commentary to avoid calling the quad 
+//meshing method.
+//Line 461 : the code below line 461 refers to Lua and might need to be removed.
+
 void lloydAlgorithm::operator () (GFace *gf) 
 {
   std::set<MVertex*> all;
@@ -57,6 +69,26 @@ void lloydAlgorithm::operator () (GFace *gf)
     triangulator.data(i++) = (*it);
   }
  
+  triangulator.Voronoi();
+  triangulator.initialize();
+  int index,number,count,max;
+  bool recent;
+  number = 0; //This variable contains the number of points to remove.
+  count = 0;
+  max = 100000;
+  for(int i=0;i<max;i++)
+  {
+    if(count>=number) break;
+	index = (int)((triangulator.numPoints-1)*((double)rand()/(double)RAND_MAX));
+	PointRecord& pt = triangulator.points[index];
+	MVertex* v = (MVertex*)pt.data;
+	if(v->onWhat()==gf && !triangulator.onHull(index)){
+	  recent = triangulator.remove_point(index);
+	  if(recent) count++;
+	}
+  }
+  triangulator.remove_all();
+	
   // compute the Voronoi diagram
   triangulator.Voronoi();
   //printf("hullSize = %d\n",triangulator.hullSize());
@@ -81,12 +113,13 @@ void lloydAlgorithm::operator () (GFace *gf)
         triangulator.voronoiCell (i,pts); 
         double E, A;
         SPoint2 p(pt.where.h,pt.where.v);
-        if (!infiniteNorm){
-          centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E, A, NULL); //backgroundMesh::current());       
-        }
-        else {
-          centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A);          
+		for(int k=0;k<pts.size();k++){
+		  rotate(p,pts[k]);
 		}
+		SPoint2 point = fixed_point(pts); //Fixed point for quad mesh generation
+		unrotate(p,point);
+		cgs(i,0) = point.x();
+		cgs(i,1) = point.y();
         ENERGY += E;
 	double d = sqrt((p.x()-cgs(i,0))*(p.x()-cgs(i,0))+
 			(p.y()-cgs(i,1))*(p.y()-cgs(i,1)));
@@ -119,6 +152,7 @@ void lloydAlgorithm::operator () (GFace *gf)
   }
 	
   // now create the vertices
+  test = 0;
   std::vector<MVertex*> mesh_vertices;
   for (int i=0; i<triangulator.numPoints;i++){
     // get the ith vertex
@@ -128,6 +162,7 @@ void lloydAlgorithm::operator () (GFace *gf)
       GPoint gp = gf->point (pt.where.h,pt.where.v);
       MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
 	  mesh_vertices.push_back(v);
+	  test++;
     }
   }
 
@@ -149,3 +184,288 @@ void lloydAlgorithm::operator () (GFace *gf)
   // clear the list of additional vertices
   gf->_additional_vertices.clear();
 }
+
+double lloydAlgorithm::get_angle(double x,double y){
+  double angle;
+  angle = M_PI/4.0; //This variable specifies the orientation of the mesh.
+  return angle;
+}
+
+//This function returns the a, b, c and d components of the rotation matrix.
+void lloydAlgorithm::get_rotation(double x,double y,double& a,double& b,double& c,double& d){
+  double angle;
+  angle = get_angle(x,y);
+  a = cos(angle);
+  b = -sin(angle);
+  c = sin(angle);
+  d = cos(angle);
+}
+
+//This function returns the a, b, c and d components of the inverse rotation matrix.
+void lloydAlgorithm::get_unrotation(double x,double y,double& _a,double& _b,double& _c,double& _d){
+  double a,b,c,d;
+  double det;
+  get_rotation(x,y,a,b,c,d);
+  det = 1.0/(a*d-b*c);
+  _a = d/det;
+  _b = -b/det;
+  _c = -c/det;
+  _d = a/det;
+}
+
+//This function rotates a vertice (point) belonging to a Voronoi cell.
+//The angle is specified by the location of the generator, because we want all vertices
+//of the Voronoi cell to be rotated by the same amount in case the orientation of the
+//mesh is not constant.
+//The generator itself is not rotated.
+void lloydAlgorithm::rotate(SPoint2 generator,SPoint2& point){
+  double x,y;
+  double new_x,new_y;
+  double a,b,c,d;
+  get_rotation(generator.x(),generator.y(),a,b,c,d);
+  x = point.x();
+  y = point.y();
+  new_x = a*x + b*y;
+  new_y = c*x + d*y;
+  point.setPosition(new_x,new_y);
+}
+
+//This function rotates a vertice (point) belonging to a Voronoi cell in the reverse direction.
+//The angle is specified by the location of the generator, because we want all vertices
+//of the Voronoi cell to be rotated by the same amount in case the orientation of the
+//mesh is not constant.
+//The generator itself is not rotated.
+void lloydAlgorithm::unrotate(SPoint2 generator,SPoint2& point){
+  double x,y;
+  double new_x,new_y;
+  double a,b,c,d;
+  get_unrotation(generator.x(),generator.y(),a,b,c,d);
+  x = point.x();
+  y = point.y();
+  new_x = a*x + b*y;
+  new_y = c*x + d*y;
+  point.setPosition(new_x,new_y);
+}
+
+//This function returns the fixed point of the voronoi cell specified in argument (polygon).
+//It finds the bounding box of the voronoi cell and then it calls the function range in x and y.
+SPoint2 lloydAlgorithm::fixed_point(const std::vector<SPoint2>& polygon){
+  int i,num;
+  double x,y;
+  double x_min,y_min;
+  double x_max,y_max;
+  double solution_x,solution_y;
+  SPoint2 point;
+  x_min = 100000000.0;
+  y_min = 100000000.0;
+  x_max = -100000000.0;
+  y_max = -100000000.0;
+  num = polygon.size();
+  for(i=0;i<num;i++){
+	point = polygon[i];
+    x = point.x();
+	y = point.y();
+	if(x < x_min) x_min = x;
+	if(y < y_min) y_min = y;
+	if(x > x_max) x_max = x;
+	if(y > y_max) y_max = y;
+  }
+  solution_x = range(polygon,x_min,x_max,1);
+  solution_y = range(polygon,y_min,y_max,0);
+  return SPoint2(solution_x,solution_y);
+}
+
+//This function returns the x or y median of the Voronoi cell specified in argument (polygon).
+//If direction = 1, the function will return the x median.
+//If direction = 0, the function will return the y median.
+//The variables u1 and u2 should be respectively equal to the minimum and the maximum of x or y,
+//depending of the value of direction.
+//This function uses the bisection method.
+double lloydAlgorithm::range(const std::vector<SPoint2>& polygon,double u1,double u2,bool direction){
+  int i,num;
+  double mid,val;
+  num = 16; //This variable specifies the number of iteration of the bisection method.
+  for(i=0;i<num;i++){
+	mid = 0.5*(u1 + u2);
+	val = gradient(polygon,mid,direction);
+	if(val >= 0) u2 = mid;
+	else u1 = mid;
+  }
+  return 0.5*(u1 + u2);
+}
+
+//This function returns the x or y component of the l1 gradient of the voronoi cell
+//specified in argument (polygon).
+//If direction = 1, it will return the x component.
+//If direction = 0, it will return the y component.
+//The variable generator should be equal to the x or y component of the point to which the 
+//voronoi cell belongs.
+//x component = A1-A2, A1 : area of the left part, A2 : area of the right part
+//y component = A3-A4, A3 : area of the lower part, A4 : area of the higher part
+double lloydAlgorithm::gradient(const std::vector<SPoint2>& polygon,double generator,bool direction){
+  int i,num,index1,index2;
+  bool intersection,side;
+  double areaA,areaB;
+  SPoint2 p1,p2,solution;
+  std::vector<bool> category;
+  std::vector<SPoint2> new_polygon;
+  std::vector<SPoint2> polygonA;
+  std::vector<SPoint2> polygonB;
+  num = polygon.size();
+  for(i=0;i<num;i++){
+	index1 = i;
+	index2 = (i+1)%num;
+	p1 = polygon[index1];
+	p2 = polygon[index2];
+	solution = line_intersection(p1,p2,generator,direction,intersection);
+	if(intersection){
+	  new_polygon.push_back(p1);
+	  new_polygon.push_back(solution);
+	  category.push_back(1);
+	  category.push_back(0);	
+	}
+	else{
+	  new_polygon.push_back(p1);
+	  category.push_back(1);	
+	}
+  }
+  num = new_polygon.size();
+  for(i=0;i<num;i++){
+	if(category[i]){
+	  side = point_side(new_polygon[i],generator,direction);
+	  if(side){
+	    polygonA.push_back(new_polygon[i]);
+	  }
+	  else{
+	    polygonB.push_back(new_polygon[i]);
+	  }
+	}
+	else{
+	  polygonA.push_back(new_polygon[i]);
+	  polygonB.push_back(new_polygon[i]);
+	}
+  }
+  areaA = polygon_area(polygonA);
+  areaB = polygon_area(polygonB);
+  return areaA-areaB;
+}
+
+//The point p is a vertice of a Voronoi cell.
+//The variable generator contains the x or y component of the point to which
+//the Voronoi cell belongs. If direction = 1, it will contain the x component,
+//if direction = 0, it will contain the y component.
+//If direction = 1, this function will return 1 if p is at the left of the generator
+//and it will return 0 if p is at the right of the generator.
+//If direction = 0, this function will return 1 if p is below the generator
+//and it will return 0 if p is higher than the generator.
+bool lloydAlgorithm::point_side(SPoint2 p,double generator,bool direction){
+  if(direction){
+    if(p.x()<=generator) return 1;
+	else return 0;
+  }
+  else{
+    if(p.y()<=generator) return 1;
+	else return 0;
+  }
+}
+
+//This function uses the Surveyors Formula to calculate the area of the Voronoi cell
+//specified in argument (polygon).
+double lloydAlgorithm::polygon_area(const std::vector<SPoint2>& polygon){
+  int i;
+  int num;
+  double area;
+  num = polygon.size();
+  if(num<3) return 0.0;
+  else{
+    area = 0.0;
+    for(i=0;i<num-1;i++){
+	  area = area + 0.5*(polygon[i].x()*polygon[i+1].y() - polygon[i+1].x()*polygon[i].y());
+    }
+    area = area + 0.5*(polygon[num-1].x()*polygon[0].y() - polygon[0].x()*polygon[num-1].y());
+    return fabs(area);
+  }
+}
+
+//This function determines if there is an intersection between the segment p1-p2 and the
+//infinite line passing through the generator, which is the point to which the Voronoi cell
+//belongs.
+//If there is an intersection, then the variable intersection will contain 1, otherwise it will
+//contain 0. In that case, the point that the function returns is not imporant.
+//If direction = 1, the variable generator will contain the x component, in that case,
+//the infinite line will be vertical.
+//If direction = 0, the variable generator will contain the y component, in that case,
+//the infinite line will be horizontal. 
+SPoint2 lloydAlgorithm::line_intersection(SPoint2 p1,SPoint2 p2,double generator,bool direction,bool& intersection){
+  double x1,x2,solution_x;
+  double y1,y2,solution_y;
+  double u,e;
+  x1 = p1.x();
+  y1 = p1.y();
+  x2 = p2.x();
+  y2 = p2.y();
+  e = 0.0000001;
+  intersection = 0;
+  if(direction){
+	if(fabs(x1 - generator)<e){
+	  intersection = 1;
+	  return SPoint2(x1,y1);
+	}
+	else if(fabs(x2 - generator)<e){
+	  intersection = 1;
+	  return SPoint2(x2,y2);
+	}
+	else if((x1<generator && x2>generator)||(x1>generator && x2<generator)){
+	  intersection = 1;
+	  u = (generator-x1)/(x2-x1);
+	  solution_x = generator;
+	  solution_y = y1 + u*(y2-y1);
+	  return SPoint2(solution_x,solution_y);
+	}
+  }
+  else{
+	if(fabs(y1 - generator)<e){
+	  intersection = 1;
+	  return SPoint2(x1,y1);
+	}
+	else if(fabs(y2 - generator)<e){
+	  intersection = 1;
+	  return SPoint2(x2,y2);
+	}
+	else if((y1<generator && y2>generator)||(y1>generator && y2<generator)){
+	  intersection = 1;
+	  u = (generator-y1)/(y2-y1);
+	  solution_x = x1 + u*(x2-x1);
+	  solution_y = generator;
+	  return SPoint2(solution_x,solution_y);
+	}
+  }
+  return SPoint2(0.0,0.0);
+}
+
+//This is the function called by the script.
+double lloydAlgorithm::optimize(int max,int flag){
+  GFace*face;
+  GModel*model = GModel::current();
+  GModel::fiter iterator;
+  lloydAlgorithm lloyd(max,flag);
+  for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
+  {
+    face = *iterator;
+	lloyd(face);
+	recombineIntoQuads(face,1,1);
+  }
+  return lloyd.test;
+}
+
+#include "Bindings.h"
+void lloydAlgorithm::registerBindings(binding *b){
+  classBinding*cb = b->addClass<lloydAlgorithm>("lloydAlgorithm");
+  cb->setDescription("This class is used to move the points according to lloyd algorithm.");
+  methodBinding*cm;
+  cm = cb->setConstructor<lloydAlgorithm>();
+  cm->setDescription("This is the constructor.");
+  cm = cb->addMethod("optimize",&lloydAlgorithm::optimize);
+  cm->setDescription("This function moves the points.");
+  cm->setArgNames("max","flag",NULL);
+}

Mesh/meshGFaceLloyd.h

@@ -5,16 +5,34 @@
 
 #ifndef _MESH_GFACE_LLOYD_H_
 #define _MESH_GFACE_LLOYD_H_
+#include "fullMatrix.h"
+#include "DivideAndConquer.h"
+#include "Bindings.h"
 
 class GFace;
 
 class lloydAlgorithm {
   int ITER_MAX;
   bool infiniteNorm;
+  double test;
  public :
   lloydAlgorithm (int itermax, bool infnorm = false)
   : ITER_MAX(itermax), infiniteNorm(infnorm) {}
+  lloydAlgorithm() {}
   void operator () (GFace *);
+  double get_angle(double,double);
+  void get_rotation(double,double,double&,double&,double&,double&);
+  void get_unrotation(double,double,double&,double&,double&,double&);
+  void rotate(SPoint2,SPoint2&);
+  void unrotate(SPoint2,SPoint2&);
+  SPoint2 fixed_point(const std::vector<SPoint2>&);
+  double range(const std::vector<SPoint2>&,double,double,bool);
+  double gradient(const std::vector<SPoint2>&,double,bool);
+  bool point_side(SPoint2,double,bool);
+  double polygon_area(const std::vector<SPoint2>&);
+  SPoint2 line_intersection(SPoint2,SPoint2,double,bool,bool&);
+  double optimize(int,int);
+  static void registerBindings(binding *b);
 };
 
 #endif