more fixes to voronoi periodization, some upstream errors now caught and construction is retried

1 files changed, 116 insertions, 34 deletions

diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 7564fa7..7c50ed0 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -6,6 +6,7 @@
 #define JCV_REAL_TYPE double
 #define JCV_ATAN2 atan2
 #define JCV_FLT_MAX 1.7976931348623157E+308
+#define NDEBUG
 #include <jc_voronoi.h>
 
 // actual mod for floats
@@ -72,6 +73,69 @@ class eulerException: public std::exception
   }
 } eulerex;
 
+class clippingException: public std::exception
+{
+  virtual const char* what() const throw()
+  {
+    return "An interior site has a clipped edge and the tesselation is malformed.";
+  }
+} clippingex;
+
+class triangleException: public std::exception
+{
+  virtual const char* what() const throw()
+  {
+    return "A dual-centered triangle has an impossible geometry and the tesselation is malformed.";
+  }
+} triex;
+
+unsigned int get_triangle_signature(unsigned int j1, unsigned int j2, unsigned int j3) {
+  // this yucky function takes three unsigned integers representing the
+  // location in the nine periodic copies of each corner of a delauney triangle
+  // and returns a signature for that triangle, which is an unsigned integer
+  // that uniquely labels the way the triangle crosses boundaries of the
+  // copies. This allows us to differentiate delauney triangles with identical
+  // vertices but which should be identified with different faces
+  unsigned int x1 = j1 % 3;
+  unsigned int y1 = j1 / 3;
+
+  unsigned int x2 = j2 % 3;
+  unsigned int y2 = j2 / 3;
+
+  unsigned int x3 = j3 % 3;
+  unsigned int y3 = j3 / 3;
+
+  if ((j1 == j2) && (j2 == j3)) {
+    return 0;
+  } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) || ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
+    return 1;
+  } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) || ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
+    return 2;
+  } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) || ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
+    return 3;
+  } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) || ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
+    return 4;
+  } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) || ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
+    return 5;
+  } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) || ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
+    return 6;
+  } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) || ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
+    return 7;
+  } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) || ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
+    return 8;
+  } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
+    return 9;
+  } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
+    return 10;
+  } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
+    return 11;
+  } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
+    return 12;
+  } else {
+    throw triex;
+  }
+}
+
 graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) {
   L = {Lx, Ly};
 
@@ -108,15 +172,15 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
     for (unsigned int i = 0; i < nv; i++) {
       const vertex& v = vertices[i];
 
-      points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
-      points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
-      points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
-      points[9 * i + 3] = {v.r.x + 0.0, v.r.y + L.y};
-      points[9 * i + 4] = {v.r.x + 0.0, v.r.y - L.y};
-      points[9 * i + 5] = {v.r.x + L.x, v.r.y + L.y};
+      points[9 * i + 0] = {v.r.x - L.x, v.r.y - L.y};
+      points[9 * i + 1] = {v.r.x + 0.0, v.r.y - L.y};
+      points[9 * i + 2] = {v.r.x + L.x, v.r.y - L.y};
+      points[9 * i + 3] = {v.r.x - L.x, v.r.y + 0.0};
+      points[9 * i + 4] = {v.r.x + 0.0, v.r.y + 0.0};
+      points[9 * i + 5] = {v.r.x + L.x, v.r.y + 0.0};
       points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
-      points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
-      points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
+      points[9 * i + 7] = {v.r.x + 0.0, v.r.y + L.y};
+      points[9 * i + 8] = {v.r.x + L.x, v.r.y + L.y};
     }
 
     // run voronoi
@@ -127,7 +191,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
     for (int i = 0; i < diagram.numsites; i++) {
       const jcv_site* site = &sites[i];
       unsigned int ind = site->index;
-      if (ind % 9 == 0) {
+      if (ind % 9 == 4) {
         double Cx = 0.0;
         double Cy = 0.0;
         double A = 0.0;
@@ -164,35 +228,30 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
 
   for (unsigned int i = 0; i < nv; i++) {
     const vertex& v = vertices[i];
-    points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
-    points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
-    points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
-    points[9 * i + 3] = {v.r.x + 0.0, v.r.y + L.y};
-    points[9 * i + 4] = {v.r.x + 0.0, v.r.y - L.y};
-    points[9 * i + 5] = {v.r.x + L.x, v.r.y + L.y};
+    points[9 * i + 0] = {v.r.x - L.x, v.r.y - L.y};
+    points[9 * i + 1] = {v.r.x + 0.0, v.r.y - L.y};
+    points[9 * i + 2] = {v.r.x + L.x, v.r.y - L.y};
+    points[9 * i + 3] = {v.r.x - L.x, v.r.y + 0.0};
+    points[9 * i + 4] = {v.r.x + 0.0, v.r.y + 0.0};
+    points[9 * i + 5] = {v.r.x + L.x, v.r.y + 0.0};
     points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
-    points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
-    points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
+    points[9 * i + 7] = {v.r.x + 0.0, v.r.y + L.y};
+    points[9 * i + 8] = {v.r.x + L.x, v.r.y + L.y};
   }
   jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram);
 
   const jcv_site* sites = jcv_diagram_get_sites(&diagram);
 
-  struct paircomp {
-    bool operator() (const std::array<unsigned int, 2>& lhs, const std::array<unsigned int, 2>& rhs) const
-    {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1];
-      else return rhs[0] < rhs[1];
-    }
-  };
-
-  std::unordered_map<std::array<unsigned int, 3>, unsigned int> known_vertices;
+  std::unordered_map<std::array<unsigned int, 4>, unsigned int> known_vertices;
 
   for (int i = 0; i < diagram.numsites; i++) {
     const jcv_site* site = &sites[i];
 
     // we only care about processing the cells of our original, central sites
-    if (site->index % 9 == 0) {
+    if (site->index % 9 == 4) {
+      bool self_bonded = false;
       unsigned int i1 = (unsigned int)(site->index / 9);
+      unsigned int j1 = (unsigned int)(site->index % 9);
       const jcv_graphedge* e = site->edges;
       const jcv_graphedge* ep = site->edges;
       while (ep->next) {
@@ -204,13 +263,26 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
         // sites whose index is larger than ours are given edges, so we
         // don't double up later
         const jcv_site* neighbor = e->neighbor;
+        if (neighbor == NULL) {
+          throw clippingex;
+        }
         unsigned int i2 = (unsigned int)(neighbor->index / 9);
+        unsigned int j2 = (unsigned int)(neighbor->index % 9);
+        unsigned int x2 = j2 % 3;
+        unsigned int y2 = j2 / 3;
 
         vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y});
 
+        if (ep->neighbor == NULL) {
+          throw clippingex;
+        }
         unsigned int i3p = (unsigned int)(ep->neighbor->index / 9);
-        std::array<unsigned int, 3> t1 = {i1, i2, i3p};
-        std::sort(t1.begin(), t1.end());
+        unsigned int j3p = (unsigned int)(ep->neighbor->index % 9);
+
+        unsigned int sig1 = get_triangle_signature(j1, j2, j3p);
+
+        std::array<unsigned int, 4> t1 = {i1, i2, i3p, sig1};
+        std::sort(t1.begin(), t1.begin() + 3);
 
         auto it1 = known_vertices.find(t1);
 
@@ -225,10 +297,12 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
           dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
         }
         
-        if (i1 < i2) {
-          unsigned int nn = neighbor->index % 9;
-          bool crossed_x = (nn == 1) || (nn == 2) || (nn > 4);
-          bool crossed_y = nn > 2;
+        if (i1 < i2 || (i1 == i2 && !self_bonded)) {
+          if (i1 == i2) {
+            self_bonded = true;
+          }
+          bool crossed_x = x2 != 1;
+          bool crossed_y = y2 != 1;
           edges.push_back({{i1, i2},
               {mod((site->p.x + neighbor->p.x) / 2, L.x),
                mod((site->p.y + neighbor->p.y) / 2, L.y)},
@@ -242,9 +316,17 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
             en = e->next;
           }
 
+          if (en->neighbor == NULL) {
+            throw clippingex;
+          }
+
           unsigned int i3n = (unsigned int)(en->neighbor->index / 9);
-          std::array<unsigned int, 3> t2 = {i1, i2, i3n};
-          std::sort(t2.begin(), t2.end());
+          unsigned int j3n = (unsigned int)(en->neighbor->index % 9);
+
+          unsigned int sig2 = get_triangle_signature(j1, j2, j3n);
+
+          std::array<unsigned int, 4> t2 = {i1, i2, i3n, sig2};
+          std::sort(t2.begin(), t2.begin() + 3);
 
           auto it2 = known_vertices.find(t2);