# $Id$ $Revision$
## Process this file with automake to produce Makefile.in
-SUBDIRS = dot tools gvpr lefty lneato dotty smyrna gvmap gvedit
+SUBDIRS = dot tools gvpr lefty lneato dotty smyrna gvmap gvedit mingle
EXTRA_DIST = Makefile.old
--- /dev/null
+# $Id$ $Revision$
+## Process this file with automake to produce Makefile.in
+
+pdfdir = $(pkgdatadir)/doc/pdf
+
+AM_CPPFLAGS = \
+ -I$(top_srcdir) \
+ -I$(top_srcdir)/lib/sparse \
+ -I$(top_srcdir)/lib/sfdpgen \
+ -I$(top_srcdir)/lib/mingle \
+ -I$(top_srcdir)/lib/ingraphs \
+ -I$(top_srcdir)/lib/cgraph \
+ -I$(top_srcdir)/lib/cdt
+
+bin_PROGRAMS = mingle
+
+#man_MANS = mingle.1
+#pdf_DATA = mingle.1.pdf
+man_MANS =
+pdf_DATA =
+
+mingle_SOURCES = minglemain.c
+mingle_CPPFLAGS = $(AM_CPPFLAGS)
+mingle_LDADD = \
+ $(top_builddir)/lib/mingle/libmingle_C.la \
+ $(top_builddir)/lib/ingraphs/libingraphs_C.la \
+ $(top_builddir)/lib/sfdpgen/libsfdpgen_C.la \
+ $(top_builddir)/lib/neatogen/libneatogen_C.la \
+ $(top_builddir)/lib/sparse/libsparse_C.la \
+ $(top_builddir)/lib/rbtree/librbtree_C.la \
+ $(top_builddir)/lib/common/libcommon_C.la \
+ $(top_builddir)/lib/cgraph/libcgraph.la \
+ $(ANN_LIBS) -lstdc++ -lm
+
+mingle.1.pdf: $(srcdir)/mingle.1
+ - @GROFF@ -Tps -man -t $(srcdir)/mingle.1 | @PS2PDF@ - - >mingle.1.pdf
+
+EXTRA_DIST = $(man_MANS) $(pdf_DATA) mingle.vcproj
+
+DISTCLEANFILES = $(pdf_DATA)
--- /dev/null
+<?xml version="1.0" encoding="Windows-1252"?>\r
+<VisualStudioProject\r
+ ProjectType="Visual C++"\r
+ Version="9.00"\r
+ Name="mingle"\r
+ ProjectGUID="{A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}"\r
+ RootNamespace="mingle"\r
+ Keyword="Win32Proj"\r
+ TargetFrameworkVersion="196613"\r
+ >\r
+ <Platforms>\r
+ <Platform\r
+ Name="Win32"\r
+ />\r
+ </Platforms>\r
+ <ToolFiles>\r
+ </ToolFiles>\r
+ <Configurations>\r
+ <Configuration\r
+ Name="Debug|Win32"\r
+ OutputDirectory="Debug"\r
+ IntermediateDirectory="Debug"\r
+ ConfigurationType="1"\r
+ CharacterSet="1"\r
+ >\r
+ <Tool\r
+ Name="VCPreBuildEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCCustomBuildTool"\r
+ />\r
+ <Tool\r
+ Name="VCXMLDataGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCWebServiceProxyGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCMIDLTool"\r
+ />\r
+ <Tool\r
+ Name="VCCLCompilerTool"\r
+ Optimization="0"\r
+ AdditionalIncludeDirectories="../../lib/ingraphs;../../lib/common;../../lib/cdt;../../lib/cgraph;../../lib/sparse;../../lib/sfdpgen;../../;../../lib/mingle"\r
+ PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+ MinimalRebuild="true"\r
+ BasicRuntimeChecks="3"\r
+ RuntimeLibrary="3"\r
+ UsePrecompiledHeader="0"\r
+ WarningLevel="3"\r
+ DebugInformationFormat="4"\r
+ />\r
+ <Tool\r
+ Name="VCManagedResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCPreLinkEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCLinkerTool"\r
+ OutputFile="c:/graphviz-ms/bin/$(ProjectName).exe"\r
+ LinkIncremental="2"\r
+ AdditionalLibraryDirectories=""C:\graphviz-ms\bin";C:\gtk\lib"\r
+ GenerateDebugInformation="true"\r
+ SubSystem="1"\r
+ TargetMachine="1"\r
+ />\r
+ <Tool\r
+ Name="VCALinkTool"\r
+ />\r
+ <Tool\r
+ Name="VCManifestTool"\r
+ />\r
+ <Tool\r
+ Name="VCXDCMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCBscMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCFxCopTool"\r
+ />\r
+ <Tool\r
+ Name="VCAppVerifierTool"\r
+ />\r
+ <Tool\r
+ Name="VCPostBuildEventTool"\r
+ />\r
+ </Configuration>\r
+ <Configuration\r
+ Name="Release|Win32"\r
+ OutputDirectory="Release"\r
+ IntermediateDirectory="Release"\r
+ ConfigurationType="1"\r
+ CharacterSet="1"\r
+ WholeProgramOptimization="1"\r
+ >\r
+ <Tool\r
+ Name="VCPreBuildEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCCustomBuildTool"\r
+ />\r
+ <Tool\r
+ Name="VCXMLDataGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCWebServiceProxyGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCMIDLTool"\r
+ />\r
+ <Tool\r
+ Name="VCCLCompilerTool"\r
+ Optimization="2"\r
+ EnableIntrinsicFunctions="true"\r
+ AdditionalIncludeDirectories="../../lib/ingraphs;../../lib/common;../../lib/cdt;../../lib/cgraph;../../lib/sparse;../../lib/sfdpgen;../../;../../lib/mingle"\r
+ PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+ RuntimeLibrary="2"\r
+ EnableFunctionLevelLinking="true"\r
+ UsePrecompiledHeader="0"\r
+ WarningLevel="3"\r
+ DebugInformationFormat="3"\r
+ />\r
+ <Tool\r
+ Name="VCManagedResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCPreLinkEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCLinkerTool"\r
+ OutputFile="c:/graphviz-ms/bin/$(ProjectName).exe"\r
+ LinkIncremental="2"\r
+ AdditionalLibraryDirectories=""C:\graphviz-ms\bin";C:\gtk\lib"\r
+ GenerateDebugInformation="true"\r
+ SubSystem="1"\r
+ OptimizeReferences="2"\r
+ EnableCOMDATFolding="2"\r
+ TargetMachine="1"\r
+ />\r
+ <Tool\r
+ Name="VCALinkTool"\r
+ />\r
+ <Tool\r
+ Name="VCManifestTool"\r
+ />\r
+ <Tool\r
+ Name="VCXDCMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCBscMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCFxCopTool"\r
+ />\r
+ <Tool\r
+ Name="VCAppVerifierTool"\r
+ />\r
+ <Tool\r
+ Name="VCPostBuildEventTool"\r
+ />\r
+ </Configuration>\r
+ </Configurations>\r
+ <References>\r
+ </References>\r
+ <Files>\r
+ <Filter\r
+ Name="Source Files"\r
+ Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"\r
+ UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"\r
+ >\r
+ <File\r
+ RelativePath=".\minglemain.c"\r
+ >\r
+ </File>\r
+ </Filter>\r
+ <Filter\r
+ Name="Header Files"\r
+ Filter="h;hpp;hxx;hm;inl;inc;xsd"\r
+ UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"\r
+ >\r
+ </Filter>\r
+ <Filter\r
+ Name="Resource Files"\r
+ Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav"\r
+ UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"\r
+ >\r
+ </Filter>\r
+ <File\r
+ RelativePath=".\ReadMe.txt"\r
+ >\r
+ </File>\r
+ </Files>\r
+ <Globals>\r
+ </Globals>\r
+</VisualStudioProject>\r
--- /dev/null
+
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ *************************************************************************/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef WIN32 /*dependencies*/
+ #pragma comment( lib, "cgraph.lib" )
+ #pragma comment( lib, "ingraphs.lib" )
+#endif
+
+#include <cgraph.h>
+#include <ingraphs.h>
+#ifdef HAVE_GETOPT_H
+#include <getopt.h>
+#else
+#include "compat_getopt.h"
+#endif
+
+#include "DotIO.h"
+#include "edge_bundling.h"
+#include "nearest_neighbor_graph.h"
+
+typedef struct {
+ int outer_iter;
+ int method;
+ int compatibility_method;
+ real K;
+ char* fmt;
+ int nneighbors;
+ int max_recursion;
+ real angle_param;
+ real angle;
+} opts_t;
+
+static char *fname;
+static FILE *outfile;
+
+static char* use_msg =
+"Usage: mingle <options> <file>\n\
+ -a t - max. turning angle [0-180] (40)\n\
+ -c i - compatability measure; 0 : distance, 1: full (default)\n\
+ -i iter: number of outer iterations/subdivisions (4)\n\
+ -k k - number of neighbors in the nearest neighbor graph of edges (10)\n\
+ -K k - the force constant\n\
+ -m method - method used. 0 (force directed), 1 (agglomerative ink saving, default), 2 (cluster+ink saving)\n\
+ -o fname - write output to file fname (stdout)\n\
+ -p t - balance for avoiding sharp angles\n\
+ The larger the t, the more sharp angles are allowed\n\
+ -r R - max. recursion level with agglomerative ink saving method (100)\n\
+ -T t - output format\n\
+ -v - verbose\n";
+
+static void
+usage (int eval)
+{
+ fputs (use_msg, stderr);
+ if (eval >= 0) exit (eval);
+}
+
+static void init(int argc, char *argv[], opts_t* opts)
+{
+ unsigned int c;
+ char* CmdName = argv[0];
+
+ opts->outer_iter = 4;
+#ifdef HAVE_ANN
+ opts->method = METHOD_INK_AGGLOMERATE;
+#else
+ opts->method = METHOD_FD;
+#endif
+ opts->compatibility_method = COMPATIBILITY_FULL;
+ opts->K = -1;
+ opts->fmt = "gv";
+ opts->nneighbors = 10;
+ opts->max_recursion = 100;
+ opts->angle_param = -1;
+ opts->angle = 40/180*M_PI;
+
+ while ((c = getopt(argc, argv, ":a:c:i:k:K:m:o:p:r:T:v")) != -1) {
+ switch (c) {
+ case 'v':
+ Verbose = 1;
+ break;
+ case ':':
+ fprintf(stderr, "%s: option -%c missing argument\n", CmdName, optopt);
+ usage(1);
+ break;
+ case '?':
+ if (optopt == '?')
+ usage(0);
+ else
+ fprintf(stderr, "%s: option -%c unrecognized - ignored\n", CmdName, optopt);
+ break;
+ }
+ }
+ argv += optind;
+ argc -= optind;
+
+ if (argc > 0)
+ Files = argv;
+ outfile = stdout;
+}
+
+static int
+bundle (Agraph_t* g, opts_t* opts)
+{
+ real *x = NULL;
+ real *label_sizes = NULL;
+ int n_edge_label_nodes;
+ int dim = 2;
+ SparseMatrix A;
+ SparseMatrix B;
+ pedge* edges;
+ real *xx, eps = 0.;
+ int nz = 0;
+ int *ia, *ja, i, j;
+ int rv = 0;
+
+ initDotIO(g);
+ A = SparseMatrix_import_dot(g, dim, &label_sizes, &x, &n_edge_label_nodes, NULL, FORMAT_CSR, NULL);
+ if (!A){
+ fprintf (stderr, "Error: could not convert graph %s (%s) into matrix\n", agnameof(g), fname);
+ return 1;
+ }
+
+ A = SparseMatrix_symmetrize(A, TRUE);
+ ia = A->ia; ja = A->ja;
+ nz = A->nz;
+ xx = MALLOC(sizeof(real)*nz*4);
+ nz = 0;
+ dim = 4;
+ for (i = 0; i < A->m; i++){
+ for (j = ia[i]; j < ia[i+1]; j++){
+ if (ja[j] > i){
+ xx[nz*dim] = x[i*2];
+ xx[nz*dim+1] = x[i*2+1];
+ xx[nz*dim+2] = x[ja[j]*2];
+ xx[nz*dim+3] = x[ja[j]*2+1];
+ nz++;
+ }
+ }
+ }
+ if (Verbose)
+ fprintf(stderr,"n = %d nz = %d\n",A->m, nz);
+
+ B = nearest_neighbor_graph(nz, MIN(opts->nneighbors, nz), dim, xx, eps);
+
+ SparseMatrix_delete(A);
+ A = B;
+ FREE(x);
+ x = xx;
+
+ dim = 2;
+
+ edges = edge_bundling(A, 2, x, opts->outer_iter, opts->K, opts->method, opts->nneighbors, opts->compatibility_method, opts->max_recursion, opts->angle_param, opts->angle, 0);
+
+ pedge_export_gv(stdout, A->m, edges);
+ return rv;
+}
+
+static Agraph_t *gread(FILE * fp)
+{
+ return agread(fp, (Agdisc_t *) 0);
+}
+
+int main(int argc, char *argv[])
+{
+ Agraph_t *g;
+ Agraph_t *prev = NULL;
+ ingraph_state ig;
+ int rv = 0;
+ opts_t opts;
+
+ init(argc, argv, &opts);
+ newIngraph(&ig, Files, gread);
+
+ while ((g = nextGraph(&ig)) != 0) {
+ if (prev)
+ agclose(prev);
+ prev = g;
+ fname = fileName(&ig);
+ if (Verbose)
+ fprintf(stderr, "Process graph %s in file %s\n", agnameof(g),
+ fname);
+ rv |= bundle (g, &opts);
+ }
+
+ return rv;
+}
fi
AM_CONDITIONAL(WITH_GTS, [test "x$use_gts" == "xYes"])
+dnl -----------------------------------
+dnl INCLUDES and LIBS for ANN.
+
+AC_ARG_WITH(ann,
+ [AS_HELP_STRING([--with-ann=yes],[ANN library])],
+ [], [with_ann=yes])
+
+if test "x$with_ann" != "xyes"; then
+ use_ann="No (disabled)"
+else
+ PKG_CHECK_MODULES(ANN, [ann],[
+ use_ann="Yes"
+ AC_DEFINE_UNQUOTED(HAVE_ANN,1,
+ [Define if you have the ann library])
+ AC_SUBST([ANN_CFLAGS])
+ AC_SUBST([ANN_LIBS])
+ ],[
+ use_ann="No (ANN library not available)"
+ ])
+fi
+AM_CONDITIONAL(WITH_ANN, [test "x$use_ann" == "xYes"])
+
dnl -----------------------------------
dnl INCLUDES and LIBS for GLADE.
lib/neatogen/Makefile
lib/fdpgen/Makefile
lib/sparse/Makefile
+ lib/mingle/Makefile
lib/label/Makefile
lib/sfdpgen/Makefile
lib/osage/Makefile
cmd/dotty/Makefile
cmd/smyrna/Makefile
cmd/gvmap/Makefile
+ cmd/mingle/Makefile
cmd/gvedit/Makefile
cmd/gvedit/gvedit.pro
cmd/gvedit/ui/Makefile
echo " fontconfig: $use_fontconfig"
echo " freetype: $use_freetype"
echo " glut: $use_glut"
+echo " ann: $use_ann"
echo " gts: $use_gts"
echo " ipsepcola: $use_ipsepcola"
echo " ltdl: $use_ltdl"
Microsoft Visual Studio Solution File, Format Version 10.00\r
-#
# Visual Studio 2008\r
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "Pathplan", "lib\pathplan\Pathplan.vcproj", "{BD347753-A09D-48B4-8752-F1D8D9CF235D}"\r
ProjectSection(ProjectDependencies) = postProject\r
{15229511-9F6C-48A5-9194-660CA6492563} = {15229511-9F6C-48A5-9194-660CA6492563}\r
EndProjectSection\r
EndProject\r
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "minglelib", "lib\mingle\minglelib.vcproj", "{70575BD2-A532-41B8-9A17-B9876E992A84}"\r
+EndProject\r
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "mingle", "cmd\mingle\mingle.vcproj", "{A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}"\r
+ ProjectSection(ProjectDependencies) = postProject\r
+ {C0663A08-F276-4DD6-B17C-E501EE066F7C} = {C0663A08-F276-4DD6-B17C-E501EE066F7C}\r
+ {15229511-9F6C-48A5-9194-660CA6492563} = {15229511-9F6C-48A5-9194-660CA6492563}\r
+ {D6CEB142-BF8E-471C-AE16-4300F2D5DEDA} = {D6CEB142-BF8E-471C-AE16-4300F2D5DEDA}\r
+ {B76BCE8C-63CC-4A99-88B5-D621D563E699} = {B76BCE8C-63CC-4A99-88B5-D621D563E699}\r
+ {C5CEB09E-79AF-4091-ACCF-D28EC3D7D90F} = {C5CEB09E-79AF-4091-ACCF-D28EC3D7D90F}\r
+ {443EB1A7-C634-4292-9F2D-C752BB2BF40F} = {443EB1A7-C634-4292-9F2D-C752BB2BF40F}\r
+ {70575BD2-A532-41B8-9A17-B9876E992A84} = {70575BD2-A532-41B8-9A17-B9876E992A84}\r
+ {D6FD0DE5-5305-458E-8CA5-FCA4B8E05B04} = {D6FD0DE5-5305-458E-8CA5-FCA4B8E05B04}\r
+ EndProjectSection\r
+EndProject\r
Global\r
GlobalSection(SolutionConfigurationPlatforms) = preSolution\r
Debug|Win32 = Debug|Win32\r
{33E97266-F213-3639-BCAB-2F3F95E15B16}.Debug|Win32.Build.0 = Debug|Win32\r
{33E97266-F213-3639-BCAB-2F3F95E15B16}.Release|Win32.ActiveCfg = Release|Win32\r
{33E97266-F213-3639-BCAB-2F3F95E15B16}.Release|Win32.Build.0 = Release|Win32\r
+ {70575BD2-A532-41B8-9A17-B9876E992A84}.Debug|Win32.ActiveCfg = Debug|Win32\r
+ {70575BD2-A532-41B8-9A17-B9876E992A84}.Debug|Win32.Build.0 = Debug|Win32\r
+ {70575BD2-A532-41B8-9A17-B9876E992A84}.Release|Win32.ActiveCfg = Release|Win32\r
+ {70575BD2-A532-41B8-9A17-B9876E992A84}.Release|Win32.Build.0 = Release|Win32\r
+ {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Debug|Win32.ActiveCfg = Debug|Win32\r
+ {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Debug|Win32.Build.0 = Debug|Win32\r
+ {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Release|Win32.ActiveCfg = Release|Win32\r
+ {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Release|Win32.Build.0 = Release|Win32\r
EndGlobalSection\r
GlobalSection(SolutionProperties) = preSolution\r
HideSolutionNode = FALSE\r
SUBDIRS = cdt graph cgraph pathplan sfio vmalloc ast \
vpsc rbtree ortho sparse patchwork expr common \
- pack xdot label gvc ingraphs topfish glcomp \
+ pack xdot label gvc ingraphs topfish glcomp mingle \
circogen dotgen dotgen2 fdpgen neatogen twopigen sfdpgen osage gvpr
EXTRA_DIST = Makefile.old gvc.vcproj gvc.def
--- /dev/null
+# $Id$Revision:
+## Process this file with automake to produce Makefile.in
+
+AM_CPPFLAGS = \
+ -I$(top_srcdir) \
+ -I$(top_srcdir)/lib/common \
+ -I$(top_srcdir)/lib/sparse \
+ -I$(top_srcdir)/lib/gvc \
+ -I$(top_srcdir)/lib/pathplan \
+ -I$(top_srcdir)/lib/sfdpgen \
+ -I$(top_srcdir)/lib/cgraph \
+ -I$(top_srcdir)/lib/cdt $(ANN_CFLAGS)
+
+noinst_HEADERS = edge_bundling.h ink.h agglomerative_bundling.h nearest_neighbor_graph.h nearest_neighbor_graph_ann.h
+
+noinst_LTLIBRARIES = libmingle_C.la
+
+libmingle_C_la_SOURCES = edge_bundling.c ink.c agglomerative_bundling.c nearest_neighbor_graph.c nearest_neighbor_graph_ann.cpp
+
+EXTRA_DIST = minglelib.vcproj
+
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#include "types.h"
+#include "globals.h"
+#include "general.h"
+#include "time.h"
+#include "SparseMatrix.h"
+#include "vector.h"
+#include "edge_bundling.h"
+#include "ink.h"
+#include "agglomerative_bundling.h"
+#include "nearest_neighbor_graph.h"
+
+#if OPENGL
+#include "gl.h"
+extern pedge *edges_global;
+extern int nedges_global;
+#endif
+
+enum {DEBUG=0};
+
+static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_init(SparseMatrix A, pedge *edges, int level){
+ Agglomerative_Ink_Bundling grid;
+ int n = A->n, i;
+
+ assert(SparseMatrix_is_symmetric(A, TRUE));
+
+ if (!A) return NULL;
+ assert(A->m == n);
+ grid = MALLOC(sizeof(struct Agglomerative_Ink_Bundling_struct));
+ grid->level = level;
+ grid->n = n;
+ grid->A = A;
+ grid->P = NULL;
+ grid->R0 = NULL;
+ grid->R = NULL;
+ grid->next = NULL;
+ grid->prev = NULL;
+ grid->inks = MALLOC(sizeof(real)*(A->m));
+ grid->edges = edges;
+ grid->delete_top_level_A = 0;
+ grid->total_ink = -1;
+ if (level == 0){
+ real total_ink = 0;
+ for (i = 0; i < n; i++) {
+ (grid->inks)[i] = ink1(edges[i]);
+ total_ink += (grid->inks)[i];
+ }
+ grid->total_ink = total_ink;
+ }
+ return grid;
+}
+
+static void Agglomerative_Ink_Bundling_delete(Agglomerative_Ink_Bundling grid){
+ if (!grid) return;
+ if (grid->A){
+ if (grid->level == 0) {
+ if (grid->delete_top_level_A) SparseMatrix_delete(grid->A);
+ } else {
+ SparseMatrix_delete(grid->A);
+ }
+ }
+ SparseMatrix_delete(grid->P);
+ /* on level 0, R0 = NULL, on level 1, R0 = R */
+ if (grid->level > 1) SparseMatrix_delete(grid->R0);
+ SparseMatrix_delete(grid->R);
+ FREE(grid->inks);
+
+ Agglomerative_Ink_Bundling_delete(grid->next);
+ FREE(grid);
+}
+
+static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){
+ /* pick is a work array of dimension n, with n the total number of original edges */
+ int *matching;
+ SparseMatrix A = grid->A;
+ int n = grid->n, level = grid->level, nc = 0;
+ int *ia = A->ia, *ja = A->ja;
+ // real *a;
+ int i, j, k, jj, jc, jmax, ni, nj, npicks;
+ int *mask;
+ pedge *edges = grid->edges;
+ real *inks = grid->inks, *cinks, inki, inkj;
+ real gain, maxgain, minink, total_gain = 0;
+ int *ip = NULL, *jp = NULL, ie;
+ Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid.
+ cedges[i] contain the list of origonal edges that make up the bundle i in the next level */
+ real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
+ point_t meet1, meet2;
+
+ if (Verbose) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n);
+ cedges = MALLOC(sizeof(Vector)*n);
+ cinks = MALLOC(sizeof(real)*n);
+ for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL);
+
+ if (grid->level > 0){
+ ip = grid->R0->ia;
+ jp = grid->R0->ja;
+ }
+
+ matching = MALLOC(sizeof(int)*n);
+ mask = MALLOC(sizeof(real)*n);
+ for (i = 0; i < n; i++) mask[i] = -1;
+
+ assert(n == A->n);
+ for (i = 0; i < n; i++) matching[i] = UNMATCHED;
+
+ // a = (real*) A->a;
+ for (i = 0; i < n; i++){
+ if (matching[i] != UNMATCHED) continue;
+
+ /* find the best matching in ink saving */
+ maxgain = 0;
+ minink = -1;
+ jmax = -1;
+ for (j = ia[i]; j < ia[i+1]; j++){
+ jj = ja[j];
+ if (jj == i) continue;
+
+ /* ink saving of merging i and j */
+ if ((jc=matching[jj]) == UNMATCHED){
+ /* neither i nor jj are matched */
+ inki = inks[i]; inkj = inks[jj];
+ if (ip && jp){/* not the first level */
+ ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
+ nj = (ip[jj+1] - ip[jj]);/* number of edges represented by jj */
+ MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
+ MEMCPY(pick+ni, &(jp[ip[jj]]), sizeof(int)*nj);
+ } else {/* first level */
+ pick[0] = i; pick[1] = jj;
+ ni = nj = 1;
+ }
+ if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj);
+ } else {
+ /* j is already matched. Its content is on cedges[jc] */
+ inki = inks[i]; inkj = cinks[jc];
+ if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc, inkj);
+ if (ip) {
+ ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
+ MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
+ } else {
+ ni = 1; pick[0] = i;
+ }
+ nj = Vector_get_length(cedges[jc]);
+ npicks = ni;
+ for (k = 0; k < nj; k++) {
+ pick[npicks++] = *((int*) Vector_get(cedges[jc], k));
+ }
+ }
+
+ npicks = ni + nj;
+ ink1 = ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle);
+ if (Verbose && DEBUG) {
+ fprintf(stderr,", if merging {");
+ for (k = 0; k < npicks; k++) fprintf(stderr,"%d,", pick[k]);
+ fprintf(stderr,"}, ");
+ fprintf(stderr, " ink0=%f, ink1=%f", inki+inkj, ink1);
+ }
+
+ gain = inki + inkj - ink1;
+ if (Verbose && DEBUG) fprintf(stderr, " gain=%f", gain);
+ if (gain > maxgain){
+ maxgain = gain;
+ minink = ink1;
+ jmax = jj;
+ if (Verbose && DEBUG) fprintf(stderr, "maxgain=%f", maxgain);
+ }
+ if (Verbose && DEBUG) fprintf(stderr, "\n");
+
+
+
+ }
+
+
+ /* now merge i and jmax */
+ if (maxgain > 0){
+ /* a good bundling of i and another edge jmax is found */
+ total_gain += maxgain;
+ jc = matching[jmax];
+ if (jc == UNMATCHED){/* i and j both unmatched. Add j in the table first */
+ 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);
+ matching[i] = matching[jmax] = nc;
+ if (ip){
+ for (k = ip[jmax]; k < ip[jmax+1]; k++) {
+ ie = jp[k];
+ Vector_add(cedges[nc], (void*) (&ie));
+ }
+ } else {
+ Vector_add(cedges[nc], (void*) (&jmax));
+ }
+ jc = nc;
+ nc++;
+ } else {/*j is already matched */
+ if (Verbose && DEBUG) printf("maxgain=%f, merge %d with existing cluster %d\n",maxgain, i, jc);
+ matching[i] = jc;
+ 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 */
+ }
+ } else {/*i can not match/bundle successfully */
+ if (Verbose && DEBUG) fprintf(stderr, "no gain in bundling node %d\n",i);
+ assert(maxgain <= 0);
+ matching[i] = nc;
+ jc = nc;
+ minink = inks[i];
+ nc++;
+ }
+
+
+ /* add i to the appropriate table */
+ if (ip){
+ for (k = ip[i]; k < ip[i+1]; k++) {
+ ie = jp[k];
+ Vector_add(cedges[jc], (void*) (&ie));
+ }
+ } else {
+ Vector_add(cedges[jc], (void*) (&i));
+ }
+ cinks[jc] = minink;
+ grand_total_ink += minink;
+ grand_total_gain += maxgain;
+
+ if (Verbose && DEBUG){
+ fprintf(stderr," coarse edge[%d]={",jc);
+ for (k = 0; k < Vector_get_length(cedges[jc]); k++) {
+ fprintf(stderr,"%d,", *((int*) Vector_get(cedges[jc], k)));
+ }
+ fprintf(stderr,"}, grand_total_gain=%f\n",grand_total_gain);
+ }
+
+ }
+
+ if (nc >= 1 && total_gain > 0){
+ /* now set up restriction and prolongation operator */
+ SparseMatrix P, R, R1, R0, B, cA;
+ real one = 1.;
+ Agglomerative_Ink_Bundling cgrid;
+
+ R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
+ for (i = 0; i < n; i++){
+ jj = matching[i];
+ SparseMatrix_coordinate_form_add_entries(R1, 1, &jj, &i, &one);
+ }
+ R = SparseMatrix_from_coordinate_format(R1);
+ SparseMatrix_delete(R1);
+ P = SparseMatrix_transpose(R);
+ B = SparseMatrix_multiply(R, A);
+ if (!B) goto RETURN;
+ cA = SparseMatrix_multiply(B, P);
+ if (!cA) goto RETURN;
+ SparseMatrix_delete(B);
+ grid->P = P;
+ grid->R = R;
+
+ level++;
+ cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level);
+
+
+ /* set up R0!!! */
+ if (grid->R0){
+ R0 = SparseMatrix_multiply(R, grid->R0);
+ } else {
+ assert(grid->level == 0);
+ R0 = R;
+ }
+ cgrid->R0 = R0;
+ cgrid->inks = cinks;
+ cgrid->total_ink = grand_total_ink;
+
+ if (Verbose) fprintf(stderr,"level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n", grid->level, cgrid->level, grid->n,
+ cgrid->n, grid->total_ink, grand_total_ink, grid->total_ink - grand_total_ink, grand_total_gain);
+ assert(ABS(grid->total_ink - cgrid->total_ink - grand_total_gain) <= 0.0001*grid->total_ink);
+
+ cgrid = Agglomerative_Ink_Bundling_establish(cgrid, pick, angle_param, angle);
+ grid->next = cgrid;
+ cgrid->prev = grid;
+
+ } else {
+ if (Verbose) fprintf(stderr,"no more improvement, orig ink = %f, gain = %f, stop and final bundling found\n", grand_total_ink, grand_total_gain);
+ /* no more improvement, stop and final bundling found */
+ for (i = 0; i < n; i++) matching[i] = i;
+ }
+
+ RETURN:
+ FREE(matching);
+ for (i = 0; i < n; i++) Vector_delete(cedges[i]);
+ FREE(cedges);
+ FREE(mask);
+ return grid;
+}
+
+
+
+static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_aggressive_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){
+ /* this does a true single-linkage clustering: find the edge that gives the best saving, merge, find again.
+ 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.
+ 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
+ 59% to 62%
+ */
+ /* pick is a work array of dimension n, with n the total number of original edges */
+ int *matching;
+ SparseMatrix A = grid->A;
+ int n = grid->n, level = grid->level, nc = 0;
+ int *ia = A->ia, *ja = A->ja;
+ // real *a;
+ int i, j, k, jj, jc = -1, jmax, imax, ni, nj, npicks;
+ int *mask;
+ pedge *edges = grid->edges;
+ real *inks = grid->inks, *cinks, inki, inkj;
+ real gain, maxgain, minink, total_gain = 0;
+ int *ip = NULL, *jp = NULL, ie;
+ Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid.
+ cedges[i] contain the list of origonal edges that make up the bundle i in the next level */
+ real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
+ point_t meet1, meet2;
+
+ if (Verbose) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n);
+ cedges = MALLOC(sizeof(Vector)*n);
+ cinks = MALLOC(sizeof(real)*n);
+ for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL);
+
+ if (grid->level > 0){
+ ip = grid->R0->ia;
+ jp = grid->R0->ja;
+ }
+
+ matching = MALLOC(sizeof(int)*n);
+ mask = MALLOC(sizeof(real)*n);
+ for (i = 0; i < n; i++) mask[i] = -1;
+
+ assert(n == A->n);
+ for (i = 0; i < n; i++) matching[i] = UNMATCHED;
+
+ //a = (real*) A->a;
+
+ do {
+ maxgain = 0;
+ imax = -1;
+ jmax = -1;
+ minink = -1;
+ for (i = 0; i < n; i++){
+ if (matching[i] != UNMATCHED) continue;
+
+ /* find the best matching in ink saving */
+ for (j = ia[i]; j < ia[i+1]; j++){
+ jj = ja[j];
+ if (jj == i) continue;
+
+ /* ink saving of merging i and j */
+ if ((jc=matching[jj]) == UNMATCHED){
+ /* neither i nor jj are matched */
+ inki = inks[i]; inkj = inks[jj];
+ if (ip && jp){/* not the first level */
+ ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
+ nj = (ip[jj+1] - ip[jj]);/* number of edges represented by jj */
+ MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
+ MEMCPY(pick+ni, &(jp[ip[jj]]), sizeof(int)*nj);
+ } else {/* first level */
+ pick[0] = i; pick[1] = jj;
+ ni = nj = 1;
+ }
+ if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj);
+ } else {
+ /* j is already matched. Its content is on cedges[jc] */
+ inki = inks[i]; inkj = cinks[jc];
+ if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc, inkj);
+ if (ip) {
+ ni = (ip[i+1] - ip[i]);/* number of edges represented by i */
+ MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni);
+ } else {
+ ni = 1; pick[0] = i;
+ }
+ nj = Vector_get_length(cedges[jc]);
+ npicks = ni;
+ for (k = 0; k < nj; k++) {
+ pick[npicks++] = *((int*) Vector_get(cedges[jc], k));
+ }
+ }
+
+ npicks = ni + nj;
+ ink1 = ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle);
+ if (Verbose && DEBUG) {
+ fprintf(stderr,", if merging {");
+ for (k = 0; k < npicks; k++) fprintf(stderr,"%d,", pick[k]);
+ fprintf(stderr,"}, ");
+ fprintf(stderr, " ink0=%f, ink1=%f", inki+inkj, ink1);
+ }
+
+ gain = inki + inkj - ink1;
+ if (Verbose && DEBUG) fprintf(stderr, " gain=%f", gain);
+ if (gain > maxgain){
+ maxgain = gain;
+ minink = ink1;
+ jmax = jj;
+ imax = i;
+ if (Verbose && DEBUG) fprintf(stderr, "maxgain=%f", maxgain);
+ }
+ if (Verbose && DEBUG) fprintf(stderr, "\n");
+
+
+
+ }
+ }
+
+ /* now merge i and jmax */
+ if (maxgain > 0){
+ /* a good bundling of i and another edge jmax is found */
+ total_gain += maxgain;
+ jc = matching[jmax];
+ if (jc == UNMATCHED){/* i and j both unmatched. Add j in the table first */
+ 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);
+ matching[imax] = matching[jmax] = nc;
+ if (ip){
+ for (k = ip[jmax]; k < ip[jmax+1]; k++) {
+ ie = jp[k];
+ Vector_add(cedges[nc], (void*) (&ie));
+ }
+ } else {
+ Vector_add(cedges[nc], (void*) (&jmax));
+ }
+ jc = nc;
+ nc++;
+ } else {/*j is already matched */
+ if (Verbose && DEBUG) printf("maxgain=%f, merge %d with existing cluster %d\n",maxgain, i, jc);
+ matching[imax] = jc;
+ 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 */
+ }
+ /* add i to the appropriate table */
+ if (ip){
+ for (k = ip[imax]; k < ip[imax+1]; k++) {
+ ie = jp[k];
+ Vector_add(cedges[jc], (void*) (&ie));
+ }
+ } else {
+ Vector_add(cedges[jc], (void*) (&imax));
+ }
+ cinks[jc] = minink;
+ grand_total_ink += minink;
+ grand_total_gain += maxgain;
+ } else {/*i can not match/bundle successfully */
+ if (Verbose && DEBUG) fprintf(stderr, "no gain in bundling node %d\n",i);
+ for (i = 0; i < n; i++){
+ assert(maxgain <= 0);
+ if (matching[i] == UNMATCHED){
+ imax = i;
+ matching[imax] = nc;
+ jc = nc;
+ minink = inks[imax];
+ nc++;
+
+ if (ip){
+ for (k = ip[imax]; k < ip[imax+1]; k++) {
+ ie = jp[k];
+ Vector_add(cedges[jc], (void*) (&ie));
+ }
+ } else {
+ Vector_add(cedges[jc], (void*) (&imax));
+ }
+ cinks[jc] = minink;
+ grand_total_ink += minink;
+ grand_total_gain += maxgain;
+
+
+
+
+ }
+ }
+ }
+
+ } while (maxgain > 0);
+
+
+ if (Verbose && DEBUG){
+ fprintf(stderr," coarse edge[%d]={",jc);
+ for (k = 0; k < Vector_get_length(cedges[jc]); k++) {
+ fprintf(stderr,"%d,", *((int*) Vector_get(cedges[jc], k)));
+ }
+ fprintf(stderr,"}\n");
+ }
+
+ if (nc >= 1 && total_gain > 0){
+ /* now set up restriction and prolongation operator */
+ SparseMatrix P, R, R1, R0, B, cA;
+ real one = 1.;
+ Agglomerative_Ink_Bundling cgrid;
+
+ R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
+ for (i = 0; i < n; i++){
+ jj = matching[i];
+ SparseMatrix_coordinate_form_add_entries(R1, 1, &jj, &i, &one);
+ }
+ R = SparseMatrix_from_coordinate_format(R1);
+ SparseMatrix_delete(R1);
+ P = SparseMatrix_transpose(R);
+ B = SparseMatrix_multiply(R, A);
+ if (!B) goto RETURN;
+ cA = SparseMatrix_multiply(B, P);
+ if (!cA) goto RETURN;
+ SparseMatrix_delete(B);
+ grid->P = P;
+ grid->R = R;
+
+ level++;
+ cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level);
+
+
+ /* set up R0!!! */
+ if (grid->R0){
+ R0 = SparseMatrix_multiply(R, grid->R0);
+ } else {
+ assert(grid->level == 0);
+ R0 = R;
+ }
+ cgrid->R0 = R0;
+ cgrid->inks = cinks;
+ cgrid->total_ink = grand_total_ink;
+
+ if (Verbose) fprintf(stderr,"level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n", grid->level, cgrid->level, grid->n,
+ cgrid->n, grid->total_ink, grand_total_ink, grid->total_ink - grand_total_ink, grand_total_gain);
+ assert(ABS(grid->total_ink - cgrid->total_ink - grand_total_gain) <= 0.0001*grid->total_ink);
+
+ cgrid = Agglomerative_Ink_Bundling_aggressive_establish(cgrid, pick, angle_param, angle);
+ grid->next = cgrid;
+ cgrid->prev = grid;
+
+ } else {
+ if (Verbose) fprintf(stderr,"no more improvement, orig ink = %f, gain = %f, stop and final bundling found\n", grand_total_ink, grand_total_gain);
+ /* no more improvement, stop and final bundling found */
+ for (i = 0; i < n; i++) matching[i] = i;
+ }
+
+ RETURN:
+ FREE(matching);
+ for (i = 0; i < n; i++) Vector_delete(cedges[i]);
+ FREE(mask);
+ return grid;
+}
+
+static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_new(SparseMatrix A0, pedge *edges, real angle_param, real angle){
+ /* give a link of edges and their nearest neighbor graph, return a multilevel of edge bundling based on ink saving */
+ Agglomerative_Ink_Bundling grid;
+ int *pick;
+ SparseMatrix A = A0;
+
+ if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){
+ A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
+ }
+ grid = Agglomerative_Ink_Bundling_init(A, edges, 0);
+
+ pick = MALLOC(sizeof(int)*A0->m);
+
+ //grid = Agglomerative_Ink_Bundling_aggressive_establish(grid, pick, angle_param, angle);
+ grid = Agglomerative_Ink_Bundling_establish(grid, pick, angle_param, angle);
+ FREE(pick);
+
+ if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */
+
+ return grid;
+}
+
+static 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){
+
+ int i, j, jj, k;
+ int *ia, *ja;
+ int *pick;
+ Agglomerative_Ink_Bundling grid, cgrid;
+ SparseMatrix R;
+ real ink0, ink1;
+ point_t meet1, meet2;
+ pedge e;
+ real TOL = 0.0001, wgt_all;
+ clock_t start;
+
+ (*recurse_level)++;
+ if (Verbose) fprintf(stderr, "agglomerative_ink_bundling_internal, recurse level ------- %d\n",*recurse_level);
+
+ assert(A->m == A->n);
+
+ *flag = 0;
+
+ start = clock();
+ grid = Agglomerative_Ink_Bundling_new(A, edges, angle_param, angle);
+ if (Verbose)
+ fprintf(stderr, "CPU for agglomerative bundling %f\n", ((real) (clock() - start))/CLOCKS_PER_SEC);
+ ink0 = grid->total_ink;
+
+ /* find coarsest */
+ cgrid = grid;
+ while (cgrid->next){
+ cgrid = cgrid->next;
+ }
+ ink1 = cgrid->total_ink;
+
+ if (*current_ink < 0){
+ *current_ink = *ink00 = ink0;
+ if (Verbose)
+ fprintf(stderr,"initial total ink = %f\n",*current_ink);
+ }
+ if (ink1 < ink0){
+ *current_ink -= ink0 - ink1;
+ }
+
+ if (Verbose) 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));
+
+ /* if no meaningful improvement (0.0001%), out, else rebundle the middle section */
+ if ((ink0-ink1)/(*ink00) < 0.000001 || *recurse_level > MAX_RECURSE_LEVEL) {
+ /* project bundles up */
+ R = cgrid->R0;
+ if (R){
+ ia = R->ia;
+ ja = R->ja;
+ for (i = 0; i < R->m; i++){
+ pick = &(ja[ia[i]]);
+
+ if (Verbose && DEBUG) fprintf(stderr,"calling ink2...\n");
+ ink1 = ink(edges, ia[i+1]-ia[i], pick, &ink0, &meet1, &meet2, angle_param, angle);
+ if (Verbose && DEBUG) fprintf(stderr,"finish calling ink2...\n");
+ assert(ABS(ink1 - cgrid->inks[i])<=MAX(TOL, TOL*ink1) && ink1 - ink0 <= TOL);
+ wgt_all = 0.;
+ if (ia[i+1]-ia[i] > 1){
+ for (j = ia[i]; j < ia[i+1]; j++){
+ /* make this edge 4 points, insert two meeting points at 1 and 2, make 3 the last point */
+ jj = ja[j];
+ 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
+ improved later */
+ e = edges[jj] = pedge_double(edges[jj]);
+
+ e->wgts = REALLOC(e->wgts, sizeof(real)*4);
+ e->x[1*dim] = meet1.x;
+ e->x[1*dim+1] = meet1.y;
+ e->x[2*dim] = meet2.x;
+ e->x[2*dim+1] = meet2.y;
+ e->x[3*dim] = e->x[4*dim];
+ e->x[3*dim+1] = e->x[4*dim+1];
+ e->npoints = 4;
+ for (k = 0; k < 3; k++) e->wgts[k] = e->wgt;
+ wgt_all += e->wgt;
+
+ }
+ for (j = ia[i]; j < ia[i+1]; j++){
+ e = edges[ja[j]];
+ e->wgts[1] = wgt_all;
+ }
+ }
+
+ }
+ }
+ } else {
+ pedge *mid_edges, midedge;/* middle section of edges that will be bundled again */
+ real *xx;
+ int ne, npp, l;
+ SparseMatrix A_mid;
+ real eps = 0., wgt, total_wgt = 0;
+
+ /* make new edges using meet1 and meet2.
+
+ call Agglomerative_Ink_Bundling_new
+
+ inherit new edges to old edges
+ */
+
+ R = cgrid->R0;
+ assert(R && cgrid != grid);/* if ink improved, we should have gone at leat 1 level down! */
+ ia = R->ia;
+ ja = R->ja;
+ ne = R->m;
+ mid_edges = MALLOC(sizeof(pedge)*ne);
+ xx = MALLOC(sizeof(real)*4*ne);
+ for (i = 0; i < R->m; i++){
+ pick = &(ja[ia[i]]);
+ wgt = 0.;
+ for (j = ia[i]; j < ia[i+1]; j++) wgt += edges[j]->wgt;
+ total_wgt += wgt;
+ if (Verbose && DEBUG) fprintf(stderr,"calling ink3...\n");
+ ink1 = ink(edges, ia[i+1]-ia[i], pick, &ink0, &meet1, &meet2, angle_param, angle);
+ if (Verbose && DEBUG) fprintf(stderr,"done calling ink3...\n");
+ assert(ABS(ink1 - cgrid->inks[i])<=MAX(TOL, TOL*ink1) && ink1 - ink0 <= TOL);
+ xx[i*4 + 0] = meet1.x;
+ xx[i*4 + 1] = meet1.y;
+ xx[i*4 + 2] = meet2.x;
+ xx[i*4 + 3] = meet2.y;
+ mid_edges[i] = pedge_wgt_new(2, dim, &(xx[i*4]), wgt);
+ //mid_edges[i] = pedge_wgt_new(2, dim, &(xx[i*4]), 1.);
+
+ }
+
+ A_mid = nearest_neighbor_graph(ne, MIN(nneighbors, ne), 4, xx, eps);
+
+ agglomerative_ink_bundling_internal(dim, A_mid, mid_edges, nneighbors, recurse_level, MAX_RECURSE_LEVEL, angle_param, angle, open_gl, current_ink, ink00, flag);
+ SparseMatrix_delete(A_mid);
+ FREE(xx);
+
+ /* patching edges with the new mid-section */
+ for (i = 0; i < R->m; i++){
+ pick = &(ja[ia[i]]);
+ midedge = mid_edges[i];
+ npp = midedge->npoints + 2;
+ for (j = ia[i]; j < ia[i+1]; j++){
+ jj = ja[j];
+ e = edges[jj] = pedge_wgts_realloc(edges[jj], npp);
+
+ assert(e->npoints = 2);
+ for (l = 0; l < dim; l++){/* move the second point to the last */
+ e->x[(npp - 1)*dim+l] = e->x[1*dim+l];
+ }
+
+ for (k = 0; k < midedge->npoints; k++){
+ for (l = 0; l < dim; l++){
+ e->x[(k+1)*dim+l] = midedge->x[k*dim+l];
+ }
+ if (k < midedge->npoints - 1){
+ if (midedge->wgts){
+ e->wgts[(k+1)] = midedge->wgts[k];
+ } else {
+ e->wgts[(k+1)] = midedge->wgt;
+ }
+ }
+ }
+ e->wgts[npp - 2] = e->wgts[0];/* the last interval take from the 1st interval */
+
+
+ e->npoints = npp;
+ }
+#ifdef OPENGL
+ if (open_gl & FALSE){
+ nedges_global = grid->n;
+ edges_global = edges;
+ drawScene();
+ // waitie(5./R->m);
+ }
+#endif
+ }
+
+ for (i = 0; i < ne; i++) pedge_delete(mid_edges[i]);
+
+ }
+
+
+#ifdef OPENGL
+ if (open_gl){
+ nedges_global = grid->n;
+ edges_global = edges;
+ drawScene();
+ // waitie(5./R->m);
+ }
+#endif
+
+ Agglomerative_Ink_Bundling_delete(grid);
+
+ return edges;
+}
+
+
+pedge* 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){
+ int recurse_level = 0;
+ real current_ink = -1, ink0;
+ pedge *edges2;
+
+ ink_count = 0;
+ edges2 = agglomerative_ink_bundling_internal(dim, A, edges, nneighbor, &recurse_level, MAX_RECURSE_LEVEL, angle_param, angle, open_gl, ¤t_ink, &ink0, flag);
+
+
+ if (Verbose)
+ 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);
+ return edges2;
+}
+
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifndef AGGLOMERATIVE_BUNDLING_H
+#define AGGLOMERATIVE_BUNDLING_H
+
+typedef struct Agglomerative_Ink_Bundling_struct *Agglomerative_Ink_Bundling;
+
+struct Agglomerative_Ink_Bundling_struct {
+ int level;/* 0, 1, ... */
+ int n;
+ SparseMatrix A; /* n x n matrix, where n is the number of edges/bundles in this level */
+ SparseMatrix P; /* prolongation matrix from level + 1 to level */
+ SparseMatrix R0;/* this is basically R[level - 1].R[level - 2]...R[0], which gives the map of bundling i to the original edges: first row of R0 gives
+ the nodes on the finest grid corresponding to the coarsest node 1, etc */
+ SparseMatrix R;/* striction mtrix from level to level + 1*/
+ Agglomerative_Ink_Bundling next;
+ Agglomerative_Ink_Bundling prev;
+ real *inks; /* amount of ink needed to draw this edge/bundle. Dimension n. */
+ real total_ink; /* amount of ink needed to draw this edge/bundle. Dimension n. */
+ pedge* edges; /* the original edge info. This does not vary level to level and is of dimenion n0, where n0 is the number of original edges */
+ int delete_top_level_A;/*whether the top level matrix should be deleted on garbage collecting the grid */
+};
+
+pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int max_recursion, real angle_param, real angle, int open_gl, int *flag);
+
+#endif /* AGGLOMERATIVE_BUNDLING_H */
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "types.h"
+#include "globals.h"
+#include "general.h"
+#include "SparseMatrix.h"
+#include "edge_bundling.h"
+#include <time.h>
+#include "clustering.h"
+#include "ink.h"
+#include "agglomerative_bundling.h"
+
+#define SMALL 1.e-10
+
+#ifdef OPENGL
+#include "gl.h"
+extern pedge *edges_global;
+extern int *clusters_global;
+#endif
+
+static real norm(int n, real *x){
+ real res = 0;
+ int i;
+
+ for (i = 0; i < n; i++) res += x[i]*x[i];
+ return sqrt(res);
+
+}
+static real sqr_dist(int dim, real *x, real *y){
+ int i;
+ real res = 0;
+ for (i = 0; i < dim; i++) res += (x[i] - y[i])*(x[i] - y[i]);
+ return res;
+}
+static real dist(int dim, real *x, real *y){
+ return sqrt(sqr_dist(dim,x,y));
+}
+
+pedge pedge_new(int np, int dim, real *x){
+ pedge e;
+
+ e = MALLOC(sizeof(struct pedge_struct));
+ e->npoints = np;
+ e->dim = dim;
+ e->len = np;
+ e->x = MALLOC(dim*(e->len)*sizeof(real));
+ MEMCPY(e->x, x, dim*(e->len)*sizeof(real));
+ e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim]));
+ e->wgt = 1.;
+ e->wgts = NULL;
+ return e;
+
+}
+pedge pedge_wgt_new(int np, int dim, real *x, real wgt){
+ pedge e;
+ int i;
+
+ e = MALLOC(sizeof(struct pedge_struct));
+ e->npoints = np;
+ e->dim = dim;
+ e->len = np;
+ e->x = MALLOC(dim*(e->len)*sizeof(real));
+ MEMCPY(e->x, x, dim*(e->len)*sizeof(real));
+ e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim]));
+ e->wgt = wgt;
+ e->wgts = MALLOC(sizeof(real)*(np - 1));
+ for (i = 0; i < np - 1; i++) e->wgts[i] = wgt;
+ return e;
+
+}
+void pedge_delete(pedge e){
+ FREE(e->x);
+ FREE(e);
+}
+
+pedge pedge_flip(pedge e){
+ /* flip the polyline so that last point becomes the first, second last the second, etc*/
+ real *y;
+ real *x = e->x;
+ int i, dim = e->dim;
+ int n = e->npoints;
+
+ y = MALLOC(sizeof(real)*e->dim);
+ for (i = 0; i < (e->npoints)/2; i++){
+ MEMCPY(y, &(x[i*dim]), sizeof(real)*dim);
+ MEMCPY(&(x[(n-1-i)*dim]), &(x[i*dim]), sizeof(real)*dim);
+ MEMCPY(&(x[i*dim]), y, sizeof(real)*dim);
+ }
+ FREE(y);
+ return e;
+}
+
+static real edge_compatibility(pedge e1, pedge e2){
+ /* two edges are u1->v1, u2->v2.
+ return 1 if two edges are exactly the same, 0 if they are very different.
+ */
+ real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2;
+ int dim = e1->dim, flipped = FALSE;
+
+ u1 = e1->x;
+ v1 = (e1->x)+(e1->npoints)*dim-dim;
+ u2 = e2->x;
+ v2 = (e2->x)+(e2->npoints)*dim-dim;
+ dist1 = sqr_dist(dim, u1, u2) + sqr_dist(dim, v1, v2);
+ dist2 = sqr_dist(dim, u1, v2) + sqr_dist(dim, v1, u2);
+ if (dist1 > dist2){
+ u = u2;
+ u2 = v2;
+ v2 = u;
+ dist1 = dist2;
+ flipped = TRUE;
+ }
+ len1 = dist(dim, u1, v1);
+ len2 = dist(dim, u2, v2);
+ //dist1 = MAX(0.1, dist1/(len1+len2+dist1));
+ dist1 = MAX(0.1, dist1/(len1+len2+0.0001*dist1));
+ if (flipped){
+ return -1/dist1;
+ } else {
+ return 1/dist1;
+ }
+}
+
+static real edge_compatibility_full(pedge e1, pedge e2){
+ /* two edges are u1->v1, u2->v2.
+ return 1 if two edges are exactly the same, 0 if they are very different.
+ This is based on Holten and van Wijk's paper
+ */
+ real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2, len;
+ real tmp, ca, cp, cs;
+ int dim = e1->dim, flipped = FALSE, i;
+
+ u1 = e1->x;
+ v1 = (e1->x)+(e1->npoints)*dim-dim;
+ u2 = e2->x;
+ v2 = (e2->x)+(e2->npoints)*dim-dim;
+ dist1 = sqr_dist(dim, u1, u2) + sqr_dist(dim, v1, v2);
+ dist2 = sqr_dist(dim, u1, v2) + sqr_dist(dim, v1, u2);
+ if (dist1 > dist2){
+ u = u2;
+ u2 = v2;
+ v2 = u;
+ dist1 = dist2;
+ flipped = TRUE;
+ }
+ len1 = MAX(dist(dim, u1, v1), SMALL);
+ len2 = MAX(dist(dim, u2, v2), SMALL);
+ len = 0.5*(len1+len2);
+
+ /* angle compatability */
+ ca = 0;
+ for (i = 0; i < dim; i++)
+ ca += (v1[i]-u1[i])*(v2[i]-u2[i]);
+ ca = ABS(ca/(len1*len2));
+ assert(ca > -0.001);
+
+ /* scale compatability */
+ //cs = 2/(len1/len2+len2/len1);
+ cs = 2/(MAX(len1,len2)/len + len/MIN(len1, len2));
+ assert(cs > -0.001 && cs < 1.001);
+
+ /* position compatability */
+ cp = 0;
+ for (i = 0; i < dim; i++) {
+ tmp = .5*(v1[i]+u1[i])-.5*(v2[i]+u2[i]);
+ cp += tmp*tmp;
+ }
+ cp = sqrt(cp);
+ cp = len/(len + cp);
+ assert(cp > -0.001 && cp < 1.001);
+
+ /* visibility compatability */
+
+ //dist1 = MAX(0.1, dist1/(len1+len2+dist1));
+ dist1 = cp*ca*cs;
+ if (flipped){
+ return -dist1;
+ } else {
+ return dist1;
+ }
+}
+
+static void fprint_rgb(FILE* fp, int r, int g, int b, int alpha){
+ fprintf(fp,"#");
+ if (r >= 16) {
+ fprintf(fp,"%2x",r);
+ } else {
+ fprintf(fp,"0%1x",r);
+ }
+ if (g >= 16) {
+ fprintf(fp,"%2x",g);
+ } else {
+ fprintf(fp,"0%1x",g);
+ }
+ if (b >= 16) {
+ fprintf(fp,"%2x",b);
+ } else {
+ fprintf(fp,"0%1x",b);
+ }
+ if (alpha >= 16) {
+ fprintf(fp,"%2x",alpha);
+ } else {
+ fprintf(fp,"0%1x",alpha);
+ }
+
+}
+
+void pedge_export_gv(FILE *fp, int ne, pedge *edges){
+ pedge edge;
+ real *x, t;
+ int i, j, k, kk, dim, sta, sto;
+ real maxwgt = 0, len, len_total, len_total0;
+ int r, g, b;
+ // real tt1[3]={0.25,0.5,0.75};
+ // real tt2[4]={0.2,0.4,0.6,0.8};
+ real tt1[3]={0.15,0.5,0.85};
+ real tt2[4]={0.15,0.4,0.6,0.85};
+ real *tt;
+
+ fprintf(fp,"strict graph{\n");
+ /* points */
+ for (i = 0; i < ne; i++){
+ edge = edges[i];
+ x = edge->x;
+ dim = edge->dim;
+ sta = 0; sto = edge->npoints - 1;
+
+ fprintf(fp, "%d [pos=\"", i);
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(fp, ",");
+ fprintf(fp, "%f", x[sta*dim+k]);
+ }
+ fprintf(fp, "\"];\n");
+
+ fprintf(fp, "%d [pos=\"", i + ne);
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(fp, ",");
+ fprintf(fp, "%f", x[sto*dim+k]);
+ }
+ fprintf(fp, "\"];\n");
+
+ }
+
+ /* figure out max number of bundled origional edges in a pedge */
+ for (i = 0; i < ne; i++){
+ edge = edges[i];
+ if (edge->wgts){
+ for (j = 0; j < edge->npoints - 1; j++){
+ maxwgt = MAX(maxwgt, edge->wgts[j]);
+ }
+ }
+ }
+
+ /* spline and colors */
+ for (i = 0; i < ne; i++){
+ fprintf(fp,"%d -- %d [pos=\"", i, i + ne);
+ edge = edges[i];
+ x = edge->x;
+ dim = edge->dim;
+ /* splines */
+ for (j = 0; j < edge->npoints; j++){
+ if (j != 0) {
+ int mm = 3;
+ fprintf(fp," ");
+ /* there are ninterval+1 points, add 3*ninterval+2 points, get rid of internal ninternal-1 points,
+ make into 3*ninterval+4 points so that gviz spline rendering can work */
+ if (j == 1 || j == edge->npoints - 1) {
+ /* every interval gets 3 points inmserted except the first and last one */
+ tt = tt2;
+ mm = 4;
+ } else {
+ tt = tt1;
+ }
+ for (kk = 1; kk <= mm; kk++){
+ //t = kk/(real) (mm+1);
+ t = tt[kk-1];
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(fp,",");
+ fprintf(fp, "%f", (x[(j-1)*dim+k]*(1-t)+x[j*dim+k]*(t)));
+ }
+ fprintf(fp," ");
+ }
+ }
+ if (j == 0 || j == edge->npoints - 1){
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(fp,",");
+ fprintf(fp, "%f", x[j*dim+k]);
+ }
+ }
+ }
+ /* colors based on how much bundling */
+ if (edge->wgts){
+ fprintf(fp, "\", wgts=\"");
+ for (j = 0; j < edge->npoints - 1; j++){
+ if (j != 0) fprintf(fp,",");
+ fprintf(fp, "%f", edge->wgts[j]);
+ }
+
+ len_total = 0;
+ len_total0 = 0;
+ fprintf(fp, "\", color=\"");
+ for (j = 0; j < edge->npoints - 1; j++){
+ len = 0;
+ for (k = 0; k < dim; k++){
+ len += (edge->x[dim*j+k] - edge->x[dim*(j+1)+k])*(edge->x[dim*j+k] - edge->x[dim*(j+1)+k]);
+ }
+ len = sqrt(len/k);
+ len_total0 += len;
+ }
+ for (j = 0; j < edge->npoints - 1; j++){
+ len = 0;
+ for (k = 0; k < dim; k++){
+ len += (edge->x[dim*j+k] - edge->x[dim*(j+1)+k])*(edge->x[dim*j+k] - edge->x[dim*(j+1)+k]);
+ }
+ len = sqrt(len/k);
+ len_total += len;
+ t = edge->wgts[j]/maxwgt;
+ /* interpotate between red (t = 1) to blue (t = 0) */
+ r = 255*t; g = 0; b = 255*(1-t); b = 255*(1-t);
+ if (j != 0) fprintf(fp,":");
+ fprint_rgb(fp, r, g, b, 85);
+ fprintf(fp,",%f",len_total/len_total0);
+ }
+
+ }
+ fprintf(fp, "\"];\n");
+
+ }
+ fprintf(fp,"}\n");
+}
+
+void pedge_export_mma(FILE *fp, int ne, pedge *edges){
+ pedge edge;
+ real *x;
+ int i, j, k, dim;
+
+ fprintf(fp,"Graphics[{");
+ /* points */
+ fprintf(fp,"{Red, ");
+ for (i = 0; i < ne; i++){
+ if (i != 0) fprintf(fp,",");
+ fprintf(fp,"Point[");
+ edge = edges[i];
+ x = edge->x;
+ dim = edge->dim;
+ fprintf(fp, "{");
+ for (j = 0; j < edge->npoints; j+= edge->npoints - 1){
+ if (j != 0) fprintf(fp,",");
+ fprintf(fp, "{");
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(fp,",");
+ fprintf(fp, "%f", x[j*dim+k]);
+ }
+ fprintf(fp, "}");
+ }
+ fprintf(fp, "}");
+ fprintf(fp, "]");
+ }
+ fprintf(fp,"},\n{GrayLevel[0.5,0.2], ");
+
+ /* spline */
+ for (i = 0; i < ne; i++){
+ if (i != 0) fprintf(fp,",");
+ fprintf(fp,"Spline[");
+ edge = edges[i];
+ x = edge->x;
+ dim = edge->dim;
+ fprintf(fp, "{");
+ for (j = 0; j < edge->npoints; j++){
+ if (j != 0) fprintf(fp,",");
+ fprintf(fp, "{");
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(fp,",");
+ fprintf(fp, "%f", x[j*dim+k]);
+ }
+ fprintf(fp, "}");
+ }
+ fprintf(fp, "}");
+ fprintf(fp, ", Bezier]");
+ }
+ fprintf(fp,"}");
+
+ fprintf(fp,"}]\n");
+}
+
+void pedge_print(char *comments, pedge e){
+ int i, j, dim;
+ dim = e->dim;
+ fprintf(stderr,"%s", comments);
+ for (i = 0; i < e->npoints; i++){
+ if (i > 0) fprintf(stderr,",");
+ fprintf(stderr,"{");
+ for (j = 0; j < dim; j++){
+ if (j > 0) fprintf(stderr,",");
+ fprintf(stderr,"%f",e->x[dim*i+j]);
+ }
+ fprintf(stderr,"}");
+ }
+ fprintf(stderr,"\n");
+}
+
+pedge pedge_realloc(pedge e, int n){
+ if (n <= e->npoints) return e;
+ e->x = REALLOC(e->x, e->dim*n*sizeof(real));
+ if (e->wgts) e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
+ e->len = n;
+ return e;
+}
+pedge pedge_wgts_realloc(pedge e, int n){
+ /* diff from pedge_alloc: allocate wgts if do not exist and initialize to wgt */
+ int i;
+ if (n <= e->npoints) return e;
+ e->x = REALLOC(e->x, e->dim*n*sizeof(real));
+ if (!(e->wgts)){
+ e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
+ for (i = 0; i < e->npoints; i++) e->wgts[i] = e->wgt;
+ } else {
+ e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
+ }
+ e->len = n;
+ return e;
+}
+
+
+pedge pedge_double(pedge e){
+ /* double the number of points (more precisely, add a point between two points in the polyline */
+ int npoints = e->npoints, len = e->len, i, dim = e->dim;
+ real *x;
+ int j, ii, ii2, np;
+
+ assert(npoints >= 2);
+ // pedge_print("before doubling edge = ", e);
+ if (npoints*2-1 > len){
+ len = 3*npoints;
+ e->x = REALLOC(e->x, dim*len*sizeof(real));
+ }
+ x = e->x;
+
+ for (i = npoints - 1; i >= 0; i--){
+ ii = 2*i;
+ for (j = 0; j < dim; j++){
+ x[dim*ii + j] = x[dim*i + j];
+ }
+ }
+
+ for (i = 0; i < npoints - 1; i++){
+ ii = 2*i;/* left and right interpolant of a new point */
+ ii2 = 2*(i+1);
+ for (j = 0; j < dim; j++){
+ x[dim*(2*i + 1) + j] = 0.5*(x[dim*ii + j] + x[dim*ii2 + j]);
+ }
+ }
+ e->len = len;
+ np = e->npoints = 2*(e->npoints) - 1;
+ e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim]));
+
+ //pedge_print("after doubling edge = ", e);
+
+ return e;
+}
+
+static void edge_tension_force(real *force, pedge e){
+ real *x = e->x;
+ int dim = e->dim;
+ int np = e->npoints;
+ int i, left, right, j;
+ //real dist_left, dist_right;
+ real s;
+
+
+ /* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless
+ */
+ //s = (np-1)*(np-1)/MAX(SMALL, e->edge_length);
+ s = (np-1)/MAX(SMALL, e->edge_length);
+ for (i = 1; i <= np - 2; i++){
+ left = i - 1;
+ right = i + 1;
+ // dist_left = dist(dim, &(x[i*dim]), &(x[left*dim]));
+ // dist_right = dist(dim, &(x[i*dim]), &(x[right*dim]));
+ for (j = 0; j < dim; j++) force[i*dim + j] += s*(x[left*dim + j] - x[i*dim + j]);
+ for (j = 0; j < dim; j++) force[i*dim + j] += s*(x[right*dim + j] - x[i*dim + j]);
+ }
+}
+
+
+#if 0
+static void edge_tension_force2(real *force, pedge e){
+ /* in this version each node is pulled towards its original position on the line */
+ real *x = e->x;
+ int dim = e->dim;
+ int np = e->npoints;
+ int i, j;
+ //int left, right;
+ //real dist_left, dist_right;
+ real s;
+
+
+ /* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless
+ */
+ s = .1/MAX(SMALL, e->edge_length);
+ for (i = 1; i <= np - 2; i++){
+ //left = i - 1;
+ // right = i + 1;
+ // dist_left = dist(dim, &(x[i*dim]), &(x[left*dim]));
+ // dist_right = dist(dim, &(x[i*dim]), &(x[right*dim]));
+ for (j = 0; j < dim; j++) force[i*dim + j] += s*((i*x[0*dim + j]+(np-1-i)*x[(np-1)*dim + j])/(np-1) - x[i*dim + j]);
+ }
+}
+#endif
+
+static void edge_attraction_force(real similarity, pedge e1, pedge e2, real *force){
+ /* attrractive force from x2 applied to x1 */
+ real *x1 = e1->x, *x2 = e2->x;
+ int dim = e1->dim;
+ int np = e1->npoints;
+ int i, j;
+ real dist, s, ss;
+ real edge_length = e1->edge_length;
+
+
+ assert(e1->npoints == e2->npoints);
+
+ /* attractive force = 1/d where d = D/||e1|| is the relative distance, D is the distance between e1 and e2.
+ so the force is norminal and unitless
+ */
+ if (similarity > 0){
+ s = edge_length;
+ s = similarity*edge_length;
+ for (i = 1; i <= np - 2; i++){
+ dist = sqr_dist(dim, &(x1[i*dim]), &(x2[i*dim]));
+ if (dist < SMALL) dist = SMALL;
+ ss = s/(dist+0.1*edge_length*sqrt(dist));
+ for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[i*dim + j] - x1[i*dim + j]);
+ }
+ } else {/* clip e2 */
+ s = -edge_length;
+ s = -similarity*edge_length;
+ for (i = 1; i <= np - 2; i++){
+ dist = sqr_dist(dim, &(x1[i*dim]), &(x2[(np - 1 - i)*dim]));
+ if (dist < SMALL) dist = SMALL;
+ ss = s/(dist+0.1*edge_length*sqrt(dist));
+ for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[(np - 1 - i)*dim + j] - x1[i*dim + j]);
+ }
+ }
+
+}
+
+static pedge* force_directed_edge_bundling(SparseMatrix A, pedge* edges, int maxit, real step0, real K, int open_gl){
+ int i, j, ne = A->n, k;
+ int *ia = A->ia, *ja = A->ja, iter = 0;
+ real *a = (real*) A->a;
+ pedge e1, e2;
+ real *force_t, *force_a;
+ int np = edges[0]->npoints, dim = edges[0]->dim;
+ real *x;
+ real step = step0;
+ real fnorm_a, fnorm_t, edge_length, start;
+
+ if (Verbose)
+ fprintf(stderr, "total interaction pairs = %d out of %d, avg neighbors per edge = %f\n",A->nz, A->m*A->m, A->nz/(real) A->m);
+
+ force_t = MALLOC(sizeof(real)*dim*(np));
+ force_a = MALLOC(sizeof(real)*dim*(np));
+ while (step > 0.001 && iter < maxit){
+ start = clock();
+ iter++;
+ for (i = 0; i < ne; i++){
+ for (j = 0; j < dim*np; j++) {
+ force_t[j] = 0.;
+ force_a[j] = 0.;
+ }
+ e1 = edges[i];
+ x = e1->x;
+ //edge_tension_force2(force_t, e1);
+ edge_tension_force(force_t, e1);
+ for (j = ia[i]; j < ia[i+1]; j++){
+ e2 = edges[ja[j]];
+ edge_attraction_force(a[j], e1, e2, force_a);
+ }
+ fnorm_t = MAX(SMALL, norm(dim*(np - 2), &(force_t[1*dim])));
+ fnorm_a = MAX(SMALL, norm(dim*(np - 2), &(force_a[1*dim])));
+ edge_length = e1->edge_length;
+
+ // fprintf(stderr,"tension force norm = %f att force norm = %f step = %f edge_length = %f\n",fnorm_t, fnorm_a, step, edge_length);
+ for (j = 1; j <= np - 2; j++){
+ for (k = 0; k < dim; k++) {
+ x[j*dim + k] += step*edge_length*(force_t[j*dim+k] + K*force_a[j*dim+k])/(sqrt(fnorm_t*fnorm_t + K*K*fnorm_a*fnorm_a));
+ }
+ }
+
+ /*
+ fprintf(stderr,"edge %d",i);
+ for (j = 0; j < np; j++){
+ if (j != 0) fprintf(stderr,",");
+ fprintf(stderr,"{");
+ for (k = 0; k < dim; k++) {
+ if (k != 0) fprintf(stderr,",");
+ fprintf(stderr,"%f", x[j*dim+k]);
+ }
+ fprintf(stderr,"}");
+ }
+ fprintf(stderr,"\n");
+ */
+
+ }
+ step = step*0.9;
+ fprintf(stderr, "iter ==== %d cpu = %f npoints = %d\n",iter, ((real) (clock() - start))/CLOCKS_PER_SEC, np - 2);
+
+#ifdef OPENGL
+ if (open_gl){
+ edges_global = edges;
+ drawScene();
+ }
+#endif
+
+ }
+
+ FREE(force_t);
+ FREE(force_a);
+ return edges;
+}
+
+static real absfun(real x){
+ return ABS(x);
+}
+
+
+
+
+pedge* modularity_ink_bundling(int dim, int ne, SparseMatrix B, pedge* edges, real angle_param, real angle){
+ int *assignment = NULL, flag, nclusters;
+ real modularity;
+ int *clusterp, *clusters;
+ SparseMatrix D, C;
+ point_t meet1, meet2;
+ real ink0, ink1;
+ pedge e;
+ int i, j, jj;
+ int use_value_for_clustering = TRUE;
+
+ //int use_value_for_clustering = FALSE;
+
+ SparseMatrix BB;
+
+ /* B may contain negative entries */
+ BB = SparseMatrix_copy(B);
+ BB = SparseMatrix_apply_fun(BB, absfun);
+ modularity_clustering(BB, TRUE, 0, use_value_for_clustering, &nclusters, &assignment, &modularity, &flag);
+ SparseMatrix_delete(BB);
+
+#ifdef OPENGL
+ clusters_global = assignment;
+#endif
+
+ assert(!flag);
+ if (Verbose) fprintf(stderr, "there are %d clusters, modularity = %f\n",nclusters, modularity);
+
+ C = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_PATTERN, FORMAT_COORD);
+
+ for (i = 0; i < ne; i++){
+ jj = assignment[i];
+ SparseMatrix_coordinate_form_add_entries(C, 1, &jj, &i, NULL);
+ }
+
+ D = SparseMatrix_from_coordinate_format(C);
+ SparseMatrix_delete(C);
+ clusterp = D->ia;
+ clusters = D->ja;
+ for (i = 0; i < nclusters; i++){
+ ink1 = ink(edges, clusterp[i+1] - clusterp[i], &(clusters[clusterp[i]]), &ink0, &meet1, &meet2, angle_param, angle);
+ fprintf(stderr,"nedges = %d ink0 = %f, ink1 = %f\n",clusterp[i+1] - clusterp[i], ink0, ink1);
+ if (ink1 < ink0){
+ for (j = clusterp[i]; j < clusterp[i+1]; j++){
+ /* make this edge 5 points, insert two meeting points at 1 and 2, make 3 the last point */
+ edges[clusters[j]] = pedge_double(edges[clusters[j]]);
+ e = edges[clusters[j]] = pedge_double(edges[clusters[j]]);
+ e->x[1*dim] = meet1.x;
+ e->x[1*dim+1] = meet1.y;
+ e->x[2*dim] = meet2.x;
+ e->x[2*dim+1] = meet2.y;
+ e->x[3*dim] = e->x[4*dim];
+ e->x[3*dim+1] = e->x[4*dim+1];
+ e->npoints = 4;
+ }
+#ifdef OPENGL
+ edges_global = edges;
+ drawScene();
+#endif
+ }
+ }
+ SparseMatrix_delete(D);
+ return edges;
+}
+
+static SparseMatrix check_compatibility(SparseMatrix A, int ne, pedge *edges, int compatibility_method, real tol){
+ /* go through the links and make sure edges are compatable */
+ SparseMatrix B, C;
+ int *ia, *ja, i, j, jj;
+ real start;
+ real dist;
+
+ B = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
+ ia = A->ia; ja = A->ja;
+ //SparseMatrix_print("A",A);
+ start = clock();
+ for (i = 0; i < ne; i++){
+ for (j = ia[i]; j < ia[i+1]; j++){
+ jj = ja[j];
+ if (i == jj) continue;
+ if (compatibility_method == COMPATIBILITY_DIST){
+ dist = edge_compatibility_full(edges[i], edges[jj]);
+ } else if (compatibility_method == COMPATIBILITY_FULL){
+ dist = edge_compatibility(edges[i], edges[jj]);
+ }
+ /*
+ fprintf(stderr,"edge %d = {",i);
+ pedge_print("", edges[i]);
+ fprintf(stderr,"edge %d = {",jj);
+ pedge_print("", edges[jj]);
+ fprintf(stderr,"compatibility=%f\n",dist);
+ */
+
+ if (ABS(dist) > tol){
+ B = SparseMatrix_coordinate_form_add_entries(B, 1, &i, &jj, &dist);
+ B = SparseMatrix_coordinate_form_add_entries(B, 1, &jj, &i, &dist);
+ }
+ }
+ }
+ C = SparseMatrix_from_coordinate_format(B);
+ //SparseMatrix_print("C",C);
+ SparseMatrix_delete(B);
+ B = C;
+ fprintf(stderr, "edge compatibilitu time = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC);
+ return B;
+}
+
+pedge* edge_bundling(SparseMatrix A0, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method,
+ int max_recursion, real angle_param, real angle, int open_gl){
+ /* bundle edges.
+ A: edge graph
+ x: edge i is at {p,q},
+ . where p = x[2*dim*i : 2*dim*i+dim-1]
+ . and q = x[2*dim*i+dim : 2*dim*i+2*dim-1]
+ maxit_outer: max outer iteration for force directed bundling. Every outer iter subdivide each edge segment into 2.
+ K: norminal edge length in force directed bundling
+ method: which method to use.
+ nneighbor: number of neighbors to be used in forming nearest neighbor graph. Used only in agglomerative method
+ compatibility_method: which method to use to calculate compatibility. Used only in force directed.
+ max_recursion: used only in agglomerative method. Specify how many level of recursion to do to bundle bundled edges again
+ open_gl: whether to plot in X.
+
+ */
+ int ne = A0->m;
+ pedge *edges;
+ SparseMatrix A = A0, B = NULL;
+ int i;
+ real tol = 0.001;
+ int k;
+ real step0 = 0.1, start;
+ int maxit = 10;
+ int flag;
+
+ assert(A->n == ne);
+ edges = MALLOC(sizeof(pedge)*ne);
+
+ for (i = 0; i < ne; i++){
+ edges[i] = pedge_new(2, dim, &(x[dim*2*i]));
+ }
+
+ A = SparseMatrix_symmetrize(A0, TRUE);
+
+
+
+ if (Verbose) start = clock();
+#ifndef HAVE_ANN
+ method = METHOD_INK;
+#endif
+ if (method == METHOD_INK){
+
+ /* go through the links and make sure edges are compatable */
+ B = check_compatibility(A, ne, edges, compatibility_method, tol);
+
+ edges = modularity_ink_bundling(dim, ne, B, edges, angle_param, angle);
+
+ } else if (method == METHOD_INK_AGGLOMERATE){
+#ifdef HAVE_ANN
+ /* plan: merge a node with its neighbors if doing so improve. Form coarsening graph, repeat until no more ink saving */
+ edges = agglomerative_ink_bundling(dim, A, edges, nneighbor, max_recursion, angle_param, angle, open_gl, &flag);
+ assert(!flag);
+#else
+ agerr (AGERR, "Graphviz built without approximate nearest neighbor library ANN; agglomerative inking not available\n");
+ edges = edges;
+#endif
+ } else if (method == METHOD_FD){/* FD method */
+
+ /* go through the links and make sure edges are compatable */
+ B = check_compatibility(A, ne, edges, compatibility_method, tol);
+
+
+ for (k = 0; k < maxit_outer; k++){
+ for (i = 0; i < ne; i++){
+ edges[i] = pedge_double(edges[i]);
+ }
+ step0 = step0/2;
+ edges = force_directed_edge_bundling(B, edges, maxit, step0, K, open_gl);
+ }
+
+ } else if (method == METHOD_NONE){
+ edges = edges;
+ } else {
+ assert(0);
+ }
+ if (Verbose)
+ fprintf(stderr, "total edge bundling cpu = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC);
+ if (B != A) SparseMatrix_delete(B);
+ if (A != A0) SparseMatrix_delete(A);
+
+ return edges;
+}
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+#ifndef EDGE_BUNDLING_H
+#define EDGE_BUNDLING_H
+
+#include <SparseMatrix.h>
+
+struct pedge_struct {
+ real wgt; /* weight, telling how many original edges this edge represent. If this edge consists of multiple sections of different weights then this is a lower bound. This only applied for agglomerative bundling */
+ int npoints;/* number of poly points */
+ int len;/* length of arra x. len >= npoints */
+ int dim;/* dim >= 2. Point i is stored from x[i*dim]*/
+ real edge_length;
+ real *x;/* coordinates of the npoints poly points. Dimension dim*npoints */
+ real *wgts;/* number of original edges each section represnet. Dimension npoint - 1. This only applied for agglomerative bundling Null for other methods */
+};
+
+typedef struct pedge_struct* pedge;
+
+pedge* edge_bundling(SparseMatrix A, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method, int max_recursion, real angle_param, real angle, int open_gl);
+void pedge_delete(pedge e);
+pedge pedge_realloc(pedge e, int np);
+pedge pedge_wgts_realloc(pedge e, int n);
+void pedge_export_mma(FILE *fp, int ne, pedge *edges);
+void pedge_export_gv(FILE *fp, int ne, pedge *edges);
+enum {METHOD_NONE = -1, METHOD_FD, METHOD_INK_AGGLOMERATE, METHOD_INK};
+enum {COMPATIBILITY_DIST = 0, COMPATIBILITY_FULL};
+pedge pedge_new(int np, int dim, real *x);
+pedge pedge_wgt_new(int np, int dim, real *x, real wgt);
+pedge pedge_double(pedge e);
+
+/* flip the polyline so that last point becomes the first, second last the second, etc*/
+pedge pedge_flip(pedge e);
+
+
+#endif /* EDGE_BUNDLING_H */
+
+
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include "types.h"
+#include "globals.h"
+#include "general.h"
+#include "SparseMatrix.h"
+#include "edge_bundling.h"
+#include "ink.h"
+
+double ink_count;
+
+point_t addPoint (point_t a, point_t b)
+{
+ a.x += b.x;
+ a.y += b.y;
+ return a;
+}
+
+point_t subPoint (point_t a, point_t b)
+{
+ a.x -= b.x;
+ a.y -= b.y;
+ return a;
+}
+
+point_t scalePoint (point_t a, double d)
+{
+ a.x *= d;
+ a.y *= d;
+ return a;
+}
+
+double dotPoint(point_t a, point_t b){
+ return a.x*b.x + a.y*b.y;
+}
+
+
+point_t Origin;
+
+/* sumLengths:
+ */
+double sumLengths_avoid_bad_angle(point_t* points, int npoints, point_t end, point_t meeting, real angle_param)
+{
+ /* avoid sharp turns, we want cos_theta to be as close to -1 as possible */
+ int i;
+ double len0, len, sum = 0;
+ double diff_x, diff_y, diff_x0, diff_y0;
+ double cos_theta, cos_max = -10;
+
+ diff_x0 = end.x-meeting.x;
+ diff_y0 = end.y-meeting.y;
+ len0 = sum = sqrt(diff_x0*diff_x0+diff_y0*diff_y0);
+
+ // distance form each of 'points' till 'meeting'
+ for (i=0; i<npoints; i++) {
+ diff_x = points[i].x-meeting.x;
+ diff_y = points[i].y-meeting.y;
+ len = sqrt(diff_x*diff_x+diff_y*diff_y);
+ sum += len;
+ cos_theta = (diff_x0*diff_x + diff_y0*diff_y)/MAX((len*len0),0.00001);
+ cos_max = MAX(cos_max, cos_theta);
+ }
+
+ // distance of single line from 'meeting' to 'end'
+ return sum*(cos_max + angle_param);/* straight line gives angle_param - 1, turning angle of 180 degree gives angle_param + 1 */
+}
+double sumLengths(point_t* points, int npoints, point_t end, point_t meeting)
+{
+ int i;
+ double sum = 0;
+ double diff_x, diff_y;
+
+ // distance form each of 'points' till 'meeting'
+ for (i=0; i<npoints; i++) {
+ diff_x = points[i].x-meeting.x;
+ diff_y = points[i].y-meeting.y;
+ sum += sqrt(diff_x*diff_x+diff_y*diff_y);
+ }
+ // distance of single line from 'meeting' to 'end'
+ diff_x = end.x-meeting.x;
+ diff_y = end.y-meeting.y;
+ sum += sqrt(diff_x*diff_x+diff_y*diff_y);
+ return sum;
+}
+
+/* bestInk:
+ */
+double bestInk(point_t* points, int npoints, point_t begin, point_t end, double prec, point_t *meet, real angle_param)
+{
+ point_t first, second, third, fourth, diff, meeting;
+ double value1, value2, value3, value4;
+
+ first = begin;
+ fourth = end;
+
+ do {
+ diff = subPoint(fourth,first);
+ second = addPoint(first,scalePoint(diff,1.0/3.0));
+ third = addPoint(first,scalePoint(diff,2.0/3.0));
+
+ if (angle_param < 1){
+ value1 = sumLengths(points, npoints, end, first);
+ value2 = sumLengths(points, npoints, end, second);
+ value3 = sumLengths(points, npoints, end, third);
+ value4 = sumLengths(points, npoints, end, fourth);
+ } else {
+ value1 = sumLengths_avoid_bad_angle(points, npoints, end, first, angle_param);
+ value2 = sumLengths_avoid_bad_angle(points, npoints, end, second, angle_param);
+ value3 = sumLengths_avoid_bad_angle(points, npoints, end, third, angle_param);
+ value4 = sumLengths_avoid_bad_angle(points, npoints, end, fourth, angle_param);
+ }
+
+ if (value1<value2) {
+ if (value1<value3) {
+ if (value1<value4) {
+ // first is smallest
+ fourth = second;
+ }
+ else {
+ // fourth is smallest
+ first = third;
+ }
+ }
+ else {
+ if (value3<value4) {
+ // third is smallest
+ first = second;
+ }
+ else {
+ // fourth is smallest
+ first = third;
+ }
+ }
+ }
+ else {
+ if (value2<value3) {
+ if (value2<value4) {
+ // second is smallest
+ fourth = third;
+ }
+ else {
+ // fourth is smallest
+ first = third;
+ }
+ }
+ else {
+ if (value3<value4) {
+ // third is smallest
+ first = second;
+ }
+ else {
+ // fourth is smallest
+ first = third;
+ }
+ }
+ }
+ } while (fabs(value1-value4)/(MIN(value1,value4)+1e-10)>prec && dotPoint(diff, diff) > 1.e-20);
+
+ meeting = scalePoint(addPoint(first,fourth),0.5);
+ *meet = meeting;
+
+ return sumLengths(points, npoints, end, meeting);
+
+}
+
+static double project_to_line(point_t pt, point_t left, point_t right, real angle){
+ /* pt
+ ^ ^
+ . \ \
+ . \ \
+ d . a\ \
+ . \ \
+ . \ \
+ . c \ alpha \ b
+ .<------left:0 ----------------------------> right:1. Find the projection of pt on the left--right line. If the turning angle is small,
+ | |
+ |<-------f---------
+ we should get a negative number. Let a := left->pt, b := left->right, then we are calculating:
+ c = |a| cos(a,b)/|b| b
+ d = a - c
+ f = -ctan(alpha)*|d|/|b| b
+ and the project of alpha degree on the left->right line is
+ c-f = |a| cos(a,b)/|b| b - -ctan(alpha)*|d|/|b| b
+ = (|a| a.b/(|a||b|) + ctan(alpha)|a-c|)/|b| b
+ = (a.b/|b| + ctan(alpha)|a-c|)/|b| b
+ the dimentionless projection is:
+ a.b/|b|^2 + ctan(alpha)|a-c|/|b|
+ = a.b/|b|^2 + ctan(alpha)|d|/|b|
+ */
+
+
+ point_t b, a;
+ real bnorm, dnorm;
+ real alpha, ccord;
+
+ if (angle <=0 || angle >= M_PI) return 2;/* return outside of the interval which should be handled as a sign of infeasible turning angle */
+ alpha = angle;
+
+ assert(alpha > 0 && alpha < M_PI);
+ b = subPoint(right, left);
+ a = subPoint(pt, left);
+ bnorm = MAX(1.e-10, dotPoint(b, b));
+ ccord = dotPoint(b, a)/bnorm;
+ dnorm = dotPoint(a,a)/bnorm - ccord*ccord;
+ if (alpha == M_PI/2){
+ return ccord;
+ } else {
+ // assert(dnorm >= MIN(-1.e-5, -1.e-5*bnorm));
+ return ccord + sqrt(MAX(0, dnorm))/tan(alpha);
+ }
+}
+
+
+
+
+
+
+
+
+
+/* ink:
+ * Compute minimal ink used the input edges are bundled.
+ * Assumes tails all occur on one side and heads on the other.
+ */
+double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle)
+{
+ int i;
+ point_t begin, end, mid, diff;
+ pedge e;
+ real *x;
+ point_t* sources = N_NEW(numEdges, point_t);
+ point_t* targets = N_NEW(numEdges, point_t);
+ double inkUsed;
+ //double eps = 1.0e-2;
+ double eps = 1.0e-2;
+ double cend = 0, cbegin = 0;
+ double wgt = 0;
+
+ // fprintf(stderr,"in ink code ========\n");
+ ink_count += numEdges;
+
+ *ink0 = 0;
+
+ /* canonicalize so that edges 1,2,3 and 3,2,1 gives the same optimal ink */
+ if (pick) vector_sort_int(numEdges, pick, TRUE);
+
+ begin = end = Origin;
+ for (i = 0; i < numEdges; i++) {
+ if (pick) {
+ e = edges[pick[i]];
+ } else {
+ e = edges[i];
+ }
+ x = e->x;
+ sources[i].x = x[0];
+ sources[i].y = x[1];
+ targets[i].x = x[e->dim*e->npoints - e->dim];
+ targets[i].y = x[e->dim*e->npoints - e->dim + 1];
+ (*ink0) += sqrt((sources[i].x - targets[i].x)*(sources[i].x - targets[i].x) + (sources[i].y - targets[i].y)*(sources[i].y - targets[i].y));
+ begin = addPoint (begin, scalePoint(sources[i], e->wgt));
+ end = addPoint (end, scalePoint(targets[i], e->wgt));
+ wgt += e->wgt;
+ //fprintf(stderr,"source={%f,%f}, target = {%f,%f}\n",sources[i].x, sources[i].y,
+ //targets[i].x, targets[i].y);
+ }
+
+ begin = scalePoint (begin, 1.0/wgt);
+ end = scalePoint (end, 1.0/wgt);
+
+
+ if (numEdges == 1){
+ *meet1 = begin;
+ *meet2 = end;
+ //fprintf(stderr,"ink used = %f\n",*ink0);
+ free (sources);
+ free (targets);
+ return *ink0;
+ }
+
+ /* shift the begin and end point to avoid sharp turns */
+ for (i = 0; i < numEdges; i++) {
+ if (pick) {
+ e = edges[pick[i]];
+ } else {
+ e = edges[i];
+ }
+ x = e->x;
+ sources[i].x = x[0];
+ sources[i].y = x[1];
+ targets[i].x = x[e->dim*e->npoints - e->dim];
+ targets[i].y = x[e->dim*e->npoints - e->dim + 1];
+ /* begin(1) ----------- mid(0) */
+ if (i == 0){
+ cbegin = project_to_line(sources[i], begin, end, angle);
+ cend = project_to_line(targets[i], end, begin, angle);
+ } else {
+ cbegin = MAX(cbegin, project_to_line(sources[i], begin, end, angle));
+ cend = MAX(cend, project_to_line(targets[i], end, begin, angle));
+ }
+ }
+
+ if (angle > 0 && angle < M_PI){
+ if (cbegin + cend > 1 || cbegin > 1 || cend > 1){
+ /* no point can be found that satisfies the angular constraints, so we give up and set ink to a large value */
+ if (Verbose && 0) fprintf(stderr,"no point satisfying any angle constraints can be found. cbeg=%f cend=%f\n",cbegin,cend);
+ inkUsed = 1000*(*ink0);
+ *meet1 = *meet2 = mid;
+ free (sources);
+ free (targets);
+ return inkUsed;
+ }
+ /* make sure the turning angle is no more than alpha degree */
+ cbegin = MAX(0, cbegin);/* make sure the new adjusted point is with in [begin,end] internal */
+ diff = subPoint(end, begin);
+ begin = addPoint(begin, scalePoint(diff, cbegin));
+
+ cend = MAX(0, cend);/* make sure the new adjusted point is with in [end,begin] internal */
+ end = subPoint(end, scalePoint(diff, cend));
+ }
+ mid = scalePoint (addPoint(begin,end),0.5);
+
+ inkUsed = (bestInk (sources, numEdges, begin, mid, eps, meet1, angle_param)
+ + bestInk (targets, numEdges, end, mid, eps, meet2, angle_param));
+ //fprintf(stderr,"beg={%f,%f}, meet1={%f,%f}, meet2={%f,%f}, mid={%f,%f}, end={%f,%f}\n",begin.x, begin.y, meet1->x, meet1->y, meet2->x, meet2->y,
+ //mid.x, mid.y, end.x, end.y);
+
+ //fprintf(stderr,"ink used = %f\n",inkUsed);
+ // fprintf(stderr,"{cb,ce}={%f, %f} end={%f,%f}, meet={%f,%f}, mid={%f,%f}\n",cbegin, cend, end.x, end.y, meet2->x, meet2->y, mid.x, mid.y);
+ free (sources);
+ free (targets);
+ return inkUsed;
+}
+
+
+double ink1(pedge e){
+
+
+ real *x, xx, yy;
+
+ real ink0 = 0;
+
+ x = e->x;
+ xx = x[0] - x[e->dim*e->npoints - e->dim];
+ yy = x[1] - x[e->dim*e->npoints - e->dim + 1];
+ ink0 += sqrt(xx*xx + yy*yy);
+ return ink0;
+}
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifndef INK_H
+#define INK_H
+
+#include <edge_bundling.h>
+
+typedef struct {
+ double x, y;
+} point_t;
+
+/* given a list of edges, find the best ink bundling by making them meet at 2 points
+ \ /
+ -meet1 ---------- meet2 -
+ / \
+ edges: list of edges
+ numEdges: number of edges
+ pick: if not NULL, consider edges pick[0], pick[1], ..., pick[numedges-1],
+ . othetwise consider edges 0, 1, ..., numEdge-1
+ ink0: ink needed if no bundling
+ meet1, meet2: meeting point
+ return: best ink needed if bundled.
+*/
+double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle);
+double ink1(pedge e);
+
+extern double ink_count;
+
+#endif /* INK_H */
--- /dev/null
+<?xml version="1.0" encoding="Windows-1252"?>\r
+<VisualStudioProject\r
+ ProjectType="Visual C++"\r
+ Version="9.00"\r
+ Name="minglelib"\r
+ ProjectGUID="{70575BD2-A532-41B8-9A17-B9876E992A84}"\r
+ RootNamespace="minglelib"\r
+ Keyword="Win32Proj"\r
+ TargetFrameworkVersion="196613"\r
+ >\r
+ <Platforms>\r
+ <Platform\r
+ Name="Win32"\r
+ />\r
+ </Platforms>\r
+ <ToolFiles>\r
+ </ToolFiles>\r
+ <Configurations>\r
+ <Configuration\r
+ Name="Debug|Win32"\r
+ OutputDirectory="Debug"\r
+ IntermediateDirectory="Debug"\r
+ ConfigurationType="4"\r
+ CharacterSet="1"\r
+ >\r
+ <Tool\r
+ Name="VCPreBuildEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCCustomBuildTool"\r
+ />\r
+ <Tool\r
+ Name="VCXMLDataGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCWebServiceProxyGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCMIDLTool"\r
+ />\r
+ <Tool\r
+ Name="VCCLCompilerTool"\r
+ Optimization="0"\r
+ AdditionalIncludeDirectories=""$(SolutionDir)/lib/sparse";"$(SolutionDir)/lib/sfpdpgen";"$(SolutionDir)/lib/common";"$(SolutionDir)/lib/gvc";"$(SolutionDir)/lib/cdt";"$(SolutionDir)/lib/pathplan";"$(SolutionDir)/lib/cgraph";"$(SolutionDir)";"C:\graphviz-ms\graphviz2\lib\mingle""\r
+ PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+ MinimalRebuild="true"\r
+ BasicRuntimeChecks="3"\r
+ RuntimeLibrary="3"\r
+ UsePrecompiledHeader="0"\r
+ WarningLevel="3"\r
+ DebugInformationFormat="4"\r
+ />\r
+ <Tool\r
+ Name="VCManagedResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCPreLinkEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCLibrarianTool"\r
+ />\r
+ <Tool\r
+ Name="VCALinkTool"\r
+ />\r
+ <Tool\r
+ Name="VCXDCMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCBscMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCFxCopTool"\r
+ />\r
+ <Tool\r
+ Name="VCPostBuildEventTool"\r
+ />\r
+ </Configuration>\r
+ <Configuration\r
+ Name="Release|Win32"\r
+ OutputDirectory="Release"\r
+ IntermediateDirectory="Release"\r
+ ConfigurationType="4"\r
+ CharacterSet="1"\r
+ WholeProgramOptimization="1"\r
+ >\r
+ <Tool\r
+ Name="VCPreBuildEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCCustomBuildTool"\r
+ />\r
+ <Tool\r
+ Name="VCXMLDataGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCWebServiceProxyGeneratorTool"\r
+ />\r
+ <Tool\r
+ Name="VCMIDLTool"\r
+ />\r
+ <Tool\r
+ Name="VCCLCompilerTool"\r
+ Optimization="2"\r
+ EnableIntrinsicFunctions="true"\r
+ AdditionalIncludeDirectories=""$(SolutionDir)/lib/sparse";"$(SolutionDir)/lib/sfpdpgen";"$(SolutionDir)/lib/common";"$(SolutionDir)/lib/gvc";"$(SolutionDir)/lib/cdt";"$(SolutionDir)/lib/pathplan";"$(SolutionDir)/lib/cgraph";"$(SolutionDir)";"C:\graphviz-ms\graphviz2\lib\mingle""\r
+ PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+ RuntimeLibrary="2"\r
+ EnableFunctionLevelLinking="true"\r
+ UsePrecompiledHeader="0"\r
+ WarningLevel="3"\r
+ DebugInformationFormat="3"\r
+ />\r
+ <Tool\r
+ Name="VCManagedResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCResourceCompilerTool"\r
+ />\r
+ <Tool\r
+ Name="VCPreLinkEventTool"\r
+ />\r
+ <Tool\r
+ Name="VCLibrarianTool"\r
+ />\r
+ <Tool\r
+ Name="VCALinkTool"\r
+ />\r
+ <Tool\r
+ Name="VCXDCMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCBscMakeTool"\r
+ />\r
+ <Tool\r
+ Name="VCFxCopTool"\r
+ />\r
+ <Tool\r
+ Name="VCPostBuildEventTool"\r
+ />\r
+ </Configuration>\r
+ </Configurations>\r
+ <References>\r
+ </References>\r
+ <Files>\r
+ <Filter\r
+ Name="Source Files"\r
+ Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"\r
+ UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"\r
+ >\r
+ <File\r
+ RelativePath=".\agglomerative_bundling.c"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\edge_bundling.c"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\ink.c"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\nearest_neighbor_graph.c"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\nearest_neighbor_graph_ann.cpp"\r
+ >\r
+ </File>\r
+ </Filter>\r
+ <Filter\r
+ Name="Header Files"\r
+ Filter="h;hpp;hxx;hm;inl;inc;xsd"\r
+ UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"\r
+ >\r
+ <File\r
+ RelativePath=".\agglomerative_bundling.h"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\edge_bundling.h"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\ink.h"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\nearest_neighbor_graph.h"\r
+ >\r
+ </File>\r
+ <File\r
+ RelativePath=".\nearest_neighbor_graph_ann.h"\r
+ >\r
+ </File>\r
+ </Filter>\r
+ <Filter\r
+ Name="Resource Files"\r
+ Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav"\r
+ UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"\r
+ >\r
+ </Filter>\r
+ <File\r
+ RelativePath=".\ReadMe.txt"\r
+ >\r
+ </File>\r
+ </Files>\r
+ <Globals>\r
+ </Globals>\r
+</VisualStudioProject>\r
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "general.h"
+#include "SparseMatrix.h"
+#include "nearest_neighbor_graph_ann.h"
+#include "nearest_neighbor_graph.h"
+
+SparseMatrix nearest_neighbor_graph(int nPts, int num_neigbors, int dim, double *x, double eps){
+ /* Gives a nearest neighbor graph of a list of dim-dimendional points. The result is a sparse matrix
+ of nPts x nPts, with num_neigbors entries per row.
+
+ nPts: number of points
+ num_neigbors: number of neighbors needed
+ dim: dimension
+ eps: error tolerance
+ x: nPts*dim vector. The i-th point is x[i*dim : i*dim + dim - 1]
+
+ */
+ int *irn = NULL, *jcn = NULL, nz;
+ real *val = NULL;
+ SparseMatrix A;
+ int k = num_neigbors;
+
+ /* need to *2 as we do two sweeps of neighbors, so could have repeats */
+ irn = MALLOC(sizeof(int)*nPts*k*2);
+ jcn = MALLOC(sizeof(int)*nPts*k*2);
+ val = MALLOC(sizeof(double)*nPts*k*2);
+
+#ifdef HAVE_ANN
+ nearest_neighbor_graph_ann(nPts, dim, num_neigbors, eps, x, &nz, &irn, &jcn, &val);
+
+ A = SparseMatrix_from_coordinate_arrays(nz, nPts, nPts, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
+#else
+ A = NULL;
+#endif
+
+ FREE(irn);
+ FREE(jcn);
+ FREE(val);
+
+ return A;
+
+
+}
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifndef NEAREST_NEIGHBOR_GRAPH_H
+#define NEAREST_NEIGHBOR_GRAPH_H
+
+SparseMatrix nearest_neighbor_graph(int nPts, int num_neigbors, int dim, double *x, double eps);
+
+#endif /* NEAREST_NEIGHBOR_GRAPH_H */
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_ANN
+//----------------------------------------------------------------------
+// File: get_nearest_neighb_graph.cpp
+//----------------------------------------------------------------------
+
+#include <cstdlib> // C standard library
+#include <cstdio> // C I/O (for sscanf)
+#include <cstring> // string manipulation
+#include <fstream> // file I/O
+#include <ANN/ANN.h> // ANN declarations
+
+using namespace std; // make std:: accessible
+int dim = 4; // dimension
+
+
+static void printPt(ostream &out, ANNpoint p) // print point
+{
+ out << "" << p[0];
+ for (int i = 1; i < dim; i++) {
+ out << "," << p[i];
+ }
+ out << "";
+}
+
+static void sortPtsX(int n, ANNpointArray pts){
+ /* sort so that edges always go from left to right in x-doordinate */
+ ANNpoint p;
+ ANNcoord x, y;
+ int i, j;
+ for (i = 0; i < n; i++){
+ for (j = 0; j < dim; j++){
+ p = pts[i];
+ if (p[0] < p[2] || (p[0] == p[2] && p[1] < p[3])) continue;
+ x = p[0]; y = p[1];
+ p[0] = p[2];
+ p[1] = p[3];
+ p[2] = x;
+ p[3] = y;
+ }
+ }
+}
+
+static void sortPtsY(int n, ANNpointArray pts){
+ /* sort so that edges always go from left to right in x-doordinate */
+ ANNpoint p;
+ ANNcoord x, y;
+ int i, j;
+ for (i = 0; i < n; i++){
+ for (j = 0; j < dim; j++){
+ p = pts[i];
+ if (p[1] < p[3] || (p[1] == p[3] && p[0] < p[2])) continue;
+ x = p[0]; y = p[1];
+ p[0] = p[2];
+ p[1] = p[3];
+ p[2] = x;
+ p[3] = y;
+ }
+ }
+}
+
+
+extern "C" void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0);
+
+void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0){
+
+ /* Gives a nearest neighbor graph is a list of dim-dimendional points. The connectivity is in irn/jcn, and the distance in val.
+
+ nPts: number of points
+ dim: dimension
+ k: number of neighbors needed
+ eps: error tolerance
+ x: nPts*dim vector. The i-th point is x[i*dim : i*dim + dim - 1]
+ nz: number of entries in the connectivity matrix irn/jcn/val
+ irn, jcn: the connectivity
+ val: the distance
+
+ note that there could be repeates
+ */
+
+ ANNpointArray dataPts; // data points
+ ANNidxArray nnIdx; // near neighbor indices
+ ANNdistArray dists; // near neighbor distances
+ ANNkd_tree* kdTree; // search structure
+
+ double *xx;
+ int *irn, *jcn;
+ double *val;
+ int nz;
+
+ irn = *irn0;
+ jcn = *jcn0;
+ val = *val0;
+
+
+ dataPts = annAllocPts(nPts, dim); // allocate data points
+ nnIdx = new ANNidx[k]; // allocate near neigh indices
+ dists = new ANNdist[k]; // allocate near neighbor dists
+
+ for (int i = 0; i < nPts; i++){
+ xx = dataPts[i];
+ for (int j = 0; j < dim; j++) xx[j] = x[i*dim + j];
+ }
+
+ //========= graph when sort based on x ========
+ nz = 0;
+ sortPtsX(nPts, dataPts);
+ kdTree = new ANNkd_tree( // build search structure
+ dataPts, // the data points
+ nPts, // number of points
+ dim); // dimension of space
+ for (int ip = 0; ip < nPts; ip++){
+ kdTree->annkSearch( // search
+ dataPts[ip], // query point
+ k, // number of near neighbors
+ nnIdx, // nearest neighbors (returned)
+ dists, // distance (returned)
+ eps); // error bound
+
+ for (int i = 0; i < k; i++) { // print summary
+ if (nnIdx[i] == ip) continue;
+ val[nz] = dists[i];
+ irn[nz] = ip;
+ jcn[nz++] = nnIdx[i];
+ //cout << ip << "--" << nnIdx[i] << " [len = " << dists[i]<< ", weight = \"" << 1./(dists[i]) << "\", wt = \"" << 1./(dists[i]) << "\"]\n";
+ //printPt(cout, dataPts[ip]);
+ //cout << "--";
+ //printPt(cout, dataPts[nnIdx[i]]);
+ //cout << "\n";
+ }
+ }
+
+
+ //========= graph when sort based on y ========
+ sortPtsY(nPts, dataPts);
+ kdTree = new ANNkd_tree( // build search structure
+ dataPts, // the data points
+ nPts, // number of points
+ dim); // dimension of space
+ for (int ip = 0; ip < nPts; ip++){
+ kdTree->annkSearch( // search
+ dataPts[ip], // query point
+ k, // number of near neighbors
+ nnIdx, // nearest neighbors (returned)
+ dists, // distance (returned)
+ eps); // error bound
+
+ for (int i = 0; i < k; i++) { // print summary
+ if (nnIdx[i] == ip) continue;
+ val[nz] = dists[i];
+ irn[nz] = ip;
+ jcn[nz++] = nnIdx[i];
+ // cout << ip << "--" << nnIdx[i] << " [len = " << dists[i]<< ", weight = \"" << 1./(dists[i]) << "\", wt = \"" << 1./(dists[i]) << "\"]\n";
+ }
+ }
+
+ delete [] nnIdx; // clean things up
+ delete [] dists;
+ delete kdTree;
+
+ *nz0 = nz;
+
+ annClose(); // done with ANN
+
+
+
+}
+
+#endif /* HAVE_ANN */
--- /dev/null
+/*************************************************************************
+ * Copyright (c) 2011 AT&T Intellectual Property
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * which accompanies this distribution, and is available at
+ * http://www.eclipse.org/legal/epl-v10.html
+ *
+ * Contributors: See CVS logs. Details at http://www.graphviz.org/
+ *************************************************************************/
+
+#ifndef NEAREST_NEIGHBOR_GRAPH_ANN_H
+#define NEAREST_NEIGHBOR_GRAPH_ANN_H
+
+void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0);
+
+#endif /* NEAREST_NEIGHBOR_GRAPH_ANN_H */