diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/src/graph.cpp | 39 |
1 files changed, 27 insertions, 12 deletions
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index 3d5647f..b9d4f03 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -6,8 +6,12 @@ #define JCV_REAL_TYPE double #define JCV_ATAN2 atan2 #define JCV_FLT_MAX 1.7976931348623157E+308 +#undef assert +#define assert_error_report_helper(cond) "assertion failed: " +#define assert(cond) {if(!(cond)) { std::cerr << assert_error_report_helper(cond) "\n"; throw assert_error_report_helper(cond); } } #include <jc_voronoi.h> +// actual mod for floats double mod(double a, double b) { if (a >= 0) { return fmod(a, b); @@ -95,26 +99,30 @@ class eulerException: public std::exception graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) { L = {Lx, Ly}; + // randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with + // probability proportional to the distance from each std::uniform_real_distribution<double> d(0.0, 1.0); unsigned int N = round(Lx * Ly / 2 + d(rng) - 0.5); unsigned int nv = N; - vertices.resize(nv); // the coordinates of the lattice, from which a delaunay triangulation - // and corresponding voronoi tessellation will be built + // and corresponding voronoi tessellation will be built. Everyone is in the + // rectangle (0, 0) (Lx, Ly) for (unsigned int i = 0; i < nv; i++) { vertices[i] = {{L.x * d(rng), L.y * d(rng)}}; } - double max_difference = std::numeric_limits<double>::max(); + // set up the voronoi objects jcv_diagram diagram; memset(&diagram, 0, sizeof(jcv_diagram)); + jcv_rect bounds = {{-Lx, -Ly}, {2 * Lx, 2 * Ly}}; std::vector<jcv_point> points(9 * nv); double rstep = sqrt(step); + double max_difference = std::numeric_limits<double>::max(); while (max_difference > pow(relax, 2) * 2 * N) { double cur_diff = 0; // to make the resulting tessellation periodic, we tile eight (!) copies of @@ -123,6 +131,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) // 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}; @@ -134,27 +143,33 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y}; } - jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); + // run voronoi + jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram); + + // relax the sites by moving the vertices towards their centroids 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; + double A = 0.0; const jcv_graphedge* e = site->edges; while (e) { - ne++; - Cx += e->pos[0].x; - Cy += e->pos[0].y; + double dA = e->pos[0].x * e->pos[1].y - e->pos[1].x * e->pos[0].y; + + A += dA; + Cx += (e->pos[0].x + e->pos[1].x) * dA; + Cy += (e->pos[0].y + e->pos[1].y) * dA; + e = e->next; } - Cx /= ne; - Cy /= ne; + A /= 2; + Cx /= 6 * A; + Cy /= 6 * A; double dx = Cx - vertices[ind / 9].r.x; double dy = Cy - vertices[ind / 9].r.y; @@ -183,7 +198,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) 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); + jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram); const jcv_site* sites = jcv_diagram_get_sites(&diagram); |