From e4f9f49b06bb0505f3e1671a2ba9f5059233ab15 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 2 Jan 2019 17:43:19 -0500 Subject: more fixes to voronoi periodization, some upstream errors now caught and construction is retried --- lib/src/graph.cpp | 150 +++++++++++++++++++++++++++++++++++++++++------------- 1 file 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 // 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& lhs, const std::array& rhs) const - {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1]; - else return rhs[0] < rhs[1]; - } - }; - - std::unordered_map, unsigned int> known_vertices; + std::unordered_map, 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 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 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 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 t2 = {i1, i2, i3n, sig2}; + std::sort(t2.begin(), t2.begin() + 3); auto it2 = known_vertices.find(t2); -- cgit v1.2.3-70-g09d2