lloyds algorithm in l1 metric

Numeric/DivideAndConquer.cpp

@@ -824,7 +824,7 @@ DocRecord::DocRecord(int n)
     numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
 {
   if(numPoints)
-    points = new PointRecord[numPoints];
+    points = new PointRecord[numPoints+1000];
 }
 
 DocRecord::~DocRecord()
@@ -864,6 +864,48 @@ void DocRecord::setPoints(fullMatrix<double> *p)
   }
 } 
 
+void DocRecord::initialize(){
+  int i;
+  for(i=0;i<numPoints;i++){
+	points[i].flag = 0;
+  }
+}
+
+bool DocRecord::remove_point(int index){
+  if(points[index].flag == 0){
+    points[index].flag = 1;
+    return 1;
+  }
+  else return 0;
+}
+
+void DocRecord::remove_all(){
+  int i;
+  int index;
+  int numPoints2;
+  PointRecord* points2;
+  numPoints2 = 0;
+  for(i=0;i<numPoints;i++){
+    if(points[i].flag==0){
+      numPoints2++;
+    }
+  }
+  points2 = new PointRecord[numPoints2];
+  index = 0;
+  for(i=0;i<numPoints;i++){
+	if(points[i].flag==0){
+	  points2[index].where.h = points[i].where.h;
+	  points2[index].where.v = points[i].where.v;
+	  points2[index].data = points[i].data;
+	  points2[index].flag = points[i].flag;
+	  index++;
+	}
+  }
+  delete [] points;
+  points = points2;
+  numPoints = numPoints2;
+}
+
 #include "Bindings.h"
 
 void DocRecord::registerBindings(binding *b)

Numeric/DivideAndConquer.h

@@ -29,6 +29,7 @@ struct PointRecord {
   DPoint where;
   DListPeek adjacent;
   void *data;
+  int flag; //0:to be kept, 1:to be removed
   PointRecord() : adjacent(0), data (0) {}
 };
 
@@ -60,7 +61,6 @@ class DocRecord{
  private:
   int _hullSize;
   PointNumero *_hull;
-  STriangle *_adjacencies;
   PointNumero Predecessor(PointNumero a, PointNumero b);
   PointNumero Successor(PointNumero a, PointNumero b);
   int FixFirst(PointNumero x, PointNumero f);
@@ -76,7 +76,6 @@ class DocRecord{
   int Insert(PointNumero a, PointNumero b);
   int DListDelete(DListPeek *dlist, PointNumero oldPoint);
   int Delete(PointNumero a, PointNumero b);
-  PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n);
   int ConvertDListToTriangles();
   void ConvertDListToVoronoiData();
   void RemoveAllDList();
@@ -84,6 +83,7 @@ class DocRecord{
   int CountPointsOnHull();
   void ConvexHull();
  public:
+  STriangle *_adjacencies;
   int numPoints;        // number of points
   PointRecord *points;  // points to triangulate
   int numTriangles;     // number of triangles
@@ -102,6 +102,10 @@ class DocRecord{
   void printMedialAxis(Octree *_octree, std::string, GFace *gf=NULL, GEdge *ge=NULL);
   double Lloyd (int);
   void voronoiCell (PointNumero pt, std::vector<SPoint2> &pts) const;
+  void initialize();
+  bool remove_point(int);
+  void remove_all();
+  PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n);
   static void registerBindings(binding *b);
 };