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 * http://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: See CVS logs. Details at http://www.graphviz.org/
9 *************************************************************************/
10
11#include "types.h"
12#include "globals.h"
13#include "general.h"
14#include "time.h"
15#include "SparseMatrix.h"
16#include "vector.h"
17#include "edge_bundling.h"
18#include "ink.h"
19#include "agglomerative_bundling.h"
20#include "nearest_neighbor_graph.h"
21
22#if OPENGL
23#include "gl.h"
24extern pedge *edges_global;
25extern int nedges_global;
26#endif
27
28enum {DEBUG=0};
29
30static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_init(SparseMatrix A, pedge *edges, int level){
31 Agglomerative_Ink_Bundling grid;
32 int n = A->n, i;
33
34 assert(SparseMatrix_is_symmetric(A, TRUE));
35
36 if (!A) return NULL;
37 assert(A->m == n);
38 grid = MALLOC(sizeof(struct Agglomerative_Ink_Bundling_struct));
39 grid->level = level;
40 grid->n = n;
41 grid->A = A;
42 grid->P = NULL;
43 grid->R0 = NULL;
44 grid->R = NULL;
45 grid->next = NULL;
46 grid->prev = NULL;
47 grid->inks = MALLOC(sizeof(real)*(A->m));
48 grid->edges = edges;
49 grid->delete_top_level_A = 0;
50 grid->total_ink = -1;
51 if (level == 0){
52 real total_ink = 0;
53 for (i = 0; i < n; i++) {
54 (grid->inks)[i] = ink1(edges[i]);
55 total_ink += (grid->inks)[i];
56 }
57 grid->total_ink = total_ink;
58 }
59 return grid;
60}
61
62static void Agglomerative_Ink_Bundling_delete(Agglomerative_Ink_Bundling grid){
63 if (!grid) return;
64 if (grid->A){
65 if (grid->level == 0) {
66 if (grid->delete_top_level_A) SparseMatrix_delete(grid->A);
67 } else {
68 SparseMatrix_delete(grid->A);
69 }
70 }
71 SparseMatrix_delete(grid->P);
72 /* on level 0, R0 = NULL, on level 1, R0 = R */
73 if (grid->level > 1) SparseMatrix_delete(grid->R0);
74 SparseMatrix_delete(grid->R);
75 FREE(grid->inks);
76
77 Agglomerative_Ink_Bundling_delete(grid->next);
78 FREE(grid);
79}
80
81static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){
82 /* pick is a work array of dimension n, with n the total number of original edges */
83 int *matching;
84 SparseMatrix A = grid->A;
85 int n = grid->n, level = grid->level, nc = 0;
86 int *ia = A->ia, *ja = A->ja;
87 // real *a;
88 int i, j, k, jj, jc, jmax, ni, nj, npicks;
89 int *mask;
90 pedge *edges = grid->edges;
91 real *inks = grid->inks, *cinks, inki, inkj;
92 real gain, maxgain, minink, total_gain = 0;
93 int *ip = NULL, *jp = NULL, ie;
94 Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid.
95 cedges[i] contain the list of origonal edges that make up the bundle i in the next level */
96 real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
97 point_t meet1, meet2;
98
99 if (Verbose > 1) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n);
100 cedges = MALLOC(sizeof(Vector)*n);
101 cinks = MALLOC(sizeof(real)*n);
102 for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL);
103
104 if (grid->level > 0){
105 ip = grid->R0->ia;
106 jp = grid->R0->ja;
107 }
108
109 matching = MALLOC(sizeof(int)*n);
110 mask = MALLOC(sizeof(real)*n);
111 for (i = 0; i < n; i++) mask[i] = -1;
112
113 assert(n == A->n);
114 for (i = 0; i < n; i++) matching[i] = UNMATCHED;
115
116 // a = (real*) A->a;
117 for (i = 0; i < n; i++){
118 if (matching[i] != UNMATCHED) continue;
119
120 /* find the best matching in ink saving */
121 maxgain = 0;
122 minink = -1;
123 jmax = -1;
124 for (j = ia[i]; j < ia[i+1]; j++){
125 jj = ja[j];
126 if (jj == i) continue;
127
128 /* ink saving of merging i and j */
129 if ((jc=matching[jj]) == UNMATCHED){
130 /* neither i nor jj are matched */
131 inki = inks[i]; inkj = inks[jj];
132 if (ip && jp){/* not the first level */
133 ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
134 nj = (ip[jj+1] - ip[jj]);/* number of edges represented by jj */
135 MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
136 MEMCPY(pick+ni, &(jp[ip[jj]]), sizeof(int)*nj);
137 } else {/* first level */
138 pick[0] = i; pick[1] = jj;
139 ni = nj = 1;
140 }
141 if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj);
142 } else {
143 /* j is already matched. Its content is on cedges[jc] */
144 inki = inks[i]; inkj = cinks[jc];
145 if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc, inkj);
146 if (ip) {
147 ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
148 MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
149 } else {
150 ni = 1; pick[0] = i;
151 }
152 nj = Vector_get_length(cedges[jc]);
153 npicks = ni;
154 for (k = 0; k < nj; k++) {
155 pick[npicks++] = *((int*) Vector_get(cedges[jc], k));
156 }
157 }
158
159 npicks = ni + nj;
160 ink1 = ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle);
161 if (Verbose && DEBUG) {
162 fprintf(stderr,", if merging {");
163 for (k = 0; k < npicks; k++) fprintf(stderr,"%d,", pick[k]);
164 fprintf(stderr,"}, ");
165 fprintf(stderr, " ink0=%f, ink1=%f", inki+inkj, ink1);
166 }
167
168 gain = inki + inkj - ink1;
169 if (Verbose && DEBUG) fprintf(stderr, " gain=%f", gain);
170 if (gain > maxgain){
171 maxgain = gain;
172 minink = ink1;
173 jmax = jj;
174 if (Verbose && DEBUG) fprintf(stderr, "maxgain=%f", maxgain);
175 }
176 if (Verbose && DEBUG) fprintf(stderr, "\n");
177
178
179
180 }
181
182
183 /* now merge i and jmax */
184 if (maxgain > 0){
185 /* a good bundling of i and another edge jmax is found */
186 total_gain += maxgain;
187 jc = matching[jmax];
188 if (jc == UNMATCHED){/* i and j both unmatched. Add j in the table first */
189 if (Verbose && DEBUG) printf("maxgain=%f, merge %d with best edge: %d to form coarsen edge %d. Ink=%f\n",maxgain, i, jmax, nc, minink);
190 matching[i] = matching[jmax] = nc;
191 if (ip){
192 for (k = ip[jmax]; k < ip[jmax+1]; k++) {
193 ie = jp[k];
194 Vector_add(cedges[nc], (void*) (&ie));
195 }
196 } else {
197 Vector_add(cedges[nc], (void*) (&jmax));
198 }
199 jc = nc;
200 nc++;
201 } else {/*j is already matched */
202 if (Verbose && DEBUG) printf("maxgain=%f, merge %d with existing cluster %d\n",maxgain, i, jc);
203 matching[i] = jc;
204 grand_total_ink -= cinks[jc];/* ink of cluster jc was already added, and will be added again as part of a larger cluster with i, so dicount */
205 }
206 } else {/*i can not match/bundle successfully */
207 if (Verbose && DEBUG) fprintf(stderr, "no gain in bundling node %d\n",i);
208 assert(maxgain <= 0);
209 matching[i] = nc;
210 jc = nc;
211 minink = inks[i];
212 nc++;
213 }
214
215
216 /* add i to the appropriate table */
217 if (ip){
218 for (k = ip[i]; k < ip[i+1]; k++) {
219 ie = jp[k];
220 Vector_add(cedges[jc], (void*) (&ie));
221 }
222 } else {
223 Vector_add(cedges[jc], (void*) (&i));
224 }
225 cinks[jc] = minink;
226 grand_total_ink += minink;
227 grand_total_gain += maxgain;
228
229 if (Verbose && DEBUG){
230 fprintf(stderr," coarse edge[%d]={",jc);
231 for (k = 0; k < Vector_get_length(cedges[jc]); k++) {
232 fprintf(stderr,"%d,", *((int*) Vector_get(cedges[jc], k)));
233 }
234 fprintf(stderr,"}, grand_total_gain=%f\n",grand_total_gain);
235 }
236
237 }
238
239 if (nc >= 1 && total_gain > 0){
240 /* now set up restriction and prolongation operator */
241 SparseMatrix P, R, R1, R0, B, cA;
242 real one = 1.;
243 Agglomerative_Ink_Bundling cgrid;
244
245 R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
246 for (i = 0; i < n; i++){
247 jj = matching[i];
248 SparseMatrix_coordinate_form_add_entries(R1, 1, &jj, &i, &one);
249 }
250 R = SparseMatrix_from_coordinate_format(R1);
251 SparseMatrix_delete(R1);
252 P = SparseMatrix_transpose(R);
253 B = SparseMatrix_multiply(R, A);
254 if (!B) goto RETURN;
255 cA = SparseMatrix_multiply(B, P);
256 if (!cA) goto RETURN;
257 SparseMatrix_delete(B);
258 grid->P = P;
259 grid->R = R;
260
261 level++;
262 cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level);
263
264
265 /* set up R0!!! */
266 if (grid->R0){
267 R0 = SparseMatrix_multiply(R, grid->R0);
268 } else {
269 assert(grid->level == 0);
270 R0 = R;
271 }
272 cgrid->R0 = R0;
273 cgrid->inks = cinks;
274 cgrid->total_ink = grand_total_ink;
275
276 if (Verbose > 1) fprintf(stderr,"level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n", grid->level, cgrid->level, grid->n,
277 cgrid->n, grid->total_ink, grand_total_ink, grid->total_ink - grand_total_ink, grand_total_gain);
278 assert(ABS(grid->total_ink - cgrid->total_ink - grand_total_gain) <= 0.0001*grid->total_ink);
279
280 cgrid = Agglomerative_Ink_Bundling_establish(cgrid, pick, angle_param, angle);
281 grid->next = cgrid;
282 cgrid->prev = grid;
283
284 } else {
285 if (Verbose > 1) fprintf(stderr,"no more improvement, orig ink = %f, gain = %f, stop and final bundling found\n", grand_total_ink, grand_total_gain);
286 /* no more improvement, stop and final bundling found */
287 for (i = 0; i < n; i++) matching[i] = i;
288 }
289
290 RETURN:
291 FREE(matching);
292 for (i = 0; i < n; i++) Vector_delete(cedges[i]);
293 FREE(cedges);
294 FREE(mask);
295 return grid;
296}
297
298
299#ifndef NOTUSED
300static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_aggressive_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){
301 /* this does a true single-linkage clustering: find the edge that gives the best saving, merge, find again.
302 As oppose to: find the edge of node i that gives the best ink saving, merge, then do the same for node i+1, etc etc.
303 Right now it is implemented as a quick hack to check whether it is worth while: it saves about 3% extra ink on airlines: from
304 59% to 62%
305 */
306 /* pick is a work array of dimension n, with n the total number of original edges */
307 int *matching;
308 SparseMatrix A = grid->A;
309 int n = grid->n, level = grid->level, nc = 0;
310 int *ia = A->ia, *ja = A->ja;
311 // real *a;
312 int i, j, k, jj, jc = -1, jmax, imax, ni, nj, npicks;
313 int *mask;
314 pedge *edges = grid->edges;
315 real *inks = grid->inks, *cinks, inki, inkj;
316 real gain, maxgain, minink, total_gain = 0;
317 int *ip = NULL, *jp = NULL, ie;
318 Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid.
319 cedges[i] contain the list of origonal edges that make up the bundle i in the next level */
320 real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
321 point_t meet1, meet2;
322
323 if (Verbose > 1) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n);
324 cedges = MALLOC(sizeof(Vector)*n);
325 cinks = MALLOC(sizeof(real)*n);
326 for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL);
327
328 if (grid->level > 0){
329 ip = grid->R0->ia;
330 jp = grid->R0->ja;
331 }
332
333 matching = MALLOC(sizeof(int)*n);
334 mask = MALLOC(sizeof(real)*n);
335 for (i = 0; i < n; i++) mask[i] = -1;
336
337 assert(n == A->n);
338 for (i = 0; i < n; i++) matching[i] = UNMATCHED;
339
340 //a = (real*) A->a;
341
342 do {
343 maxgain = 0;
344 imax = -1;
345 jmax = -1;
346 minink = -1;
347 for (i = 0; i < n; i++){
348 if (matching[i] != UNMATCHED) continue;
349
350 /* find the best matching in ink saving */
351 for (j = ia[i]; j < ia[i+1]; j++){
352 jj = ja[j];
353 if (jj == i) continue;
354
355 /* ink saving of merging i and j */
356 if ((jc=matching[jj]) == UNMATCHED){
357 /* neither i nor jj are matched */
358 inki = inks[i]; inkj = inks[jj];
359 if (ip && jp){/* not the first level */
360 ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
361 nj = (ip[jj+1] - ip[jj]);/* number of edges represented by jj */
362 MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
363 MEMCPY(pick+ni, &(jp[ip[jj]]), sizeof(int)*nj);
364 } else {/* first level */
365 pick[0] = i; pick[1] = jj;
366 ni = nj = 1;
367 }
368 if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj);
369 } else {
370 /* j is already matched. Its content is on cedges[jc] */
371 inki = inks[i]; inkj = cinks[jc];
372 if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc, inkj);
373 if (ip) {
374 ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
375 MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
376 } else {
377 ni = 1; pick[0] = i;
378 }
379 nj = Vector_get_length(cedges[jc]);
380 npicks = ni;
381 for (k = 0; k < nj; k++) {
382 pick[npicks++] = *((int*) Vector_get(cedges[jc], k));
383 }
384 }
385
386 npicks = ni + nj;
387 ink1 = ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle);
388 if (Verbose && DEBUG) {
389 fprintf(stderr,", if merging {");
390 for (k = 0; k < npicks; k++) fprintf(stderr,"%d,", pick[k]);
391 fprintf(stderr,"}, ");
392 fprintf(stderr, " ink0=%f, ink1=%f", inki+inkj, ink1);
393 }
394
395 gain = inki + inkj - ink1;
396 if (Verbose && DEBUG) fprintf(stderr, " gain=%f", gain);
397 if (gain > maxgain){
398 maxgain = gain;
399 minink = ink1;
400 jmax = jj;
401 imax = i;
402 if (Verbose && DEBUG) fprintf(stderr, "maxgain=%f", maxgain);
403 }
404 if (Verbose && DEBUG) fprintf(stderr, "\n");
405
406
407
408 }
409 }
410
411 /* now merge i and jmax */
412 if (maxgain > 0){
413 /* a good bundling of i and another edge jmax is found */
414 total_gain += maxgain;
415 jc = matching[jmax];
416 if (jc == UNMATCHED){/* i and j both unmatched. Add j in the table first */
417 if (Verbose && DEBUG) printf("maxgain=%f, merge %d with best edge: %d to form coarsen edge %d. Ink=%f\n",maxgain, imax, jmax, nc, minink);
418 matching[imax] = matching[jmax] = nc;
419 if (ip){
420 for (k = ip[jmax]; k < ip[jmax+1]; k++) {
421 ie = jp[k];
422 Vector_add(cedges[nc], (void*) (&ie));
423 }
424 } else {
425 Vector_add(cedges[nc], (void*) (&jmax));
426 }
427 jc = nc;
428 nc++;
429 } else {/*j is already matched */
430 if (Verbose && DEBUG) printf("maxgain=%f, merge %d with existing cluster %d\n",maxgain, i, jc);
431 matching[imax] = jc;
432 grand_total_ink -= cinks[jc];/* ink of cluster jc was already added, and will be added again as part of a larger cluster with i, so dicount */
433 }
434 /* add i to the appropriate table */
435 if (ip){
436 for (k = ip[imax]; k < ip[imax+1]; k++) {
437 ie = jp[k];
438 Vector_add(cedges[jc], (void*) (&ie));
439 }
440 } else {
441 Vector_add(cedges[jc], (void*) (&imax));
442 }
443 cinks[jc] = minink;
444 grand_total_ink += minink;
445 grand_total_gain += maxgain;
446 } else {/*i can not match/bundle successfully */
447 if (Verbose && DEBUG) fprintf(stderr, "no gain in bundling node %d\n",i);
448 for (i = 0; i < n; i++){
449 assert(maxgain <= 0);
450 if (matching[i] == UNMATCHED){
451 imax = i;
452 matching[imax] = nc;
453 jc = nc;
454 minink = inks[imax];
455 nc++;
456
457 if (ip){
458 for (k = ip[imax]; k < ip[imax+1]; k++) {
459 ie = jp[k];
460 Vector_add(cedges[jc], (void*) (&ie));
461 }
462 } else {
463 Vector_add(cedges[jc], (void*) (&imax));
464 }
465 cinks[jc] = minink;
466 grand_total_ink += minink;
467 grand_total_gain += maxgain;
468
469
470
471
472 }
473 }
474 }
475
476 } while (maxgain > 0);
477
478
479 if (Verbose && DEBUG){
480 fprintf(stderr," coarse edge[%d]={",jc);
481 for (k = 0; k < Vector_get_length(cedges[jc]); k++) {
482 fprintf(stderr,"%d,", *((int*) Vector_get(cedges[jc], k)));
483 }
484 fprintf(stderr,"}\n");
485 }
486
487 if (nc >= 1 && total_gain > 0){
488 /* now set up restriction and prolongation operator */
489 SparseMatrix P, R, R1, R0, B, cA;
490 real one = 1.;
491 Agglomerative_Ink_Bundling cgrid;
492
493 R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
494 for (i = 0; i < n; i++){
495 jj = matching[i];
496 SparseMatrix_coordinate_form_add_entries(R1, 1, &jj, &i, &one);
497 }
498 R = SparseMatrix_from_coordinate_format(R1);
499 SparseMatrix_delete(R1);
500 P = SparseMatrix_transpose(R);
501 B = SparseMatrix_multiply(R, A);
502 if (!B) goto RETURN;
503 cA = SparseMatrix_multiply(B, P);
504 if (!cA) goto RETURN;
505 SparseMatrix_delete(B);
506 grid->P = P;
507 grid->R = R;
508
509 level++;
510 cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level);
511
512
513 /* set up R0!!! */
514 if (grid->R0){
515 R0 = SparseMatrix_multiply(R, grid->R0);
516 } else {
517 assert(grid->level == 0);
518 R0 = R;
519 }
520 cgrid->R0 = R0;
521 cgrid->inks = cinks;
522 cgrid->total_ink = grand_total_ink;
523
524 if (Verbose > 1) fprintf(stderr,"level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n", grid->level, cgrid->level, grid->n,
525 cgrid->n, grid->total_ink, grand_total_ink, grid->total_ink - grand_total_ink, grand_total_gain);
526 assert(ABS(grid->total_ink - cgrid->total_ink - grand_total_gain) <= 0.0001*grid->total_ink);
527
528 cgrid = Agglomerative_Ink_Bundling_aggressive_establish(cgrid, pick, angle_param, angle);
529 grid->next = cgrid;
530 cgrid->prev = grid;
531
532 } else {
533 if (Verbose > 1) fprintf(stderr,"no more improvement, orig ink = %f, gain = %f, stop and final bundling found\n", grand_total_ink, grand_total_gain);
534 /* no more improvement, stop and final bundling found */
535 for (i = 0; i < n; i++) matching[i] = i;
536 }
537
538 RETURN:
539 FREE(matching);
540 for (i = 0; i < n; i++) Vector_delete(cedges[i]);
541 FREE(mask);
542 return grid;
543}
544#endif
545
546static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_new(SparseMatrix A0, pedge *edges, real angle_param, real angle){
547 /* give a link of edges and their nearest neighbor graph, return a multilevel of edge bundling based on ink saving */
548 Agglomerative_Ink_Bundling grid;
549 int *pick;
550 SparseMatrix A = A0;
551
552 if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){
553 A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
554 }
555 grid = Agglomerative_Ink_Bundling_init(A, edges, 0);
556
557 pick = MALLOC(sizeof(int)*A0->m);
558
559 //grid = Agglomerative_Ink_Bundling_aggressive_establish(grid, pick, angle_param, angle);
560 grid = Agglomerative_Ink_Bundling_establish(grid, pick, angle_param, angle);
561 FREE(pick);
562
563 if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */
564
565 return grid;
566}
567
568static pedge* agglomerative_ink_bundling_internal(int dim, SparseMatrix A, pedge* edges, int nneighbors, int *recurse_level, int MAX_RECURSE_LEVEL, real angle_param, real angle, int open_gl, real *current_ink, real *ink00, int *flag){
569
570 int i, j, jj, k;
571 int *ia, *ja;
572 int *pick;
573 Agglomerative_Ink_Bundling grid, cgrid;
574 SparseMatrix R;
575 real ink0, ink1;
576 point_t meet1, meet2;
577 pedge e;
578 real TOL = 0.0001, wgt_all;
579 clock_t start;
580
581 (*recurse_level)++;
582 if (Verbose > 1) fprintf(stderr, "agglomerative_ink_bundling_internal, recurse level ------- %d\n",*recurse_level);
583
584 assert(A->m == A->n);
585
586 *flag = 0;
587
588 start = clock();
589 grid = Agglomerative_Ink_Bundling_new(A, edges, angle_param, angle);
590 if (Verbose > 1)
591 fprintf(stderr, "CPU for agglomerative bundling %f\n", ((real) (clock() - start))/CLOCKS_PER_SEC);
592 ink0 = grid->total_ink;
593
594 /* find coarsest */
595 cgrid = grid;
596 while (cgrid->next){
597 cgrid = cgrid->next;
598 }
599 ink1 = cgrid->total_ink;
600
601 if (*current_ink < 0){
602 *current_ink = *ink00 = ink0;
603 if (Verbose > 1)
604 fprintf(stderr,"initial total ink = %f\n",*current_ink);
605 }
606 if (ink1 < ink0){
607 *current_ink -= ink0 - ink1;
608 }
609
610 if (Verbose > 1) fprintf(stderr,"ink: %f->%f, edges: %d->%d, current ink = %f, percentage gain over original = %f\n", ink0, ink1, grid->n, cgrid->n, *current_ink, (ink0-ink1)/(*ink00));
611
612 /* if no meaningful improvement (0.0001%), out, else rebundle the middle section */
613 if ((ink0-ink1)/(*ink00) < 0.000001 || *recurse_level > MAX_RECURSE_LEVEL) {
614 /* project bundles up */
615 R = cgrid->R0;
616 if (R){
617 ia = R->ia;
618 ja = R->ja;
619 for (i = 0; i < R->m; i++){
620 pick = &(ja[ia[i]]);
621
622 if (Verbose && DEBUG) fprintf(stderr,"calling ink2...\n");
623 ink1 = ink(edges, ia[i+1]-ia[i], pick, &ink0, &meet1, &meet2, angle_param, angle);
624 if (Verbose && DEBUG) fprintf(stderr,"finish calling ink2...\n");
625 assert(ABS(ink1 - cgrid->inks[i])<=MAX(TOL, TOL*ink1) && ink1 - ink0 <= TOL);
626 wgt_all = 0.;
627 if (ia[i+1]-ia[i] > 1){
628 for (j = ia[i]; j < ia[i+1]; j++){
629 /* make this edge 4 points, insert two meeting points at 1 and 2, make 3 the last point */
630 jj = ja[j];
631 edges[jj] = pedge_double(edges[jj]);/* has to call pedge_double twice: from 2 points to 3 points to 5 points. The last point not used, may be
632 improved later */
633 e = edges[jj] = pedge_double(edges[jj]);
634
635 e->wgts = REALLOC(e->wgts, sizeof(real)*4);
636 e->x[1*dim] = meet1.x;
637 e->x[1*dim+1] = meet1.y;
638 e->x[2*dim] = meet2.x;
639 e->x[2*dim+1] = meet2.y;
640 e->x[3*dim] = e->x[4*dim];
641 e->x[3*dim+1] = e->x[4*dim+1];
642 e->npoints = 4;
643 for (k = 0; k < 3; k++) e->wgts[k] = e->wgt;
644 wgt_all += e->wgt;
645
646 }
647 for (j = ia[i]; j < ia[i+1]; j++){
648 e = edges[ja[j]];
649 e->wgts[1] = wgt_all;
650 }
651 }
652
653 }
654 }
655 } else {
656 pedge *mid_edges, midedge;/* middle section of edges that will be bundled again */
657 real *xx;
658 int ne, npp, l;
659 SparseMatrix A_mid;
660 real eps = 0., wgt, total_wgt = 0;
661
662 /* make new edges using meet1 and meet2.
663
664 call Agglomerative_Ink_Bundling_new
665
666 inherit new edges to old edges
667 */
668
669 R = cgrid->R0;
670 assert(R && cgrid != grid);/* if ink improved, we should have gone at leat 1 level down! */
671 ia = R->ia;
672 ja = R->ja;
673 ne = R->m;
674 mid_edges = MALLOC(sizeof(pedge)*ne);
675 xx = MALLOC(sizeof(real)*4*ne);
676 for (i = 0; i < R->m; i++){
677 pick = &(ja[ia[i]]);
678 wgt = 0.;
679 for (j = ia[i]; j < ia[i+1]; j++) wgt += edges[j]->wgt;
680 total_wgt += wgt;
681 if (Verbose && DEBUG) fprintf(stderr,"calling ink3...\n");
682 ink1 = ink(edges, ia[i+1]-ia[i], pick, &ink0, &meet1, &meet2, angle_param, angle);
683 if (Verbose && DEBUG) fprintf(stderr,"done calling ink3...\n");
684 assert(ABS(ink1 - cgrid->inks[i])<=MAX(TOL, TOL*ink1) && ink1 - ink0 <= TOL);
685 xx[i*4 + 0] = meet1.x;
686 xx[i*4 + 1] = meet1.y;
687 xx[i*4 + 2] = meet2.x;
688 xx[i*4 + 3] = meet2.y;
689 mid_edges[i] = pedge_wgt_new(2, dim, &(xx[i*4]), wgt);
690 //mid_edges[i] = pedge_wgt_new(2, dim, &(xx[i*4]), 1.);
691
692 }
693
694 A_mid = nearest_neighbor_graph(ne, MIN(nneighbors, ne), 4, xx, eps);
695
696 agglomerative_ink_bundling_internal(dim, A_mid, mid_edges, nneighbors, recurse_level, MAX_RECURSE_LEVEL, angle_param, angle, open_gl, current_ink, ink00, flag);
697 SparseMatrix_delete(A_mid);
698 FREE(xx);
699
700 /* patching edges with the new mid-section */
701 for (i = 0; i < R->m; i++){
702 pick = &(ja[ia[i]]);
703 midedge = mid_edges[i];
704 npp = midedge->npoints + 2;
705 for (j = ia[i]; j < ia[i+1]; j++){
706 jj = ja[j];
707 e = edges[jj] = pedge_wgts_realloc(edges[jj], npp);
708
709 assert(e->npoints = 2);
710 for (l = 0; l < dim; l++){/* move the second point to the last */
711 e->x[(npp - 1)*dim+l] = e->x[1*dim+l];
712 }
713
714 for (k = 0; k < midedge->npoints; k++){
715 for (l = 0; l < dim; l++){
716 e->x[(k+1)*dim+l] = midedge->x[k*dim+l];
717 }
718 if (k < midedge->npoints - 1){
719 if (midedge->wgts){
720 e->wgts[(k+1)] = midedge->wgts[k];
721 } else {
722 e->wgts[(k+1)] = midedge->wgt;
723 }
724 }
725 }
726 e->wgts[npp - 2] = e->wgts[0];/* the last interval take from the 1st interval */
727
728
729 e->npoints = npp;
730 }
731#ifdef OPENGL
732 if (open_gl & FALSE){
733 nedges_global = grid->n;
734 edges_global = edges;
735 drawScene();
736 // waitie(5./R->m);
737 }
738#endif
739 }
740
741 for (i = 0; i < ne; i++) pedge_delete(mid_edges[i]);
742
743 }
744
745
746#ifdef OPENGL
747 if (open_gl){
748 nedges_global = grid->n;
749 edges_global = edges;
750 drawScene();
751 // waitie(5./R->m);
752 }
753#endif
754
755 Agglomerative_Ink_Bundling_delete(grid);
756
757 return edges;
758}
759
760
761pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int MAX_RECURSE_LEVEL, real angle_param, real angle, int open_gl, int *flag){
762 int recurse_level = 0;
763 real current_ink = -1, ink0;
764 pedge *edges2;
765
766 ink_count = 0;
767 edges2 = agglomerative_ink_bundling_internal(dim, A, edges, nneighbor, &recurse_level, MAX_RECURSE_LEVEL, angle_param, angle, open_gl, &current_ink, &ink0, flag);
768
769
770 if (Verbose > 1)
771 fprintf(stderr,"initial total ink = %f, final total ink = %f, inksaving = %f percent, total ink_calc = %f, avg ink_calc per edge = %f\n", ink0, current_ink, (ink0-current_ink)/ink0, ink_count, ink_count/(real) A->m);
772 return edges2;
773}
774
775