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