1 files changed, 125 insertions, 0 deletions
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index f2e8065..f25cf87 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,5 +1,13 @@
 
 #include "graph.hpp"
+#include <cstring>
+#include <iostream>
+
+#define JC_VORONOI_IMPLEMENTATION
+#define JCV_REAL_TYPE double
+#define JCV_ATAN2 atan2
+#define JCV_FLT_MAX 1.7976931348623157E+308
+#include <jc_voronoi.h>
 
 graph::graph(unsigned int Nx, unsigned int Ny) {
   L = {(double)Nx, (double)Ny};
@@ -36,3 +44,120 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
   }
 }
 
+graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
+  L = {(double)Nx, (double)Ny};
+
+  unsigned int dnv = Nx * Ny / 2;
+
+  dual_vertices.resize(dnv);
+  
+  std::normal_distribution<double> d(0.0, spread);
+
+  // the coordinates of the dual lattice, from which a delaunay triangulation
+  // and corresponding voronoi tesselation will be built
+  for (unsigned int i = 0; i < dnv; i++) {
+    double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng);
+    double ry = (double)(i / (Nx / 2)) + d(rng);
+    dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
+  }
+
+  // to make the resulting tesselation periodic, we tile eight (!) copies of
+  // the original dual points for a total of nine. note that the index of each
+  // point quotient 9 is equal to the original index (we'll use this to
+  // translate back)
+  std::vector<jcv_point> points;
+  points.reserve(9 * dnv);
+  for (const vertex &v : dual_vertices) {
+    points.push_back({v.r.x, v.r.y});
+    points.push_back({v.r.x + L.x, v.r.y + 0.0});
+    points.push_back({v.r.x - L.x, v.r.y + 0.0});
+    points.push_back({v.r.x + 0.0, v.r.y + L.y});
+    points.push_back({v.r.x + 0.0, v.r.y - L.y});
+    points.push_back({v.r.x + L.x, v.r.y + L.y});
+    points.push_back({v.r.x - L.x, v.r.y + L.y});
+    points.push_back({v.r.x + L.x, v.r.y - L.y});
+    points.push_back({v.r.x - L.x, v.r.y - L.y});
+  }
+
+  jcv_diagram diagram;
+  memset(&diagram, 0, sizeof(jcv_diagram));
+
+  jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram);
+
+  const jcv_site* sites = jcv_diagram_get_sites(&diagram);
+
+  struct coorcomp {
+    bool operator() (const coordinate& lhs, const coordinate& rhs) const
+    {return lhs.x<rhs.x;}
+  };
+
+  std::map<coordinate, unsigned int, coorcomp> 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) {
+      const jcv_graphedge* e = site->edges;
+      // for each edge on the site's cell
+      while(e) {
+        // assess whether the dual edge needs to be added. only neighboring
+        // dual sites whose index is larger than ours are given edges, so we
+        // don't double up later
+        const jcv_site* neighbor = e->neighbor;
+        unsigned int index = (unsigned int)(site->index / 9);
+        unsigned int real_index = (unsigned int)(neighbor->index / 9);
+        
+        if (index < real_index) {
+          dual_edges.push_back({index, real_index});
+          coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))};
+          coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))};
+
+          std::map<coordinate, unsigned int>::iterator it1 = known_vertices.find(r1);
+          std::map<coordinate, unsigned int>::iterator it2 = known_vertices.find(r2);
+
+          unsigned int vi1, vi2;
+
+          if (it1 == known_vertices.end()) {
+            vi1 = vertices.size();
+            vertices.push_back({r1});
+            known_vertices[r1] = vi1;
+          } else {
+            vi1 = it1->second;
+          }
+          if (it2 == known_vertices.end()) {
+            vi2 = vertices.size();
+            vertices.push_back({r2});
+            known_vertices[r2] = vi2;
+          } else {
+            vi2 = it2->second;
+          }
+
+          edges.push_back({vi1, vi2});
+        }
+
+        e = e->next;
+      }
+    }
+  }
+
+  for (edge &e : edges) { 
+    std::cout << e[0] << " " << e[1] << " ";
+  }
+  std::cout << "\n";
+  for (vertex &v : vertices) { 
+    std::cout << v.r.x << " " << v.r.y << " ";
+  }
+  std::cout << "\n";
+  for (edge &e : dual_edges) { 
+    std::cout << e[0] << " " << e[1] << " ";
+  }
+  std::cout << "\n";
+  for (vertex &v : dual_vertices) { 
+    std::cout << v.r.x << " " << v.r.y << " ";
+  }
+  std::cout << "\n";
+
+  jcv_diagram_free(&diagram);
+}
+