-
Notifications
You must be signed in to change notification settings - Fork 68
/
par_bubbles.h
1451 lines (1325 loc) · 50.2 KB
/
par_bubbles.h
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
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// BUBBLES :: https://github.com/prideout/par
// Simple C library for packing circles into hierarchical (or flat) diagrams.
//
// Implements "Visualization of Large Hierarchical Data by Circle Packing" from
// Wang et al (2006).
//
// Also contains an implementation of Emo Welzl's "Smallest enclosing disks"
// algorithm (1991) for enclosing points with circles, and a related method
// from Mike Bostock for enclosing circles with the smallest possible circle.
//
// The API is divided into four sections:
//
// - Enclosing. Compute the smallest bounding circle for points or circles.
// - Packing. Pack circles together, or into other circles.
// - Queries. Given a touch point, pick a circle from a hierarchy, etc.
// - Deep Zoom. Uses relative coordinate systems to support arbitrary depth.
//
// In addition to the comment block above each function declaration, the API
// has informal documentation here:
//
// https://prideout.net/bubbles
//
// Distributed under the MIT License, see bottom of file.
#ifndef PAR_BUBBLES_H
#define PAR_BUBBLES_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stdbool.h>
#include <stdint.h>
// This can be any signed integer type.
#ifndef PAR_BUBBLES_INT
#define PAR_BUBBLES_INT int32_t
#endif
// This must be "float" or "double" or "long double". Note that you should not
// need high precision if you use the relative coordinate systems API.
#ifndef PAR_BUBBLES_FLT
#define PAR_BUBBLES_FLT double
#endif
// Enclosing / Touching --------------------------------------------------------
// Read an array of (x,y) coordinates, write a single 3-tuple (x,y,radius).
void par_bubbles_enclose_points(PAR_BUBBLES_FLT const* xy, PAR_BUBBLES_INT npts,
PAR_BUBBLES_FLT* result);
// Read an array of 3-tuples (x,y,radius), write a 3-tuple (x,y,radius).
void par_bubbles_enclose_disks(PAR_BUBBLES_FLT const* xyr,
PAR_BUBBLES_INT ndisks, PAR_BUBBLES_FLT* result);
// Find the circle (x,y,radius) that is tangent to 3 points (x,y).
void par_bubbles_touch_three_points(PAR_BUBBLES_FLT const* xy,
PAR_BUBBLES_FLT* result);
// Find a position for disk "c" that makes it tangent to "a" and "b".
// Note that the ordering of a and b can affect where c will land.
// All three arguments are pointers to three-tuples (x,y,radius).
void par_bubbles_touch_two_disks(PAR_BUBBLES_FLT* c, PAR_BUBBLES_FLT const* a,
PAR_BUBBLES_FLT const* b);
// Returns the smallest circle that intersects the three specified circles.
// This is the problem of Problem of Apollonius!
void par_bubbles_touch_three_disks(PAR_BUBBLES_FLT const* xyr1,
PAR_BUBBLES_FLT const* xyr2, PAR_BUBBLES_FLT const* xyr3,
PAR_BUBBLES_FLT* result);
// Packing ---------------------------------------------------------------------
// Tiny POD structure returned by all packing functions. Private data is
// attached after the public fields, so clients should call the provided
// free function rather than freeing the memory manually.
typedef struct {
PAR_BUBBLES_FLT* xyr; // array of 3-tuples (x y radius) in input order
PAR_BUBBLES_INT count; // number of 3-tuples in "xyr"
PAR_BUBBLES_INT* ids; // populated by par_bubbles_cull
} par_bubbles_t;
void par_bubbles_free_result(par_bubbles_t*);
// Entry point for unbounded non-hierarchical packing. Takes a list of radii.
par_bubbles_t* par_bubbles_pack(PAR_BUBBLES_FLT const* radiuses,
PAR_BUBBLES_INT nradiuses);
// Consume a hierarchy defined by a list of integers. Each integer is an index
// to its parent. The root node is its own parent, and it must be the first node
// in the list. Clients do not have control over individual radiuses, only the
// radius of the outermost enclosing disk.
par_bubbles_t* par_bubbles_hpack_circle(PAR_BUBBLES_INT* nodes,
PAR_BUBBLES_INT nnodes, PAR_BUBBLES_FLT radius);
// Queries ---------------------------------------------------------------------
// Find the node at the given position. Children are on top of their parents.
// If the result is -1, there is no node at the given pick coordinate.
PAR_BUBBLES_INT par_bubbles_pick(par_bubbles_t const*, PAR_BUBBLES_FLT x,
PAR_BUBBLES_FLT y);
// Get bounding box; take a pointer to 4 floats and set them to min xy, max xy.
void par_bubbles_compute_aabb(par_bubbles_t const*, PAR_BUBBLES_FLT* aabb);
// Check if the given circle (3-tuple) intersects the given aabb (4-tuple).
bool par_bubbles_check_aabb(PAR_BUBBLES_FLT const* disk,
PAR_BUBBLES_FLT const* aabb);
// Clip the bubble diagram to the given AABB (4-tuple of left,bottom,right,top)
// and return the result. Circles smaller than the given world-space
// "minradius" are removed. Optionally, an existing diagram (dst) can be passed
// in to receive the culled dataset, which reduces the number of memory allocs
// when calling this function frequently. Pass null to "dst" to create a new
// culled diagram.
par_bubbles_t* par_bubbles_cull(par_bubbles_t const* src,
PAR_BUBBLES_FLT const* aabb, PAR_BUBBLES_FLT minradius, par_bubbles_t* dst);
// Dump out a SVG file for diagnostic purposes.
void par_bubbles_export(par_bubbles_t const* bubbles, char const* filename);
// Returns a pointer to a list of children nodes.
void par_bubbles_get_children(par_bubbles_t const* bubbles, PAR_BUBBLES_INT idx,
PAR_BUBBLES_INT** pchildren, PAR_BUBBLES_INT* nchildren);
// Returns the given node's parent, or 0 if it's the root.
PAR_BUBBLES_INT par_bubbles_get_parent(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT idx);
// Finds the height of the tree and returns one of its deepest leaves.
void par_bubbles_get_maxdepth(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT* maxdepth, PAR_BUBBLES_INT* leaf);
// Finds the height of the tree at a certain node.
PAR_BUBBLES_INT par_bubbles_get_depth(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT node);
// Returns a 4-tuple (min xy, max xy) for the given node.
void par_bubbles_compute_aabb_for_node(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT node, PAR_BUBBLES_FLT* aabb);
// Find the deepest node that is an ancestor of both A and B. Classic!
PAR_BUBBLES_INT par_bubbles_lowest_common_ancestor(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT node_a, PAR_BUBBLES_INT node_b);
// Deep Zoom -------------------------------------------------------------------
// Similar to hpack, but maintains precision by storing disk positions within
// the local coordinate system of their parent. After calling this function,
// clients can use cull_local to flatten the coordinate systems.
par_bubbles_t* par_bubbles_hpack_local(PAR_BUBBLES_INT* nodes,
PAR_BUBBLES_INT nnodes);
// Similar to par_bubbles_cull, but takes a root node rather than an AABB,
// and returns a result within the local coordinate system of the new root.
// In other words, the new root will have radius 1, centered at (0,0). The
// minradius is also expressed in this coordinate system.
par_bubbles_t* par_bubbles_cull_local(par_bubbles_t const* src,
PAR_BUBBLES_FLT const* aabb, PAR_BUBBLES_FLT minradius,
PAR_BUBBLES_INT root, par_bubbles_t* dst);
// Finds the smallest node in the given bubble diagram that completely encloses
// the given axis-aligned bounding box (min xy, max xy). The AABB coordinates
// are expressed in the local coordinate system of the given root node.
PAR_BUBBLES_INT par_bubbles_find_local(par_bubbles_t const* src,
PAR_BUBBLES_FLT const* aabb, PAR_BUBBLES_INT root);
// Similar to pick, but expects (x,y) to be in the coordinate system of the
// given root node.
PAR_BUBBLES_INT par_bubbles_pick_local(par_bubbles_t const*, PAR_BUBBLES_FLT x,
PAR_BUBBLES_FLT y, PAR_BUBBLES_INT root, PAR_BUBBLES_FLT minradius);
// Obtains the scale and translation (which should be applied in that order)
// that can move a point from the node0 coord system to the node1 coord system.
// The "xform" argument should point to three floats, which will be populated
// with: x translation, y translation, and scale.
bool par_bubbles_transform_local(par_bubbles_t const* bubbles,
PAR_BUBBLES_FLT* xform, PAR_BUBBLES_INT node0, PAR_BUBBLES_INT node1);
// Dump out a SVG file for diagnostic purposes.
void par_bubbles_export_local(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT idx, char const* filename);
typedef enum {
PAR_BUBBLES_FILTER_DEFAULT,
PAR_BUBBLES_FILTER_DISCARD_LAST_CHILD,
PAR_BUBBLES_FILTER_KEEP_ONLY_LAST_CHILD
} par_bubbles_filter;
// Special-case function that affects the behavior of subsequent calls to
// cull_local. Allows clients to filter the children list of each non-leaf
// node, which is especially useful when using placeholder bubbles for labels.
void par_bubbles_set_filter(par_bubbles_t* src, par_bubbles_filter f);
typedef enum {
PAR_BUBBLES_HORIZONTAL,
PAR_BUBBLES_VERTICAL
} par_bubbles_orientation;
// Sets some global state that affect subsequent calls to hpack. The first two
// children can either be placed horizontally (default) or vertically. The
// effect of this is subtle, since overall layout is obviously circular.
void par_bubbles_set_orientation(par_bubbles_orientation );
#ifndef PAR_PI
#define PAR_PI (3.14159265359)
#define PAR_MIN(a, b) (a > b ? b : a)
#define PAR_MAX(a, b) (a > b ? a : b)
#define PAR_CLAMP(v, lo, hi) PAR_MAX(lo, PAR_MIN(hi, v))
#define PAR_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
#define PAR_SQR(a) ((a) * (a))
#endif
#ifndef PAR_MALLOC
#define PAR_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
#define PAR_CALLOC(T, N) ((T*) calloc(N * sizeof(T), 1))
#define PAR_REALLOC(T, BUF, N) ((T*) realloc(BUF, sizeof(T) * (N)))
#define PAR_FREE(BUF) free(BUF)
#endif
#ifdef __cplusplus
}
#endif
// -----------------------------------------------------------------------------
// END PUBLIC API
// -----------------------------------------------------------------------------
#ifdef PAR_BUBBLES_IMPLEMENTATION
#define PARINT PAR_BUBBLES_INT
#define PARFLT PAR_BUBBLES_FLT
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <assert.h>
static par_bubbles_orientation par_bubbles__ostate = PAR_BUBBLES_HORIZONTAL;
typedef struct {
PARINT prev;
PARINT next;
} par_bubbles__node;
typedef struct {
PARFLT* xyr; // results array
PARINT count; // client-provided count
PARINT* ids; // populated by par_bubbles_cull
PARFLT const* radiuses; // client-provided radius list
par_bubbles__node* chain; // counterclockwise enveloping chain
PARINT const* graph_parents; // client-provided parent indices
PARINT* graph_children; // flat list of children indices
PARINT* graph_heads; // list of "pointers" to first child
PARINT* graph_tails; // list of "pointers" to one-past-last child
PARINT npacked;
PARINT maxwidth;
PARINT capacity;
par_bubbles_filter filter;
} par_bubbles__t;
static PARFLT par_bubbles__len2(PARFLT const* a)
{
return a[0] * a[0] + a[1] * a[1];
}
static void par_bubbles__initgraph(par_bubbles__t* bubbles)
{
PARINT const* parents = bubbles->graph_parents;
PARINT* nchildren = PAR_CALLOC(PARINT, bubbles->count);
for (PARINT i = 0; i < bubbles->count; i++) {
nchildren[parents[i]]++;
}
PARINT c = 0;
bubbles->graph_heads = PAR_CALLOC(PARINT, bubbles->count * 2);
bubbles->graph_tails = bubbles->graph_heads + bubbles->count;
for (PARINT i = 0; i < bubbles->count; i++) {
bubbles->maxwidth = PAR_MAX(bubbles->maxwidth, nchildren[i]);
bubbles->graph_heads[i] = bubbles->graph_tails[i] = c;
c += nchildren[i];
}
bubbles->graph_heads[0] = bubbles->graph_tails[0] = 1;
bubbles->graph_children = PAR_MALLOC(PARINT, c);
for (PARINT i = 1; i < bubbles->count; i++) {
PARINT parent = parents[i];
bubbles->graph_children[bubbles->graph_tails[parent]++] = i;
}
PAR_FREE(nchildren);
}
static void par_bubbles__initflat(par_bubbles__t* bubbles)
{
PARFLT* xyr = bubbles->xyr;
PARFLT const* radii = bubbles->radiuses;
par_bubbles__node* chain = bubbles->chain;
PARFLT x0, y0, x1, y1;
if (par_bubbles__ostate == PAR_BUBBLES_HORIZONTAL) {
x0 = -radii[0];
y0 = 0;
x1 = radii[1];
y1 = 0;
} else {
x0 = 0;
y0 = -radii[0];
x1 = 0;
y1 = radii[1];
}
*xyr++ = x0;
*xyr++ = y0;
*xyr++ = *radii++;
if (bubbles->count == ++bubbles->npacked) {
return;
}
*xyr++ = x1;
*xyr++ = y1;
*xyr++ = *radii++;
if (bubbles->count == ++bubbles->npacked) {
return;
}
xyr[2] = *radii;
par_bubbles_touch_two_disks(xyr, xyr - 6, xyr - 3);
if (bubbles->count == ++bubbles->npacked) {
return;
}
chain[0].prev = 2;
chain[0].next = 1;
chain[1].prev = 0;
chain[1].next = 2;
chain[2].prev = 1;
chain[2].next = 0;
}
// March forward or backward along the enveloping chain, starting with the
// node at "cn" and testing for collision against the node at "ci".
static PARINT par_bubbles__collide(par_bubbles__t* bubbles, PARINT ci,
PARINT cn, PARINT* cj, PARINT direction)
{
PARFLT const* ci_xyr = bubbles->xyr + ci * 3;
par_bubbles__node* chain = bubbles->chain;
PARINT nsteps = 1;
if (direction > 0) {
for (PARINT i = chain[cn].next; i != cn; i = chain[i].next, ++nsteps) {
PARFLT const* i_xyr = bubbles->xyr + i * 3;
PARFLT dx = i_xyr[0] - ci_xyr[0];
PARFLT dy = i_xyr[1] - ci_xyr[1];
PARFLT dr = i_xyr[2] + ci_xyr[2];
if (0.999 * dr * dr > dx * dx + dy * dy) {
*cj = i;
return nsteps;
}
}
return 0;
}
for (PARINT i = chain[cn].prev; i != cn; i = chain[i].prev, ++nsteps) {
PARFLT const* i_xyr = bubbles->xyr + i * 3;
PARFLT dx = i_xyr[0] - ci_xyr[0];
PARFLT dy = i_xyr[1] - ci_xyr[1];
PARFLT dr = i_xyr[2] + ci_xyr[2];
if (0.999 * dr * dr > dx * dx + dy * dy) {
*cj = i;
return nsteps;
}
}
return 0;
}
static void par_bubbles__packflat(par_bubbles__t* bubbles)
{
PARFLT const* radii = bubbles->radiuses;
PARFLT* xyr = bubbles->xyr;
par_bubbles__node* chain = bubbles->chain;
// Find the circle closest to the origin, known as "Cm" in the paper.
PARINT cm = 0;
PARFLT mindist = par_bubbles__len2(xyr + 0);
PARFLT dist = par_bubbles__len2(xyr + 3);
if (dist > mindist) {
cm = 1;
}
dist = par_bubbles__len2(xyr + 6);
if (dist > mindist) {
cm = 2;
}
// In the paper, "Cn" is always the node that follows "Cm".
PARINT ci, cn = chain[cm].next;
for (ci = bubbles->npacked; ci < bubbles->count; ) {
PARFLT* ci_xyr = xyr + ci * 3;
ci_xyr[2] = radii[ci];
PARFLT* cm_xyr = xyr + cm * 3;
PARFLT* cn_xyr = xyr + cn * 3;
par_bubbles_touch_two_disks(ci_xyr, cn_xyr, cm_xyr);
// Check for a collision. In the paper, "Cj" is the intersecting node.
PARINT cj_f;
PARINT nfsteps = par_bubbles__collide(bubbles, ci, cn, &cj_f, +1);
if (!nfsteps) {
chain[cm].next = ci;
chain[ci].prev = cm;
chain[ci].next = cn;
chain[cn].prev = ci;
cm = ci++;
continue;
}
// Search backwards for a collision, in case it is closer.
PARINT cj_b;
PARINT nbsteps = par_bubbles__collide(bubbles, ci, cm, &cj_b, -1);
// Intersection occurred after Cn.
if (nfsteps <= nbsteps) {
cn = cj_f;
chain[cm].next = cn;
chain[cn].prev = cm;
continue;
}
// Intersection occurred before Cm.
cm = cj_b;
chain[cm].next = cn;
chain[cn].prev = cm;
}
bubbles->npacked = bubbles->count;
}
static void par__disk_from_two_points(PARFLT const* xy1, PARFLT const* xy2,
PARFLT* result)
{
PARFLT dx = xy1[0] - xy2[0];
PARFLT dy = xy1[1] - xy2[1];
result[0] = 0.5 * (xy1[0] + xy2[0]);
result[1] = 0.5 * (xy1[1] + xy2[1]);
result[2] = sqrt(dx * dx + dy * dy) / 2.0;
}
static PARINT par__disk_contains_point(PARFLT const* xyr, PARFLT const* xy)
{
PARFLT dx = xyr[0] - xy[0];
PARFLT dy = xyr[1] - xy[1];
return dx * dx + dy * dy <= PAR_SQR(xyr[2]);
}
static void par__easydisk_from_points(PARFLT* disk, PARFLT const* edgepts,
PARINT nedgepts)
{
if (nedgepts == 0) {
disk[0] = 0;
disk[1] = 0;
disk[2] = 0;
return;
}
if (nedgepts == 1) {
disk[0] = edgepts[0];
disk[1] = edgepts[1];
disk[2] = 0;
return;
}
par__disk_from_two_points(edgepts, edgepts + 2, disk);
if (nedgepts == 2 || par__disk_contains_point(disk, edgepts + 4)) {
return;
}
par__disk_from_two_points(edgepts, edgepts + 4, disk);
if (par__disk_contains_point(disk, edgepts + 2)) {
return;
}
par__disk_from_two_points(edgepts + 2, edgepts + 4, disk);
if (par__disk_contains_point(disk, edgepts)) {
return;
}
par_bubbles_touch_three_points(edgepts, disk);
}
static void par__minidisk_points(PARFLT* disk, PARFLT const* pts, PARINT npts,
PARFLT const* edgepts, PARINT nedgepts)
{
if (npts == 0 || nedgepts == 3) {
par__easydisk_from_points(disk, edgepts, nedgepts);
return;
}
PARFLT const* pt = pts + (--npts) * 2;
par__minidisk_points(disk, pts, npts, edgepts, nedgepts);
if (!par__disk_contains_point(disk, pt)) {
PARFLT edgepts1[6];
for (PARINT i = 0; i < nedgepts * 2; i += 2) {
edgepts1[i] = edgepts[i];
edgepts1[i + 1] = edgepts[i + 1];
}
edgepts1[2 * nedgepts] = pt[0];
edgepts1[2 * nedgepts + 1] = pt[1];
par__minidisk_points(disk, pts, npts, edgepts1, ++nedgepts);
}
}
// Returns true if the specified circle1 contains the specified circle2.
static bool par__disk_contains_disk(PARFLT const* xyr1, PARFLT const* xyr2)
{
PARFLT xc0 = xyr1[0] - xyr2[0];
PARFLT yc0 = xyr1[1] - xyr2[1];
return sqrt(xc0 * xc0 + yc0 * yc0) < xyr1[2] - xyr2[2] + 1e-6;
}
// Returns the smallest circle that intersects the two specified circles.
static void par__disk_from_two_disks(PARFLT const* xyr1, PARFLT const* xyr2,
PARFLT* result)
{
PARFLT x1 = xyr1[0], y1 = xyr1[1], r1 = xyr1[2];
PARFLT x2 = xyr2[0], y2 = xyr2[1], r2 = xyr2[2];
PARFLT x12 = x2 - x1, y12 = y2 - y1, r12 = r2 - r1;
PARFLT l = sqrt(x12 * x12 + y12 * y12);
result[0] = (x1 + x2 + x12 / l * r12) / 2;
result[1] = (y1 + y2 + y12 / l * r12) / 2;
result[2] = (l + r1 + r2) / 2;
}
static void par__easydisk_from_disks(PARFLT* disk, PARFLT const* edgedisks,
PARINT nedgedisks)
{
assert(nedgedisks <= 3);
if (nedgedisks == 0) {
disk[0] = 0;
disk[1] = 0;
disk[2] = 0;
return;
}
if (nedgedisks == 1) {
disk[0] = edgedisks[0];
disk[1] = edgedisks[1];
disk[2] = edgedisks[2];
return;
}
if (nedgedisks == 2) {
par__disk_from_two_disks(edgedisks, edgedisks + 3, disk);
return;
}
par_bubbles_touch_three_disks(edgedisks, edgedisks + 3, edgedisks + 6,
disk);
}
static void par__minidisk_disks(PARFLT* result, PARFLT const* disks,
PARINT ndisks, PARFLT const* edgedisks, PARINT nedgedisks)
{
if (ndisks == 0 || nedgedisks == 3) {
par__easydisk_from_disks(result, edgedisks, nedgedisks);
return;
}
PARFLT const* disk = disks + (--ndisks) * 3;
par__minidisk_disks(result, disks, ndisks, edgedisks, nedgedisks);
if (!par__disk_contains_disk(result, disk)) {
PARFLT edgedisks1[9];
for (PARINT i = 0; i < nedgedisks * 3; i += 3) {
edgedisks1[i] = edgedisks[i];
edgedisks1[i + 1] = edgedisks[i + 1];
edgedisks1[i + 2] = edgedisks[i + 2];
}
edgedisks1[3 * nedgedisks] = disk[0];
edgedisks1[3 * nedgedisks + 1] = disk[1];
edgedisks1[3 * nedgedisks + 2] = disk[2];
par__minidisk_disks(result, disks, ndisks, edgedisks1, ++nedgedisks);
}
}
static void par_bubbles__copy_disk(par_bubbles__t const* src,
par_bubbles__t* dst, PARINT parent)
{
PARINT i = dst->count++;
if (dst->capacity < dst->count) {
dst->capacity = PAR_MAX(16, dst->capacity) * 2;
dst->xyr = PAR_REALLOC(PARFLT, dst->xyr, 3 * dst->capacity);
dst->ids = PAR_REALLOC(PARINT, dst->ids, dst->capacity);
}
PARFLT const* xyr = src->xyr + parent * 3;
dst->xyr[i * 3] = xyr[0];
dst->xyr[i * 3 + 1] = xyr[1];
dst->xyr[i * 3 + 2] = xyr[2];
dst->ids[i] = parent;
}
void par_bubbles_enclose_points(PARFLT const* xy, PARINT npts, PARFLT* result)
{
par__minidisk_points(result, xy, npts, 0, 0);
}
void par_bubbles_enclose_disks(PARFLT const* xyr, PARINT ndisks, PARFLT* result)
{
par__minidisk_disks(result, xyr, ndisks, 0, 0);
}
void par_bubbles_touch_three_points(PARFLT const* xy, PARFLT* xyr)
{
// Many thanks to Stephen Schmitts:
// http://www.abecedarical.com/zenosamples/zs_circle3pts.html
PARFLT p1x = xy[0], p1y = xy[1];
PARFLT p2x = xy[2], p2y = xy[3];
PARFLT p3x = xy[4], p3y = xy[5];
PARFLT a = p2x - p1x, b = p2y - p1y;
PARFLT c = p3x - p1x, d = p3y - p1y;
PARFLT e = a * (p2x + p1x) * 0.5 + b * (p2y + p1y) * 0.5;
PARFLT f = c * (p3x + p1x) * 0.5 + d * (p3y + p1y) * 0.5;
PARFLT det = a*d - b*c;
PARFLT cx = xyr[0] = (d*e - b*f) / det;
PARFLT cy = xyr[1] = (-c*e + a*f) / det;
xyr[2] = sqrt((p1x - cx)*(p1x - cx) + (p1y - cy)*(p1y - cy));
}
void par_bubbles_touch_two_disks(PARFLT* c, PARFLT const* a, PARFLT const* b)
{
PARFLT db = a[2] + c[2], dx = b[0] - a[0], dy = b[1] - a[1];
if (db && (dx || dy)) {
PARFLT da = b[2] + c[2], dc = dx * dx + dy * dy;
da *= da;
db *= db;
PARFLT x = 0.5 + (db - da) / (2 * dc);
PARFLT db1 = db - dc;
PARFLT y0 = PAR_MAX(0, 2 * da * (db + dc) - db1 * db1 - da * da);
PARFLT y = sqrt(y0) / (2 * dc);
c[0] = a[0] + x * dx + y * dy;
c[1] = a[1] + x * dy - y * dx;
} else {
c[0] = a[0] + db;
c[1] = a[1];
}
}
void par_bubbles_touch_three_disks(PARFLT const* xyr1, PARFLT const* xyr2,
PARFLT const* xyr3, PARFLT* result)
{
PARFLT x1 = xyr1[0], y1 = xyr1[1], r1 = xyr1[2],
x2 = xyr2[0], y2 = xyr2[1], r2 = xyr2[2],
x3 = xyr3[0], y3 = xyr3[1], r3 = xyr3[2],
a2 = 2 * (x1 - x2),
b2 = 2 * (y1 - y2),
c2 = 2 * (r2 - r1),
d2 = x1 * x1 + y1 * y1 - r1 * r1 - x2 * x2 - y2 * y2 + r2 * r2,
a3 = 2 * (x1 - x3),
b3 = 2 * (y1 - y3),
c3 = 2 * (r3 - r1),
d3 = x1 * x1 + y1 * y1 - r1 * r1 - x3 * x3 - y3 * y3 + r3 * r3,
ab = a3 * b2 - a2 * b3,
xa = (b2 * d3 - b3 * d2) / ab - x1,
xb = (b3 * c2 - b2 * c3) / ab,
ya = (a3 * d2 - a2 * d3) / ab - y1,
yb = (a2 * c3 - a3 * c2) / ab,
A = xb * xb + yb * yb - 1,
B = 2 * (xa * xb + ya * yb + r1),
C = xa * xa + ya * ya - r1 * r1,
r = (-B - sqrt(B * B - 4 * A * C)) / (2 * A);
result[0] = xa + xb * r + x1;
result[1] = ya + yb * r + y1;
result[2] = r;
}
void par_bubbles_free_result(par_bubbles_t* pubbub)
{
par_bubbles__t* bubbles = (par_bubbles__t*) pubbub;
PAR_FREE(bubbles->graph_children);
PAR_FREE(bubbles->graph_heads);
PAR_FREE(bubbles->chain);
PAR_FREE(bubbles->xyr);
PAR_FREE(bubbles->ids);
PAR_FREE(bubbles);
}
par_bubbles_t* par_bubbles_pack(PARFLT const* radiuses, PARINT nradiuses)
{
par_bubbles__t* bubbles = PAR_CALLOC(par_bubbles__t, 1);
if (nradiuses > 0) {
bubbles->radiuses = radiuses;
bubbles->count = nradiuses;
bubbles->chain = PAR_MALLOC(par_bubbles__node, nradiuses);
bubbles->xyr = PAR_MALLOC(PARFLT, 3 * nradiuses);
par_bubbles__initflat(bubbles);
par_bubbles__packflat(bubbles);
}
return (par_bubbles_t*) bubbles;
}
// Assigns a radius to every node according to its number of descendants.
void par_bubbles__generate_radii(par_bubbles__t* bubbles,
par_bubbles__t* worker, PARINT parent)
{
PARINT head = bubbles->graph_heads[parent];
PARINT tail = bubbles->graph_tails[parent];
PARINT nchildren = tail - head;
PARINT pr = parent * 3 + 2;
bubbles->xyr[pr] = 1;
if (nchildren == 0) {
return;
}
for (PARINT cindex = head; cindex != tail; cindex++) {
PARINT child = bubbles->graph_children[cindex];
par_bubbles__generate_radii(bubbles, worker, child);
bubbles->xyr[pr] += bubbles->xyr[child * 3 + 2];
}
// The following square root seems to produce a nicer, more space-filling,
// distribution of radiuses in randomly-generated trees.
bubbles->xyr[pr] = sqrtf(bubbles->xyr[pr]);
}
void par_bubbles__hpack(par_bubbles__t* bubbles, par_bubbles__t* worker,
PARINT parent, bool local)
{
PARINT head = bubbles->graph_heads[parent];
PARINT tail = bubbles->graph_tails[parent];
PARINT nchildren = tail - head;
if (nchildren == 0) {
return;
}
// Cast away const because we're using the worker as a cache to avoid
// a kazillion malloc / free calls.
PARFLT* radiuses = (PARFLT*) worker->radiuses;
// We perform flat layout twice: once without padding (to determine scale)
// and then again with scaled padding.
PARFLT enclosure[3];
PARFLT px = bubbles->xyr[parent * 3 + 0];
PARFLT py = bubbles->xyr[parent * 3 + 1];
PARFLT pr = bubbles->xyr[parent * 3 + 2];
const PARFLT PAR_HPACK_PADDING1 = 0.15;
const PARFLT PAR_HPACK_PADDING2 = 0.025;
PARFLT scaled_padding = 0.0;
while (1) {
worker->npacked = 0;
worker->count = nchildren;
PARINT c = 0;
for (PARINT cindex = head; cindex != tail; cindex++) {
PARINT child = bubbles->graph_children[cindex];
radiuses[c++] = bubbles->xyr[child * 3 + 2] + scaled_padding;
}
par_bubbles__initflat(worker);
par_bubbles__packflat(worker);
// Using Welzl's algorithm instead of a simple AABB enclosure is
// slightly slower and doesn't yield much aesthetic improvement.
#if PAR_BUBBLES_HPACK_WELZL
par_bubbles_enclose_disks(worker->xyr, nchildren, enclosure);
#else
PARFLT aabb[6];
par_bubbles_compute_aabb((par_bubbles_t const*) worker, aabb);
enclosure[0] = 0.5 * (aabb[0] + aabb[2]);
enclosure[1] = 0.5 * (aabb[1] + aabb[3]);
enclosure[2] = 0;
for (PARINT c = 0; c < nchildren; c++) {
PARFLT x = worker->xyr[c * 3 + 0] - enclosure[0];
PARFLT y = worker->xyr[c * 3 + 1] - enclosure[1];
PARFLT r = worker->xyr[c * 3 + 2];
enclosure[2] = PAR_MAX(enclosure[2], r + sqrtf(x * x + y * y));
}
#endif
if (scaled_padding || !PAR_HPACK_PADDING1) {
break;
} else {
scaled_padding = PAR_HPACK_PADDING1 / enclosure[2];
}
}
PARFLT cx = enclosure[0], cy = enclosure[1], cr = enclosure[2];
scaled_padding *= cr;
cr += PAR_HPACK_PADDING2 * cr;
// Transform the children to fit nicely into either (a) the unit circle,
// or (b) their parent. The former is used if "local" is true.
PARFLT scale, tx, ty;
if (local) {
scale = 1.0 / cr;
tx = 0;
ty = 0;
} else {
scale = pr / cr;
tx = px;
ty = py;
}
PARFLT const* src = worker->xyr;
for (PARINT cindex = head; cindex != tail; cindex++, src += 3) {
PARFLT* dst = bubbles->xyr + 3 * bubbles->graph_children[cindex];
dst[0] = tx + scale * (src[0] - cx);
dst[1] = ty + scale * (src[1] - cy);
dst[2] = scale * (src[2] - scaled_padding);
}
// Recursion. TODO: It might be better to use our own stack here.
for (PARINT cindex = head; cindex != tail; cindex++) {
par_bubbles__hpack(bubbles, worker, bubbles->graph_children[cindex],
local);
}
}
par_bubbles_t* par_bubbles_hpack_circle(PARINT* nodes, PARINT nnodes,
PARFLT radius)
{
par_bubbles__t* bubbles = PAR_CALLOC(par_bubbles__t, 1);
if (nnodes > 0) {
bubbles->graph_parents = nodes;
bubbles->count = nnodes;
bubbles->chain = PAR_MALLOC(par_bubbles__node, nnodes);
bubbles->xyr = PAR_MALLOC(PARFLT, 3 * nnodes);
par_bubbles__initgraph(bubbles);
par_bubbles__t* worker = PAR_CALLOC(par_bubbles__t, 1);
worker->radiuses = PAR_MALLOC(PARFLT, bubbles->maxwidth);
worker->chain = PAR_MALLOC(par_bubbles__node, bubbles->maxwidth);
worker->xyr = PAR_MALLOC(PARFLT, 3 * bubbles->maxwidth);
par_bubbles__generate_radii(bubbles, worker, 0);
bubbles->xyr[0] = 0;
bubbles->xyr[1] = 0;
bubbles->xyr[2] = radius;
par_bubbles__hpack(bubbles, worker, 0, false);
par_bubbles_free_result((par_bubbles_t*) worker);
}
return (par_bubbles_t*) bubbles;
}
// TODO: use a stack instead of recursion
static PARINT par_bubbles__pick(par_bubbles__t const* bubbles, PARINT parent,
PARFLT x, PARFLT y)
{
PARFLT const* xyr = bubbles->xyr + parent * 3;
PARFLT d2 = PAR_SQR(x - xyr[0]) + PAR_SQR(y - xyr[1]);
if (d2 > PAR_SQR(xyr[2])) {
return -1;
}
PARINT head = bubbles->graph_heads[parent];
PARINT tail = bubbles->graph_tails[parent];
for (PARINT cindex = head; cindex != tail; cindex++) {
PARINT child = bubbles->graph_children[cindex];
PARINT result = par_bubbles__pick(bubbles, child, x, y);
if (result > -1) {
return result;
}
}
return parent;
}
PARINT par_bubbles_pick(par_bubbles_t const* cbubbles, PARFLT x, PARFLT y)
{
par_bubbles__t const* bubbles = (par_bubbles__t const*) cbubbles;
if (bubbles->count == 0) {
return -1;
}
return par_bubbles__pick(bubbles, 0, x, y);
}
void par_bubbles_compute_aabb(par_bubbles_t const* bubbles, PARFLT* aabb)
{
if (bubbles->count == 0) {
return;
}
PARFLT const* xyr = bubbles->xyr;
aabb[0] = aabb[2] = xyr[0];
aabb[1] = aabb[3] = xyr[1];
for (PARINT i = 0; i < bubbles->count; i++, xyr += 3) {
aabb[0] = PAR_MIN(xyr[0] - xyr[2], aabb[0]);
aabb[1] = PAR_MIN(xyr[1] - xyr[2], aabb[1]);
aabb[2] = PAR_MAX(xyr[0] + xyr[2], aabb[2]);
aabb[3] = PAR_MAX(xyr[1] + xyr[2], aabb[3]);
}
}
bool par_bubbles_check_aabb(PARFLT const* disk, PARFLT const* aabb)
{
PARFLT cx = PAR_CLAMP(disk[0], aabb[0], aabb[2]);
PARFLT cy = PAR_CLAMP(disk[1], aabb[1], aabb[3]);
PARFLT dx = disk[0] - cx;
PARFLT dy = disk[1] - cy;
PARFLT d2 = dx * dx + dy * dy;
return d2 < (disk[2] * disk[2]);
}
static void par_bubbles__cull(par_bubbles__t const* src, PARFLT const* aabb,
PARFLT minradius, par_bubbles__t* dst, PARINT parent)
{
PARFLT const* xyr = src->xyr + parent * 3;
if (xyr[2] < minradius || !par_bubbles_check_aabb(xyr, aabb)) {
return;
}
par_bubbles__copy_disk(src, dst, parent);
PARINT head = src->graph_heads[parent];
PARINT tail = src->graph_tails[parent];
for (PARINT cindex = head; cindex != tail; cindex++) {
PARINT child = src->graph_children[cindex];
par_bubbles__cull(src, aabb, minradius, dst, child);
}
}
par_bubbles_t* par_bubbles_cull(par_bubbles_t const* psrc,
PARFLT const* aabb, PARFLT minradius, par_bubbles_t* pdst)
{
par_bubbles__t const* src = (par_bubbles__t const*) psrc;
par_bubbles__t* dst = (par_bubbles__t*) pdst;
if (!dst) {
dst = PAR_CALLOC(par_bubbles__t, 1);
pdst = (par_bubbles_t*) dst;
} else {
dst->count = 0;
}
if (src->count == 0) {
return pdst;
}
par_bubbles__cull(src, aabb, minradius, dst, 0);
return pdst;
}
void par_bubbles_export(par_bubbles_t const* bubbles, char const* filename)
{
PARFLT aabb[4];
par_bubbles_compute_aabb(bubbles, aabb);
PARFLT maxextent = PAR_MAX(aabb[2] - aabb[0], aabb[3] - aabb[1]);
PARFLT padding = 0.05 * maxextent;
FILE* svgfile = fopen(filename, "wt");
if (!svgfile) {
fprintf(stderr, "Unable to open %s for writing.\n", filename);
return;
}
fprintf(svgfile,
"<svg viewBox='%f %f %f %f' width='640px' height='640px' "
"version='1.1' "
"xmlns='http://www.w3.org/2000/svg'>\n"
"<g stroke-width='0.5' stroke-opacity='0.5' stroke='black' "
"fill-opacity='0.2' fill='#2A8BB6'>\n"
"<rect fill-opacity='0.1' stroke='none' fill='#2A8BB6' x='%f' y='%f' "
"width='100%%' height='100%%'/>\n",
aabb[0] - padding, aabb[1] - padding,
aabb[2] - aabb[0] + 2 * padding, aabb[3] - aabb[1] + 2 * padding,
aabb[0] - padding, aabb[1] - padding);
PARFLT const* xyr = bubbles->xyr;
for (PARINT i = 0; i < bubbles->count; i++, xyr += 3) {
fprintf(svgfile, "<circle stroke-width='%f' cx='%f' cy='%f' r='%f'/>\n",
xyr[2] * 0.01, xyr[0], xyr[1], xyr[2]);
fprintf(svgfile, "<text text-anchor='middle' stroke='none' "
"x='%f' y='%f' font-size='%f'>%d</text>\n",
xyr[0], xyr[1] + xyr[2] * 0.125, xyr[2] * 0.5, (int) i);
}
fputs("</g>\n</svg>", svgfile);
fclose(svgfile);
}
void par_bubbles_get_children(par_bubbles_t const* pbubbles, PARINT node,
PARINT** pchildren, PARINT* nchildren)
{
par_bubbles__t const* bubbles = (par_bubbles__t const*) pbubbles;
*pchildren = bubbles->graph_children + bubbles->graph_heads[node];
*nchildren = bubbles->graph_tails[node] - bubbles->graph_heads[node];
}
PARINT par_bubbles_get_parent(par_bubbles_t const* pbubbles, PARINT node)
{
par_bubbles__t const* bubbles = (par_bubbles__t const*) pbubbles;
return bubbles->graph_parents[node];
}
void par_bubbles__get_maxdepth(par_bubbles__t const* bubbles, PARINT* maxdepth,
PARINT* leaf, PARINT parent, PARINT depth)
{
if (depth > *maxdepth) {
*leaf = parent;
*maxdepth = depth;
}
PARINT* children;
PARINT nchildren;
par_bubbles_t const* pbubbles = (par_bubbles_t const*) bubbles;
par_bubbles_get_children(pbubbles, parent, &children, &nchildren);
for (PARINT c = 0; c < nchildren; c++) {
par_bubbles__get_maxdepth(bubbles, maxdepth, leaf, children[c],
depth + 1);
}
}
void par_bubbles_get_maxdepth(par_bubbles_t const* pbubbles, PARINT* maxdepth,
PARINT* leaf)
{
par_bubbles__t const* bubbles = (par_bubbles__t const*) pbubbles;
*maxdepth = -1;
*leaf = -1;
return par_bubbles__get_maxdepth(bubbles, maxdepth, leaf, 0, 0);
}
PARINT par_bubbles_get_depth(par_bubbles_t const* pbubbles, PARINT node)
{
par_bubbles__t const* bubbles = (par_bubbles__t const*) pbubbles;
PARINT const* parents = bubbles->graph_parents;
PARINT depth = 0;
while (node) {
node = parents[node];
depth++;
}
return depth;
}
void par_bubbles_compute_aabb_for_node(par_bubbles_t const* bubbles,
PAR_BUBBLES_INT node, PAR_BUBBLES_FLT* aabb)
{
PARFLT const* xyr = bubbles->xyr + 3 * node;
aabb[0] = aabb[2] = xyr[0];
aabb[1] = aabb[3] = xyr[1];
aabb[0] = PAR_MIN(xyr[0] - xyr[2], aabb[0]);
aabb[1] = PAR_MIN(xyr[1] - xyr[2], aabb[1]);
aabb[2] = PAR_MAX(xyr[0] + xyr[2], aabb[2]);
aabb[3] = PAR_MAX(xyr[1] + xyr[2], aabb[3]);