29static double norm(
int n,
const double *x) {
33 for (i = 0; i < n; i++) res += x[i]*x[i];
38static double sqr_dist(
int dim,
const double *x,
const double *y) {
41 for (i = 0; i < dim; i++) res += (x[i] - y[i])*(x[i] - y[i]);
45static double dist(
int dim,
const double *x,
const double *y) {
55 e.
x.assign(x, &x[dim * e.
len]);
68 e.
x.assign(x, &x[dim * e.
len]);
71 e.
wgts = std::vector<double>(np - 1, wgt);
81 double dist1,
dist2, len1, len2;
82 const int dim = e1.
dim;
85 const double *u1 = e1.
x.data();
86 const double *v1 = e1.
x.data() + e1.
npoints * dim - dim;
87 const double *u2 = e2.
x.data();
88 const double *v2 = e2.
x.data() + e2.
npoints * dim - dim;
96 len1 =
dist(dim, u1, v1);
97 len2 =
dist(dim, u2, v2);
98 dist1 = std::max(0.1, dist1 / (len1 + len2 + 0.0001 * dist1));
111 double dist1,
dist2, len1, len2,
len;
112 double tmp, ca, cp, cs;
114 bool flipped =
false;
116 const double *u1 = e1.
x.data();
117 const double *v1 = e1.
x.data() + e1.
npoints * dim - dim;
118 const double *u2 = e2.
x.data();
119 const double *v2 = e2.
x.data() + e2.
npoints * dim - dim;
127 len1 = std::max(
dist(dim, u1, v1),
SMALL);
128 len2 = std::max(
dist(dim, u2, v2),
SMALL);
129 len = 0.5*(len1+len2);
133 for (i = 0; i < dim; i++)
134 ca += (v1[i]-u1[i])*(v2[i]-u2[i]);
135 ca = fabs(ca/(len1*len2));
139 cs = 2 / (std::max(len1, len2) /
len +
len / std::min(len1, len2));
140 assert(cs > -0.001 && cs < 1.001);
144 for (i = 0; i < dim; i++) {
145 tmp = .5*(v1[i]+u1[i])-.5*(v2[i]+u2[i]);
150 assert(cp > -0.001 && cp < 1.001);
163 fprintf(fp,
"#%02x%02x%02x%02x", r, g, b,
alpha);
168 int i, j, k, kk, dim, sta, sto;
169 double maxwgt = 0,
len, len_total0;
171 double tt1[3]={0.15,0.5,0.85};
172 double tt2[4]={0.15,0.4,0.6,0.85};
175 fprintf(fp,
"strict graph{\n");
177 for (i = 0; i < ne; i++){
179 const std::vector<double> &x =
edge.
x;
182 sto =
edge.npoints - 1;
184 fprintf(fp,
"%d [pos=\"", i);
185 for (k = 0; k < dim; k++) {
186 if (k != 0) fprintf(fp,
",");
187 fprintf(fp,
"%f", x[sta*dim+k]);
189 fprintf(fp,
"\"];\n");
191 fprintf(fp,
"%d [pos=\"", i + ne);
192 for (k = 0; k < dim; k++) {
193 if (k != 0) fprintf(fp,
",");
194 fprintf(fp,
"%f", x[sto*dim+k]);
196 fprintf(fp,
"\"];\n");
201 for (i = 0; i < ne; i++){
203 if (!
edge.wgts.empty()) {
204 for (j = 0; j <
edge.npoints - 1; j++) {
205 maxwgt = std::max(maxwgt,
edge.wgts[j]);
211 for (i = 0; i < ne; i++){
212 fprintf(fp,
"%d -- %d [pos=\"", i, i + ne);
214 const std::vector<double> &x =
edge.
x;
217 for (j = 0; j <
edge.npoints; j++) {
223 if (j == 1 || j ==
edge.npoints - 1) {
230 for (kk = 1; kk <= mm; kk++){
232 for (k = 0; k < dim; k++) {
233 if (k != 0) fprintf(fp,
",");
234 fprintf(fp,
"%f", (x[(j-1)*dim+k]*(1-t)+x[j*dim+k]*(t)));
239 if (j == 0 || j ==
edge.npoints - 1){
240 for (k = 0; k < dim; k++) {
241 if (k != 0) fprintf(fp,
",");
242 fprintf(fp,
"%f", x[j*dim+k]);
247 if (!
edge.wgts.empty()) {
248 fprintf(fp,
"\", wgts=\"");
249 for (j = 0; j <
edge.npoints - 1; j++){
250 if (j != 0) fprintf(fp,
",");
251 fprintf(fp,
"%f",
edge.wgts[j]);
255 fprintf(fp,
"\", color=\"");
256 for (j = 0; j <
edge.npoints - 1; j++){
258 for (k = 0; k < dim; k++){
259 len += (
edge.x[dim * j + k] -
edge.x[dim * (j + 1) + k]) * (
edge.x[dim * j + k] -
edge.x[dim * (j + 1) + k]);
264 for (j = 0; j <
edge.npoints - 1; j++) {
266 for (k = 0; k < dim; k++){
267 len += (
edge.x[dim * j + k] -
edge.x[dim * (j + 1) + k]) * (
edge.x[dim * j + k] -
edge.x[dim * (j + 1) + k]);
270 t =
edge.wgts[j] / maxwgt;
272 r = 255*t; g = 0; b = 255*(1-t); b = 255*(1-t);
273 if (j != 0) fprintf(fp,
":");
275 if (j <
edge.npoints - 2) fprintf(fp,
";%f",
len / len_total0);
279 fprintf(fp,
"\"];\n");
286static void pedge_print(
char *comments,
const pedge &e) {
289 fprintf(stderr,
"%s", comments);
290 for (i = 0; i < e.
npoints; i++){
291 if (i > 0) fprintf(stderr,
",");
293 for (j = 0; j < dim; j++){
294 if (j > 0) fprintf(stderr,
",");
295 fprintf(stderr,
"%f", e.
x[dim * i + j]);
299 fprintf(stderr,
"\n");
307 e.
x.resize(e.
dim * n, 0);
308 if (e.
wgts.empty()) {
309 e.
wgts.resize(n - 1);
310 for (i = 0; i < e.
npoints; i++)
313 e.
wgts.resize(n - 1);
323 assert(npoints >= 2);
324 if (npoints*2-1 >
len){
326 e.
x.resize(dim *
len, 0);
328 std::vector<double> &x = e.
x;
330 for (i = npoints - 1; i >= 0; i--){
332 for (j = 0; j < dim; j++){
333 x[dim*ii + j] = x[dim*i + j];
337 for (i = 0; i < npoints - 1; i++){
340 for (j = 0; j < dim; j++){
341 x[dim*(2*i + 1) + j] = 0.5*(x[dim*ii + j] + x[dim*ii2 + j]);
346 e.
edge_length =
dist(dim, &x.data()[0 * dim], &x.data()[(np - 1) * dim]);
350 const std::vector<double> &x = e.
x;
351 const int dim = e.
dim;
359 for (i = 1; i <= np - 2; i++){
362 for (j = 0; j < dim; j++) force[i*dim + j] +=
s*(x[
left*dim + j] - x[i*dim + j]);
363 for (j = 0; j < dim; j++) force[i*dim + j] +=
s*(x[
right*dim + j] - x[i*dim + j]);
368 const pedge &e2, std::vector<double> &force) {
370 const std::vector<double> &x1 = e1.
x, &x2 = e2.
x;
371 const int dim = e1.
dim;
384 s = similarity*edge_length;
385 for (i = 1; i <= np - 2; i++){
386 dist =
sqr_dist(dim, &x1.data()[i*dim], &x2.data()[i*dim]);
388 ss =
s/(
dist+0.1*edge_length*sqrt(
dist));
389 for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[i*dim + j] - x1[i*dim + j]);
393 s = -similarity*edge_length;
394 for (i = 1; i <= np - 2; i++){
395 dist =
sqr_dist(dim, &x1.data()[i*dim], &x2.data()[(np - 1 - i)*dim]);
397 ss =
s/(
dist+0.1*edge_length*sqrt(
dist));
398 for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[(np - 1 - i)*dim + j] - x1[i*dim + j]);
405 std::vector<pedge> &edges,
int maxit,
406 double step0,
double K) {
407 int i, j, ne =
A->
n, k;
408 int *ia =
A->ia, *ja =
A->ja, iter = 0;
409 double *a = (
double*)
A->a;
410 const int np = edges[0].npoints, dim = edges[0].dim;
412 double fnorm_a, fnorm_t, edge_length, start;
415 fprintf(stderr,
"total interaction pairs = %d out of %d, avg neighbors per edge = %f\n",
A->nz,
A->m*
A->m,
A->nz/(
double)
A->m);
417 std::vector<double> force_t(dim * np);
418 std::vector<double> force_a(dim * np);
419 while (step > 0.001 && iter <
maxit){
422 for (i = 0; i < ne; i++){
423 for (j = 0; j < dim*np; j++) {
427 pedge &e1 = edges[i];
428 std::vector<double> &x = e1.
x;
430 for (j = ia[i]; j < ia[i+1]; j++){
431 const pedge &e2 = edges[ja[j]];
434 fnorm_t = std::max(
SMALL,
norm(dim * (np - 2), &force_t.data()[dim]));
435 fnorm_a = std::max(
SMALL,
norm(dim * (np - 2), &force_a.data()[dim]));
438 for (j = 1; j <= np - 2; j++){
439 for (k = 0; k < dim; k++) {
440 x[j * dim + k] += step * edge_length
441 * (force_t[j * dim + k] + K * force_a[j * dim+k])
442 / hypot(fnorm_t, K * fnorm_a);
449 fprintf(stderr,
"iter ==== %d cpu = %f npoints = %d\n",iter, ((
double) (clock() - start))/CLOCKS_PER_SEC, np - 2);
454 std::vector<pedge> &edges,
455 double angle_param,
double angle) {
456 int *assignment =
NULL, nclusters;
458 int *clusterp, *clusters;
472 if (
Verbose > 1) fprintf(stderr,
"there are %d clusters, modularity = %f\n",nclusters, modularity);
476 for (i = 0; i < ne; i++){
485 for (i = 0; i < nclusters; i++) {
486 ink1 =
ink(edges, clusterp[i + 1] - clusterp[i], &clusters[clusterp[i]],
487 &ink0, &meet1, &meet2, angle_param, angle);
489 fprintf(stderr,
"nedges = %d ink0 = %f, ink1 = %f\n",clusterp[i+1] - clusterp[i], ink0,
ink1);
491 for (j = clusterp[i]; j < clusterp[i+1]; j++){
495 pedge &e = edges[clusters[j]];
496 e.
x[1 * dim] = meet1.
x;
497 e.
x[1 * dim + 1] = meet1.
y;
498 e.
x[2 * dim] = meet2.
x;
499 e.
x[2 * dim + 1] = meet2.
y;
500 e.
x[3 * dim] = e.
x[4 * dim];
501 e.
x[3 * dim + 1] = e.
x[4 * dim + 1];
510 const std::vector<pedge> &edges,
511 int compatibility_method,
double tol) {
514 int *ia, *ja, i, j, jj;
519 ia =
A->ia; ja =
A->ja;
521 for (i = 0; i < ne; i++){
522 for (j = ia[i]; j < ia[i+1]; j++){
524 if (i == jj)
continue;
541 fprintf(stderr,
"edge compatibilitu time = %f\n",((
double) (clock() - start))/CLOCKS_PER_SEC);
546 const std::vector<double> &x,
int maxit_outer,
547 double K,
int method,
int nneighbor,
548 int compatibility_method,
int max_recursion,
549 double angle_param,
double angle) {
568 double step0 = 0.1, start = 0.0;
572 std::vector<pedge> edges;
575 for (i = 0; i < ne; i++){
576 edges.emplace_back(
pedge_new(2, dim, &x.data()[dim * 2 * i]));
599 agerrorf(
"Graphviz built without approximate nearest neighbor library ANN; agglomerative inking not available\n");
607 for (k = 0; k < maxit_outer; k++){
608 for (i = 0; i < ne; i++){
621 fprintf(stderr,
"total edge bundling cpu = %f\n",((
double) (clock() - start))/CLOCKS_PER_SEC);
SparseMatrix SparseMatrix_from_coordinate_format(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_delete(SparseMatrix A)
SparseMatrix SparseMatrix_copy(SparseMatrix A)
SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double(*fun)(double x))
SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format)
void agglomerative_ink_bundling(int dim, SparseMatrix A, std::vector< pedge > &edges, int nneighbor, int MAX_RECURSE_LEVEL, double angle_param, double angle)
void modularity_clustering(SparseMatrix A, bool inplace, int ncluster_target, int *nclusters, int **assignment, double *modularity)
static SparseMatrix check_compatibility(SparseMatrix A, int ne, const std::vector< pedge > &edges, int compatibility_method, double tol)
pedge pedge_wgt_new(int np, int dim, double *x, double wgt)
static void force_directed_edge_bundling(SparseMatrix A, std::vector< pedge > &edges, int maxit, double step0, double K)
void pedge_wgts_realloc(pedge &e, int n)
static void edge_attraction_force(double similarity, const pedge &e1, const pedge &e2, std::vector< double > &force)
static void fprint_rgb(FILE *fp, int r, int g, int b, int alpha)
static void modularity_ink_bundling(int dim, int ne, SparseMatrix B, std::vector< pedge > &edges, double angle_param, double angle)
std::vector< pedge > edge_bundling(SparseMatrix A0, int dim, const std::vector< double > &x, int maxit_outer, double K, int method, int nneighbor, int compatibility_method, int max_recursion, double angle_param, double angle)
static double norm(int n, const double *x)
void pedge_export_gv(FILE *fp, int ne, const std::vector< pedge > &edges)
static void edge_tension_force(std::vector< double > &force, const pedge &e)
static double edge_compatibility(const pedge &e1, const pedge &e2)
void pedge_double(pedge &e)
static double edge_compatibility_full(const pedge &e1, const pedge &e2)
static double sqr_dist(int dim, const double *x, const double *y)
static double dist(int dim, const double *x, const double *y)
void pedge_delete(pedge &)
static pedge pedge_new(int np, int dim, const double *x)
static double len(glCompPoint p)
void agerrorf(const char *fmt,...)
double ink(const std::vector< pedge > &edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, double angle_param, double angle)
double ink1(const pedge &e)
PATHUTIL_API COORD dist2(Ppoint_t, Ppoint_t)
std::vector< double > x
coordinates of the npoints poly points. Dimension dim*npoints
std::vector< double > wgts
graphs, nodes and edges info: Agraphinfo_t, Agnodeinfo_t and Agedgeinfo_t