extern void printvis(vconfig_t * cp);
extern int in_poly(Ppoly_t argpoly, Ppoint_t q);
+#ifdef WITH_CGRAPH
+static boolean spline_merge(node_t * n)
+{
+ return FALSE;
+}
+
+static boolean swap_ends_p(edge_t * e)
+{
+ return FALSE;
+}
+
+static splineInfo sinfo = { swap_ends_p, spline_merge };
+
+static void
+make_barriers(Ppoly_t ** poly, int npoly, int pp, int qp,
+ Pedge_t ** barriers, int *n_barriers)
+{
+ int i, j, k, n, b;
+ Pedge_t *bar;
+
+ n = 0;
+ for (i = 0; i < npoly; i++) {
+ if (i == pp)
+ continue;
+ if (i == qp)
+ continue;
+ n = n + poly[i]->pn;
+ }
+ bar = N_GNEW(n, Pedge_t);
+ b = 0;
+ for (i = 0; i < npoly; i++) {
+ if (i == pp)
+ continue;
+ if (i == qp)
+ continue;
+ for (j = 0; j < poly[i]->pn; j++) {
+ k = j + 1;
+ if (k >= poly[i]->pn)
+ k = 0;
+ bar[b].a = poly[i]->ps[j];
+ bar[b].b = poly[i]->ps[k];
+ b++;
+ }
+ }
+ assert(b == n);
+ *barriers = bar;
+ *n_barriers = n;
+}
+
+/* genPt:
+ */
+static Ppoint_t genPt(double x, double y, point c)
+{
+ Ppoint_t p;
+
+ p.x = x + c.x;
+ p.y = y + c.y;
+ return p;
+}
+
+
+/* recPt:
+ */
+static Ppoint_t recPt(double x, double y, point c, expand_t* m)
+{
+ Ppoint_t p;
+
+ p.x = (x * m->x) + c.x;
+ p.y = (y * m->y) + c.y;
+ return p;
+}
+static void makePortLabels(edge_t * e)
+{
+ if (ED_head_label(e) && !ED_head_label(e)->set) {
+ place_portlabel(e, TRUE);
+ updateBB(agraphof(agtail(e)), ED_head_label(e));
+ }
+ if (ED_tail_label(e) && !ED_tail_label(e)->set) {
+ place_portlabel(e, FALSE);
+ updateBB(agraphof(agtail(e)), ED_tail_label(e));
+ }
+}
+
+/* endPoints:
+ * Extract the actual end points of the spline, where
+ * they touch the node.
+ */
+static void endPoints(splines * spl, point * p, point * q)
+{
+ bezier bz;
+
+ bz = spl->list[0];
+ if (bz.sflag)
+ *p = bz.sp;
+ else
+ *p = bz.list[0];
+ bz = spl->list[spl->size - 1];
+ if (bz.eflag)
+ *q = bz.ep;
+ else
+ *q = bz.list[bz.size - 1];
+}
+
+/* polylineMidpoint;
+ * Find midpoint of polyline.
+ * pp and pq are set to the endpoints of the line segment containing it.
+ */
+static point
+polylineMidpoint (splines* spl, point* pp, point* pq)
+{
+ bezier bz;
+ int i, j, k;
+ double d, dist = 0;
+ point m;
+ pointf pf, qf, mf;
+
+ for (i = 0; i < spl->size; i++) {
+ bz = spl->list[i];
+ for (j = 0, k=3; k < bz.size; j+=3,k+=3) {
+ P2PF (bz.list[j],pf);
+ P2PF (bz.list[k],qf);
+ dist += DIST(pf,qf);
+ }
+ }
+ dist /= 2;
+ for (i = 0; i < spl->size; i++) {
+ bz = spl->list[i];
+ for (j = 0, k=3; k < bz.size; j+=3,k+=3) {
+ P2PF (bz.list[j],pf);
+ P2PF (bz.list[k],qf);
+ d = DIST(pf,qf);
+ if (d >= dist) {
+ *pp = bz.list[j];
+ *pq = bz.list[k];
+ mf.x = ((qf.x*dist) + (pf.x*(d-dist)))/d;
+ mf.y = ((qf.y*dist) + (pf.y*(d-dist)))/d;
+ PF2P(mf, m);
+ return m;
+ }
+ else
+ dist -= d;
+ }
+ }
+ assert (FALSE); /* should never get here */
+ return m;
+}
+
+#define LEFTOF(a,b,c) (((a.y - b.y)*(c.x - b.x) - (c.y - b.y)*(a.x - b.x)) > 0)
+#define MAXLABELWD (POINTS_PER_INCH/2.0)
+
+/* addEdgeLabels:
+ * rp and rq are the port points of the tail and head node.
+ * Adds label, headlabel and taillabel.
+ * The use of 2 and 4 in computing ld.x and ld.y are fudge factors, to
+ * introduce a bit of spacing.
+ * Updates bounding box.
+ * We try to use the actual endpoints of the spline, as they may differ
+ * significantly from rp and rq, but if the spline is degenerate (e.g.,
+ * the nodes overlap), we use rp and rq.
+ */
+void addEdgeLabels(edge_t * e, point rp, point rq)
+{
+ int et = EDGE_TYPE (agraphof(aghead(e))->root);
+ point p, q;
+ point d; /* midpoint of segment p-q */
+ point ld;
+ point sp;
+ point del;
+ pointf spf;
+ double f, ht, wd, dist2;
+ int leftOf;
+
+ if (ED_label(e) && !ED_label(e)->set) {
+ endPoints(ED_spl(e), &p, &q);
+ if ((p.x == q.x) && (p.y == q.y)) { /* degenerate spline */
+ p = rp;
+ q = rq;
+ sp = p;
+ }
+ else if (et == ET_SPLINE) {
+ d.x = (q.x + p.x) / 2;
+ d.y = (p.y + q.y) / 2;
+ sp = dotneato_closest(ED_spl(e), d);
+ }
+ else { /* ET_PLINE or ET_LINE */
+ sp = polylineMidpoint (ED_spl(e), &p, &q);
+ }
+ del.x = q.x - p.x;
+ del.y = q.y - p.y;
+ dist2 = del.x*del.x + del.y*del.y;
+ ht = (ED_label(e)->dimen.y + 2)/2.0;
+ spf.x = sp.x;
+ spf.y = sp.y;
+ if (dist2) {
+ wd = (MIN(ED_label(e)->dimen.x + 2, MAXLABELWD))/2.0;
+ leftOf = LEFTOF(p, q, sp);
+ if ((leftOf && (del.y >= 0)) || (!leftOf && (del.y < 0))) {
+ if (del.x*del.y >= 0)
+ ht *= -1;
+ }
+ else {
+ wd *= -1;
+ if (del.x*del.y < 0)
+ ht *= -1;
+ }
+ f = (del.y*wd - del.x*ht)/dist2;
+ ld.x = -f*del.y;
+ ld.y = f*del.x;
+ }
+ else { /* end points the same */
+ ld.x = 0;
+ ld.y = -ht;
+ }
+
+ ED_label(e)->p.x = spf.x + ld.x;
+ ED_label(e)->p.y = spf.y + ld.y;
+ ED_label(e)->set = TRUE;
+ updateBB(agraphof(agtail(e)), ED_label(e));
+ }
+ makePortLabels(e);
+}
+
+typedef struct {
+ node_t *n1;
+ point p1;
+ node_t *n2;
+ point p2;
+} edgeinfo;
+typedef struct {
+ Dtlink_t link;
+ edgeinfo id;
+ edge_t *e;
+} edgeitem;
+
+static void *newitem(Dt_t * d, edgeitem * obj, Dtdisc_t * disc)
+{
+ edgeitem *newp;
+
+ NOTUSED(disc);
+ newp = NEW(edgeitem);
+ newp->id = obj->id;
+ newp->e = obj->e;
+ ED_count(newp->e) = 1;
+
+ return newp;
+}
+
+static void freeitem(Dt_t * d, edgeitem * obj, Dtdisc_t * disc)
+{
+ free(obj);
+}
+
+static int
+cmpitems(Dt_t * d, edgeinfo * key1, edgeinfo * key2, Dtdisc_t * disc)
+{
+ int x;
+
+ if (key1->n1 > key2->n1)
+ return 1;
+ if (key1->n1 < key2->n1)
+ return -1;
+ if (key1->n2 > key2->n2)
+ return 1;
+ if (key1->n2 < key2->n2)
+ return -1;
+
+ if ((x = key1->p1.x - key2->p1.x))
+ return x;
+ if ((x = key1->p1.y - key2->p1.y))
+ return x;
+ if ((x = key1->p2.x - key2->p2.x))
+ return x;
+ return (key1->p2.y - key2->p2.y);
+}
+
+Dtdisc_t edgeItemDisc = {
+ offsetof(edgeitem, id),
+ sizeof(edgeinfo),
+ offsetof(edgeitem, link),
+ (Dtmake_f) newitem,
+ (Dtfree_f) freeitem,
+ (Dtcompar_f) cmpitems,
+ 0,
+ 0,
+ 0
+};
+
+/* equivEdge:
+ * See if we have already encountered an edge between the same
+ * node:port pairs. If so, return the earlier edge. If not,
+ * this edge is added to map and returned.
+ * We first have to canonicalize the key fields using a lexicographic
+ * ordering.
+ */
+static edge_t *equivEdge(Dt_t * map, edge_t * e)
+{
+ edgeinfo test;
+ edgeitem dummy;
+ edgeitem *ip;
+
+ if (agtail(e) < aghead(e)) {
+ test.n1 = agtail(e);
+ test.p1 = ED_tail_port(e).p;
+ test.n2 = aghead(e);
+ test.p2 = ED_head_port(e).p;
+ } else if (agtail(e) > aghead(e)) {
+ test.n2 = agtail(e);
+ test.p2 = ED_tail_port(e).p;
+ test.n1 = aghead(e);
+ test.p1 = ED_head_port(e).p;
+ } else {
+ point hp = ED_head_port(e).p;
+ point tp = ED_tail_port(e).p;
+ if (tp.x < hp.x) {
+ test.p1 = tp;
+ test.p2 = hp;
+ } else if (tp.x > hp.x) {
+ test.p1 = hp;
+ test.p2 = tp;
+ } else if (tp.y < hp.y) {
+ test.p1 = tp;
+ test.p2 = hp;
+ } else if (tp.y > hp.y) {
+ test.p1 = hp;
+ test.p2 = tp;
+ } else {
+ test.p1 = test.p2 = tp;
+ }
+ test.n2 = test.n1 = agtail(e);
+ }
+ dummy.id = test;
+ dummy.e = e;
+ ip = dtinsert(map, &dummy);
+ return ip->e;
+}
+
+
+/* makeSelfArcs:
+ * Generate loops. We use the library routine makeSelfEdge
+ * which also places the labels.
+ * We have to handle port labels here.
+ * as well as update the bbox from edge labels.
+ */
+void makeSelfArcs(path * P, edge_t * e, int stepx)
+{
+ int cnt = ED_count(e);
+
+ if (cnt == 1) {
+ edge_t *edges1[1];
+ edges1[0] = e;
+ makeSelfEdge(P, edges1, 0, 1, stepx, stepx, &sinfo);
+ if (ED_label(e))
+ updateBB(agraphof(agtail(e)), ED_label(e));
+ makePortLabels(e);
+ } else {
+ int i;
+ edge_t **edges = N_GNEW(cnt, edge_t *);
+ for (i = 0; i < cnt; i++) {
+ edges[i] = e;
+ e = ED_to_virt(e);
+ }
+ makeSelfEdge(P, edges, 0, cnt, stepx, stepx, &sinfo);
+ for (i = 0; i < cnt; i++) {
+ e = edges[i];
+ if (ED_label(e))
+ updateBB(agraphof(agtail(e)), ED_label(e));
+ makePortLabels(e);
+ }
+ free(edges);
+ }
+}
+
+/* makeStraightEdge:
+ *
+ * FIX: handle ports on boundary?
+ */
+void
+makeStraightEdge(graph_t * g, edge_t * e, int doPolyline)
+{
+ point dumb[4];
+ node_t *n = agtail(e);
+ node_t *head = aghead(e);
+ int e_cnt = ED_count(e);
+ pointf perp;
+ point del;
+ edge_t *e0;
+ int i, j, xstep, dx;
+ double l_perp;
+ point dumber[4];
+ point p, q;
+
+ p = dumb[1] = dumb[0] = add_points(ND_coord_i(n), ED_tail_port(e).p);
+ q = dumb[2] = dumb[3] =
+ add_points(ND_coord_i(head), ED_head_port(e).p);
+
+ if (e_cnt == 1) {
+ clip_and_install(e, aghead(e), dumb, 4, &sinfo);
+ addEdgeLabels(e, p, q);
+ return;
+ }
+
+ e0 = e;
+ perp.x = dumb[0].y - dumb[3].y;
+ perp.y = dumb[3].x - dumb[0].x;
+ if ((perp.x == 0) && (perp.y == 0)) {
+ /* degenerate case */
+ dumb[1] = dumb[0];
+ dumb[2] = dumb[3];
+ del.x = 0;
+ del.y = 0;
+ }
+ else {
+ l_perp = sqrt(perp.x * perp.x + perp.y * perp.y);
+ xstep = GD_nodesep(g);
+ dx = xstep * (e_cnt - 1) / 2;
+ dumb[1].x = dumb[0].x + (dx * perp.x) / l_perp;
+ dumb[1].y = dumb[0].y + (dx * perp.y) / l_perp;
+ dumb[2].x = dumb[3].x + (dx * perp.x) / l_perp;
+ dumb[2].y = dumb[3].y + (dx * perp.y) / l_perp;
+ del.x = -xstep * perp.x / l_perp;
+ del.y = -xstep * perp.y / l_perp;
+ }
+
+ for (i = 0; i < e_cnt; i++) {
+ if (aghead(e0) == head) {
+ p = dumb[0];
+ q = dumb[3];
+ for (j = 0; j < 4; j++) {
+ dumber[j] = dumb[j];
+ }
+ } else {
+ p = dumb[3];
+ q = dumb[0];
+ for (j = 0; j < 4; j++) {
+ dumber[3 - j] = dumb[j];
+ }
+ }
+ if (doPolyline) {
+ Ppoint_t pts[4];
+ Ppolyline_t spl, line;
+ point* ispline;
+
+ line.pn = 4;
+ line.ps = pts;
+ for (i=0; i < 4; i++)
+ P2PF (dumber[i], pts[i]);
+ make_polyline (line, &spl);
+ ispline = N_GNEW(spl.pn, point);
+ for (i=0; i < spl.pn; i++)
+ PF2P (spl.ps[i], ispline[i]);
+ clip_and_install(e0, aghead(e0), ispline, spl.pn, &sinfo);
+ free(ispline);
+ }
+ else
+ clip_and_install(e0, aghead(e0), dumber, 4, &sinfo);
+
+ addEdgeLabels(e0, p, q);
+ e0 = ED_to_virt(e0);
+ dumb[1].x += del.x;
+ dumb[1].y += del.y;
+ dumb[2].x += del.x;
+ dumb[2].y += del.y;
+ }
+}
+
+#define LEN(x,y) sqrt((x)*(x)+(y)*(y))
+
+/* makeObstacle:
+ * Given a node, return an obstacle reflecting the
+ * node's geometry. pmargin specifies how much space to allow
+ * around the polygon.
+ * Returns the constructed polygon on success, NULL on failure.
+ * Failure means the node shape is not supported.
+ *
+ * The polygon has its vertices in CW order.
+ *
+ */
+Ppoly_t *makeObstacle(node_t * n, expand_t* pmargin)
+{
+ Ppoly_t *obs;
+ polygon_t *poly;
+ double adj = 0.0;
+ int j, sides;
+ pointf polyp;
+ box b;
+ point pt;
+ field_t *fld;
+ epsf_t *desc;
+
+ switch (shapeOf(n)) {
+ case SH_POLY:
+ case SH_POINT:
+ obs = NEW(Ppoly_t);
+ poly = (polygon_t *) ND_shape_info(n);
+ if (poly->sides >= 3) {
+ sides = poly->sides;
+ } else { /* ellipse */
+ sides = 8;
+ adj = drand48() * .01;
+ }
+ obs->pn = sides;
+ obs->ps = N_NEW(sides, Ppoint_t);
+ /* assuming polys are in CCW order, and pathplan needs CW */
+ for (j = 0; j < sides; j++) {
+ double xmargin = 0.0, ymargin = 0.0;
+ if (poly->sides >= 3) {
+ if (pmargin->doAdd) {
+ if (poly->sides == 4) {
+ switch (j) {
+ case 0 :
+ xmargin = pmargin->x;
+ ymargin = pmargin->y;
+ break;
+ case 1 :
+ xmargin = -pmargin->x;
+ ymargin = pmargin->y;
+ break;
+ case 2 :
+ xmargin = -pmargin->x;
+ ymargin = -pmargin->y;
+ break;
+ case 3 :
+ xmargin = pmargin->x;
+ ymargin = -pmargin->y;
+ break;
+ }
+ polyp.x = poly->vertices[j].x + xmargin;
+ polyp.y = poly->vertices[j].y + ymargin;
+ }
+ else {
+ double h = LEN(poly->vertices[j].x,poly->vertices[j].y);
+ polyp.x = poly->vertices[j].x * (1.0 + pmargin->x/h);
+ polyp.y = poly->vertices[j].y * (1.0 + pmargin->y/h);
+ }
+ }
+ else {
+ polyp.x = poly->vertices[j].x * pmargin->x;
+ polyp.y = poly->vertices[j].y * pmargin->y;
+ }
+ } else {
+ double c, s;
+ c = cos(2.0 * M_PI * j / sides + adj);
+ s = sin(2.0 * M_PI * j / sides + adj);
+ if (pmargin->doAdd) {
+ polyp.x = c*(ND_lw_i(n)+ND_rw_i(n)+pmargin->x) / 2.0;
+ polyp.y = s*(ND_ht_i(n)+pmargin->y) / 2.0;
+ }
+ else {
+ polyp.x = pmargin->x * c * (ND_lw_i(n) + ND_rw_i(n)) / 2.0;
+ polyp.y = pmargin->y * s * ND_ht_i(n) / 2.0;
+ }
+ }
+ obs->ps[sides - j - 1].x = polyp.x + ND_coord_i(n).x;
+ obs->ps[sides - j - 1].y = polyp.y + ND_coord_i(n).y;
+ }
+ break;
+ case SH_RECORD:
+ fld = (field_t *) ND_shape_info(n);
+ b = fld->b;
+ obs = NEW(Ppoly_t);
+ obs->pn = 4;
+ obs->ps = N_NEW(4, Ppoint_t);
+ /* CW order */
+ pt = ND_coord_i(n);
+ if (pmargin->doAdd) {
+ obs->ps[0] = genPt(b.LL.x-pmargin->x, b.LL.y-pmargin->y, pt);
+ obs->ps[1] = genPt(b.LL.x-pmargin->x, b.UR.y+pmargin->y, pt);
+ obs->ps[2] = genPt(b.UR.x+pmargin->x, b.UR.y+pmargin->y, pt);
+ obs->ps[3] = genPt(b.UR.x+pmargin->x, b.LL.y-pmargin->y, pt);
+ }
+ else {
+ obs->ps[0] = recPt(b.LL.x, b.LL.y, pt, pmargin);
+ obs->ps[1] = recPt(b.LL.x, b.UR.y, pt, pmargin);
+ obs->ps[2] = recPt(b.UR.x, b.UR.y, pt, pmargin);
+ obs->ps[3] = recPt(b.UR.x, b.LL.y, pt, pmargin);
+ }
+ break;
+ case SH_EPSF:
+ desc = (epsf_t *) (ND_shape_info(n));
+ obs = NEW(Ppoly_t);
+ obs->pn = 4;
+ obs->ps = N_NEW(4, Ppoint_t);
+ /* CW order */
+ pt = ND_coord_i(n);
+ if (pmargin->doAdd) {
+ obs->ps[0] = genPt(-ND_lw_i(n)-pmargin->x, -ND_ht_i(n)-pmargin->y,pt);
+ obs->ps[1] = genPt(-ND_lw_i(n)-pmargin->x, ND_ht_i(n)+pmargin->y,pt);
+ obs->ps[2] = genPt(ND_rw_i(n)+pmargin->x, ND_ht_i(n)+pmargin->y,pt);
+ obs->ps[3] = genPt(ND_rw_i(n)+pmargin->x, -ND_ht_i(n)-pmargin->y,pt);
+ }
+ else {
+ obs->ps[0] = recPt(-ND_lw_i(n), -ND_ht_i(n), pt, pmargin);
+ obs->ps[1] = recPt(-ND_lw_i(n), ND_ht_i(n), pt, pmargin);
+ obs->ps[2] = recPt(ND_rw_i(n), ND_ht_i(n), pt, pmargin);
+ obs->ps[3] = recPt(ND_rw_i(n), -ND_ht_i(n), pt, pmargin);
+ }
+ break;
+ default:
+ obs = NULL;
+ break;
+ }
+ return obs;
+}
+
+/* getPath
+ * Construct the shortest path from one endpoint of e to the other.
+ * The obstacles and their number are given by obs and npoly.
+ * vconfig is a precomputed data structure to help in the computation.
+ * If chkPts is true, the function finds the polygons, if any, containing
+ * the endpoints and tells the shortest path computation to ignore them.
+ * Assumes this info is set in ND_lim, usually from _spline_edges.
+ * Returns the shortest path.
+ */
+Ppolyline_t
+getPath(edge_t * e, vconfig_t * vconfig, int chkPts, Ppoly_t ** obs,
+ int npoly)
+{
+ Ppolyline_t line;
+ int pp, qp;
+ Ppoint_t p, q;
+ point p1, q1;
+
+ p1 = add_points(ND_coord_i(agtail(e)), ED_tail_port(e).p);
+ q1 = add_points(ND_coord_i(aghead(e)), ED_head_port(e).p);
+ P2PF(p1, p);
+ P2PF(q1, q);
+
+ /* determine the polygons (if any) that contain the endpoints */
+ pp = qp = POLYID_NONE;
+ if (chkPts) {
+ pp = ND_lim(agtail(e));
+ qp = ND_lim(aghead(e));
+/*
+ for (i = 0; i < npoly; i++) {
+ if ((pp == POLYID_NONE) && in_poly(*obs[i], p))
+ pp = i;
+ if ((qp == POLYID_NONE) && in_poly(*obs[i], q))
+ qp = i;
+ }
+*/
+ }
+ Pobspath(vconfig, p, pp, q, qp, &line);
+ return line;
+}
+
+/* makePolyline:
+ */
+static void
+makePolyline(edge_t * e)
+{
+ int i;
+ Ppolyline_t spl, line = ED_path(e);
+ point* ispline;
+ point p1, q1;
+ Ppoint_t p0, q0;
+
+ p0 = line.ps[0];
+ q0 = line.ps[line.pn - 1];
+ make_polyline (line, &spl);
+ ispline = N_GNEW(spl.pn, point);
+ for (i=0; i < spl.pn; i++) {
+ PF2P (spl.ps[i], ispline[i]);
+ }
+
+ if (Verbose > 1)
+ fprintf(stderr, "polyline %s %s\n", agnameof(agtail(e)), agnameof(aghead(e)));
+ clip_and_install(e, aghead(e), ispline, spl.pn, &sinfo);
+ free(ispline);
+ PF2P(p0, p1);
+ PF2P(q0, q1);
+ addEdgeLabels(e, p1, q1);
+}
+
+/* makeSpline:
+ * Construct a spline connecting the endpoints of e, avoiding the npoly
+ * obstacles obs.
+ * The resultant spline is attached to the edge, the positions of any
+ * edge labels are computed, and the graph's bounding box is recomputed.
+ *
+ * If chkPts is true, the function checks if one or both of the endpoints
+ * is on or inside one of the obstacles and, if so, tells the shortest path
+ * computation to ignore them.
+ */
+void makeSpline(edge_t * e, Ppoly_t ** obs, int npoly, boolean chkPts)
+{
+ Ppolyline_t line, spline;
+ Pvector_t slopes[2];
+ int i, n_barriers;
+ int pp, qp;
+ Ppoint_t p, q;
+ point *ispline;
+ Pedge_t *barriers;
+ point p1, q1;
+
+ line = ED_path(e);
+ p = line.ps[0];
+ q = line.ps[line.pn - 1];
+ /* determine the polygons (if any) that contain the endpoints */
+ pp = qp = POLYID_NONE;
+ if (chkPts)
+ for (i = 0; i < npoly; i++) {
+ if ((pp == POLYID_NONE) && in_poly(*obs[i], p))
+ pp = i;
+ if ((qp == POLYID_NONE) && in_poly(*obs[i], q))
+ qp = i;
+ }
+
+ make_barriers(obs, npoly, pp, qp, &barriers, &n_barriers);
+ slopes[0].x = slopes[0].y = 0.0;
+ slopes[1].x = slopes[1].y = 0.0;
+ Proutespline(barriers, n_barriers, line, slopes, &spline);
+
+ /* north why did you ever use int coords */
+ ispline = N_GNEW(spline.pn, point);
+ for (i = 0; i < spline.pn; i++) {
+ ispline[i].x = ROUND(spline.ps[i].x);
+ ispline[i].y = ROUND(spline.ps[i].y);
+ }
+ if (Verbose > 1)
+ fprintf(stderr, "spline %s %s\n", agnameof(agtail(e)),agnameof(aghead(e)));
+ clip_and_install(e, aghead(e), ispline, spline.pn, &sinfo);
+ free(ispline);
+ free(barriers);
+ PF2P(p, p1);
+ PF2P(q, q1);
+ addEdgeLabels(e, p1, q1);
+}
+
+ /* True if either head or tail has a port on its boundary */
+#define BOUNDARY_PORT(e) ((ED_tail_port(e).side)||(ED_head_port(e).side))
+
+/* _spline_edges:
+ * Basic default routine for creating edges.
+ * If splines are requested, we construct the obstacles.
+ * If not, or nodes overlap, the function reverts to line segments.
+ * NOTE: intra-cluster edges are not constrained to
+ * remain in the cluster's bounding box and, conversely, a cluster's box
+ * is not altered to reflect intra-cluster edges.
+ * If Nop > 1 and the spline exists, it is just copied.
+ */
+static int _spline_edges(graph_t * g, expand_t* pmargin, int edgetype)
+{
+ node_t *n;
+ edge_t *e;
+ edge_t *e0;
+ Ppoly_t **obs = 0;
+ Ppoly_t *obp;
+ int cnt, i = 0, npoly;
+ vconfig_t *vconfig = 0;
+ path *P = NULL;
+ int useEdges = (Nop > 1);
+#ifdef ORTHO
+ extern void orthoEdges (Agraph_t* g, int useLbls, splineInfo* sinfo);
+#endif
+ router_t* rtr = 0;
+
+ /* build configuration */
+ if (edgetype != ET_LINE) {
+ obs = N_NEW(agnnodes(g), Ppoly_t *);
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ obp = makeObstacle(n, pmargin);
+ if (obp) {
+ ND_lim(n) = i;
+ obs[i++] = obp;
+ }
+ else
+ ND_lim(n) = POLYID_NONE;
+ }
+ } else {
+ obs = 0;
+ }
+ npoly = i;
+ if (obs) {
+ if (Plegal_arrangement(obs, npoly)) {
+ if (edgetype != ET_ORTHO) vconfig = Pobsopen(obs, npoly);
+ }
+ else if (Verbose)
+ fprintf(stderr,
+ "nodes touch - falling back to straight line edges\n");
+ }
+
+ /* route edges */
+ if (Verbose)
+ fprintf(stderr, "Creating edges using %s\n",
+ (edgetype == ET_ORTHO) ? "orthogonal lines" :
+ (vconfig ? (edgetype == ET_SPLINE ? "splines" : "polylines") :
+ "line segments"));
+ if (vconfig) {
+ /* path-finding pass */
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ for (e = agfstout(g, n); e; e = agnxtout(g, e)) {
+ ED_path(e) = getPath(e, vconfig, TRUE, obs, npoly);
+ }
+ }
+ }
+#ifdef ORTHO
+ else if (edgetype == ET_ORTHO) {
+ orthoEdges (g, 0, &sinfo);
+ useEdges = 1;
+ }
+#endif
+
+ /* spline-drawing pass */
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ for (e = agfstout(g, n); e; e = agnxtout(g, e)) {
+/* fprintf (stderr, "%s -- %s %d\n", e->tail->name, e->head->name, ED_count(e)); */
+ node_t *head = aghead(e);
+ if (useEdges && ED_spl(e)) {
+ addEdgeLabels(e,
+ add_points(ND_coord_i(n), ED_tail_port(e).p),
+ add_points(ND_coord_i(head),
+ ED_head_port(e).p));
+ }
+ else if (ED_count(e) == 0) continue; /* only do representative */
+ else if (n == head) { /* self arc */
+ if (!P) {
+ P = NEW(path);
+ P->boxes = N_NEW(agnnodes(g) + 20 * 2 * 9, box);
+ }
+ makeSelfArcs(P, e, GD_nodesep(g));
+ } else if (vconfig) { /* ET_SPLINE or ET_POLYLINE */
+#ifdef HAVE_GTS
+ if ((ED_count(e) > 1) || BOUNDARY_PORT(e)) {
+ int fail = 0;
+ if ((ED_path(e).pn == 2) && !BOUNDARY_PORT(e))
+ /* if a straight line can connect the ends */
+ makeStraightEdge(g, e, edgetype == ET_PLINE);
+ else {
+ if (!rtr) rtr = mkRouter (obs, npoly);
+ fail = makeMultiSpline(e, rtr, edgetype == ET_PLINE);
+ }
+ if (!fail) continue;
+ }
+ /* We can probably remove this branch and just use
+ * makeMultiSpline. It can also catch the makeStraightEdge
+ * case. We could then eliminate all of the vconfig stuff.
+ */
+#endif
+ cnt = ED_count(e);
+ e0 = e;
+ for (i = 0; i < cnt; i++) {
+ if (edgetype == ET_SPLINE)
+ makeSpline(e0, obs, npoly, TRUE);
+ else
+ makePolyline(e0);
+ e0 = ED_to_virt(e0);
+ }
+ } else {
+ makeStraightEdge(g, e, 0);
+ }
+ }
+ }
+
+#ifdef HAVE_GTS
+ if (rtr)
+ freeRouter (rtr);
+#endif
+
+ if (vconfig)
+ Pobsclose (vconfig);
+ if (P) {
+ free(P->boxes);
+ free(P);
+ }
+ if (obs) {
+ for (i=0; i < npoly; i++)
+ free (obs[i]);
+ free (obs);
+ }
+ return 0;
+}
+
+/* splineEdges:
+ * Main wrapper code for generating edges.
+ * Sets desired separation.
+ * Coalesces equivalent edges (edges * with the same endpoints).
+ * It then calls the edge generating function, and marks the
+ * spline phase complete.
+ * Returns 0 on success.
+ *
+ * The edge function is given the graph, the separation to be added
+ * around obstacles, and the type of edge. It must guarantee
+ * that all bounding boxes are current; in particular, the bounding box of
+ * g must reflect the addition of the edges.
+ */
+int
+splineEdges(graph_t * g, int (*edgefn) (graph_t *, expand_t*, int),
+ int edgetype)
+{
+ node_t *n;
+ edge_t *e;
+ expand_t margin;
+ Dt_t *map;
+
+ margin = esepFactor (g);
+
+ /* find equivalent edges */
+ map = dtopen(&edgeItemDisc, Dtoset);
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ for (e = agfstout(g, n); e; e = agnxtout(g, e)) {
+ edge_t *leader = equivEdge(map, e);
+ if (leader != e) {
+ ED_count(leader)++;
+ ED_to_virt(e) = ED_to_virt(leader);
+ ED_to_virt(leader) = e;
+ }
+ }
+ }
+ dtclose(map);
+
+ if (edgefn(g, &margin, edgetype))
+ return 1;
+
+ State = GVSPLINES;
+ return 0;
+}
+
+/* spline_edges1:
+ * Construct edges using default algorithm and given splines value.
+ * Return 0 on success.
+ */
+int spline_edges1(graph_t * g, int edgetype)
+{
+ return splineEdges(g, _spline_edges, edgetype);
+}
+
+/* spline_edges0:
+ * Sets the graph's aspect ratio.
+ * Check splines attribute and construct edges using default algorithm.
+ * If the splines attribute is defined but equal to "", skip edge routing.
+ *
+ * Assumes u.bb for has been computed for g and all clusters
+ * (not just top-level clusters), and that GD_bb(g).LL is at the origin.
+ *
+ * This last criterion is, I believe, mainly to simplify the code
+ * in neato_set_aspect. It would be good to remove this constraint,
+ * as this would allow nodes pinned on input to have the same coordinates
+ * when output in dot or plain format.
+ *
+ */
+void spline_edges0(graph_t * g)
+{
+ int et = EDGE_TYPE (agroot(g));
+ neato_set_aspect(g);
+ if (et == ET_NONE) return;
+#ifndef ORTHO
+ if (et == ET_ORTHO) {
+ agerr (AGWARN, "Orthogonal edges not yet supported\n");
+ et = ET_PLINE;
+ GD_flags(agroot(g)) &= ~ET_ORTHO;
+ GD_flags(agroot(g)) |= ET_PLINE;
+ }
+#endif
+ spline_edges1(g, et);
+}
+
+/* spline_edges:
+ * Compute bounding box, translate graph to origin,
+ * then construct all edges. We assume the graph
+ * has no clusters, and only nodes have been
+ * positioned.
+ */
+void spline_edges(graph_t * g)
+{
+ node_t *n;
+ pointf offset;
+
+ compute_bb(g);
+ offset = cvt2ptf(GD_bb(g).LL);
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ ND_pos(n)[0] -= offset.x;
+ ND_pos(n)[1] -= offset.y;
+ }
+ GD_bb(g).UR.x -= GD_bb(g).LL.x;
+ GD_bb(g).UR.y -= GD_bb(g).LL.y;
+ GD_bb(g).LL.x = 0;
+ GD_bb(g).LL.y = 0;
+ spline_edges0(g);
+}
+
+/* scaleEdge:
+ * Scale edge by given factor.
+ * Assume ED_spl != NULL.
+ */
+static void scaleEdge(edge_t * e, double xf, double yf)
+{
+ int i, j;
+ point *pt;
+ bezier *bez;
+ point delh, delt;
+
+ delh.x = POINTS(ND_pos(aghead(e))[0] * (xf - 1.0));
+ delh.y = POINTS(ND_pos(aghead(e))[1] * (yf - 1.0));
+ delt.x = POINTS(ND_pos(agtail(e))[0] * (xf - 1.0));
+ delt.y = POINTS(ND_pos(agtail(e))[1] * (yf - 1.0));
+
+ bez = ED_spl(e)->list;
+ for (i = 0; i < ED_spl(e)->size; i++) {
+ pt = bez->list;
+ for (j = 0; j < bez->size; j++) {
+ if ((i == 0) && (j == 0)) {
+ pt->x += delt.x;
+ pt->y += delt.y;
+ }
+ else if ((i == ED_spl(e)->size-1) && (j == bez->size-1)) {
+ pt->x += delh.x;
+ pt->y += delh.y;
+ }
+ else {
+ pt->x *= xf;
+ pt->y *= yf;
+ }
+ pt++;
+ }
+ if (bez->sflag) {
+ bez->sp.x += delt.x;
+ bez->sp.y += delt.y;
+ }
+ if (bez->eflag) {
+ bez->ep.x += delh.x;
+ bez->ep.y += delh.y;
+ }
+ bez++;
+ }
+
+ if (ED_label(e) && ED_label(e)->set) {
+ ED_label(e)->p.x *= xf;
+ ED_label(e)->p.y *= yf;
+ }
+ if (ED_head_label(e) && ED_head_label(e)->set) {
+ ED_head_label(e)->p.x += delh.x;
+ ED_head_label(e)->p.y += delh.y;
+ }
+ if (ED_tail_label(e) && ED_tail_label(e)->set) {
+ ED_tail_label(e)->p.x += delt.x;
+ ED_tail_label(e)->p.y += delt.y;
+ }
+}
+
+/* scaleBB:
+ * Scale bounding box of clusters of g by given factors.
+ */
+static void scaleBB(graph_t * g, double xf, double yf)
+{
+ int i;
+
+ GD_bb(g).UR.x *= xf;
+ GD_bb(g).UR.y *= yf;
+ GD_bb(g).LL.x *= xf;
+ GD_bb(g).LL.y *= yf;
+
+ if (GD_label(g) && GD_label(g)->set) {
+ GD_label(g)->p.x *= xf;
+ GD_label(g)->p.y *= yf;
+ }
+
+ for (i = 1; i <= GD_n_cluster(g); i++)
+ scaleBB(GD_clust(g)[i], xf, yf);
+}
+
+/* _neato_set_aspect;
+ * Assume all bounding boxes are correct and
+ * that GD_bb(g).LL is at origin.
+ */
+static void _neato_set_aspect(graph_t * g)
+{
+ /* int i; */
+ double xf, yf, actual, desired;
+ node_t *n;
+
+ /* compute_bb(g); */
+ if (GD_drawing(g)->ratio_kind) {
+ /* normalize */
+ assert(GD_bb(g).LL.x == 0);
+ assert(GD_bb(g).LL.y == 0);
+ if (GD_flip(g)) {
+ int t = GD_bb(g).UR.x;
+ GD_bb(g).UR.x = GD_bb(g).UR.y;
+ GD_bb(g).UR.y = t;
+ }
+ if (GD_drawing(g)->ratio_kind == R_FILL) {
+ /* fill is weird because both X and Y can stretch */
+ if (GD_drawing(g)->size.x <= 0)
+ return;
+ xf = (double) GD_drawing(g)->size.x / (double) GD_bb(g).UR.x;
+ yf = (double) GD_drawing(g)->size.y / (double) GD_bb(g).UR.y;
+ /* handle case where one or more dimensions is too big */
+ if ((xf < 1.0) || (yf < 1.0)) {
+ if (xf < yf) {
+ yf = yf / xf;
+ xf = 1.0;
+ } else {
+ xf = xf / yf;
+ yf = 1.0;
+ }
+ }
+ } else if (GD_drawing(g)->ratio_kind == R_EXPAND) {
+ if (GD_drawing(g)->size.x <= 0)
+ return;
+ xf = (double) GD_drawing(g)->size.x / (double) GD_bb(g).UR.x;
+ yf = (double) GD_drawing(g)->size.y / (double) GD_bb(g).UR.y;
+ if ((xf > 1.0) && (yf > 1.0)) {
+ double scale = MIN(xf, yf);
+ xf = yf = scale;
+ } else
+ return;
+ } else if (GD_drawing(g)->ratio_kind == R_VALUE) {
+ desired = GD_drawing(g)->ratio;
+ actual = ((double) GD_bb(g).UR.y) / ((double) GD_bb(g).UR.x);
+ if (actual < desired) {
+ yf = desired / actual;
+ xf = 1.0;
+ } else {
+ xf = actual / desired;
+ yf = 1.0;
+ }
+ } else
+ return;
+ if (GD_flip(g)) {
+ double t = xf;
+ xf = yf;
+ yf = t;
+ }
+
+ if (Nop > 1) {
+ edge_t *e;
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ for (e = agfstout(g, n); e; e = agnxtout(g, e))
+ if (ED_spl(e))
+ scaleEdge(e, xf, yf);
+ }
+ }
+ /* Not relying on neato_nlist here allows us not to have to
+ * allocate it in the root graph and the connected components.
+ */
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ ND_pos(n)[0] = ND_pos(n)[0] * xf;
+ ND_pos(n)[1] = ND_pos(n)[1] * yf;
+ }
+ scaleBB(g, xf, yf);
+ }
+}
+
+/* neato_set_aspect:
+ * Sets aspect ratio if necessary; real work done in _neato_set_aspect;
+ * This also copies the internal layout coordinates (ND_pos) to the
+ * external ones (ND_coord_i).
+ */
+void neato_set_aspect(graph_t * g)
+{
+ node_t *n;
+
+ _neato_set_aspect(g);
+ for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
+ ND_coord_i(n).x = POINTS(ND_pos(n)[0]);
+ ND_coord_i(n).y = POINTS(ND_pos(n)[1]);
+ }
+}
+#else
static boolean spline_merge(node_t * n)
{
return FALSE;
ND_coord_i(n).y = POINTS(ND_pos(n)[1]);
}
}
-
+#endif
}
}
+#ifdef WITH_CGRAPH
+static double doubleattr(void *obj, Agsym_t* index, double defval)
+{
+ double val;
+ if (index < 0)
+ return defval;
+ if (sscanf(agxget(obj, index), "%lf", &val) < 1)
+ return defval;
+ return val;
+}
+
+/* degreeKind;
+ * Returns degree of n ignoring loops and multiedges.
+ * Returns 0, 1 or many (2)
+ * For case of 1, returns other endpoint of edge.
+ */
+static int degreeKind(graph_t * g, node_t * n, node_t ** op)
+{
+ edge_t *ep;
+ int deg = 0;
+ node_t *other = NULL;
+
+ for (ep = agfstedge(g, n); ep; ep = agnxtedge(g, ep, n)) {
+ if (aghead(ep) == aghead(ep))
+ continue; /* ignore loops */
+ if (deg == 1) {
+ if (((agtail(ep) == n) && (aghead(ep) == other)) || /* ignore multiedge */
+ ((agtail(ep) == other) && (aghead(ep) == n)))
+ continue;
+ return 2;
+ } else { /* deg == 0 */
+ if (agtail(ep) == n)
+ other = aghead(ep);
+ else
+ other = agtail(ep);
+ *op = other;
+ deg++;
+ }
+ }
+ return deg;
+}
+
+
+/* prune:
+ * np is at end of a chain. If its degree is 0, remove it.
+ * If its degree is 1, remove it and recurse.
+ * If its degree > 1, stop.
+ * The node next is the next node to be visited during iteration.
+ * If it is equal to a node being deleted, set it to next one.
+ * Delete from root graph, since G may be a connected component subgraph.
+ * Return next.
+ */
+static node_t *prune(graph_t * G, node_t * np, node_t * next)
+{
+ node_t *other;
+ int deg;
+
+ while (np) {
+ deg = degreeKind(G, np, &other);
+ if (deg == 0) {
+ if (next == np)
+ next = agnxtnode(G, np);
+ agdelete(agroot(G), np);
+ np = 0;
+ } else if (deg == 1) {
+ if (next == np)
+ next = agnxtnode(G, np);
+ agdelete(agroot(G), np);
+ np = other;
+ } else
+ np = 0;
+
+ }
+ return next;
+}
+
+static double setEdgeLen(graph_t * G, node_t * np, Agsym_t* lenx)
+{
+ edge_t *ep;
+ double total_len = 0.0;
+ double len;
+ for (ep = agfstout(G, np); ep; ep = agnxtout(G, ep)) {
+ len = doubleattr(ep, lenx, 1.0);
+ if (len <= 0) {
+ agerr(AGERR, "bad edge len %f in %s ignored\n", len,agnameof(G));
+ len = 1.0;
+ }
+ ED_dist(ep) = len;
+ total_len += len;
+ }
+ return total_len;
+}
+
+/* scan_graph_mode:
+ * Prepare the graph and data structures depending on the layout mode.
+ * If Reduce is true, eliminate singletons and trees. Since G may be a
+ * subgraph, we remove the nodes from the root graph.
+ * Return the number of nodes in the reduced graph.
+ */
+int scan_graph_mode(graph_t * G, int mode)
+{
+ int i, nV, nE, deg;
+ Agsym_t* lenx;
+ char *str;
+ node_t *np, *xp, *other;
+ double total_len = 0.0;
+
+ if (Verbose)
+ fprintf(stderr, "Scanning graph %s, %d nodes\n", agnameof(G),
+ agnnodes(G));
+
+ /* Eliminate singletons and trees */
+ if (Reduce) {
+ for (np = agfstnode(G); np; np = xp) {
+ xp = agnxtnode(G, np);
+ deg = degreeKind(G, np, &other);
+ if (deg == 0) { /* singleton node */
+ agdelete(agroot(G), np);
+ } else if (deg == 1) {
+ agdelete(agroot(G), np);
+ xp = prune(G, other, xp);
+ }
+ }
+ }
+ nV = agnnodes(G);
+ nE = agnedges(G);
+
+ lenx = agattrsym(agroot(G),"len");
+ if (mode == MODE_KK) {
+ Epsilon = .0001 * nV;
+ getdouble(G, "epsilon", &Epsilon);
+ if ((str = agget(agroot(G), "Damping")))
+ Damping = atof(str);
+ else
+ Damping = .99;
+ GD_neato_nlist(G) = N_NEW(nV + 1, node_t *);
+ for (i = 0, np = agfstnode(G); np; np = agnxtnode(G, np)) {
+ GD_neato_nlist(G)[i] = np;
+ ND_id(np) = i++;
+ ND_heapindex(np) = -1;
+ total_len += setEdgeLen(G, np, lenx);
+ }
+ } else {
+ Epsilon = DFLT_TOLERANCE;
+ getdouble(G, "epsilon", &Epsilon);
+ for (i = 0, np = agfstnode(G); np; np = agnxtnode(G, np)) {
+ ND_id(np) = i++;
+ total_len += setEdgeLen(G, np, lenx);
+ }
+ }
+
+ str = agget(G, "defaultdist");
+ if (str && str[0])
+ Initial_dist = MAX(Epsilon, atof(str));
+ else
+ Initial_dist = total_len / (nE > 0 ? nE : 1) * sqrt(nV) + 1;
+
+ if (!Nop && (mode == MODE_KK)) {
+ GD_dist(G) = new_array(nV, nV, Initial_dist);
+ GD_spring(G) = new_array(nV, nV, 1.0);
+ GD_sum_t(G) = new_array(nV, Ndim, 1.0);
+ GD_t(G) = new_3array(nV, nV, Ndim, 0.0);
+ }
+
+ return nV;
+}
+
+int scan_graph(graph_t * g)
+{
+ return scan_graph_mode(g, MODE_KK);
+}
+
+void free_scan_graph(graph_t * g)
+{
+ /* int nG; */
+ free(GD_neato_nlist(g));
+ if (!Nop) {
+ free_array(GD_dist(g));
+ free_array(GD_spring(g));
+ free_array(GD_sum_t(g));
+ /* nG = agnnodes(g); */
+ free_3array(GD_t(g));
+ GD_t(g) = NULL;
+ }
+}
+
+void jitter_d(node_t * np, int nG, int n)
+{
+ int k;
+ for (k = n; k < Ndim; k++)
+ ND_pos(np)[k] = nG * drand48();
+}
+
+void jitter3d(node_t * np, int nG)
+{
+ jitter_d(np, nG, 2);
+}
+
+void randompos(node_t * np, int nG)
+{
+ ND_pos(np)[0] = nG * drand48();
+ ND_pos(np)[1] = nG * drand48();
+ if (Ndim > 2)
+ jitter3d(np, nG);
+}
+
+void initial_positions(graph_t * G, int nG)
+{
+ int init, i;
+ node_t *np;
+ static int once = 0;
+
+ if (Verbose)
+ fprintf(stderr, "Setting initial positions\n");
+
+ init = checkStart(G, nG, INIT_RANDOM);
+ if (init == INIT_REGULAR)
+ return;
+ if ((init == INIT_SELF) && (once == 0)) {
+ agerr(AGWARN, "start=%s not supported with mode=self - ignored\n");
+ once = 1;
+ }
+
+ for (i = 0; (np = GD_neato_nlist(G)[i]); i++) {
+ if (hasPos(np))
+ continue;
+ randompos(np, 1);
+ }
+}
+
+void diffeq_model(graph_t * G, int nG)
+{
+ int i, j, k;
+ double dist, **D, **K, del[MAXDIM], f;
+ node_t *vi, *vj;
+ edge_t *e;
+
+ if (Verbose) {
+ fprintf(stderr, "Setting up spring model: ");
+ start_timer();
+ }
+ /* init springs */
+ K = GD_spring(G);
+ D = GD_dist(G);
+ for (i = 0; i < nG; i++) {
+ for (j = 0; j < i; j++) {
+ f = Spring_coeff / (D[i][j] * D[i][j]);
+ if ((e =
+ agfindedge(G, GD_neato_nlist(G)[i],
+ GD_neato_nlist(G)[j])))
+ f = f * ED_factor(e);
+ K[i][j] = K[j][i] = f;
+ }
+ }
+
+ /* init differential equation solver */
+ for (i = 0; i < nG; i++)
+ for (k = 0; k < Ndim; k++)
+ GD_sum_t(G)[i][k] = 0.0;
+
+ for (i = 0; (vi = GD_neato_nlist(G)[i]); i++) {
+ for (j = 0; j < nG; j++) {
+ if (i == j)
+ continue;
+ vj = GD_neato_nlist(G)[j];
+ dist = distvec(ND_pos(vi), ND_pos(vj), del);
+ for (k = 0; k < Ndim; k++) {
+ GD_t(G)[i][j][k] =
+ GD_spring(G)[i][j] * (del[k] -
+ GD_dist(G)[i][j] * del[k] /
+ dist);
+ GD_sum_t(G)[i][k] += GD_t(G)[i][j][k];
+ }
+ }
+ }
+ if (Verbose) {
+ fprintf(stderr, "%.2f sec\n", elapsed_sec());
+ }
+}
+
+/* total_e:
+ * Return 2*energy of system.
+ */
+static double total_e(graph_t * G, int nG)
+{
+ int i, j, d;
+ double e = 0.0; /* 2*energy */
+ double t0; /* distance squared */
+ double t1;
+ node_t *ip, *jp;
+
+ for (i = 0; i < nG - 1; i++) {
+ ip = GD_neato_nlist(G)[i];
+ for (j = i + 1; j < nG; j++) {
+ jp = GD_neato_nlist(G)[j];
+ for (t0 = 0.0, d = 0; d < Ndim; d++) {
+ t1 = (ND_pos(ip)[d] - ND_pos(jp)[d]);
+ t0 += t1 * t1;
+ }
+ e = e + GD_spring(G)[i][j] *
+ (t0 + GD_dist(G)[i][j] * GD_dist(G)[i][j]
+ - 2.0 * GD_dist(G)[i][j] * sqrt(t0));
+ }
+ }
+ return e;
+}
+
+void solve_model(graph_t * G, int nG)
+{
+ node_t *np;
+
+ Epsilon2 = Epsilon * Epsilon;
+
+ while ((np = choose_node(G, nG))) {
+ move_node(G, nG, np);
+ }
+ if (Verbose) {
+ fprintf(stderr, "\nfinal e = %f", total_e(G, nG));
+ fprintf(stderr, " %d%s iterations %.2f sec\n",
+ GD_move(G), (GD_move(G) == MaxIter ? "!" : ""),
+ elapsed_sec());
+ }
+ if (GD_move(G) == MaxIter)
+ agerr(AGWARN, "Max. iterations (%d) reached on graph %s\n",
+ MaxIter, agnameof(G));
+}
+
+void update_arrays(graph_t * G, int nG, int i)
+{
+ int j, k;
+ double del[MAXDIM], dist, old;
+ node_t *vi, *vj;
+
+ vi = GD_neato_nlist(G)[i];
+ for (k = 0; k < Ndim; k++)
+ GD_sum_t(G)[i][k] = 0.0;
+ for (j = 0; j < nG; j++) {
+ if (i == j)
+ continue;
+ vj = GD_neato_nlist(G)[j];
+ dist = distvec(ND_pos(vi), ND_pos(vj), del);
+ for (k = 0; k < Ndim; k++) {
+ old = GD_t(G)[i][j][k];
+ GD_t(G)[i][j][k] =
+ GD_spring(G)[i][j] * (del[k] -
+ GD_dist(G)[i][j] * del[k] / dist);
+ GD_sum_t(G)[i][k] += GD_t(G)[i][j][k];
+ old = GD_t(G)[j][i][k];
+ GD_t(G)[j][i][k] = -GD_t(G)[i][j][k];
+ GD_sum_t(G)[j][k] += (GD_t(G)[j][i][k] - old);
+ }
+ }
+}
+
+#define Msub(i,j) M[(i)*Ndim+(j)]
+void D2E(graph_t * G, int nG, int n, double *M)
+{
+ int i, l, k;
+ node_t *vi, *vn;
+ double scale, sq, t[MAXDIM];
+ double **K = GD_spring(G);
+ double **D = GD_dist(G);
+
+ vn = GD_neato_nlist(G)[n];
+ for (l = 0; l < Ndim; l++)
+ for (k = 0; k < Ndim; k++)
+ Msub(l, k) = 0.0;
+ for (i = 0; i < nG; i++) {
+ if (n == i)
+ continue;
+ vi = GD_neato_nlist(G)[i];
+ sq = 0.0;
+ for (k = 0; k < Ndim; k++) {
+ t[k] = ND_pos(vn)[k] - ND_pos(vi)[k];
+ sq += (t[k] * t[k]);
+ }
+ scale = 1 / fpow32(sq);
+ for (k = 0; k < Ndim; k++) {
+ for (l = 0; l < k; l++)
+ Msub(l, k) += K[n][i] * D[n][i] * t[k] * t[l] * scale;
+ Msub(k, k) +=
+ K[n][i] * (1.0 - D[n][i] * (sq - (t[k] * t[k])) * scale);
+ }
+ }
+ for (k = 1; k < Ndim; k++)
+ for (l = 0; l < k; l++)
+ Msub(k, l) = Msub(l, k);
+}
+
+void final_energy(graph_t * G, int nG)
+{
+ fprintf(stderr, "iterations = %d final e = %f\n", GD_move(G),
+ total_e(G, nG));
+}
+
+node_t *choose_node(graph_t * G, int nG)
+{
+ int i, k;
+ double m, max;
+ node_t *choice, *np;
+ static int cnt = 0;
+#if 0
+ double e;
+ static double save_e = MAXDOUBLE;
+#endif
+
+ cnt++;
+ if (GD_move(G) >= MaxIter)
+ return NULL;
+#if 0
+ if ((cnt % 100) == 0) {
+ e = total_e(G, nG);
+ if (e - save_e > 0)
+ return NULL;
+ save_e = e;
+ }
+#endif
+ max = 0.0;
+ choice = NULL;
+ for (i = 0; i < nG; i++) {
+ np = GD_neato_nlist(G)[i];
+ if (ND_pinned(np) > P_SET)
+ continue;
+ for (m = 0.0, k = 0; k < Ndim; k++)
+ m += (GD_sum_t(G)[i][k] * GD_sum_t(G)[i][k]);
+ /* could set the color=energy of the node here */
+ if (m > max) {
+ choice = np;
+ max = m;
+ }
+ }
+ if (max < Epsilon2)
+ choice = NULL;
+ else {
+ if (Verbose && (cnt % 100 == 0)) {
+ fprintf(stderr, "%.3f ", sqrt(max));
+ if (cnt % 1000 == 0)
+ fprintf(stderr, "\n");
+ }
+#if 0
+ e = total_e(G, nG);
+ if (fabs((e - save_e) / save_e) < 1e-5) {
+ choice = NULL;
+ }
+#endif
+ }
+ return choice;
+}
+
+void move_node(graph_t * G, int nG, node_t * n)
+{
+ int i, m;
+ static double *a, b[MAXDIM], c[MAXDIM];
+
+ m = ND_id(n);
+ a = ALLOC(Ndim * Ndim, a, double);
+ D2E(G, nG, m, a);
+ for (i = 0; i < Ndim; i++)
+ c[i] = -GD_sum_t(G)[m][i];
+ solve(a, b, c, Ndim);
+ for (i = 0; i < Ndim; i++) {
+ b[i] = (Damping + 2 * (1 - Damping) * drand48()) * b[i];
+ ND_pos(n)[i] += b[i];
+ }
+ GD_move(G)++;
+ update_arrays(G, nG, m);
+ if (test_toggle()) {
+ double sum = 0;
+ for (i = 0; i < Ndim; i++) {
+ sum += fabs(b[i]);
+ } /* Why not squared? */
+ sum = sqrt(sum);
+ fprintf(stderr, "%s %.3f\n",agnameof(n), sum);
+ }
+}
+
+static node_t **Heap;
+static int Heapsize;
+static node_t *Src;
+
+void heapup(node_t * v)
+{
+ int i, par;
+ node_t *u;
+
+ for (i = ND_heapindex(v); i > 0; i = par) {
+ par = (i - 1) / 2;
+ u = Heap[par];
+ if (ND_dist(u) <= ND_dist(v))
+ break;
+ Heap[par] = v;
+ ND_heapindex(v) = par;
+ Heap[i] = u;
+ ND_heapindex(u) = i;
+ }
+}
+
+void heapdown(node_t * v)
+{
+ int i, left, right, c;
+ node_t *u;
+
+ i = ND_heapindex(v);
+ while ((left = 2 * i + 1) < Heapsize) {
+ right = left + 1;
+
+ if ((right < Heapsize)
+ && (ND_dist(Heap[right]) < ND_dist(Heap[left])))
+ c = right;
+ else
+ c = left;
+ u = Heap[c];
+ if (ND_dist(v) <= ND_dist(u))
+ break;
+ Heap[c] = v;
+ ND_heapindex(v) = c;
+ Heap[i] = u;
+ ND_heapindex(u) = i;
+ i = c;
+ }
+}
+
+void neato_enqueue(node_t * v)
+{
+ int i;
+
+ assert(ND_heapindex(v) < 0);
+ i = Heapsize++;
+ ND_heapindex(v) = i;
+ Heap[i] = v;
+ if (i > 0)
+ heapup(v);
+}
+
+node_t *neato_dequeue(void)
+{
+ int i;
+ node_t *rv, *v;
+
+ if (Heapsize == 0)
+ return NULL;
+ rv = Heap[0];
+ i = --Heapsize;
+ v = Heap[i];
+ Heap[0] = v;
+ ND_heapindex(v) = 0;
+ if (i > 1)
+ heapdown(v);
+ ND_heapindex(rv) = -1;
+ return rv;
+}
+
+void shortest_path(graph_t * G, int nG)
+{
+ node_t *v;
+
+ Heap = N_NEW(nG + 1, node_t *);
+ if (Verbose) {
+ fprintf(stderr, "Calculating shortest paths: ");
+ start_timer();
+ }
+ for (v = agfstnode(G); v; v = agnxtnode(G, v))
+ s1(G, v);
+ if (Verbose) {
+ fprintf(stderr, "%.2f sec\n", elapsed_sec());
+ }
+ free(Heap);
+}
+
+void s1(graph_t * G, node_t * node)
+{
+ node_t *v, *u;
+ edge_t *e;
+ int t;
+ double f;
+
+ for (t = 0; (v = GD_neato_nlist(G)[t]); t++)
+ ND_dist(v) = Initial_dist;
+ Src = node;
+ ND_dist(Src) = 0;
+ ND_hops(Src) = 0;
+ neato_enqueue(Src);
+
+ while ((v = neato_dequeue())) {
+ if (v != Src)
+ make_spring(G, Src, v, ND_dist(v));
+ for (e = agfstedge(G, v); e; e = agnxtedge(G, e, v)) {
+ if ((u = aghead(e)) == v)
+ u = agtail(e);
+ f = ND_dist(v) + ED_dist(e);
+ if (ND_dist(u) > f) {
+ ND_dist(u) = f;
+ if (ND_heapindex(u) >= 0)
+ heapup(u);
+ else {
+ ND_hops(u) = ND_hops(v) + 1;
+ neato_enqueue(u);
+ }
+ }
+ }
+ }
+}
+
+void make_spring(graph_t * G, node_t * u, node_t * v, double f)
+{
+ int i, j;
+
+ i = ND_id(u);
+ j = ND_id(v);
+ GD_dist(G)[i][j] = GD_dist(G)[j][i] = f;
+}
+
+int allow_edits(int nsec)
+{
+#ifdef INTERACTIVE
+ static int onetime = TRUE;
+ static FILE *fp;
+ static fd_set fd;
+ static struct timeval tv;
+
+ char buf[256], name[256];
+ double x, y;
+ node_t *np;
+
+ if (onetime) {
+ fp = fopen("/dev/tty", "r");
+ if (fp == NULL)
+ exit(1);
+ setbuf(fp, NULL);
+ tv.tv_usec = 0;
+ onetime = FALSE;
+ }
+ tv.tv_sec = nsec;
+ FD_ZERO(&fd);
+ FD_SET(fileno(fp), &fd);
+ if (select(32, &fd, (fd_set *) 0, (fd_set *) 0, &tv) > 0) {
+ fgets(buf, sizeof(buf), fp);
+ switch (buf[0]) {
+ case 'm': /* move node */
+ if (sscanf(buf + 1, "%s %lf%lf", name, &x, &y) == 3) {
+ np = getnode(G, name);
+ if (np) {
+ NP_pos(np)[0] = x;
+ NP_pos(np)[1] = y;
+ diffeq_model();
+ }
+ }
+ break;
+ case 'q':
+ return FALSE;
+ default:
+ agerr(AGERR, "unknown command '%s', ignored\n", buf);
+ }
+ return TRUE;
+ }
+#endif
+ return FALSE;
+}
+
+#else
double doubleattr(void *obj, int index, double defval)
{
double val;
return defval;
return val;
}
-
/* degreeKind;
* Returns degree of n ignoring loops and multiedges.
* Returns 0, 1 or many (2)
#endif
return FALSE;
}
+#endif
\ No newline at end of file
projects / graphviz / commitdiff