summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/src/graph.cpp150
1 files changed, 116 insertions, 34 deletions
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 7564fa7..7c50ed0 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -6,6 +6,7 @@
#define JCV_REAL_TYPE double
#define JCV_ATAN2 atan2
#define JCV_FLT_MAX 1.7976931348623157E+308
+#define NDEBUG
#include <jc_voronoi.h>
// actual mod for floats
@@ -72,6 +73,69 @@ class eulerException: public std::exception
}
} eulerex;
+class clippingException: public std::exception
+{
+ virtual const char* what() const throw()
+ {
+ return "An interior site has a clipped edge and the tesselation is malformed.";
+ }
+} clippingex;
+
+class triangleException: public std::exception
+{
+ virtual const char* what() const throw()
+ {
+ return "A dual-centered triangle has an impossible geometry and the tesselation is malformed.";
+ }
+} triex;
+
+unsigned int get_triangle_signature(unsigned int j1, unsigned int j2, unsigned int j3) {
+ // this yucky function takes three unsigned integers representing the
+ // location in the nine periodic copies of each corner of a delauney triangle
+ // and returns a signature for that triangle, which is an unsigned integer
+ // that uniquely labels the way the triangle crosses boundaries of the
+ // copies. This allows us to differentiate delauney triangles with identical
+ // vertices but which should be identified with different faces
+ unsigned int x1 = j1 % 3;
+ unsigned int y1 = j1 / 3;
+
+ unsigned int x2 = j2 % 3;
+ unsigned int y2 = j2 / 3;
+
+ unsigned int x3 = j3 % 3;
+ unsigned int y3 = j3 / 3;
+
+ if ((j1 == j2) && (j2 == j3)) {
+ return 0;
+ } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) || ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
+ return 1;
+ } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) || ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
+ return 2;
+ } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) || ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
+ return 3;
+ } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) || ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
+ return 4;
+ } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) || ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
+ return 5;
+ } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) || ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
+ return 6;
+ } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) || ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
+ return 7;
+ } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) || ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
+ return 8;
+ } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
+ return 9;
+ } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
+ return 10;
+ } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
+ return 11;
+ } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
+ return 12;
+ } else {
+ throw triex;
+ }
+}
+
graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) {
L = {Lx, Ly};
@@ -108,15 +172,15 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
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 + 0] = {v.r.x - L.x, v.r.y - L.y};
+ points[9 * i + 1] = {v.r.x + 0.0, v.r.y - L.y};
+ points[9 * i + 2] = {v.r.x + L.x, v.r.y - L.y};
+ points[9 * i + 3] = {v.r.x - L.x, v.r.y + 0.0};
+ points[9 * i + 4] = {v.r.x + 0.0, v.r.y + 0.0};
+ points[9 * i + 5] = {v.r.x + L.x, v.r.y + 0.0};
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};
+ points[9 * i + 7] = {v.r.x + 0.0, v.r.y + L.y};
+ points[9 * i + 8] = {v.r.x + L.x, v.r.y + L.y};
}
// run voronoi
@@ -127,7 +191,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
for (int i = 0; i < diagram.numsites; i++) {
const jcv_site* site = &sites[i];
unsigned int ind = site->index;
- if (ind % 9 == 0) {
+ if (ind % 9 == 4) {
double Cx = 0.0;
double Cy = 0.0;
double A = 0.0;
@@ -164,35 +228,30 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
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 + 0] = {v.r.x - L.x, v.r.y - L.y};
+ points[9 * i + 1] = {v.r.x + 0.0, v.r.y - L.y};
+ points[9 * i + 2] = {v.r.x + L.x, v.r.y - L.y};
+ points[9 * i + 3] = {v.r.x - L.x, v.r.y + 0.0};
+ points[9 * i + 4] = {v.r.x + 0.0, v.r.y + 0.0};
+ points[9 * i + 5] = {v.r.x + L.x, v.r.y + 0.0};
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};
+ points[9 * i + 7] = {v.r.x + 0.0, 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(), &bounds, &diagram);
const jcv_site* sites = jcv_diagram_get_sites(&diagram);
- struct paircomp {
- bool operator() (const std::array<unsigned int, 2>& lhs, const std::array<unsigned int, 2>& rhs) const
- {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1];
- else return rhs[0] < rhs[1];
- }
- };
-
- std::unordered_map<std::array<unsigned int, 3>, unsigned int> known_vertices;
+ std::unordered_map<std::array<unsigned int, 4>, unsigned int> 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) {
+ if (site->index % 9 == 4) {
+ bool self_bonded = false;
unsigned int i1 = (unsigned int)(site->index / 9);
+ unsigned int j1 = (unsigned int)(site->index % 9);
const jcv_graphedge* e = site->edges;
const jcv_graphedge* ep = site->edges;
while (ep->next) {
@@ -204,13 +263,26 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
// sites whose index is larger than ours are given edges, so we
// don't double up later
const jcv_site* neighbor = e->neighbor;
+ if (neighbor == NULL) {
+ throw clippingex;
+ }
unsigned int i2 = (unsigned int)(neighbor->index / 9);
+ unsigned int j2 = (unsigned int)(neighbor->index % 9);
+ unsigned int x2 = j2 % 3;
+ unsigned int y2 = j2 / 3;
vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y});
+ if (ep->neighbor == NULL) {
+ throw clippingex;
+ }
unsigned int i3p = (unsigned int)(ep->neighbor->index / 9);
- std::array<unsigned int, 3> t1 = {i1, i2, i3p};
- std::sort(t1.begin(), t1.end());
+ unsigned int j3p = (unsigned int)(ep->neighbor->index % 9);
+
+ unsigned int sig1 = get_triangle_signature(j1, j2, j3p);
+
+ std::array<unsigned int, 4> t1 = {i1, i2, i3p, sig1};
+ std::sort(t1.begin(), t1.begin() + 3);
auto it1 = known_vertices.find(t1);
@@ -225,10 +297,12 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
}
- if (i1 < i2) {
- unsigned int nn = neighbor->index % 9;
- bool crossed_x = (nn == 1) || (nn == 2) || (nn > 4);
- bool crossed_y = nn > 2;
+ if (i1 < i2 || (i1 == i2 && !self_bonded)) {
+ if (i1 == i2) {
+ self_bonded = true;
+ }
+ bool crossed_x = x2 != 1;
+ bool crossed_y = y2 != 1;
edges.push_back({{i1, i2},
{mod((site->p.x + neighbor->p.x) / 2, L.x),
mod((site->p.y + neighbor->p.y) / 2, L.y)},
@@ -242,9 +316,17 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
en = e->next;
}
+ if (en->neighbor == NULL) {
+ throw clippingex;
+ }
+
unsigned int i3n = (unsigned int)(en->neighbor->index / 9);
- std::array<unsigned int, 3> t2 = {i1, i2, i3n};
- std::sort(t2.begin(), t2.end());
+ unsigned int j3n = (unsigned int)(en->neighbor->index % 9);
+
+ unsigned int sig2 = get_triangle_signature(j1, j2, j3n);
+
+ std::array<unsigned int, 4> t2 = {i1, i2, i3n, sig2};
+ std::sort(t2.begin(), t2.begin() + 3);
auto it2 = known_vertices.find(t2);