From b1b18ae49b0d22d3fbd5146eb6416c8b9e4dd62c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 12 Dec 2018 23:46:58 -0500 Subject: added support for voronoi graph generation—existing example code does not use it yet MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- lib/CMakeLists.txt | 2 +- lib/include/graph.hpp | 3 ++ lib/src/graph.cpp | 125 ++++++++++++++++++++++++++++++++++++++++++++++++++ lib/voronoi | 1 + 4 files changed, 130 insertions(+), 1 deletion(-) create mode 160000 lib/voronoi (limited to 'lib') diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index bf3c3b5..49a9602 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -1,4 +1,4 @@ add_library(frac SHARED src/graph.cpp src/network.cpp) -target_include_directories(frac PUBLIC include) +target_include_directories(frac PUBLIC include PRIVATE voronoi/src) diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp index 9a630ff..75d54f4 100644 --- a/lib/include/graph.hpp +++ b/lib/include/graph.hpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include class graph { public: @@ -27,5 +29,6 @@ class graph { std::vector dual_edges; graph(unsigned int Nx, unsigned int Ny); + graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double = 0.25); }; 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 +#include + +#define JC_VORONOI_IMPLEMENTATION +#define JCV_REAL_TYPE double +#define JCV_ATAN2 atan2 +#define JCV_FLT_MAX 1.7976931348623157E+308 +#include 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 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 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 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::iterator it1 = known_vertices.find(r1); + std::map::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); +} + diff --git a/lib/voronoi b/lib/voronoi new file mode 160000 index 0000000..0e0f6ea --- /dev/null +++ b/lib/voronoi @@ -0,0 +1 @@ +Subproject commit 0e0f6ea352087ea6ad1e0484a55c6ab27462f4f7 -- cgit v1.2.3-70-g09d2