diff options
Diffstat (limited to 'lib/src')
| -rw-r--r-- | lib/src/graph.cpp | 124 |
1 files changed, 92 insertions, 32 deletions
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index 10f62b7..3d5647f 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -8,6 +8,14 @@ #define JCV_FLT_MAX 1.7976931348623157E+308 #include <jc_voronoi.h> +double mod(double a, double b) { + if (a >= 0) { + return fmod(a, b); + } else { + return fmod(a + b * ceil(-a / b), b); + } +} + graph::graph(unsigned int Nx, unsigned int Ny) { L = {(double)Nx, (double)Ny}; @@ -84,45 +92,97 @@ class eulerException: public std::exception } } eulerex; +graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) { + L = {Lx, Ly}; -graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) { - L = {(double)Nx, (double)Ny}; + std::uniform_real_distribution<double> d(0.0, 1.0); + unsigned int N = round(Lx * Ly / 2 + d(rng) - 0.5); - unsigned int nv = Nx * Ny / 2; + unsigned int nv = N; vertices.resize(nv); - - std::normal_distribution<double> d(0.0, spread); // the coordinates of the lattice, from which a delaunay triangulation // and corresponding voronoi tessellation will be built for (unsigned int i = 0; i < nv; i++) { - double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng); - double ry = (double)(i / (Nx / 2)) + d(rng); - vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}}; - } - - // to make the resulting tessellation periodic, we tile eight (!) copies of - // the original 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 * nv); - for (const vertex &v : 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}); + vertices[i] = {{L.x * d(rng), L.y * d(rng)}}; } + double max_difference = std::numeric_limits<double>::max(); jcv_diagram diagram; memset(&diagram, 0, sizeof(jcv_diagram)); + std::vector<jcv_point> points(9 * nv); + + double rstep = sqrt(step); + + while (max_difference > pow(relax, 2) * 2 * N) { + double cur_diff = 0; + // to make the resulting tessellation periodic, we tile eight (!) copies of + // the original 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) + 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 + 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}; + } + + jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); + const jcv_site* sites = jcv_diagram_get_sites(&diagram); + + for (int i = 0; i < diagram.numsites; i++) { + const jcv_site* site = &sites[i]; + unsigned int ind = site->index; + if (ind % 9 == 0) { + double Cx = 0.0; + double Cy = 0.0; + unsigned int ne = 0; + + const jcv_graphedge* e = site->edges; + while (e) { + ne++; + Cx += e->pos[0].x; + Cy += e->pos[0].y; + e = e->next; + } + Cx /= ne; + Cy /= ne; + + double dx = Cx - vertices[ind / 9].r.x; + double dy = Cy - vertices[ind / 9].r.y; + + double dist = pow(dx, 2) + pow(dy, 2); + if (dist > cur_diff) { + cur_diff = dist; + } + vertices[ind / 9] = {{mod(vertices[ind / 9].r.x + rstep * dx, L.x), mod(vertices[ind / 9].r.y + rstep * dy, L.y)}}; + } + } + max_difference = cur_diff; + jcv_diagram_free(&diagram); + memset(&diagram, 0, sizeof(jcv_diagram)); + } + + 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 + 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}; + } jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); const jcv_site* sites = jcv_diagram_get_sites(&diagram); @@ -167,7 +227,7 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) if (it1 == known_vertices.end()) { vi1 = dual_vertices.size(); - dual_vertices.push_back({{fmod(L.x + e->pos[0].x, L.x), fmod(L.y + e->pos[0].y, L.y)},{{site->p.x, site->p.y}}}); + dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}}); known_vertices[t1] = vi1; } else { vi1 = it1->second; @@ -176,8 +236,8 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) if (i1 < i2) { edges.push_back({{i1, i2}, - {fmod(L.x + (site->p.x + neighbor->p.x) / 2, L.x), - fmod(L.y + (site->p.y + neighbor->p.y) / 2, L.y)} + {mod((site->p.x + neighbor->p.x) / 2, L.x), + mod((site->p.y + neighbor->p.y) / 2, L.y)} }); jcv_graphedge *en; @@ -197,15 +257,15 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) if (it2 == known_vertices.end()) { vi2 = dual_vertices.size(); - dual_vertices.push_back({{fmod(L.x + e->pos[1].x, L.x), fmod(L.y + e->pos[1].y, L.y)},{}}); + dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)},{}}); known_vertices[t2] = vi2; } else { vi2 = it2->second; } dual_edges.push_back({{vi1, vi2}, - {fmod(L.x + (e->pos[0].x + e->pos[1].x) / 2, L.x), - fmod(L.y + (e->pos[0].y + e->pos[1].y) / 2, L.y)} + {mod((e->pos[0].x + e->pos[1].x) / 2, L.x), + mod((e->pos[0].y + e->pos[1].y) / 2, L.y)} }); } |