1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
#include "graph.hpp"
#include <cstring>
#include <iostream>
#define JC_VORONOI_IMPLEMENTATION
#define JCV_REAL_TYPE double
#define JCV_ATAN2 atan2
#define JCV_FLT_MAX 1.7976931348623157E+308
#include <jc_voronoi.h>
graph::graph(unsigned int Nx, unsigned int Ny) {
L = {(double)Nx, (double)Ny};
unsigned int ne = Nx * Ny;
unsigned int nv = ne / 2;
vertices.resize(nv);
edges.reserve(ne);
dual_vertices.resize(nv);
dual_edges.reserve(ne);
for (unsigned int i = 0; i < nv; i++) {
vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
vertices[i].r.y = (double)(i / (Nx / 2));
dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
dual_vertices[i].r.y = (double)(i / (Nx / 2));
}
for (unsigned int x = 0; x < Ny; x++) {
for (unsigned int y = 0; y < Nx; y++) {
unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2);
unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv;
edges.push_back({v1, v2});
unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2);
unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv;
dual_edges.push_back({dv1, dv2});
}
}
}
graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
L = {(double)Nx, (double)Ny};
unsigned int dnv = Nx * Ny / 2;
dual_vertices.resize(dnv);
std::normal_distribution<double> d(0.0, spread);
// the coordinates of the dual lattice, from which a delaunay triangulation
// and corresponding voronoi tesselation will be built
for (unsigned int i = 0; i < dnv; i++) {
double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng);
double ry = (double)(i / (Nx / 2)) + d(rng);
dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
}
// to make the resulting tesselation periodic, we tile eight (!) copies of
// the original dual 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 * dnv);
for (const vertex &v : dual_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});
}
jcv_diagram diagram;
memset(&diagram, 0, sizeof(jcv_diagram));
jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram);
const jcv_site* sites = jcv_diagram_get_sites(&diagram);
struct coorcomp {
bool operator() (const coordinate& lhs, const coordinate& rhs) const
{return lhs.x<rhs.x;}
};
std::map<coordinate, unsigned int, coorcomp> 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) {
const jcv_graphedge* e = site->edges;
// for each edge on the site's cell
while(e) {
// assess whether the dual edge needs to be added. only neighboring
// dual sites whose index is larger than ours are given edges, so we
// don't double up later
const jcv_site* neighbor = e->neighbor;
unsigned int index = (unsigned int)(site->index / 9);
unsigned int real_index = (unsigned int)(neighbor->index / 9);
if (index < real_index) {
dual_edges.push_back({index, real_index});
coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))};
coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))};
std::map<coordinate, unsigned int>::iterator it1 = known_vertices.find(r1);
std::map<coordinate, unsigned int>::iterator it2 = known_vertices.find(r2);
unsigned int vi1, vi2;
if (it1 == known_vertices.end()) {
vi1 = vertices.size();
vertices.push_back({r1});
known_vertices[r1] = vi1;
} else {
vi1 = it1->second;
}
if (it2 == known_vertices.end()) {
vi2 = vertices.size();
vertices.push_back({r2});
known_vertices[r2] = vi2;
} else {
vi2 = it2->second;
}
edges.push_back({vi1, vi2});
}
e = e->next;
}
}
}
for (edge &e : edges) {
std::cout << e[0] << " " << e[1] << " ";
}
std::cout << "\n";
for (vertex &v : vertices) {
std::cout << v.r.x << " " << v.r.y << " ";
}
std::cout << "\n";
for (edge &e : dual_edges) {
std::cout << e[0] << " " << e[1] << " ";
}
std::cout << "\n";
for (vertex &v : dual_vertices) {
std::cout << v.r.x << " " << v.r.y << " ";
}
std::cout << "\n";
jcv_diagram_free(&diagram);
}
|
