summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/graph.cpp125
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);
+}
+