Graphviz 13.0.0~dev.20250424.1043
Loading...
Searching...
No Matches
make_map.c
Go to the documentation of this file.
1/*************************************************************************
2 * Copyright (c) 2011 AT&T Intellectual Property
3 * All rights reserved. This program and the accompanying materials
4 * are made available under the terms of the Eclipse Public License v1.0
5 * which accompanies this distribution, and is available at
6 * https://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: Details at https://graphviz.org
9 *************************************************************************/
10
11#define STANDALONE
12#include <assert.h>
13#include <sparse/SparseMatrix.h>
14#include <sparse/general.h>
15#include <math.h>
16#include <sparse/QuadTree.h>
17#include <stdbool.h>
18#include <stddef.h>
19#include <string.h>
20#include <cgraph/cgraph.h>
21#include "make_map.h"
24#include <sparse/colorutil.h>
25#include <neatogen/delaunay.h>
26#include <util/agxbuf.h>
27#include <util/alloc.h>
28#include <util/debug.h>
29#include <util/list.h>
30#include <util/prisize_t.h>
31
32#include <edgepaint/lab.h>
34
35void map_palette_optimal_coloring(char *color_scheme, SparseMatrix A0,
36 float **rgb_r, float **rgb_g, float **rgb_b){
37 /*
38 for a graph A, get a distinctive color of its nodes so that the color distanmce among all nodes are maximized. Here
39 color distance on a node is defined as the minimum of color differences between a node and its neighbors.
40 color_scheme: rgb, gray, lab, or one of the color palettes in color_palettes.h, or a list of hex rgb colors separaterd by comma like "#ff0000,#00ff00"
41 A: the graph of n nodes
42 cdim: dimension of the color space
43 rgb_r, rgb_g, rgb_b: float array of length A->m + 1, which contains color for each country. 1-based
44 */
45
46 /*color: On input an array of size n*cdim, if NULL, will be allocated. On exit the final color assignment for node i is [cdim*i,cdim*(i+1)), in RGB (between 0 to 1)
47 */
48 double *colors = NULL;
49 int n = A0->m, i, cdim;
50
52 bool weightedQ = true;
53
54 {double *dist = NULL;
55 A = SparseMatrix_symmetrize(A0, false);
59 free(dist);
61 SparseMatrix_export(stdout, A);
62 }
63
64 // lightness: of the form 0,70, specifying the range of lightness of LAB
65 // color. Ignored if scheme is not COLOR_LAB.
66 int lightness[] = {0, 100};
67
68 // accuracy is the threshold given so that when finding the coloring for each
69 // node, the optimal is with in "accuracy" of the true global optimal.
70 const double accuracy = 0.01;
71
72 // seed: random_seed. If negative, consider -seed as the number of random
73 // start iterations
74 const int seed = -10;
75
76 node_distinct_coloring(color_scheme, lightness, weightedQ, A, accuracy, seed,
77 &cdim, &colors);
78
79 if (A != A0){
81 }
82 *rgb_r = gv_calloc(n + 1, sizeof(float));
83 *rgb_g = gv_calloc(n + 1, sizeof(float));
84 *rgb_b = gv_calloc(n + 1, sizeof(float));
85
86 for (i = 0; i < n; i++){
87 (*rgb_r)[i+1] = (float) colors[cdim*i];
88 (*rgb_g)[i+1] = (float) colors[cdim*i + 1];
89 (*rgb_b)[i+1] = (float) colors[cdim*i + 2];
90 }
91 free(colors);
92}
93
94void map_optimal_coloring(int seed, SparseMatrix A, float *rgb_r, float *rgb_g, float *rgb_b){
95 int *p = NULL;
96 float *u = NULL;
97 int n = A->m;
98 int i;
99
101
102 rgb_r++; rgb_b++; rgb_g++;/* seems necessary, but need to better think about cases when clusters are not contiguous */
103 vector_float_take(n, rgb_r, n, p, &u);
104 for (i = 0; i < n; i++) rgb_r[i] = u[i];
105 vector_float_take(n, rgb_g, n, p, &u);
106 for (i = 0; i < n; i++) rgb_g[i] = u[i];
107 vector_float_take(n, rgb_b, n, p, &u);
108 for (i = 0; i < n; i++) rgb_b[i] = u[i];
109 free(u);
110}
111
112static int get_poly_id(int ip, SparseMatrix point_poly_map){
113 return point_poly_map->ja[point_poly_map->ia[ip]];
114}
115
116void improve_contiguity(int n, int dim, int *grouping, SparseMatrix poly_point_map, double *x, SparseMatrix graph){
117 /*
118 grouping: which group each of the vertex belongs to
119 poly_point_map: a matrix of dimension npolys x (n + nrandom), poly_point_map[i,j] != 0 if polygon i contains the point j.
120 . If j < n, it is the original point, otherwise it is artificial point (forming the rectangle around a label) or random points.
121 */
122 int i, j, *ia, *ja, u, v;
123 SparseMatrix point_poly_map, D;
124 double dist;
125 int nbad = 0, flag;
126 int maxit = 10;
127
129
130 assert(graph->m == n);
131 ia = D->ia; ja = D->ja;
132 double *a = D->a;
133
134 /* point_poly_map: each row i has only 1 entry at column j, which says that point i is in polygon j */
135 point_poly_map = SparseMatrix_transpose(poly_point_map);
136
137 for (i = 0; i < n; i++){
138 u = i;
139 for (j = ia[i]; j < ia[i+1]; j++){
140 v = ja[j];
141 dist = distance_cropped(x, dim, u, v);
142 if (grouping[u] != grouping[v]){
143 a[j] = 1.1*dist;
144 } else if (get_poly_id(u, point_poly_map) == get_poly_id(v, point_poly_map)){
145 a[j] = dist;
146 } else {
147 nbad++;
148 a[j] = 0.9*dist;
149 }
150
151 }
152 }
153
154 GV_INFO("ratio (edges among discontiguous regions vs total edges)=%f", (double)nbad / ia[n]);
155 stress_model(dim, D, &x, maxit, &flag);
156
157 assert(!flag);
158
160 SparseMatrix_delete(point_poly_map);
161}
162
163struct Triangle {
164 int vertices[3];/* 3 points */
165 double center[2]; /* center of the triangle */
166};
167
168static void normal(double v[], double normal[]){
169 if (v[0] == 0){
170 normal[0] = 1; normal[1] = 0;
171 } else {
172 normal[0] = -v[1];
173 normal[1] = v[0];
174 }
175}
176
177static void triangle_center(double x[], double y[], double z[], double c[]){
178 /* find the "center" c, which is the intersection of the 3 vectors that are normal to each
179 of the edges respectively, and which passes through the center of the edges respectively
180 center[{x_, y_, z_}] := Module[
181 {xy = 0.5*(x + y), yz = 0.5*(y + z), zx = 0.5*(z + x), nxy, nyz,
182 beta, cen},
183 nxy = normal[y - x];
184 nyz = normal[y - z];
185 beta = (y-x).(xy - yz)/(nyz.(y-x));
186 cen = yz + beta*nyz;
187 Graphics[{Line[{x, y, z, x}], Red, Point[cen], Line[{cen, xy}],
188 Line[{cen, yz}], Green, Line[{cen, zx}]}]
189
190 ]
191 */
192 double xy[2], yz[2], nxy[2], nyz[2], ymx[2], ymz[2], beta, bot;
193 int i;
194
195 for (i = 0; i < 2; i++) ymx[i] = y[i] - x[i];
196 for (i = 0; i < 2; i++) ymz[i] = y[i] - z[i];
197 for (i = 0; i < 2; i++) xy[i] = 0.5*(x[i] + y[i]);
198 for (i = 0; i < 2; i++) yz[i] = 0.5*(y[i] + z[i]);
199
200
201 normal(ymx, nxy);
202 normal(ymz, nyz);
203 bot = nyz[0]*(x[0]-y[0])+nyz[1]*(x[1]-y[1]);
204 if (bot == 0){/* xy and yz are parallel */
205 c[0] = xy[0]; c[1] = xy[1];
206 return;
207 }
208 beta = ((x[0] - y[0])*(xy[0] - yz[0])+(x[1] - y[1])*(xy[1] - yz[1]))/bot;
209 c[0] = yz[0] + beta*nyz[0];
210 c[1] = yz[1] + beta*nyz[1];
211}
212
213static SparseMatrix matrix_add_entry(SparseMatrix A, int i, int j, int val){
214 int i1 = i, j1 = j;
215 if (i < j) {
216 i1 = j; j1 = i;
217 }
219 return SparseMatrix_coordinate_form_add_entry(A, i1, j1, &val);
220}
221
222static void plot_dot_edges(FILE *f, SparseMatrix A){
223 int i, *ia, *ja, j;
224
225
226 int n = A->m;
227 ia = A->ia;
228 ja = A->ja;
229 for (i = 0; i < n; i++){
230 for (j = ia[i]; j < ia[i+1]; j++){
231 if (ja[j] == i) continue;
232 fprintf(f,"%d -- %d;\n",i,ja[j]);
233 }
234 }
235}
236
237static void plot_dot_labels(FILE *f, int n, int dim, double *x, char **labels, float *fsz){
238 int i;
239
240 for (i = 0; i < n; i++){
241 if (fsz){
242 fprintf(f, "%d [label=\"%s\", pos=\"%lf,%lf\", fontsize=%f];\n",i, labels[i], x[i*dim], x[i*dim+1], fsz[i]);
243 } else {
244 fprintf(f, "%d [label=\"%s\", pos=\"%lf,%lf\"];\n",i, labels[i], x[i*dim], x[i*dim+1]);
245 }
246 }
247
248}
249
250DEFINE_LIST(doubles, double)
251
252static void dot_polygon(agxbuf *sbuff, doubles_t xp, doubles_t yp,
253 double line_width, bool fill, const char *cstring) {
254
255 assert(doubles_size(&xp) == doubles_size(&yp));
256 if (!doubles_is_empty(&xp)){
257 if (fill) {
258 agxbprint(sbuff,
259 " c %" PRISIZE_T " -%s C %" PRISIZE_T " -%s P %" PRISIZE_T " ",
260 strlen(cstring), cstring, strlen(cstring), cstring,
261 doubles_size(&xp));
262 } else {
263 if (line_width > 0){
264 size_t len_swidth = (size_t)snprintf(NULL, 0, "%f", line_width);
265 agxbprint(sbuff, " c %" PRISIZE_T " -%s S %" PRISIZE_T
266 " -setlinewidth(%f) L %" PRISIZE_T " ", strlen(cstring), cstring,
267 len_swidth + 14, line_width, doubles_size(&xp));
268 } else {
269 agxbprint(sbuff, " c %" PRISIZE_T " -%s L %" PRISIZE_T " ", strlen(cstring),
270 cstring, doubles_size(&xp));
271 }
272 }
273 for (size_t i = 0; i < doubles_size(&xp); i++) {
274 agxbprint(sbuff, " %f %f", doubles_get(&xp, i), doubles_get(&yp, i));
275 }
276 }
277}
278
279static void plot_dot_polygons(agxbuf *sbuff, double line_width,
280 const char *line_color, SparseMatrix polys,
281 double *x_poly, int *polys_groups, float *r,
282 float *g, float *b, const char *opacity) {
283 int i, j, *ia = polys->ia, *ja = polys->ja, *a = polys->a, npolys = polys->m, nverts = polys->n, ipoly,first;
284 const bool fill = false;
285 const bool use_line = line_width >= 0;
286
287 agxbuf cstring_buffer = {0};
288 agxbput(&cstring_buffer, "#aaaaaaff");
289 char *cstring = agxbuse(&cstring_buffer);
290
291 doubles_t xp = {0};
292 doubles_t yp = {0};
293
294 GV_INFO("npolys = %d", npolys);
295 first = abs(a[0]); ipoly = first + 1;
296 for (i = 0; i < npolys; i++){
297 for (j = ia[i]; j < ia[i+1]; j++){
298 assert(ja[j] < nverts && ja[j] >= 0);
299 (void)nverts;
300 if (abs(a[j]) != ipoly){/* the first poly, or a hole */
301 ipoly = abs(a[j]);
302 if (r && g && b) {
303 rgb2hex(r[polys_groups[i]], g[polys_groups[i]], b[polys_groups[i]],
304 &cstring_buffer, opacity);
305 cstring = agxbuse(&cstring_buffer);
306 }
307 dot_polygon(sbuff, xp, yp, line_width, fill, cstring);
308 // start a new polygon
309 doubles_clear(&xp);
310 doubles_clear(&yp);
311 }
312 doubles_append(&xp, x_poly[2 * ja[j]]);
313 doubles_append(&yp, x_poly[2 * ja[j] + 1]);
314 }
315 if (use_line) {
316 dot_polygon(sbuff, xp, yp, line_width, fill, line_color);
317 } else {
318 /* why set fill to polys_groups[i]?*/
319 dot_polygon(sbuff, xp, yp, -1, true, cstring);
320 }
321 }
322 agxbfree(&cstring_buffer);
323 doubles_free(&xp);
324 doubles_free(&yp);
325}
326
327void plot_dot_map(Agraph_t* gr, int n, int dim, double *x, SparseMatrix polys,
328 SparseMatrix poly_lines, double line_width,
329 const char *line_color, double *x_poly, int *polys_groups,
330 char **labels, float *fsz, float *r, float *g, float *b,
331 const char* opacity, SparseMatrix A, FILE* f) {
332 /* if graph object exist, we just modify some attributes, otherwise we dump the whole graph */
333 bool plot_polyQ = true;
334 agxbuf sbuff = {0};
335
336 if (!r || !g || !b) plot_polyQ = false;
337
338 if (!gr) {
339 fprintf(f, "graph map {\n node [margin = 0 width=0.0001 height=0.00001 shape=plaintext];\n graph [outputorder=edgesfirst, bgcolor=\"#dae2ff\"]\n edge [color=\"#55555515\",fontname=\"Helvetica-Bold\"]\n");
340 } else {
341 agattr_text(gr, AGNODE, "margin", "0");
342 agattr_text(gr, AGNODE, "width", "0.0001");
343 agattr_text(gr, AGNODE, "height", "0.0001");
344 agattr_text(gr, AGNODE, "shape", "plaintext");
345 agattr_text(gr, AGNODE, "margin", "0");
346 agattr_text(gr, AGNODE, "fontname", "Helvetica-Bold");
347 agattr_text(gr, AGRAPH, "outputorder", "edgesfirst");
348 agattr_text(gr, AGRAPH, "bgcolor", "#dae2ff");
349 if (!A) agattr_text(gr, AGEDGE, "style","invis");/* do not plot edges */
350 }
351
352 /*polygons */
353 if (plot_polyQ) {
354 if (!gr) fprintf(f,"_background = \"");
355 plot_dot_polygons(&sbuff, -1., NULL, polys, x_poly, polys_groups, r, g, b, opacity);
356 }
357
358 /* polylines: line width is set here */
359 if (line_width >= 0){
360 plot_dot_polygons(&sbuff, line_width, line_color, poly_lines, x_poly, polys_groups, NULL, NULL, NULL, NULL);
361 }
362 if (!gr) {
363 fprintf(f,"%s",agxbuse(&sbuff));
364 fprintf(f,"\"\n");/* close polygons/lines */
365 } else {
366 agattr_text(gr, AGRAPH, "_background", agxbuse(&sbuff));
367 agwrite(gr, f);
368 }
369
370 /* nodes */
371 if (!gr && labels) plot_dot_labels(f, n, dim, x, labels, fsz);
372 /* edges */
373 if (!gr && A) plot_dot_edges(f, A);
374
375 /* background color + plot label?*/
376
377 if (!gr) fprintf(f, "}\n");
378
379 agxbfree(&sbuff);
380}
381
396static int get_tri(int n, int dim, double *x, int *nt, struct Triangle **T,
397 SparseMatrix *E) {
398 int i, j, i0, i1, i2, ntri;
400
401 int* trilist = get_triangles(x, n, &ntri);
402 if (trilist == NULL) {
403 return -1;
404 }
405
406 *T = gv_calloc(ntri, sizeof(struct Triangle));
407
409 for (i = 0; i < ntri; i++) {
410 for (j = 0; j < 3; j++) {
411 (*T)[i].vertices[j] = trilist[i * 3 + j];
412 }
413 i0 = (*T)[i].vertices[0]; i1 = (*T)[i].vertices[1]; i2 = (*T)[i].vertices[2];
414
415 triangle_center(&x[i0*dim], &x[i1*dim], &x[i2*dim], (*T)[i].center);
416 A = matrix_add_entry(A, i0, i1, i);
417 A = matrix_add_entry(A, i1, i2, i);
418 A = matrix_add_entry(A, i2, i0, i);
419 }
420
424 *E = B;
425
426 *nt = ntri;
427
428 free(trilist);
429 return 0;
430}
431
432static SparseMatrix get_country_graph(int n, SparseMatrix A, int *groups, int GRP_RANDOM, int GRP_BBOX){
433 /* form a graph each vertex is a group (a country), and a vertex is connected to another if the two countries shares borders.
434 since the group ID may not be contiguous (e.g., only groups 2,3,5, -1), we will return NULL if one of the group has non-positive ID! */
435 int *ia, *ja;
436 int one = 1, jj, i, j, ig1, ig2;
437 SparseMatrix B, BB;
438 int min_grp, max_grp;
439
440 min_grp = max_grp = groups[0];
441 for (i = 0; i < n; i++) {
442 max_grp = MAX(groups[i], max_grp);
443 min_grp = MIN(groups[i], min_grp);
444 }
445 if (min_grp <= 0) return NULL;
446 B = SparseMatrix_new(max_grp, max_grp, 1, MATRIX_TYPE_INTEGER, FORMAT_COORD);
447 ia = A->ia;
448 ja = A->ja;
449 for (i = 0; i < n; i++){
450 ig1 = groups[i]-1;/* add a diagonal entry */
452 for (j = ia[i]; j < ia[i+1]; j++){
453 jj = ja[j];
454 if (i != jj && groups[i] != groups[jj] && groups[jj] != GRP_RANDOM && groups[jj] != GRP_BBOX){
455 ig1 = groups[i]-1; ig2 = groups[jj]-1;
457 }
458 }
459 }
462 return BB;
463}
464
465static void conn_comp(int n, SparseMatrix A, int *groups, SparseMatrix *poly_point_map){
466 /* form a graph where only vertices that are connected as well as in the same group are connected */
467 int *ia, *ja;
468 int one = 1, jj, i, j;
469 SparseMatrix B, BB;
470 int ncomps, *comps = NULL;
471
473 ia = A->ia;
474 ja = A->ja;
475 for (i = 0; i < n; i++){
476 for (j = ia[i]; j < ia[i+1]; j++){
477 jj = ja[j];
478 if (i != jj && groups[i] == groups[jj]){
480 }
481 }
482 }
484
485 int *comps_ptr = SparseMatrix_weakly_connected_components(BB, &ncomps, &comps);
488 *poly_point_map = SparseMatrix_new(ncomps, n, n, MATRIX_TYPE_PATTERN, FORMAT_CSR);
489 free((*poly_point_map)->ia);
490 free((*poly_point_map)->ja);
491 (*poly_point_map)->ia = comps_ptr;
492 (*poly_point_map)->ja = comps;
493 (*poly_point_map)->nz = n;
494
495}
496
497static void get_poly_lines(int nt, SparseMatrix E, int ncomps, int *comps_ptr,
498 int *comps, int *groups, SparseMatrix *poly_lines,
499 int **polys_groups, int GRP_RANDOM, int GRP_BBOX) {
500 /*============================================================
501
502 polygon outlines
503
504 ============================================================*/
505 int i, *tlist, nz, ipoly, nnt, ii, jj, t1, t2, t, cur, next, nn, j, nlink, sta;
506 int *elist, edim = 3;/* a list tell which vertex a particular vertex is linked with during poly construction.
507 since the surface is a cycle, each can only link with 2 others, the 3rd position is used to record how many links
508 */
509 int *ie = E->ia, *je = E->ja, *e = E->a;
511
512 int *mask = gv_calloc(nt, sizeof(int));
513 for (i = 0; i < nt; i++) mask[i] = -1;
514 /* loop over every point in each connected component */
515 elist = gv_calloc(nt * edim, sizeof(int));
516 tlist = gv_calloc(nt * 2, sizeof(int));
517 *poly_lines = SparseMatrix_new(ncomps, nt, 1, MATRIX_TYPE_INTEGER, FORMAT_COORD);
518 *polys_groups = gv_calloc(ncomps, sizeof(int));
519
520 for (i = 0; i < nt; i++) elist[i*edim + 2] = 0;
521 nz = ie[E->m] - ie[0];
522
523 ipoly = 1;
524
525 for (i = 0; i < ncomps; i++){
526 nnt = 0;
527 for (j = comps_ptr[i]; j < comps_ptr[i+1]; j++){
528 ii = comps[j];
529
530 (*polys_groups)[i] = groups[ii];/* assign the grouping of each poly */
531
532 /* skip the country formed by random points */
533 if (groups[ii] == GRP_RANDOM || groups[ii] == GRP_BBOX) continue;
534
535 for (jj = ie[ii]; jj < ie[ii+1]; jj++){
536 if (groups[je[jj]] != groups[ii] && jj < nz - 1 && je[jj] == je[jj+1]){/* an triangle edge neighboring 2 triangles and two ends not in the same groups */
537 t1 = e[jj];
538 t2 = e[jj+1];
539
540 nlink = elist[t1*edim + 2]%2;
541 elist[t1*edim + nlink] = t2;/* t1->t2*/
542 elist[t1*edim + 2]++;
543
544 nlink = elist[t2*edim + 2]%2;
545 elist[t2*edim + nlink] = t1;/* t1->t2*/
546 elist[t2*edim + 2]++;
547
548 tlist[nnt++] = t1; tlist[nnt++] = t2;
549 jj++;
550 }
551 }
552 }/* done poly edges for this component i */
553
554 /* form one or more (if there is a hole) polygon outlines for this component */
555 for (j = 0; j < nnt; j++){
556 t = tlist[j];
557 if (mask[t] != i){
558 cur = sta = t; mask[cur] = i;
559 next = neighbor(t, 1, edim, elist);
560 SparseMatrix_coordinate_form_add_entry(*poly_lines, i, cur, &ipoly);
561 while (next != sta){
562 mask[next] = i;
563
564 SparseMatrix_coordinate_form_add_entry(*poly_lines, i, next, &ipoly);
565
566 nn = neighbor(next, 0, edim, elist);
567 if (nn == cur) {
568 nn = neighbor(next, 1, edim, elist);
569 }
570 assert(nn != cur);
571
572 cur = next;
573 next = nn;
574 }
575
576 SparseMatrix_coordinate_form_add_entry(*poly_lines, i, sta, &ipoly);/* complete a cycle by adding starting point */
577
578 ipoly++;
579 }
580
581 }/* found poly_lines for this comp */
582 }
583
585 SparseMatrix_delete(*poly_lines);
586 *poly_lines = A;
587
588 free(tlist);
589 free(elist);
590 free(mask);
591}
592
593static void cycle_print(int head, int *cycle, int *edge_table){
594 int cur, next;
595
596 cur = head;
597 fprintf(stderr, "cycle (edges): {");
598 while ((next = cycle_next(cur)) != head){
599 fprintf(stderr, "%d,",cur);
600 cur = next;
601 }
602 fprintf(stderr, "%d}\n",cur);
603
604 cur = head;
605 fprintf(stderr, "cycle (vertices): ");
606 while ((next = cycle_next(cur)) != head){
607 fprintf(stderr, "%d--",edge_head(cur));
608 cur = next;
609 }
610 fprintf(stderr, "%d--%d\n",edge_head(cur),edge_tail(cur));
611}
612
613static int same_edge(int ecur, int elast, int *edge_table){
614 return (edge_head(ecur) == edge_head(elast) && edge_tail(ecur) == edge_tail(elast))
615 || (edge_head(ecur) == edge_tail(elast) && edge_tail(ecur) == edge_head(elast));
616}
617
618static void get_polygon_solids(int nt, SparseMatrix E, int ncomps,
619 int *comps_ptr, int *comps, SparseMatrix *polys)
620{
621 /*============================================================
622
623 polygon solids that will be colored
624
625 ============================================================*/
626 int *edge_table;/* a table of edges of the triangle graph. If two vertex u and v are connected and are adjacent to two triangles
627 t1 and t2, then from u there are two edges to v, one denoted as t1->t2, and the other t2->t1. They are
628 numbered as e1 and e2. edge_table[e1]={t1,t2} and edge_table[e2]={t2,t1}
629 */
630 SparseMatrix half_edges;/* a graph of triangle edges. If two vertex u and v are connected and are adjacent to two triangles
631 t1 and t2, then from u there are two edges to v, one denoted as t1->t2, and the other t2->t1. They are
632 numbered as e1 and e2. Likewise from v to u there are also two edges e1 and e2.
633 */
634
635 int n = E->m, *ie = E->ia, *je = E->ja, *e = E->a, ne, i, j, t1, t2, jj, ii;
636 int *cycle, cycle_head = 0;/* a list of edges that form a cycle that describe the polygon. cycle[e][0] gives the prev edge in the cycle from e,
637 cycle[e][1] gives the next edge
638 */
639 int *edge_cycle_map, NOT_ON_CYCLE = -1;/* map an edge e to its position on cycle, unless it does not exist (NOT_ON_CYCLE) */
640 int *emask;/* whether an edge is seen this iter */
641 enum {NO_DUPLICATE = -1};
642 int *elist, edim = 3;/* a list tell which edge a particular vertex is linked with when a voro cell has been visited,
643 since the surface is a cycle, each vertex can only link with 2 edges, the 3rd position is used to record how many links
644 */
645
646 int k, duplicate, ee = 0, ecur, enext, eprev, cur, next, nn, nlink, head, elast = 0, etail, tail, ehead, efirst;
647
648 int DEBUG_CYCLE = 0;
650
651 ne = E->nz;
652 edge_table = gv_calloc(ne * 2, sizeof(int));
653
654 half_edges = SparseMatrix_new(n, n, 1, MATRIX_TYPE_INTEGER, FORMAT_COORD);
655
656 ne = 0;
657 for (i = 0; i < n; i++){
658 for (j = ie[i]; j < ie[i+1]; j++){
659 if (j < ie[n] - ie[0] - 1 && i > je[j] && je[j] == je[j+1]){/* an triangle edge neighboring 2 triangles. Since E is symmetric, we only do one edge of E*/
660 t1 = e[j];
661 t2 = e[j+1];
662 jj = je[j];
663 assert(jj < n);
664 edge_table[ne*2] = t1;/*t1->t2*/
665 edge_table[ne*2+1] = t2;
666 half_edges = SparseMatrix_coordinate_form_add_entry(half_edges, i, jj, &ne);
667 half_edges = SparseMatrix_coordinate_form_add_entry(half_edges, jj, i, &ne);
668 ne++;
669
670 edge_table[ne*2] = t2;/*t2->t1*/
671 edge_table[ne*2+1] = t1;
672 half_edges = SparseMatrix_coordinate_form_add_entry(half_edges, i, jj, &ne);
673 half_edges = SparseMatrix_coordinate_form_add_entry(half_edges, jj, i, &ne);
674
675
676 ne++;
677 j++;
678 }
679 }
680 }
681 assert(E->nz >= ne);
682
683 cycle = gv_calloc(ne * 2, sizeof(int));
685 SparseMatrix_delete(half_edges);half_edges = B;
686
687 edge_cycle_map = gv_calloc(ne, sizeof(int));
688 emask = gv_calloc(ne, sizeof(int));
689 for (i = 0; i < ne; i++) edge_cycle_map[i] = NOT_ON_CYCLE;
690 for (i = 0; i < ne; i++) emask[i] = -1;
691
692 ie = half_edges->ia;
693 je = half_edges->ja;
694 e = half_edges->a;
695 elist = gv_calloc(nt * 3, sizeof(int));
696 for (i = 0; i < nt; i++) elist[i*edim + 2] = 0;
697
698 *polys = SparseMatrix_new(ncomps, nt, 1, MATRIX_TYPE_INTEGER, FORMAT_COORD);
699
700 for (i = 0; i < ncomps; i++){
701 if (DEBUG_CYCLE) fprintf(stderr, "\n ============ comp %d has %d members\n",i, comps_ptr[i+1]-comps_ptr[i]);
702 for (k = comps_ptr[i]; k < comps_ptr[i+1]; k++){
703 ii = comps[k];
704 duplicate = NO_DUPLICATE;
705 if (DEBUG_CYCLE) fprintf(stderr,"member = %d has %d neighbors\n",ii, ie[ii+1]-ie[ii]);
706 for (j = ie[ii]; j < ie[ii+1]; j++){
707 jj = je[j];
708 ee = e[j];
709 t1 = edge_head(ee);
710 if (DEBUG_CYCLE) fprintf(stderr," linked with %d using half-edge %d, {head,tail} of the edge = {%d, %d}\n",jj, ee, t1, edge_tail(ee));
711 nlink = elist[t1*edim + 2]%2;
712 elist[t1*edim + nlink] = ee;/* t1->t2*/
713 elist[t1*edim + 2]++;
714
715 if (edge_cycle_map[ee] != NOT_ON_CYCLE) duplicate = ee;
716 emask[ee] = ii;
717 }
718
719 if (duplicate == NO_DUPLICATE){
720 /* this must be the first time the cycle is being established, a new voro cell*/
721 ecur = ee;
722 cycle_head = ecur;
723 cycle_next(ecur) = ecur;
724 cycle_prev(ecur) = ecur;
725 edge_cycle_map[ecur] = 1;
726 head = cur = edge_head(ecur);
727 next = edge_tail(ecur);
728 if (DEBUG_CYCLE) fprintf(stderr, "NEW CYCLE\n starting with edge %d, {head,tail}={%d,%d}\n", ee, head, next);
729 while (next != head){
730 enext = neighbor(next, 0, edim, elist);/* two voro edges linked with triangle "next" */
731 if ((edge_head(enext) == cur && edge_tail(enext) == next)
732 || (edge_head(enext) == next && edge_tail(enext) == cur)){/* same edge */
733 enext = neighbor(next, 1, edim, elist);
734 };
735 if (DEBUG_CYCLE) fprintf(stderr, "cur edge = %d, next edge %d, {head,tail}={%d,%d},\n",ecur, enext, edge_head(enext), edge_tail(enext));
736 nn = edge_head(enext);
737 if (nn == next) nn = edge_tail(enext);
738 cycle_next(enext) = cycle_next(ecur);
739 cycle_prev(enext) = ecur;
740 cycle_next(ecur) = enext;
741 cycle_prev(ee) = enext;
742 edge_cycle_map[enext] = 1;
743
744 ecur = enext;
745 cur = next;
746 next = nn;
747 }
748 if (DEBUG_CYCLE) cycle_print(ee, cycle,edge_table);
749 } else {
750 /* we found a duplicate edge, remove that, and all contiguous neighbors that overlap with the current voro
751 */
752 ecur = ee = duplicate;
753 while (emask[ecur] == ii){
754 /* contiguous overlapping edges, Cycling is not possible
755 since the cycle can not complete surround the new voro cell and yet
756 do not contain any other edges
757 */
758 ecur = cycle_next(ecur);
759 }
760 if (DEBUG_CYCLE) fprintf(stderr," duplicating edge = %d, starting from the a non-duplicating edge %d, search backwards\n",ee, ecur);
761
762 ecur = cycle_prev(ecur);
763 efirst = ecur;
764 while (emask[ecur] == ii){
765 if (DEBUG_CYCLE) fprintf(stderr," remove edge %d (%d--%d)\n",ecur, edge_head(ecur), edge_tail(ecur));
766 /* short this duplicating edge */
767 edge_cycle_map[ecur] = NOT_ON_CYCLE;
768 enext = cycle_next(ecur);
769 eprev = cycle_prev(ecur);
770 cycle_next(ecur) = ecur;/* isolate this edge */
771 cycle_prev(ecur) = ecur;
772 cycle_next(eprev) = enext;/* short */
773 cycle_prev(enext) = eprev;
774 elast = ecur;/* record the last removed edge */
775 ecur = eprev;
776 }
777
778 if (DEBUG_CYCLE) {
779 fprintf(stderr, "remaining (broken) cycle = ");
780 cycle_print(cycle_next(ecur), cycle,edge_table);
781 }
782
783 /* we now have a broken cycle of head = edge_tail(ecur) and tail = edge_head(cycle_next(ecur)) */
784 ehead = ecur; etail = cycle_next(ecur);
785 cycle_head = ehead;
786 head = edge_tail(ehead);
787 tail = edge_head(etail);
788
789 /* pick an edge ev from head in the voro that is a removed edge: since the removed edges form a path starting from
790 efirst, and at elast (head of elast is head), usually we just need to check that ev is not the same as elast,
791 but in the case of a voro filling in a hole, we also need to check that ev is not efirst,
792 since in this case every edge of the voro cell is removed
793 */
794 ecur = neighbor(head, 0, edim, elist);
795 if (same_edge(ecur, elast, edge_table)){
796 ecur = neighbor(head, 1, edim, elist);
797 };
798
799 if (DEBUG_CYCLE) fprintf(stderr, "forwarding now from edge %d = {%d, %d}, try to reach vtx %d, first edge from voro = %d\n",
800 ehead, edge_head(ehead), edge_tail(ehead), tail, ecur);
801
802 /* now go along voro edges till we reach the tail of the broken cycle*/
803 cycle_next(ehead) = ecur;
804 cycle_prev(ecur) = ehead;
805 cycle_prev(etail) = ecur;
806 cycle_next(ecur) = etail;
807 if (same_edge(ecur, efirst, edge_table)){
808 if (DEBUG_CYCLE) fprintf(stderr, "this voro cell fill in a hole completely!!!!\n");
809 } else {
810
811 edge_cycle_map[ecur] = 1;
812 head = cur = edge_head(ecur);
813 next = edge_tail(ecur);
814 if (DEBUG_CYCLE) fprintf(stderr, "starting with edge %d, {head,tail}={%d,%d}\n", ecur, head, next);
815 while (next != tail){
816 enext = neighbor(next, 0, edim, elist);/* two voro edges linked with triangle "next" */
817 if ((edge_head(enext) == cur && edge_tail(enext) == next)
818 || (edge_head(enext) == next && edge_tail(enext) == cur)){/* same edge */
819 enext = neighbor(next, 1, edim, elist);
820 };
821 if (DEBUG_CYCLE) fprintf(stderr, "cur edge = %d, next edge %d, {head,tail}={%d,%d},\n",ecur, enext, edge_head(enext), edge_tail(enext));
822
823
824 nn = edge_head(enext);
825 if (nn == next) nn = edge_tail(enext);
826 cycle_next(enext) = cycle_next(ecur);
827 cycle_prev(enext) = ecur;
828 cycle_next(ecur) = enext;
829 cycle_prev(etail) = enext;
830 edge_cycle_map[enext] = 1;
831
832 ecur = enext;
833 cur = next;
834 next = nn;
835 }
836 }
837
838 }
839
840 }
841 /* done this component, load to sparse matrix, unset edge_map*/
842 ecur = cycle_head;
843 while ((enext = cycle_next(ecur)) != cycle_head){
844 edge_cycle_map[ecur] = NOT_ON_CYCLE;
845 head = edge_head(ecur);
847 ecur = enext;
848 }
849 edge_cycle_map[ecur] = NOT_ON_CYCLE;
850 head = edge_head(ecur); tail = edge_tail(ecur);
852 SparseMatrix_coordinate_form_add_entry(*polys, i, tail, &i);
853
854
855 /* unset edge_map */
856 }
857
859 SparseMatrix_delete(*polys);
860 *polys = B;
861
862 SparseMatrix_delete(half_edges);
863 free(cycle);
864 free(edge_cycle_map);
865 free(elist);
866 free(emask);
867 free(edge_table);
868}
869
870static void get_polygons(int n, int nrandom, int dim, int *grouping, int nt,
871 struct Triangle *Tp, SparseMatrix E, int *nverts,
872 double **x_poly, SparseMatrix *poly_lines,
873 SparseMatrix *polys, int **polys_groups,
874 SparseMatrix *poly_point_map,
875 SparseMatrix *country_graph) {
876 int i, j;
877 int *groups;
878 int maxgrp;
879 int *comps = NULL, *comps_ptr = NULL, ncomps;
880 int GRP_RANDOM, GRP_BBOX;
881
882 assert(dim == 2);
883 *nverts = nt;
884
885 groups = gv_calloc(n + nrandom, sizeof(int));
886 maxgrp = grouping[0];
887 for (i = 0; i < n; i++) {
888 maxgrp = MAX(maxgrp, grouping[i]);
889 groups[i] = grouping[i];
890 }
891
892 GRP_RANDOM = maxgrp + 1; GRP_BBOX = maxgrp + 2;
893 for (i = n; i < n + nrandom - 4; i++) {/* all random points in the same group */
894 groups[i] = GRP_RANDOM;
895 }
896 for (i = n + nrandom - 4; i < n + nrandom; i++) {/* last 4 pts of the expanded bonding box in the same group */
897 groups[i] = GRP_BBOX;
898 }
899
900 /* finding connected components: vertices that are connected in the triangle graph, as well as in the same group */
901 conn_comp(n + nrandom, E, groups, poly_point_map);
902
903 ncomps = (*poly_point_map)->m;
904 comps = (*poly_point_map)->ja;
905 comps_ptr = (*poly_point_map)->ia;
906
907 /* connected components are such that the random points and the bounding box 4 points forms the last
908 remaining components */
909 for (i = ncomps - 1; i >= 0; i--) {
910 if (groups[comps[comps_ptr[i]]] != GRP_RANDOM &&
911 groups[comps[comps_ptr[i]]] != GRP_BBOX) break;
912 }
913 ncomps = i + 1;
914 GV_INFO("ncomps = %d", ncomps);
915
916 *x_poly = gv_calloc(dim * nt, sizeof(double));
917 for (i = 0; i < nt; i++){
918 for (j = 0; j < dim; j++){
919 (*x_poly)[i*dim+j] = Tp[i].center[j];
920 }
921 }
922
923 /*============================================================
924
925 polygon outlines
926
927 ============================================================*/
928 get_poly_lines(nt, E, ncomps, comps_ptr, comps, groups, poly_lines,
929 polys_groups, GRP_RANDOM, GRP_BBOX);
930
931 /*============================================================
932
933 polygon solids
934
935 ============================================================*/
936 get_polygon_solids(nt, E, ncomps, comps_ptr, comps, polys);
937
938 *country_graph = get_country_graph(n, E, groups, GRP_RANDOM, GRP_BBOX);
939
940 free(groups);
941}
942
943static int make_map_internal(bool include_OK_points, int n, int dim, double *x0,
944 int *grouping0, SparseMatrix graph,
945 double bounding_box_margin, int nrandom,
946 int nedgep, double shore_depth_tol, int *nverts,
947 double **x_poly, SparseMatrix *poly_lines,
948 SparseMatrix *polys, int **polys_groups,
949 SparseMatrix *poly_point_map,
950 SparseMatrix *country_graph, int highlight_cluster) {
951
952
953 double xmax[2], xmin[2], area, *x = x0;
954 int i, j;
955 QuadTree qt;
956 int dim2 = 2, nn = 0;
957 int max_qtree_level = 10;
958 double ymin[2], min;
959 int imin, nzok = 0, nzok0 = 0, nt;
960 double *xran, point[2];
961 struct Triangle *Tp;
963 double boxsize[2];
964 bool INCLUDE_OK_POINTS = include_OK_points;/* OK points are random points inserted and found to be within shore_depth_tol of real/artificial points,
965 including them instead of throwing away increase realism of boundary */
966 int *grouping = grouping0;
967
968 int HIGHLIGHT_SET = highlight_cluster;
969
970 for (j = 0; j < dim2; j++) {
971 xmax[j] = x[j];
972 xmin[j] = x[j];
973 }
974
975 for (i = 0; i < n; i++){
976 for (j = 0; j < dim2; j++) {
977 xmax[j] = fmax(xmax[j], x[i*dim+j]);
978 xmin[j] = fmin(xmin[j], x[i*dim+j]);
979 }
980 }
981 boxsize[0] = xmax[0] - xmin[0];
982 boxsize[1] = xmax[1] - xmin[1];
983 area = boxsize[0]*boxsize[1];
984
985 if (nrandom == 0) {
986 nrandom = n;
987 } else if (nrandom < 0){
988 nrandom = -nrandom * n;
989 } else if (nrandom < 4) {/* by default we add 4 point on 4 corners anyway */
990 nrandom = 0;
991 } else {
992 nrandom -= 4;
993 }
994
995 if (shore_depth_tol < 0) shore_depth_tol = sqrt(area/(double) n); /* set to average distance for random distribution */
996 GV_INFO("nrandom=%d shore_depth_tol=%.08f", nrandom, shore_depth_tol);
997
998
999 /* add artificial points along each edge to avoid as much as possible
1000 two connected components be separated due to small shore depth */
1001 {
1002 int nz;
1003 double *y;
1004 int k, t, np=nedgep;
1005 if (graph && np){
1006 fprintf(stderr,"add art np = %d\n",np);
1007 nz = graph->nz;
1008 y = gv_calloc(dim * n + dim * nz * np, sizeof(double));
1009 for (i = 0; i < n*dim; i++) y[i] = x[i];
1010 grouping = gv_calloc(n + nz * np, sizeof(int));
1011 for (i = 0; i < n; i++) grouping[i] = grouping0[i];
1012 nz = n;
1013 for (i = 0; i < graph->m; i++){
1014
1015 for (j = graph->ia[i]; j < graph->ia[i+1]; j++){
1016 if (!HIGHLIGHT_SET || (grouping[i] == grouping[graph->ja[j]] && grouping[i] == HIGHLIGHT_SET)){
1017 for (t = 0; t < np; t++){
1018 for (k = 0; k < dim; k++){
1019 y[nz*dim+k] = t/((double) np)*x[i*dim+k] + (1-t/((double) np))*x[(graph->ja[j])*dim + k];
1020 }
1021 assert(n + (nz-n)*np + t < n + nz*np && n + (nz-n)*np + t >= 0);
1022 if (t/((double) np) > 0.5){
1023 grouping[nz] = grouping[i];
1024 } else {
1025 grouping[nz] = grouping[graph->ja[j]];
1026 }
1027 nz++;
1028 }
1029 }
1030 }
1031 }
1032 fprintf(stderr, "after adding edge points, n:%d->%d\n",n, nz);
1033 n = nz;
1034 x = y;
1035 qt = QuadTree_new_from_point_list(dim, nz, max_qtree_level, y);
1036 } else {
1037 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x);
1038 }
1039 }
1040
1041 /* generate random points for lake/sea effect */
1042 if (nrandom != 0){
1043 for (i = 0; i < dim2; i++) {
1044 if (bounding_box_margin > 0){
1045 xmin[i] -= bounding_box_margin;
1046 xmax[i] += bounding_box_margin;
1047 } else if (bounding_box_margin < 0) {
1048 xmin[i] -= boxsize[i]*(-bounding_box_margin);
1049 xmax[i] += boxsize[i]*(-bounding_box_margin);
1050 } else { // auto bounding box
1051 xmin[i] -= fmax(boxsize[i] * 0.2, 2.* shore_depth_tol);
1052 xmax[i] += fmax(boxsize[i] * 0.2, 2 * shore_depth_tol);
1053 }
1054 }
1055 if (Verbose) {
1056 double bbm = bounding_box_margin;
1057 if (bbm > 0)
1058 fprintf (stderr, "bounding box margin: %.06f", bbm);
1059 else if (bbm < 0)
1060 fprintf (stderr, "bounding box margin: (%.06f * %.06f)", boxsize[0], -bbm);
1061 else
1062 fprintf(stderr, "bounding box margin: %.06f",
1063 fmax(boxsize[0] * 0.2, 2 * shore_depth_tol));
1064 }
1065 if (nrandom < 0) {
1066 const double area2 = (xmax[1] - xmin[1]) * (xmax[0] - xmin[0]);
1067 const double n1 = floor(area2 / (shore_depth_tol * shore_depth_tol));
1068 const double n2 = n * floor(area2 / area);
1069 nrandom = fmax(n1, n2);
1070 }
1071 srand(123);
1072 xran = gv_calloc((nrandom + 4) * dim2, sizeof(double));
1073 int nz = 0;
1074 if (INCLUDE_OK_POINTS){
1075 nzok0 = nzok = nrandom - 1;/* points that are within tolerance of real or artificial points */
1076 if (grouping == grouping0) {
1077 int *grouping2 = gv_calloc(n + nrandom, sizeof(int));
1078 memcpy(grouping2, grouping, sizeof(int)*n);
1079 grouping = grouping2;
1080 } else {
1081 grouping = gv_recalloc(grouping, n, n + nrandom, sizeof(int));
1082 }
1083 }
1084 nn = n;
1085
1086 for (i = 0; i < nrandom; i++){
1087
1088 for (j = 0; j < dim2; j++){
1089 point[j] = xmin[j] + (xmax[j] - xmin[j])*drand();
1090 }
1091
1092 QuadTree_get_nearest(qt, point, ymin, &imin, &min);
1093
1094 if (min > shore_depth_tol){/* point not too close, accepted */
1095 for (j = 0; j < dim2; j++){
1096 xran[nz*dim2+j] = point[j];
1097 }
1098 nz++;
1099 } else if (INCLUDE_OK_POINTS && min > shore_depth_tol/10){/* avoid duplicate points */
1100 for (j = 0; j < dim2; j++){
1101 xran[nzok*dim2+j] = point[j];
1102 }
1103 grouping[nn++] = grouping[imin];
1104 nzok--;
1105
1106 }
1107
1108 }
1109 nrandom = nz;
1110 if (Verbose) fprintf(stderr, "nn nrandom=%d\n", nrandom);
1111 } else {
1112 xran = gv_calloc(4 * dim2, sizeof(double));
1113 }
1114
1115
1116
1117 /* add 4 corners even if nrandom = 0. The corners should be further away from the other points to avoid skinny triangles */
1118 for (i = 0; i < dim2; i++) xmin[i] -= 0.2*(xmax[i]-xmin[i]);
1119 for (i = 0; i < dim2; i++) xmax[i] += 0.2*(xmax[i]-xmin[i]);
1120 i = nrandom;
1121 for (j = 0; j < dim2; j++) xran[i*dim2+j] = xmin[j];
1122 i++;
1123 for (j = 0; j < dim2; j++) xran[i*dim2+j] = xmax[j];
1124 i++;
1125 xran[i*dim2] = xmin[0]; xran[i*dim2+1] = xmax[1];
1126 i++;
1127 xran[i*dim2] = xmax[0]; xran[i*dim2+1] = xmin[1];
1128 nrandom += 4;
1129
1130
1131 double *xcombined;
1132 if (INCLUDE_OK_POINTS){
1133 xcombined = gv_calloc((nn + nrandom) * dim2, sizeof(double));
1134 } else {
1135 xcombined = gv_calloc((n + nrandom) * dim2, sizeof(double));
1136 }
1137 for (i = 0; i < n; i++) {
1138 for (j = 0; j < dim2; j++) xcombined[i*dim2+j] = x[i*dim+j];
1139 }
1140 for (i = 0; i < nrandom; i++) {
1141 for (j = 0; j < dim2; j++) xcombined[(i + nn)*dim2+j] = xran[i*dim+j];
1142 }
1143
1144 if (INCLUDE_OK_POINTS){
1145 for (i = 0; i < nn - n; i++) {
1146 for (j = 0; j < dim2; j++) xcombined[(i + n)*dim2+j] = xran[(nzok0 - i)*dim+j];
1147 }
1148 n = nn;
1149 }
1150
1151
1152 {
1153 int nz, nh = 0;/* the set to highlight */
1154 if (HIGHLIGHT_SET){
1155 if (Verbose) fprintf(stderr," highlight cluster %d, n = %d\n",HIGHLIGHT_SET, n);
1156 /* shift set to the beginning */
1157 nz = 0;
1158 for (i = 0; i < n; i++){
1159 if (grouping[i] == HIGHLIGHT_SET){
1160 nh++;
1161 for (j = 0; j < dim; j++){
1162 xcombined[nz++] = x[i*dim+j];
1163 }
1164 }
1165 }
1166 for (i = 0; i < n; i++){
1167 if (grouping[i] != HIGHLIGHT_SET){
1168 for (j = 0; j < dim; j++){
1169 xcombined[nz++] = x[i*dim+j];
1170 }
1171 }
1172 }
1173 assert(nz == n*dim);
1174 for (i = 0; i < nh; i++){
1175 grouping[i] = 1;
1176 }
1177 for (i = nh; i < n; i++){
1178 grouping[i] = 2;
1179 }
1180 nrandom += n - nh;/* count everything except cluster HIGHLIGHT_SET as random */
1181 n = nh;
1182 if (Verbose) fprintf(stderr,"nh = %d\n",nh);
1183 }
1184 }
1185
1186 int rc = 0;
1187 if (get_tri(n + nrandom, dim2, xcombined, &nt, &Tp, &E) != 0) {
1188 rc = -1;
1189 goto done;
1190 }
1191 get_polygons(n, nrandom, dim2, grouping, nt, Tp, E, nverts, x_poly,
1192 poly_lines, polys, polys_groups, poly_point_map, country_graph);
1193
1195 free(Tp);
1196done:
1197 free(xcombined);
1198 free(xran);
1199 if (grouping != grouping0) free(grouping);
1200 if (x != x0) free(x);
1201 return rc;
1202}
1203
1204static void add_point(int *n, int igrp, double **x, int *nmax, double point[], int **groups){
1205
1206 if (*n >= *nmax){
1207 int old_nmax = *nmax;
1208 *nmax = 20 + *n;
1209 *x = gv_recalloc(*x, 2 * old_nmax, 2 * *nmax, sizeof(double));
1210 *groups = gv_recalloc(*groups, old_nmax, *nmax, sizeof(int));
1211 }
1212
1213 (*x)[(*n)*2] = point[0];
1214 (*x)[(*n)*2+1] = point[1];
1215 (*groups)[*n] = igrp;
1216 (*n)++;
1217}
1218
1219static void get_boundingbox(int n, int dim, double *x, double *width, double *bbox){
1220 int i;
1221 bbox[0] = bbox[1] = x[0];
1222 bbox[2] = bbox[3] = x[1];
1223
1224 for (i = 0; i < n; i++){
1225 bbox[0] = fmin(bbox[0], x[i * dim] - width[i * dim]);
1226 bbox[1] = fmax(bbox[1], x[i * dim] + width[i * dim]);
1227 bbox[2] = fmin(bbox[2], x[i * dim + 1] - width[i * dim + 1]);
1228 bbox[3] = fmax(bbox[3], x[i * dim + 1] + width[i * dim + 1]);
1229 }
1230}
1231
1232int make_map_from_rectangle_groups(bool include_OK_points,
1233 int n, int dim, double *x, double *sizes,
1234 int *grouping, SparseMatrix graph, double bounding_box_margin, int nrandom, int *nart, int nedgep,
1235 double shore_depth_tol,
1236 int *nverts, double **x_poly,
1237 SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map,
1238 SparseMatrix *country_graph, int highlight_cluster){
1239
1240 /* create a list of polygons from a list of rectangles in 2D. rectangles belong to groups. rectangles in the same group that are also close
1241 geometrically will be in the same polygon describing the outline of the group. The main difference for this function and
1242 make_map_from_point_groups is that in this function, the input are points with width/heights, and we try not to place
1243 "lakes" inside these rectangles. This is achieved approximately by adding artificial points along the perimeter of the rectangles,
1244 as well as near the center.
1245
1246 input:
1247 include_OK_points: OK points are random points inserted and found to be within shore_depth_tol of real/artificial points,
1248 . including them instead of throwing away increase realism of boundary
1249 n: number of points
1250 dim: dimension of the points. If dim > 2, only the first 2D is used.
1251 x: coordinates
1252 sizes: width and height
1253 grouping: which group each of the vertex belongs to
1254 graph: the link structure between points. If graph == NULL, this is not used. otherwise
1255 . it is assumed that matrix is symmetric and the graph is undirected
1256 bounding_box_margin: margin used to form the bounding box.
1257 . if negative, it is taken as relative. i.e., -0.5 means a margin of 0.5*box_size
1258 nrandom (input): number of random points to insert in the bounding box to figure out lakes and seas.
1259 . If nrandom = 0, no points are inserted, if nrandom < 0, the number is decided automatically.
1260 .
1261 nart: on entry, number of artificial points to be added along each side of a rectangle enclosing the labels. if < 0, auto-selected.
1262 . On exit, actual number of artificial points added.
1263 nedgep: number of artificial points are adding along edges to establish as much as possible a bright between nodes
1264 . connected by the edge, and avoid islands that are connected. k = 0 mean no points.
1265 shore_depth_tol: nrandom random points are inserted in the bounding box of the points,
1266 . such random points are then weeded out if it is within distance of shore_depth_tol from
1267 . real points. If 0, auto assigned
1268
1269 output:
1270 nverts: number of vertices in the Voronoi diagram
1271 x_poly: the 2D coordinates of these polygons, dimension nverts*2
1272 poly_lines: the sparse matrix representation of the polygon indices, as well as their identity. The matrix is of size
1273 . npolygons x nverts. The i-th polygon is formed by linking vertices with index in the i-th row of the sparse matrix.
1274 . Each row is of the form {{i,j1,m},...{i,jk,m},{i,j1,m},{i,l1,m+1},...}, where j1--j2--jk--j1 form one loop,
1275 . and l1 -- l2 -- ... form another. Each row can have more than 1 loop only when the connected region the polylines represent
1276 . has at least 1 holes.
1277 polys: the sparse matrix representation of the polygon indices, as well as their identity. The matrix is of size
1278 . npolygons x nverts. The i-th polygon is formed by linking vertices with index in the i-th row of the sparse matrix.
1279 . Unlike poly_lines, here each row represent an one stroke drawing of the SOLID polygon, vertices
1280 . along this path may repeat
1281 polys_groups: the group (color) each polygon belongs to, this include all groups of the real points,
1282 . plus the random point group and the bounding box group
1283 poly_point_map: a matrix of dimension npolys x (n + nrandom), poly_point_map[i,j] != 0 if polygon i contains the point j.
1284 . If j < n, it is the original point, otherwise it is artificial point (forming the rectangle around a label) or random points.
1285 country_graph: shows which country is a neighbor of which country.
1286 . if country i and country j are neighbor, then the {i,j} entry is the total number of vertices that
1287 . belongs to i and j, and share an edge of the triangulation. In addition, {i,i} and {j,j} have values equal
1288 . to the number of vertices in each of the countries. If the input "grouping" has negative or zero value, then
1289 . country_graph = NULL.
1290
1291
1292 */
1293 double *X;
1294 int N, nmax, i, j, k, igrp;
1295 int *groups, K = *nart;/* average number of points added per side of rectangle */
1296
1297 double avgsize[2], avgsz, h[2], p1, p0;
1298 double point[2];
1299 int nadded[2];
1300 double delta[2];
1301 double bbox[4];
1302
1303 if (K < 0){
1304 K = (int) 10/(1+n/400.);/* 0 if n > 3600*/
1305 }
1306 *nart = 0;
1307 if (Verbose){
1308 int maxgp = grouping[0];
1309 int mingp = grouping[0];
1310 for (i = 0; i < n; i++) {
1311 maxgp = MAX(maxgp, grouping[i]);
1312 mingp = MIN(mingp, grouping[i]);
1313 }
1314 fprintf(stderr, "max grouping - min grouping + 1 = %d\n",maxgp - mingp + 1);
1315 }
1316
1317 int rc = 0;
1318 if (!sizes){
1319 return make_map_internal(include_OK_points, n, dim, x, grouping, graph,
1320 bounding_box_margin, nrandom, nedgep,
1321 shore_depth_tol, nverts, x_poly, poly_lines, polys,
1322 polys_groups, poly_point_map, country_graph,
1323 highlight_cluster);
1324 } else {
1325
1326 /* add artificial node due to node sizes */
1327 avgsize[0] = 0;
1328 avgsize[1] = 0;
1329 for (i = 0; i < n; i++){
1330 for (j = 0; j < 2; j++) {
1331 avgsize[j] += sizes[i*dim+j];
1332 }
1333 }
1334 for (i = 0; i < 2; i++) avgsize[i] /= n;
1335 avgsz = 0.5*(avgsize[0] + avgsize[1]);
1336 GV_INFO("avgsize = {%f, %f}", avgsize[0], avgsize[1]);
1337
1338 nmax = 2*n;
1339 X = gv_calloc(dim * (n + nmax), sizeof(double));
1340 groups = gv_calloc(n + nmax, sizeof(int));
1341 for (i = 0; i < n; i++) {
1342 groups[i] = grouping[i];
1343 for (j = 0; j < 2; j++){
1344 X[i*2+j] = x[i*dim+j];
1345 }
1346 }
1347 N = n;
1348
1349 if (shore_depth_tol < 0) {
1350 shore_depth_tol = -(shore_depth_tol)*avgsz;
1351 } else if (shore_depth_tol == 0){
1352 get_boundingbox(n, dim, x, sizes, bbox);
1353 const double area = (bbox[1] - bbox[0]) * (bbox[3] - bbox[2]);
1354 shore_depth_tol = sqrt(area / n);
1355 GV_INFO("setting shore length ======%f", shore_depth_tol);
1356 } else {
1357 }
1358
1359 /* add artificial points in an anti-clockwise fashion */
1360
1361 if (K > 0){
1362 delta[0] = .5*avgsize[0]/K; delta[1] = .5*avgsize[1]/K;/* small perturbation to make boundary between labels looks more fractal */
1363 } else {
1364 delta[0] = delta[1] = 0.;
1365 }
1366 for (i = 0; i < n; i++){
1367 igrp = grouping[i];
1368 for (j = 0; j < 2; j++) {
1369 if (avgsz == 0){
1370 nadded[j] = 0;
1371 } else {
1372 nadded[j] = (int) K*sizes[i*dim+j]/avgsz;
1373 }
1374 }
1375
1376 /*top: left to right */
1377 if (nadded[0] > 0){
1378 h[0] = sizes[i*dim]/nadded[0];
1379 point[0] = x[i*dim] - sizes[i*dim]/2;
1380 p1 = point[1] = x[i*dim+1] + sizes[i*dim + 1]/2;
1381 add_point(&N, igrp, &X, &nmax, point, &groups);
1382 for (k = 0; k < nadded[0] - 1; k++){
1383 point[0] += h[0];
1384 point[1] = p1 + (0.5-drand())*delta[1];
1385 add_point(&N, igrp, &X, &nmax, point, &groups);
1386 }
1387
1388 /* bot: right to left */
1389 point[0] = x[i*dim] + sizes[i*dim]/2;
1390 p1 = point[1] = x[i*dim+1] - sizes[i*dim + 1]/2;
1391 add_point(&N, igrp, &X, &nmax, point, &groups);
1392 for (k = 0; k < nadded[0] - 1; k++){
1393 point[0] -= h[0];
1394 point[1] = p1 + (0.5-drand())*delta[1];
1395 add_point(&N, igrp, &X, &nmax, point, &groups);
1396 }
1397 }
1398
1399 if (nadded[1] > 0){
1400 /* left: bot to top */
1401 h[1] = sizes[i*dim + 1]/nadded[1];
1402 p0 = point[0] = x[i*dim] - sizes[i*dim]/2;
1403 point[1] = x[i*dim+1] - sizes[i*dim + 1]/2;
1404 add_point(&N, igrp, &X, &nmax, point, &groups);
1405 for (k = 0; k < nadded[1] - 1; k++){
1406 point[0] = p0 + (0.5-drand())*delta[0];
1407 point[1] += h[1];
1408 add_point(&N, igrp, &X, &nmax, point, &groups);
1409 }
1410
1411 /* right: top to bot */
1412 p0 = point[0] = x[i*dim] + sizes[i*dim]/2;
1413 point[1] = x[i*dim+1] + sizes[i*dim + 1]/2;
1414 add_point(&N, igrp, &X, &nmax, point, &groups);
1415 for (k = 0; k < nadded[1] - 1; k++){
1416 point[0] = p0 + (0.5-drand())*delta[0];
1417 point[1] -= h[1];
1418 add_point(&N, igrp, &X, &nmax, point, &groups);
1419 }
1420 }
1421 *nart = N - n;
1422
1423 }/* done adding artificial points due to node size*/
1424
1425 rc = make_map_internal(include_OK_points, N, dim, X, groups, graph,
1426 bounding_box_margin, nrandom, nedgep,
1427 shore_depth_tol, nverts, x_poly, poly_lines, polys,
1428 polys_groups, poly_point_map, country_graph,
1429 highlight_cluster);
1430 free(groups);
1431 free(X);
1432 }
1433 return rc;
1434}
void QuadTree_get_nearest(QuadTree qt, double *x, double *ymin, int *imin, double *min)
Definition QuadTree.c:682
QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, double *coord)
Definition QuadTree.c:309
void SparseMatrix_distance_matrix(SparseMatrix D0, double **dist0)
SparseMatrix SparseMatrix_from_coordinate_format(SparseMatrix A)
SparseMatrix SparseMatrix_transpose(SparseMatrix A)
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, bool pattern_symmetric_only)
SparseMatrix SparseMatrix_coordinate_form_add_entry(SparseMatrix A, int irn, int jcn, const void *val)
void SparseMatrix_export(FILE *f, SparseMatrix A)
void SparseMatrix_delete(SparseMatrix A)
SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A)
SparseMatrix SparseMatrix_sort(SparseMatrix A)
int * SparseMatrix_weakly_connected_components(SparseMatrix A0, int *ncomp, int **comps)
SparseMatrix SparseMatrix_from_coordinate_format_not_compacted(SparseMatrix A)
SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A)
SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format)
SparseMatrix SparseMatrix_from_dense(int m, int n, double *x)
@ FORMAT_COORD
@ FORMAT_CSR
@ MATRIX_TYPE_PATTERN
@ MATRIX_TYPE_INTEGER
static void agxbfree(agxbuf *xb)
free any malloced resources
Definition agxbuf.h:78
static int agxbprint(agxbuf *xb, const char *fmt,...)
Printf-style output to an agxbuf.
Definition agxbuf.h:234
static WUR char * agxbuse(agxbuf *xb)
Definition agxbuf.h:307
Memory allocation wrappers that exit on failure.
static void * gv_recalloc(void *ptr, size_t old_nmemb, size_t new_nmemb, size_t size)
Definition alloc.h:73
static void * gv_calloc(size_t nmemb, size_t size)
Definition alloc.h:26
#define MIN(a, b)
Definition arith.h:28
#define N(n)
Definition bcomps.c:58
abstract graph C library, Cgraph API
void rgb2hex(float r, float g, float b, agxbuf *cstring, const char *opacity)
Definition colorutil.c:21
void country_graph_coloring(int seed, SparseMatrix A, int **p)
helpers for verbose/debug printing
#define GV_INFO(...)
Definition debug.h:15
int * get_triangles(double *x, int n, int *tris)
Definition delaunay.c:697
#define head
Definition dthdr.h:15
static long seed
Definition exeval.c:1043
#define A(n, t)
Definition expr.h:76
static double dist(int dim, double *x, double *y)
#define E
Definition gdefs.h:6
#define X(prefix, name, str, type, subtype,...)
Definition gdefs.h:14
double drand(void)
Definition general.c:22
void vector_float_take(int n, float *v, int m, int *p, float **u)
Definition general.c:53
double distance_cropped(double *x, int dim, int i, int j)
Definition general.c:116
double xmax
Definition geometry.c:15
double ymin
Definition geometry.c:15
double xmin
Definition geometry.c:15
static bool Verbose
Definition gml2gv.c:23
void free(void *)
node NULL
Definition grammar.y:163
Agsym_t * agattr_text(Agraph_t *g, int kind, char *name, const char *value)
creates or looks up text attributes of a graph
Definition attr.c:357
int agwrite(Agraph_t *g, void *chan)
Return 0 on success, EOF on failure.
Definition write.c:693
@ AGEDGE
Definition cgraph.h:207
@ AGNODE
Definition cgraph.h:207
@ AGRAPH
Definition cgraph.h:207
Agraph_t * graph(char *name)
Definition gv.cpp:30
static int imin(int a, int b)
minimum of two integers
Definition gv_math.h:28
agxbput(xb, staging)
static int z
#define B
Definition hierarchy.c:118
#define D
Definition hierarchy.c:120
#define DEFINE_LIST(name, type)
Definition list.h:22
static int make_map_internal(bool include_OK_points, int n, int dim, double *x0, int *grouping0, SparseMatrix graph, double bounding_box_margin, int nrandom, int nedgep, double shore_depth_tol, int *nverts, double **x_poly, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map, SparseMatrix *country_graph, int highlight_cluster)
Definition make_map.c:943
static void triangle_center(double x[], double y[], double z[], double c[])
Definition make_map.c:177
void map_palette_optimal_coloring(char *color_scheme, SparseMatrix A0, float **rgb_r, float **rgb_g, float **rgb_b)
Definition make_map.c:35
static void plot_dot_edges(FILE *f, SparseMatrix A)
Definition make_map.c:222
static void get_polygons(int n, int nrandom, int dim, int *grouping, int nt, struct Triangle *Tp, SparseMatrix E, int *nverts, double **x_poly, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map, SparseMatrix *country_graph)
Definition make_map.c:870
static void add_point(int *n, int igrp, double **x, int *nmax, double point[], int **groups)
Definition make_map.c:1204
static SparseMatrix get_country_graph(int n, SparseMatrix A, int *groups, int GRP_RANDOM, int GRP_BBOX)
Definition make_map.c:432
static int same_edge(int ecur, int elast, int *edge_table)
Definition make_map.c:613
static void get_boundingbox(int n, int dim, double *x, double *width, double *bbox)
Definition make_map.c:1219
int make_map_from_rectangle_groups(bool include_OK_points, int n, int dim, double *x, double *sizes, int *grouping, SparseMatrix graph, double bounding_box_margin, int nrandom, int *nart, int nedgep, double shore_depth_tol, int *nverts, double **x_poly, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map, SparseMatrix *country_graph, int highlight_cluster)
Definition make_map.c:1232
static void get_poly_lines(int nt, SparseMatrix E, int ncomps, int *comps_ptr, int *comps, int *groups, SparseMatrix *poly_lines, int **polys_groups, int GRP_RANDOM, int GRP_BBOX)
Definition make_map.c:497
static int get_tri(int n, int dim, double *x, int *nt, struct Triangle **T, SparseMatrix *E)
Definition make_map.c:396
static SparseMatrix matrix_add_entry(SparseMatrix A, int i, int j, int val)
Definition make_map.c:213
static void normal(double v[], double normal[])
Definition make_map.c:168
static void cycle_print(int head, int *cycle, int *edge_table)
Definition make_map.c:593
void map_optimal_coloring(int seed, SparseMatrix A, float *rgb_r, float *rgb_g, float *rgb_b)
Definition make_map.c:94
static void plot_dot_labels(FILE *f, int n, int dim, double *x, char **labels, float *fsz)
Definition make_map.c:237
static void dot_polygon(agxbuf *sbuff, doubles_t xp, doubles_t yp, double line_width, bool fill, const char *cstring)
Definition make_map.c:252
static int get_poly_id(int ip, SparseMatrix point_poly_map)
Definition make_map.c:112
static void get_polygon_solids(int nt, SparseMatrix E, int ncomps, int *comps_ptr, int *comps, SparseMatrix *polys)
Definition make_map.c:618
static void plot_dot_polygons(agxbuf *sbuff, double line_width, const char *line_color, SparseMatrix polys, double *x_poly, int *polys_groups, float *r, float *g, float *b, const char *opacity)
Definition make_map.c:279
void plot_dot_map(Agraph_t *gr, int n, int dim, double *x, SparseMatrix polys, SparseMatrix poly_lines, double line_width, const char *line_color, double *x_poly, int *polys_groups, char **labels, float *fsz, float *r, float *g, float *b, const char *opacity, SparseMatrix A, FILE *f)
Definition make_map.c:327
void improve_contiguity(int n, int dim, int *grouping, SparseMatrix poly_point_map, double *x, SparseMatrix graph)
Definition make_map.c:116
static void conn_comp(int n, SparseMatrix A, int *groups, SparseMatrix *poly_point_map)
Definition make_map.c:465
#define cycle_prev(e)
Definition make_map.h:44
#define cycle_next(e)
Definition make_map.h:45
#define edge_tail(e)
Definition make_map.h:43
#define neighbor(t, i, edim, elist)
Definition make_map.h:41
#define edge_head(e)
Definition make_map.h:42
#define delta
Definition maze.c:136
static boxf bbox(Ppoly_t **obsp, int npoly, int *np)
static const int dim
int node_distinct_coloring(const char *color_scheme, int *lightness, bool weightedQ, SparseMatrix A0, double accuracy, int seed, int *cdim0, double **colors)
PATHUTIL_API COORD area2(Ppoint_t, Ppoint_t, Ppoint_t)
Definition visibility.c:44
static const int maxit
Definition power.c:16
#define PRISIZE_T
Definition prisize_t.h:25
void stress_model(int dim, SparseMatrix B, double **x, int maxit_sm, int *flag)
Definition stress_model.c:9
graph or subgraph
Definition cgraph.h:424
double center[2]
Definition make_map.c:165
int vertices[3]
Definition make_map.c:164
Definition types.h:251
Definition geom.h:27
static point center(point vertex[], size_t n)
static mytime_t T
Definition timing.c:41
#define MAX(a, b)
Definition write.c:31