← Back to branch summary

~ubuntu-branches/ubuntu/precise/graphviz/precise-security

« back to all changes in this revision

Viewing changes to lib/neatogen/multispline.c

  • Committer: Bazaar Package Importer
  • Author(s): David Claughton
  • Date: 2010-03-24 22:45:18 UTC
  • mfrom: (1.2.7 upstream) (6.1.7 sid)
  • Revision ID: james.westby@ubuntu.com-20100324224518-do441tthbqjaqjzd
Tags: 2.26.3-4
http://bugs.debian.org/575255
Add patch to fix segfault in circo. Backported from upstream snapshot
release.  Thanks to Francis Russell for his work on this.
(Closes: #575255)

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/neatogen/multispline.c
 
1
/* $Id: multispline.c,v 1.18 2009/08/11 21:21:20 erg Exp $Revision: */
 
2
/* vim:set shiftwidth=4 ts=8: */
 
3
 
 
4
/**********************************************************
 
5
*      This software is part of the graphviz package      *
 
6
*                http://www.graphviz.org/                 *
 
7
*                                                         *
 
8
*            Copyright (c) 1994-2004 AT&T Corp.           *
 
9
*                and is licensed under the                *
 
10
*            Common Public License, Version 1.0           *
 
11
*                      by AT&T Corp.                      *
 
12
*                                                         *
 
13
*        Information and Software Systems Research        *
 
14
*              AT&T Research, Florham Park NJ             *
 
15
**********************************************************/
 
16
 
 
17
#include <multispline.h>
 
18
#include <delaunay.h>
 
19
#include <neatoprocs.h>
 
20
#include <math.h>
 
21
 
 
22
 
 
23
static boolean spline_merge(node_t * n)
 
24
{
 
25
    return FALSE;
 
26
}
 
27
 
 
28
static boolean swap_ends_p(edge_t * e)
 
29
{
 
30
    return FALSE;
 
31
}
 
32
 
 
33
static splineInfo sinfo = { swap_ends_p, spline_merge };
 
34
 
 
35
typedef struct {
 
36
    int i, j;
 
37
} ipair;
 
38
 
 
39
typedef struct _tri {
 
40
    ipair v;
 
41
    struct _tri *nxttri;
 
42
} tri;
 
43
 
 
44
typedef struct {
 
45
    Ppoly_t poly;       /* base polygon used for routing an edge */
 
46
    tri **triMap;       /* triMap[j] is list of all opposite sides of 
 
47
                           triangles containing vertex j, represented
 
48
                           as the indices of the two points in poly */
 
49
} tripoly_t;
 
50
 
 
51
/*
 
52
 * Support for map from I x I -> I
 
53
 */
 
54
typedef struct {
 
55
    Dtlink_t link;              /* cdt data */
 
56
    int a[2];                   /* key */
 
57
    int t;
 
58
} item;
 
59
 
 
60
static int cmpItem(Dt_t * d, int p1[], int p2[], Dtdisc_t * disc)
 
61
{
 
62
    NOTUSED(d);
 
63
    NOTUSED(disc);
 
64
 
 
65
    if (p1[0] < p2[0])
 
66
        return -1;
 
67
    else if (p1[0] > p2[0])
 
68
        return 1;
 
69
    else if (p1[1] < p2[1])
 
70
        return -1;
 
71
    else if (p1[1] > p2[1])
 
72
        return 1;
 
73
    else
 
74
        return 0;
 
75
}
 
76
 
 
77
/* newItem:
 
78
 */
 
79
static void *newItem(Dt_t * d, item * objp, Dtdisc_t * disc)
 
80
{
 
81
    item *newp = NEW(item);
 
82
 
 
83
    NOTUSED(disc);
 
84
    newp->a[0] = objp->a[0];
 
85
    newp->a[1] = objp->a[1];
 
86
    newp->t = objp->t;
 
87
 
 
88
    return newp;
 
89
}
 
90
 
 
91
static void freeItem(Dt_t * d, item * obj, Dtdisc_t * disc)
 
92
{
 
93
    free(obj);
 
94
}
 
95
 
 
96
static Dtdisc_t itemdisc = {
 
97
    offsetof(item, a),
 
98
    2 * sizeof(int),
 
99
    offsetof(item, link),
 
100
    (Dtmake_f) newItem,
 
101
    (Dtfree_f) freeItem,
 
102
    (Dtcompar_f) cmpItem,
 
103
    NIL(Dthash_f),
 
104
    NIL(Dtmemory_f),
 
105
    NIL(Dtevent_f)
 
106
};
 
107
 
 
108
static void addMap(Dt_t * map, int a, int b, int t)
 
109
{
 
110
    item it;
 
111
    int tmp;
 
112
    if (a > b) {
 
113
        tmp = a;
 
114
        a = b;
 
115
        b = tmp;
 
116
    }
 
117
    it.a[0] = a;
 
118
    it.a[1] = b;
 
119
    it.t = t;
 
120
    dtinsert(map, &it);
 
121
}
 
122
 
 
123
/* mapSegToTri:
 
124
 * Create mapping from indices of side endpoints to triangle id 
 
125
 * We use a set rather than a bag because the segments used for lookup
 
126
 * will always be a side of a polygon and hence have a unique triangle.
 
127
 */
 
128
static Dt_t *mapSegToTri(surface_t * sf)
 
129
{
 
130
    Dt_t *map = dtopen(&itemdisc, Dtoset);
 
131
    int i, a, b, c;
 
132
    int *ps = sf->faces;
 
133
    for (i = 0; i < sf->nfaces; i++) {
 
134
        a = *ps++;
 
135
        b = *ps++;
 
136
        c = *ps++;
 
137
        addMap(map, a, b, i);
 
138
        addMap(map, b, c, i);
 
139
        addMap(map, c, a, i);
 
140
    }
 
141
    return map;
 
142
}
 
143
 
 
144
static int findMap(Dt_t * map, int a, int b)
 
145
{
 
146
    item it;
 
147
    item *ip;
 
148
    if (a > b) {
 
149
        int tmp = a;
 
150
        a = b;
 
151
        b = tmp;
 
152
    }
 
153
    it.a[0] = a;
 
154
    it.a[1] = b;
 
155
    ip = (item *) dtsearch(map, &it);
 
156
    assert(ip);
 
157
    return ip->t;
 
158
}
 
159
 
 
160
/*
 
161
 * Support for map from I -> I
 
162
 */
 
163
 
 
164
typedef struct {
 
165
    Dtlink_t link;              /* cdt data */
 
166
    int i;                      /* key */
 
167
    int j;
 
168
} Ipair;
 
169
 
 
170
static int cmpIpair(Dt_t * d, int *p1, int *p2, Dtdisc_t * disc)
 
171
{
 
172
    NOTUSED(d);
 
173
    NOTUSED(disc);
 
174
 
 
175
    return (*p1 - *p2);
 
176
}
 
177
 
 
178
static void *newIpair(Dt_t * d, Ipair * objp, Dtdisc_t * disc)
 
179
{
 
180
    Ipair *newp = NEW(Ipair);
 
181
 
 
182
    NOTUSED(disc);
 
183
    newp->i = objp->i;
 
184
    newp->j = objp->j;
 
185
 
 
186
    return newp;
 
187
}
 
188
 
 
189
static void freeIpair(Dt_t * d, Ipair * obj, Dtdisc_t * disc)
 
190
{
 
191
    free(obj);
 
192
}
 
193
 
 
194
static Dtdisc_t ipairdisc = {
 
195
    offsetof(Ipair, i),
 
196
    sizeof(int),
 
197
    offsetof(Ipair, link),
 
198
    (Dtmake_f) newIpair,
 
199
    (Dtfree_f) freeIpair,
 
200
    (Dtcompar_f) cmpIpair,
 
201
    NIL(Dthash_f),
 
202
    NIL(Dtmemory_f),
 
203
    NIL(Dtevent_f)
 
204
};
 
205
 
 
206
static void vmapAdd(Dt_t * map, int i, int j)
 
207
{
 
208
    Ipair obj;
 
209
    obj.i = i;
 
210
    obj.j = j;
 
211
    dtinsert(map, &obj);
 
212
}
 
213
 
 
214
static int vMap(Dt_t * map, int i)
 
215
{
 
216
    Ipair *ip;
 
217
    ip = (Ipair *) dtmatch(map, &i);
 
218
    return ip->j;
 
219
}
 
220
 
 
221
/* mapTri:
 
222
 * Map vertex indices from router_t to tripoly_t coordinates.
 
223
 */
 
224
static void mapTri(Dt_t * map, tri * tp)
 
225
{
 
226
    for (; tp; tp = tp->nxttri) {
 
227
        tp->v.i = vMap(map, tp->v.i);
 
228
        tp->v.j = vMap(map, tp->v.j);
 
229
    }
 
230
}
 
231
 
 
232
/* addTri:
 
233
 */
 
234
static tri *
 
235
addTri(int i, int j, tri * oldp)
 
236
{
 
237
    tri *tp = NEW(tri);
 
238
    tp->v.i = i;
 
239
    tp->v.j = j;
 
240
    tp->nxttri = oldp;
 
241
    return tp;
 
242
}
 
243
 
 
244
/* bisect:
 
245
 * Return the angle bisecting the angle pp--cp--np
 
246
 */
 
247
static double bisect(pointf pp, pointf cp, pointf np)
 
248
{
 
249
    double theta, phi;
 
250
    theta = atan2(np.y - cp.y, np.x - cp.x);
 
251
    phi = atan2(pp.y - cp.y, pp.x - cp.x);
 
252
    return (theta + phi) / 2.0;
 
253
}
 
254
 
 
255
/* raySeg:
 
256
 * Check if ray v->w intersects segment a--b.
 
257
 */
 
258
static int raySeg(pointf v, pointf w, pointf a, pointf b)
 
259
{
 
260
    int wa = wind(v, w, a);
 
261
    int wb = wind(v, w, b);
 
262
    if (wa == wb)
 
263
        return 0;
 
264
    if (wa == 0) {
 
265
        return (wind(v, b, w) * wind(v, b, a) >= 0);
 
266
    } else {
 
267
        return (wind(v, a, w) * wind(v, a, b) >= 0);
 
268
    }
 
269
}
 
270
 
 
271
/* raySegIntersect:
 
272
 * Find the point p where ray v->w intersects segment ai-bi, if any.
 
273
 * Return 1 on success, 0 on failure
 
274
 */
 
275
static int
 
276
raySegIntersect(pointf v, pointf w, pointf a, pointf b, pointf * p)
 
277
{
 
278
    if (raySeg(v, w, a, b))
 
279
        return line_intersect(v, w, a, b, p);
 
280
    else
 
281
        return 0;
 
282
}
 
283
 
 
284
/* triPoint:
 
285
 * Given the triangle vertex v, and point w so that v->w points
 
286
 * into the polygon, return where the ray v->w intersects the
 
287
 * polygon. The search uses all of the opposite sides of triangles
 
288
 * with v as vertex.
 
289
 * Return 0 on success; 1 on failure.
 
290
 */
 
291
static int
 
292
triPoint(tripoly_t * trip, int vx, pointf v, pointf w, pointf * ip)
 
293
{
 
294
    tri *tp;
 
295
 
 
296
    for (tp = trip->triMap[vx]; tp; tp = tp->nxttri) {
 
297
        if (raySegIntersect
 
298
            (v, w, trip->poly.ps[tp->v.i], trip->poly.ps[tp->v.j], ip))
 
299
            return 0;
 
300
    }
 
301
#ifdef DEBUG
 
302
    fprintf(stderr, "%%Failure in triPoint\n");
 
303
    fprintf(stderr, "0 0 1 setrgbcolor\n");
 
304
    drawLine(v, w);
 
305
    for (tp = trip->triMap[vx]; tp; tp = tp->nxttri) {
 
306
        drawTri(v, trip->poly.ps[tp->v.i], trip->poly.ps[tp->v.j]);
 
307
    }
 
308
#endif
 
309
 
 
310
    return 1;
 
311
}
 
312
 
 
313
/* ctrlPtIdx:
 
314
 * Find the index of v in the points polys->ps.
 
315
 * We start at 1 since the point corresponding to 0 
 
316
 * will never be used as v.
 
317
 */
 
318
static int ctrlPtIdx(pointf v, Ppoly_t * polys)
 
319
{
 
320
    pointf w;
 
321
    int i;
 
322
    for (i = 1; i < polys->pn; i++) {
 
323
        w = polys->ps[i];
 
324
        if ((w.x == v.x) && (w.y == v.y))
 
325
            return i;
 
326
    }
 
327
    return -1;
 
328
}
 
329
 
 
330
#define SEP 15
 
331
 
 
332
/* mkCtrlPts:
 
333
 * Generate mult points associated with v.
 
334
 * The points will lie on the ray bisecting the angle prev--v--nxt.
 
335
 * The first point will aways be v.
 
336
 * The rest are positioned equally spaced with maximum spacing SEP.
 
337
 * In addition, they all lie within the polygon trip->poly.
 
338
 * Parameter s gives the index after which a vertex likes on the
 
339
 * opposite side. This is necessary to get the "curvature" of the
 
340
 * path correct.
 
341
 */
 
342
static pointf *mkCtrlPts(int s, int mult, pointf prev, pointf v,
 
343
                           pointf nxt, tripoly_t * trip)
 
344
{
 
345
    pointf *ps;
 
346
    int idx = ctrlPtIdx(v, &(trip->poly));
 
347
    int i;
 
348
    double d, sep, theta, sinTheta, cosTheta;
 
349
    pointf q, w;
 
350
 
 
351
    if (idx < 0)
 
352
        return NULL;
 
353
 
 
354
    ps = N_GNEW(mult, pointf);
 
355
    theta = bisect(prev, v, nxt);
 
356
    sinTheta = sin(theta);
 
357
    cosTheta = cos(theta);
 
358
    w.x = v.x + 100 * cosTheta;
 
359
    w.y = v.y + 100 * sinTheta;
 
360
    if ((idx > s) && (wind(prev, v, w) != 1)) {
 
361
        sinTheta *= -1;
 
362
        cosTheta *= -1;
 
363
        w.x = v.x + 100 * cosTheta;
 
364
        w.y = v.y + 100 * sinTheta;
 
365
    } else if (wind(prev, v, w) != -1) {
 
366
        sinTheta *= -1;
 
367
        cosTheta *= -1;
 
368
        w.x = v.x + 100 * cosTheta;
 
369
        w.y = v.y + 100 * sinTheta;
 
370
    }
 
371
 
 
372
    if (triPoint(trip, idx, v, w, &q)) {
 
373
        return 0;
 
374
    }
 
375
 
 
376
    d = DIST(q, v);
 
377
    if (d >= mult * SEP)
 
378
        sep = SEP;
 
379
    else
 
380
        sep = d / mult;
 
381
    if (idx < s) {
 
382
        for (i = 0; i < mult; i++) {
 
383
            ps[i].x = v.x + i * sep * cosTheta;
 
384
            ps[i].y = v.y + i * sep * sinTheta;
 
385
        }
 
386
    } else {
 
387
        for (i = 0; i < mult; i++) {
 
388
            ps[mult - i - 1].x = v.x + i * sep * cosTheta;
 
389
            ps[mult - i - 1].y = v.y + i * sep * sinTheta;
 
390
        }
 
391
    }
 
392
    return ps;
 
393
}
 
394
 
 
395
/*
 
396
 * Simple graph structure for recording the triangle graph.
 
397
 */
 
398
 
 
399
typedef struct {
 
400
    int ne;         /* no. of edges. */
 
401
    int *edges;     /* indices of edges adjacent to node. */
 
402
    pointf ctr;     /* center of triangle. */
 
403
} tnode;
 
404
 
 
405
typedef struct {
 
406
    int t, h;       /* indices of head and tail nodes */
 
407
    ipair seg;      /* indices of points forming shared segment */
 
408
    double dist;    /* length of edge; usually distance between centers */
 
409
} tedge;
 
410
 
 
411
typedef struct {
 
412
    tnode *nodes;
 
413
    tedge *edges;
 
414
    int nedges;    /* no. of edges; no. of nodes is stored in router_t */
 
415
} tgraph;
 
416
 
 
417
struct router_s {
 
418
    int pn;     /* no. of points */
 
419
    pointf *ps; /* all points in configuration */
 
420
    int *obs;   /* indices in obstacle i are obs[i]...obs[i+1]-1 */
 
421
    int *tris;  /* indices of triangle i are tris[3*i]...tris[3*i+2] */
 
422
    Dt_t *trimap; /* map from obstacle side (a,b) to index of adj. triangle */
 
423
    int tn;       /* no. of nodes in tg */
 
424
    tgraph *tg;   /* graph of triangles */
 
425
};
 
426
 
 
427
/* triCenter:
 
428
 * Given an array of points and 3 integer indices,
 
429
 * compute and return the center of the triangle.
 
430
 */
 
431
static pointf triCenter(pointf * pts, int *idxs)
 
432
{
 
433
    pointf a = pts[*idxs++];
 
434
    pointf b = pts[*idxs++];
 
435
    pointf c = pts[*idxs++];
 
436
    pointf p;
 
437
    p.x = (a.x + b.x + c.x) / 3.0;
 
438
    p.y = (a.y + b.y + c.y) / 3.0;
 
439
    return p;
 
440
}
 
441
 
 
442
#define MARGIN 32
 
443
 
 
444
/* bbox:
 
445
 * Compute bounding box of polygons, and return it
 
446
 * with an added margin of MARGIN.
 
447
 * Store total number of points in *np.
 
448
 */
 
449
static boxf bbox(Ppoly_t** obsp, int npoly, int *np)
 
450
{
 
451
    boxf bb;
 
452
    int i, j, cnt = 0;
 
453
    pointf p;
 
454
    Ppoly_t* obs;
 
455
 
 
456
    bb.LL.x = bb.LL.y = MAXDOUBLE;
 
457
    bb.UR.x = bb.UR.y = -MAXDOUBLE;
 
458
 
 
459
    for (i = 0; i < npoly; i++) {
 
460
        obs = *obsp++;
 
461
        for (j = 0; j < obs->pn; j++) {
 
462
            p = obs->ps[j];
 
463
            if (p.x < bb.LL.x)
 
464
                bb.LL.x = p.x;
 
465
            if (p.x > bb.UR.x)
 
466
                bb.UR.x = p.x;
 
467
            if (p.y < bb.LL.y)
 
468
                bb.LL.y = p.y;
 
469
            if (p.y > bb.UR.y)
 
470
                bb.UR.y = p.y;
 
471
            cnt++;
 
472
        }
 
473
    }
 
474
 
 
475
    *np = cnt;
 
476
 
 
477
    bb.LL.x -= MARGIN;
 
478
    bb.LL.y -= MARGIN;
 
479
    bb.UR.x += MARGIN;
 
480
    bb.UR.y += MARGIN;
 
481
 
 
482
    return bb;
 
483
}
 
484
 
 
485
static int *mkTriIndices(surface_t * sf)
 
486
{
 
487
    int *tris = N_GNEW(3 * sf->nfaces, int);
 
488
    memcpy(tris, sf->faces, 3 * sf->nfaces * sizeof(int));
 
489
    return tris;
 
490
}
 
491
 
 
492
/* degT:
 
493
 * Given a pointer to an array of 3 integers, return
 
494
 * the number not equal to -1
 
495
 */
 
496
static int degT(int *ip)
 
497
{
 
498
    ip++;
 
499
    if (*ip++ == -1)
 
500
        return 1;
 
501
    if (*ip == -1)
 
502
        return 2;
 
503
    else
 
504
        return 3;
 
505
}
 
506
 
 
507
/* sharedEdge:
 
508
 * Returns a pair of integer (x,y), x < y, where x and y are the
 
509
 * indices of the two vertices of the shared edge.
 
510
 */
 
511
static ipair sharedEdge(int *p, int *q)
 
512
{
 
513
    ipair pt;
 
514
    int tmp, p1, p2;
 
515
    p1 = *p;
 
516
    p2 = *(p + 1);
 
517
    if (p1 == *q) {
 
518
        if ((p2 != *(q + 1)) && (p2 != *(q + 2))) {
 
519
            p2 = *(p + 2);
 
520
        }
 
521
    } else if (p1 == *(q + 1)) {
 
522
        if ((p2 != *q) && (p2 != *(q + 2))) {
 
523
            p2 = *(p + 2);
 
524
        }
 
525
    } else if (p1 == *(q + 2)) {
 
526
        if ((p2 != *q) && (p2 != *(q + 1))) {
 
527
            p2 = *(p + 2);
 
528
        }
 
529
    } else {
 
530
        p1 = *(p + 2);
 
531
    }
 
532
 
 
533
    if (p1 > p2) {
 
534
        tmp = p1;
 
535
        p1 = p2;
 
536
        p2 = tmp;
 
537
    }
 
538
    pt.i = p1;
 
539
    pt.j = p2;
 
540
    return pt;
 
541
}
 
542
 
 
543
/* addTriEdge:
 
544
 * Add an edge to g, with tail t, head h, distance d, and shared
 
545
 * segment seg.
 
546
 */
 
547
static void addTriEdge(tgraph * g, int t, int h, double d, ipair seg)
 
548
{
 
549
    tedge *ep = g->edges + g->nedges;
 
550
    tnode *tp = g->nodes + t;
 
551
    tnode *hp = g->nodes + h;
 
552
 
 
553
    ep->t = t;
 
554
    ep->h = h;
 
555
    ep->dist = DIST(tp->ctr, hp->ctr);
 
556
    ep->seg = seg;
 
557
 
 
558
    tp->edges[tp->ne++] = g->nedges;
 
559
    hp->edges[hp->ne++] = g->nedges;
 
560
 
 
561
    g->nedges++;
 
562
}
 
563
 
 
564
static void freeTriGraph(tgraph * tg)
 
565
{
 
566
    free(tg->nodes->edges);
 
567
    free(tg->nodes);
 
568
    free(tg->edges);
 
569
    free(tg);
 
570
}
 
571
 
 
572
/* mkTriGraph:
 
573
 * Generate graph with triangles as nodes and an edge iff two triangles
 
574
 * share an edge.
 
575
 */
 
576
static tgraph *mkTriGraph(surface_t * sf, int maxv, pointf * pts)
 
577
{
 
578
    tgraph *g;
 
579
    tnode *np;
 
580
    int j, i, ne = 0;
 
581
    int *edgei;
 
582
    int *jp;
 
583
 
 
584
    /* ne is twice no. of edges */
 
585
    for (i = 0; i < 3 * sf->nfaces; i++)
 
586
        if (sf->neigh[i] != -1)
 
587
            ne++;
 
588
 
 
589
    g = GNEW(tgraph);
 
590
 
 
591
    /* plus 2 for nodes added as endpoints of an edge */
 
592
    g->nodes = N_GNEW(sf->nfaces + 2, tnode);
 
593
 
 
594
    /* allow 1 possible extra edge per triangle, plus 
 
595
     * obstacles can have at most maxv triangles touching 
 
596
     */
 
597
    edgei = N_GNEW(sf->nfaces + ne + 2 * maxv, int);
 
598
    g->edges = N_GNEW(ne/2 + 2 * maxv, tedge);
 
599
    g->nedges = 0;
 
600
 
 
601
    for (i = 0; i < sf->nfaces; i++) {
 
602
        np = g->nodes + i;
 
603
        np->ne = 0;
 
604
        np->edges = edgei;
 
605
        np->ctr = triCenter(pts, sf->faces + 3 * i);
 
606
        edgei += degT(sf->neigh + 3 * i) + 1;
 
607
    }
 
608
    /* initialize variable nodes */
 
609
    np = g->nodes + i;
 
610
    np->ne = 0;
 
611
    np->edges = edgei;
 
612
    np++;
 
613
    np->ne = 0;
 
614
    np->edges = edgei + maxv;
 
615
 
 
616
    for (i = 0; i < sf->nfaces; i++) {
 
617
        np = g->nodes + i;
 
618
        jp = sf->neigh + 3 * i;
 
619
        ne = 0;
 
620
        while ((ne < 3) && ((j = *jp++) != -1)) {
 
621
            if (i < j) {
 
622
                double dist = DIST(np->ctr, (g->nodes + j)->ctr);
 
623
                ipair seg =
 
624
                    sharedEdge(sf->faces + 3 * i, sf->faces + 3 * j);
 
625
                addTriEdge(g, i, j, dist, seg);
 
626
            }
 
627
            ne++;
 
628
        }
 
629
    }
 
630
 
 
631
    return g;
 
632
}
 
633
 
 
634
void freeRouter(router_t * rtr)
 
635
{
 
636
    free(rtr->ps);
 
637
    free(rtr->obs);
 
638
    free(rtr->tris);
 
639
    dtclose(rtr->trimap);
 
640
    freeTriGraph(rtr->tg);
 
641
    free(rtr);
 
642
}
 
643
 
 
644
#ifdef DEBUG
 
645
#include <psdbg.c>
 
646
static void
 
647
prTriGraph (router_t* rtr, int n)
 
648
{
 
649
    FILE* fp = fopen ("dump","w");
 
650
    int i;
 
651
    pointf* pts = rtr->ps;
 
652
    tnode* nodes = rtr->tg->nodes;
 
653
    char buf[BUFSIZ];
 
654
    
 
655
    psInit();
 
656
    for (i=0;i < rtr->tn; i++) {
 
657
        pointf a = pts[rtr->tris[3*i]];
 
658
        pointf b = pts[rtr->tris[3*i+1]];
 
659
        pointf c = pts[rtr->tris[3*i+2]];
 
660
        psTri (a, b,c);
 
661
        sprintf (buf, "%d", i);
 
662
        psTxt (buf, nodes[i].ctr);
 
663
    }
 
664
    for (i=rtr->tn;i < n; i++) {
 
665
        psTxt (buf, nodes[i].ctr);
 
666
    }
 
667
    psColor ("1 0 0");
 
668
    for (i=0;i < rtr->tg->nedges; i++) {
 
669
        tedge* e = rtr->tg->edges+i;
 
670
        psSeg (nodes[e->t].ctr, nodes[e->h].ctr);
 
671
    }
 
672
    psOut(fp);
 
673
    fclose(fp);
 
674
}
 
675
#endif
 
676
 
 
677
router_t *mkRouter(Ppoly_t** obsp, int npoly)
 
678
{
 
679
    router_t *rtr = NEW(router_t);
 
680
    Ppoly_t* obs;
 
681
    boxf bb;
 
682
    pointf *pts;
 
683
    int npts;
 
684
    surface_t *sf;
 
685
    int *segs;
 
686
    double *x;
 
687
    double *y;
 
688
    int maxv = 4; /* default max. no. of vertices in an obstacle; set below */
 
689
    /* points in obstacle i have indices obsi[i] through obsi[i+1]-1 in pts
 
690
     */
 
691
    int *obsi = N_NEW(npoly + 1, int);
 
692
    int i, j, ix = 4, six = 0;
 
693
 
 
694
    bb = bbox(obsp, npoly, &npts);
 
695
    npts += 4;                  /* 4 points of bounding box */
 
696
    pts = N_GNEW(npts, pointf); /* all points are stored in pts */
 
697
    segs = N_GNEW(2 * npts, int);       /* indices of points forming segments */
 
698
 
 
699
    /* store bounding box in CCW order */
 
700
    pts[0] = bb.LL;
 
701
    pts[1].x = bb.UR.x;
 
702
    pts[1].y = bb.LL.y;
 
703
    pts[2] = bb.UR;
 
704
    pts[3].x = bb.LL.x;
 
705
    pts[3].y = bb.UR.y;
 
706
    for (i = 1; i <= 4; i++) {
 
707
        segs[six++] = i - 1;
 
708
        if (i < 4)
 
709
            segs[six++] = i;
 
710
        else
 
711
            segs[six++] = 0;
 
712
    }
 
713
 
 
714
    /* store obstacles in CW order and generate constraint segments */
 
715
    for (i = 0; i < npoly; i++) {
 
716
        obsi[i] = ix;
 
717
        obs = *obsp++;
 
718
        for (j = 1; j <= obs->pn; j++) {
 
719
            segs[six++] = ix;
 
720
            if (j < obs->pn)
 
721
                segs[six++] = ix + 1;
 
722
            else
 
723
                segs[six++] = obsi[i];
 
724
            pts[ix++] = obs->ps[j - 1];
 
725
        }
 
726
        if (obs->pn > maxv)
 
727
            maxv = obs->pn;
 
728
    }
 
729
    obsi[i] = ix;
 
730
 
 
731
    /* copy points into coordinate arrays */
 
732
    x = N_GNEW(npts, double);
 
733
    y = N_GNEW(npts, double);
 
734
    for (i = 0; i < npts; i++) {
 
735
        x[i] = pts[i].x;
 
736
        y[i] = pts[i].y;
 
737
    }
 
738
    sf = mkSurface(x, y, npts, segs, npts);
 
739
    free(x);
 
740
    free(y);
 
741
    free(segs);
 
742
 
 
743
    rtr->ps = pts;
 
744
    rtr->pn = npts;
 
745
    rtr->obs = obsi;
 
746
    rtr->tris = mkTriIndices(sf);
 
747
    rtr->trimap = mapSegToTri(sf);
 
748
    rtr->tn = sf->nfaces;
 
749
    rtr->tg = mkTriGraph(sf, maxv, pts);
 
750
 
 
751
    freeSurface(sf);
 
752
    return rtr;
 
753
}
 
754
 
 
755
/* finishEdge:
 
756
 * Finish edge generation, clipping to nodes and adding arrowhead
 
757
 * if necessary, and adding edge labels
 
758
 */
 
759
static void
 
760
finishEdge (graph_t* g, edge_t* e, Ppoly_t spl, int flip, pointf p, pointf q)
 
761
{
 
762
    int j;
 
763
    pointf *spline = N_GNEW(spl.pn, pointf);
 
764
    pointf p1, q1;
 
765
 
 
766
    if (flip) {
 
767
        for (j = 0; j < spl.pn; j++) {
 
768
            spline[spl.pn - 1 - j] = spl.ps[j];
 
769
        }
 
770
        p1 = q;
 
771
        q1 = p;
 
772
    }
 
773
    else {
 
774
        for (j = 0; j < spl.pn; j++) {
 
775
            spline[j] = spl.ps[j];
 
776
        }
 
777
        p1 = p;
 
778
        q1 = q;
 
779
    }
 
780
    if (Verbose > 1)
 
781
        fprintf(stderr, "spline %s %s\n", agnameof(agtail(e)), agnameof(aghead(e)));
 
782
    clip_and_install(e, aghead(e), spline, spl.pn, &sinfo);
 
783
    free(spline);
 
784
 
 
785
    addEdgeLabels(g, e, p1, q1);
 
786
}
 
787
 
 
788
#define EQPT(p,q) (((p).x==(q).x)&&((p).y==(q).y))
 
789
 
 
790
/* tweakEnd:
 
791
 * Hack because path routing doesn't know about the interiors
 
792
 * of polygons. If the first or last segment of the shortest path
 
793
 * lies along one of the polygon boundaries, the path may flip
 
794
 * inside the polygon. To avoid this, we shift the point a bit.
 
795
 *
 
796
 * If the edge p(=poly.ps[s])-q of the shortest path is also an
 
797
 * edge of the border polygon, move p slightly inside the polygon
 
798
 * and return it. If prv and nxt are the two vertices adjacent to
 
799
 * p in the polygon, let m be the midpoint of prv--nxt. We then
 
800
 * move a tiny bit along the ray p->m.
 
801
 *
 
802
 * Otherwise, return p unchanged.
 
803
 */
 
804
static Ppoint_t
 
805
tweakEnd (Ppoly_t poly, int s, Ppolyline_t pl, Ppoint_t q)
 
806
{
 
807
    Ppoint_t prv, nxt, p;
 
808
 
 
809
    p = poly.ps[s];
 
810
    nxt = poly.ps[(s + 1) % poly.pn];
 
811
    if (s == 0)
 
812
        prv = poly.ps[poly.pn-1];
 
813
    else
 
814
        prv = poly.ps[s - 1];
 
815
    if (EQPT(q, nxt) || EQPT(q, prv) ){
 
816
        Ppoint_t m;
 
817
        double d;
 
818
        m.x = (nxt.x + prv.x)/2.0 - p.x;
 
819
        m.y = (nxt.y + prv.y)/2.0 - p.y;
 
820
        d = LEN(m.x,m.y);
 
821
        p.x += 0.1*m.x/d;
 
822
        p.y += 0.1*m.y/d;
 
823
    }
 
824
    return p;
 
825
}
 
826
 
 
827
static void
 
828
tweakPath (Ppoly_t poly, int s, int t, Ppolyline_t pl)
 
829
{
 
830
    pl.ps[0] = tweakEnd (poly, s, pl, pl.ps[1]);
 
831
    pl.ps[pl.pn-1] = tweakEnd (poly, t, pl, pl.ps[pl.pn-2]);
 
832
}
 
833
 
 
834
 
 
835
/* genroute:
 
836
 * Generate splines for e and cohorts.
 
837
 * Edges go from s to t.
 
838
 * Return 0 on success.
 
839
 */
 
840
static int 
 
841
genroute(graph_t* g, tripoly_t * trip, int s, int t, edge_t * e, int doPolyline)
 
842
{
 
843
    pointf eps[2];
 
844
    Pvector_t evs[2];
 
845
    pointf **cpts;              /* lists of control points */
 
846
    Ppoly_t poly;
 
847
    Ppolyline_t pl, spl;
 
848
    int i, j;
 
849
    Ppolyline_t mmpl;
 
850
    Pedge_t *medges = N_GNEW(trip->poly.pn, Pedge_t);
 
851
    int pn;
 
852
    int mult = ED_count(e);
 
853
    node_t* head = aghead(e);
 
854
 
 
855
    eps[0].x = trip->poly.ps[s].x, eps[0].y = trip->poly.ps[s].y;
 
856
    eps[1].x = trip->poly.ps[t].x, eps[1].y = trip->poly.ps[t].y;
 
857
    Pshortestpath(&(trip->poly), eps, &pl);
 
858
 
 
859
    if (pl.pn == 2) {
 
860
        makeStraightEdge(agraphof(head), e, doPolyline);
 
861
        return 0;
 
862
    }
 
863
 
 
864
    evs[0].x = evs[0].y = 0;
 
865
    evs[1].x = evs[1].y = 0;
 
866
 
 
867
    if ((mult == 1) || Concentrate) {
 
868
        poly = trip->poly;
 
869
        for (j = 0; j < poly.pn; j++) {
 
870
            medges[j].a = poly.ps[j];
 
871
            medges[j].b = poly.ps[(j + 1) % poly.pn];
 
872
        }
 
873
        tweakPath (poly, s, t, pl);
 
874
        Proutespline(medges, poly.pn, pl, evs, &spl);
 
875
        finishEdge (g, e, spl, aghead(e) != head, eps[0], eps[1]);
 
876
        free(medges);
 
877
 
 
878
        return 0;
 
879
    }
 
880
    
 
881
    pn = 2 * (pl.pn - 1);
 
882
 
 
883
    cpts = N_NEW(pl.pn - 2, pointf *);
 
884
    for (i = 0; i < pl.pn - 2; i++) {
 
885
        cpts[i] =
 
886
            mkCtrlPts(t, mult+1, pl.ps[i], pl.ps[i + 1], pl.ps[i + 2], trip);
 
887
        if (!cpts[i]) {
 
888
            agerr(AGWARN, "Could not create control points for multiple spline for edge (%s,%s)\n", agnameof(agtail(e)), agnameof(aghead(e)));
 
889
            return 1;
 
890
        }
 
891
    }
 
892
 
 
893
    poly.ps = N_GNEW(pn, pointf);
 
894
    poly.pn = pn;
 
895
 
 
896
    for (i = 0; i < mult; i++) {
 
897
        poly.ps[0] = eps[0];
 
898
        for (j = 1; j < pl.pn - 1; j++) {
 
899
            poly.ps[j] = cpts[j - 1][i];
 
900
        }
 
901
        poly.ps[pl.pn - 1] = eps[1];
 
902
        for (j = 1; j < pl.pn - 1; j++) {
 
903
            poly.ps[pn - j] = cpts[j - 1][i + 1];
 
904
        }
 
905
        Pshortestpath(&poly, eps, &mmpl);
 
906
 
 
907
        if (doPolyline) {
 
908
            make_polyline (mmpl, &spl);
 
909
        }
 
910
        else {
 
911
            for (j = 0; j < poly.pn; j++) {
 
912
                medges[j].a = poly.ps[j];
 
913
                medges[j].b = poly.ps[(j + 1) % poly.pn];
 
914
            }
 
915
            tweakPath (poly, 0, pl.pn-1, mmpl);
 
916
            Proutespline(medges, poly.pn, mmpl, evs, &spl);
 
917
        }
 
918
        finishEdge (g, e, spl, aghead(e) != head, eps[0], eps[1]);
 
919
 
 
920
        e = ED_to_virt(e);
 
921
    }
 
922
 
 
923
    for (i = 0; i < pl.pn - 2; i++)
 
924
        free(cpts[i]);
 
925
    free(cpts);
 
926
    free(medges);
 
927
    free(poly.ps);
 
928
    return 0;
 
929
}
 
930
 
 
931
#define NSMALL -0.0000000001
 
932
 
 
933
/* inCone:
 
934
 * Returns true iff q is in the convex cone a-b-c
 
935
 */
 
936
static int
 
937
inCone (pointf a, pointf b, pointf c, pointf q)
 
938
{
 
939
    return ((area2(q,a,b) >= NSMALL) && (area2(q,b,c) >= NSMALL));
 
940
}
 
941
 
 
942
static pointf north = {0, 1};
 
943
static pointf northeast = {1, 1};
 
944
static pointf east = {1, 0};
 
945
static pointf southeast = {1, -1};
 
946
static pointf south = {0, -1};
 
947
static pointf southwest = {-1, -1};
 
948
static pointf west = {-1, 0};
 
949
static pointf northwest = {-1, 1};
 
950
 
 
951
/* addEndpoint:
 
952
 * Add node to graph representing spline end point p inside obstruction obs_id.
 
953
 * For each side of obstruction, add edge from p to corresponding triangle.
 
954
 * The node id of the new node in the graph is v_id.
 
955
 * If p lies on the side of its node (sides != 0), we limit the triangles
 
956
 * to those within 45 degrees of each side of the natural direction of p.
 
957
 */
 
958
static void addEndpoint(router_t * rtr, pointf p, node_t* v, int v_id, int sides)
 
959
{
 
960
    int obs_id = ND_lim(v);
 
961
    int starti = rtr->obs[obs_id];
 
962
    int endi = rtr->obs[obs_id + 1];
 
963
    pointf* pts = rtr->ps;
 
964
    int i, t;
 
965
    double d;
 
966
    pointf vr, v0, v1;
 
967
 
 
968
    switch (sides) {
 
969
    case TOP :
 
970
        vr = add_pointf (p, north);
 
971
        v0 = add_pointf (p, northwest);
 
972
        v1 = add_pointf (p, northeast);
 
973
        break;
 
974
    case TOP|RIGHT :
 
975
        vr = add_pointf (p, northeast);
 
976
        v0 = add_pointf (p, north);
 
977
        v1 = add_pointf (p, east);
 
978
        break;
 
979
    case RIGHT :
 
980
        vr = add_pointf (p, east);
 
981
        v0 = add_pointf (p, northeast);
 
982
        v1 = add_pointf (p, southeast);
 
983
        break;
 
984
    case BOTTOM|RIGHT :
 
985
        vr = add_pointf (p, southeast);
 
986
        v0 = add_pointf (p, east);
 
987
        v1 = add_pointf (p, south);
 
988
        break;
 
989
    case BOTTOM :
 
990
        vr = add_pointf (p, south);
 
991
        v0 = add_pointf (p, southeast);
 
992
        v1 = add_pointf (p, southwest);
 
993
        break;
 
994
    case BOTTOM|LEFT :
 
995
        vr = add_pointf (p, southwest);
 
996
        v0 = add_pointf (p, south);
 
997
        v1 = add_pointf (p, west);
 
998
        break;
 
999
    case LEFT :
 
1000
        vr = add_pointf (p, west);
 
1001
        v0 = add_pointf (p, southwest);
 
1002
        v1 = add_pointf (p, northwest);
 
1003
        break;
 
1004
    case TOP|LEFT :
 
1005
        vr = add_pointf (p, northwest);
 
1006
        v0 = add_pointf (p, west);
 
1007
        v1 = add_pointf (p, north);
 
1008
        break;
 
1009
    case 0 :
 
1010
        break;
 
1011
    default :
 
1012
        assert (0);
 
1013
        break;
 
1014
    }
 
1015
 
 
1016
    rtr->tg->nodes[v_id].ne = 0;
 
1017
    rtr->tg->nodes[v_id].ctr = p;
 
1018
    for (i = starti; i < endi; i++) {
 
1019
        ipair seg;
 
1020
        seg.i = i;
 
1021
        if (i < endi - 1)
 
1022
            seg.j = i + 1;
 
1023
        else
 
1024
            seg.j = starti;
 
1025
        t = findMap(rtr->trimap, seg.i, seg.j);
 
1026
        if (sides && !inCone (v0, p, v1, pts[seg.i]) && !inCone (v0, p, v1, pts[seg.j]) && !raySeg(p,vr,pts[seg.i],pts[seg.j]))
 
1027
            continue;
 
1028
        d = DIST(p, (rtr->tg->nodes + t)->ctr);
 
1029
        addTriEdge(rtr->tg, v_id, t, d, seg);
 
1030
    }
 
1031
}
 
1032
 
 
1033
/* edgeToSeg:
 
1034
 * Given edge from i to j, find segment associated
 
1035
 * with the edge.
 
1036
 *
 
1037
 * This lookup could be made faster by modifying the 
 
1038
 * shortest path algorithm to store the edges rather than
 
1039
 * the nodes.
 
1040
 */
 
1041
static ipair edgeToSeg(tgraph * tg, int i, int j)
 
1042
{
 
1043
    ipair ip;
 
1044
    tnode *np = tg->nodes + i;
 
1045
    tedge *ep;
 
1046
 
 
1047
    for (i = 0; i < np->ne; i++) {
 
1048
        ep = tg->edges + np->edges[i];
 
1049
        if ((ep->t == j) || (ep->h == j))
 
1050
            return (ep->seg);
 
1051
    }
 
1052
 
 
1053
    assert(0);
 
1054
    return ip;
 
1055
}
 
1056
 
 
1057
static void
 
1058
freeTripoly (tripoly_t* trip)
 
1059
{
 
1060
    int i;
 
1061
    tri* tp;
 
1062
    tri* nxt;
 
1063
 
 
1064
    free (trip->poly.ps);
 
1065
    for (i = 0; i < trip->poly.pn; i++) {
 
1066
        for (tp = trip->triMap[i]; tp; tp = nxt) {
 
1067
            nxt = tp->nxttri;
 
1068
            free (tp);
 
1069
        }
 
1070
    }
 
1071
    free (trip->triMap);
 
1072
    free (trip);
 
1073
}
 
1074
 
 
1075
/* Auxiliary data structure used to translate a path of rectangles
 
1076
 * into a polygon. Each side_t represents a vertex on one side of
 
1077
 * the polygon. v is the index of the vertex in the global router_t,
 
1078
 * and ts is a linked list of the indices of segments of sides opposite
 
1079
 * to v in some triangle on the path. These lists will be translated
 
1080
 * to polygon indices by mapTri, and stored in tripoly_t.triMap.
 
1081
 */
 
1082
typedef struct {
 
1083
    int v;
 
1084
    tri *ts;
 
1085
} side_t;
 
1086
 
 
1087
/* mkPoly:
 
1088
 * Construct simple polygon from shortest path from t to s in g.
 
1089
 * dad gives the indices of the triangles on path.
 
1090
 * sx used to store index of s in points.
 
1091
 * index of t is always 0
 
1092
 */
 
1093
static tripoly_t *mkPoly(router_t * rtr, int *dad, int s, int t,
 
1094
                         pointf p_s, pointf p_t, int *sx)
 
1095
{
 
1096
    tripoly_t *ps;
 
1097
    int nxt;
 
1098
    ipair p;
 
1099
    int nt = 0;
 
1100
    side_t *side1;
 
1101
    side_t *side2;
 
1102
    int i, idx;
 
1103
    int cnt1 = 0;
 
1104
    int cnt2 = 0;
 
1105
    pointf *pts;
 
1106
    pointf *pps;
 
1107
    /* maps vertex index used in router_t to vertex index used in tripoly */
 
1108
    Dt_t *vmap;
 
1109
    tri **trim;
 
1110
 
 
1111
    /* count number of triangles in path */
 
1112
    for (nxt = dad[t]; nxt != s; nxt = dad[nxt])
 
1113
        nt++;
 
1114
 
 
1115
    side1 = N_NEW(nt + 4, side_t);
 
1116
    side2 = N_NEW(nt + 4, side_t);
 
1117
 
 
1118
    nxt = dad[t];
 
1119
    p = edgeToSeg(rtr->tg, nxt, t);
 
1120
    side1[cnt1].ts = addTri(-1, p.j, NULL);
 
1121
    side1[cnt1++].v = p.i;
 
1122
    side2[cnt2].ts = addTri(-1, p.i, NULL);
 
1123
    side2[cnt2++].v = p.j;
 
1124
 
 
1125
    t = nxt;
 
1126
    for (nxt = dad[t]; nxt >= 0; nxt = dad[nxt]) {
 
1127
        p = edgeToSeg(rtr->tg, t, nxt);
 
1128
        if (p.i == side1[cnt1 - 1].v) {
 
1129
            side1[cnt1 - 1].ts =
 
1130
                addTri(side2[cnt2 - 1].v, p.j, side1[cnt1 - 1].ts);
 
1131
            side2[cnt2 - 1].ts =
 
1132
                addTri(side1[cnt1 - 1].v, p.j, side2[cnt2 - 1].ts);
 
1133
            side2[cnt2].ts =
 
1134
                addTri(side2[cnt2 - 1].v, side1[cnt1 - 1].v, NULL);
 
1135
            side2[cnt2++].v = p.j;
 
1136
        } else if (p.i == side2[cnt2 - 1].v) {
 
1137
            side1[cnt1 - 1].ts =
 
1138
                addTri(side2[cnt2 - 1].v, p.j, side1[cnt1 - 1].ts);
 
1139
            side2[cnt2 - 1].ts =
 
1140
                addTri(side1[cnt1 - 1].v, p.j, side2[cnt2 - 1].ts);
 
1141
            side1[cnt1].ts =
 
1142
                addTri(side2[cnt2 - 1].v, side1[cnt1 - 1].v, NULL);
 
1143
            side1[cnt1++].v = p.j;
 
1144
        } else if (p.j == side1[cnt1 - 1].v) {
 
1145
            side1[cnt1 - 1].ts =
 
1146
                addTri(side2[cnt2 - 1].v, p.i, side1[cnt1 - 1].ts);
 
1147
            side2[cnt2 - 1].ts =
 
1148
                addTri(side1[cnt1 - 1].v, p.i, side2[cnt2 - 1].ts);
 
1149
            side2[cnt2].ts =
 
1150
                addTri(side2[cnt2 - 1].v, side1[cnt1 - 1].v, NULL);
 
1151
            side2[cnt2++].v = p.i;
 
1152
        } else {
 
1153
            side1[cnt1 - 1].ts =
 
1154
                addTri(side2[cnt2 - 1].v, p.i, side1[cnt1 - 1].ts);
 
1155
            side2[cnt2 - 1].ts =
 
1156
                addTri(side1[cnt1 - 1].v, p.i, side2[cnt2 - 1].ts);
 
1157
            side1[cnt1].ts =
 
1158
                addTri(side2[cnt2 - 1].v, side1[cnt1 - 1].v, NULL);
 
1159
            side1[cnt1++].v = p.i;
 
1160
        }
 
1161
        t = nxt;
 
1162
    }
 
1163
    side1[cnt1 - 1].ts = addTri(-2, side2[cnt2 - 1].v, side1[cnt1 - 1].ts);
 
1164
    side2[cnt2 - 1].ts = addTri(-2, side1[cnt1 - 1].v, side2[cnt2 - 1].ts);
 
1165
 
 
1166
    /* store points in pts starting with t in 0, 
 
1167
     * then side1, then s, then side2 
 
1168
     */
 
1169
    vmap = dtopen(&ipairdisc, Dtoset);
 
1170
    vmapAdd(vmap, -1, 0);
 
1171
    vmapAdd(vmap, -2, cnt1 + 1);
 
1172
    pps = pts = N_GNEW(nt + 4, pointf);
 
1173
    trim = N_NEW(nt + 4, tri *);
 
1174
    *pps++ = p_t;
 
1175
    idx = 1;
 
1176
    for (i = 0; i < cnt1; i++) {
 
1177
        vmapAdd(vmap, side1[i].v, idx);
 
1178
        *pps++ = rtr->ps[side1[i].v];
 
1179
        trim[idx++] = side1[i].ts;
 
1180
    }
 
1181
    *pps++ = p_s;
 
1182
    idx++;
 
1183
    for (i = cnt2 - 1; i >= 0; i--) {
 
1184
        vmapAdd(vmap, side2[i].v, idx);
 
1185
        *pps++ = rtr->ps[side2[i].v];
 
1186
        trim[idx++] = side2[i].ts;
 
1187
    }
 
1188
 
 
1189
    for (i = 0; i < nt + 4; i++) {
 
1190
        mapTri(vmap, trim[i]);
 
1191
    }
 
1192
 
 
1193
    ps = NEW(tripoly_t);
 
1194
    ps->poly.pn = nt + 4;  /* nt triangles gives nt+2 points plus s and t */
 
1195
    ps->poly.ps = pts;
 
1196
    ps->triMap = trim;
 
1197
 
 
1198
    free (side1);
 
1199
    free (side2);
 
1200
    dtclose(vmap);
 
1201
    *sx = cnt1 + 1;             /* index of s in ps */
 
1202
    return ps;
 
1203
}
 
1204
 
 
1205
/* resetGraph:
 
1206
 * Remove edges and nodes added for current edge routing
 
1207
 */
 
1208
static void resetGraph(tgraph * g, int ncnt, int ecnt)
 
1209
{
 
1210
    int i;
 
1211
    tnode *np = g->nodes;
 
1212
    g->nedges = ecnt;
 
1213
    for (i = 0; i < ncnt; i++) {
 
1214
        if (np->edges + np->ne == (np + 1)->edges)
 
1215
            np->ne--;
 
1216
        np++;
 
1217
    }
 
1218
}
 
1219
 
 
1220
#define PQTYPE int
 
1221
#define PQVTYPE float
 
1222
 
 
1223
#define PQ_TYPES
 
1224
#include <fPQ.h>
 
1225
#undef PQ_TYPES
 
1226
 
 
1227
typedef struct {
 
1228
    PQ pq;
 
1229
    PQVTYPE *vals;
 
1230
    int *idxs;
 
1231
} PPQ;
 
1232
 
 
1233
#define N_VAL(pq,n) ((PPQ*)pq)->vals[n]
 
1234
#define N_IDX(pq,n) ((PPQ*)pq)->idxs[n]
 
1235
 
 
1236
#define PQ_CODE
 
1237
#include <fPQ.h>
 
1238
#undef PQ_CODE
 
1239
 
 
1240
#define N_DAD(n) dad[n]
 
1241
#define E_WT(e) (e->dist)
 
1242
#define UNSEEN -MAXFLOAT
 
1243
 
 
1244
/* triPath:
 
1245
 * Find the shortest path with lengths in g from
 
1246
 * v0 to v1. The returned vector (dad) encodes the
 
1247
 * shorted path from v1 to v0. That path is given by
 
1248
 * v1, dad[v1], dad[dad[v1]], ..., v0.
 
1249
 */
 
1250
static int *
 
1251
triPath(tgraph * g, int n, int v0, int v1, PQ * pq)
 
1252
{
 
1253
    int i, j, adjn;
 
1254
    double d;
 
1255
    tnode *np;
 
1256
    tedge *e;
 
1257
    int *dad = N_NEW(n, int);
 
1258
 
 
1259
    for (i = 0; i < pq->PQsize; i++)
 
1260
        N_VAL(pq, i) = UNSEEN;
 
1261
 
 
1262
    PQinit(pq);
 
1263
    N_DAD(v0) = -1;
 
1264
    N_VAL(pq, v0) = 0;
 
1265
    PQinsert(pq, v0);
 
1266
 
 
1267
    while ((i = PQremove(pq)) != -1) {
 
1268
        N_VAL(pq, i) *= -1;
 
1269
        if (i == v1)
 
1270
            break;
 
1271
        np = g->nodes + i;
 
1272
        for (j = 0; j < np->ne; j++) {
 
1273
            e = g->edges + np->edges[j];
 
1274
            if (e->t == i)
 
1275
                adjn = e->h;
 
1276
            else
 
1277
                adjn = e->t;
 
1278
            if (N_VAL(pq, adjn) < 0) {
 
1279
                d = -(N_VAL(pq, i) + E_WT(e));
 
1280
                if (N_VAL(pq, adjn) == UNSEEN) {
 
1281
                    N_VAL(pq, adjn) = d;
 
1282
                    N_DAD(adjn) = i;
 
1283
                    PQinsert(pq, adjn);
 
1284
                } else if (N_VAL(pq, adjn) < d) {
 
1285
                    PQupdate(pq, adjn, d);
 
1286
                    N_DAD(adjn) = i;
 
1287
                }
 
1288
            }
 
1289
        }
 
1290
    }
 
1291
    return dad;
 
1292
}
 
1293
 
 
1294
/* makeMultiSpline:
 
1295
 * FIX: we don't really use the shortest path provided by ED_path,
 
1296
 * so avoid in neato spline code.
 
1297
 * Return 0 on success.
 
1298
 */
 
1299
int makeMultiSpline(graph_t* g,  edge_t* e, router_t * rtr, int doPolyline)
 
1300
{
 
1301
    Ppolyline_t line = ED_path(e);
 
1302
    node_t *t = agtail(e);
 
1303
    node_t *h = aghead(e);
 
1304
    pointf t_p = line.ps[0];
 
1305
    pointf h_p = line.ps[line.pn - 1];
 
1306
    tripoly_t *poly;
 
1307
    int idx;
 
1308
    int *sp;
 
1309
    int t_id = rtr->tn;
 
1310
    int h_id = rtr->tn + 1;
 
1311
    int ecnt = rtr->tg->nedges;
 
1312
    PPQ pq;
 
1313
    PQTYPE *idxs;
 
1314
    PQVTYPE *vals;
 
1315
    int ret;
 
1316
 
 
1317
        /* Add endpoints to triangle graph */
 
1318
    addEndpoint(rtr, t_p, t, t_id, ED_tail_port(e).side);
 
1319
    addEndpoint(rtr, h_p, h, h_id, ED_head_port(e).side);
 
1320
 
 
1321
        /* Initialize priority queue */
 
1322
    PQgen(&pq.pq, rtr->tn + 2, -1);
 
1323
    idxs = N_GNEW(pq.pq.PQsize + 1, PQTYPE);
 
1324
    vals = N_GNEW(pq.pq.PQsize + 1, PQVTYPE);
 
1325
    vals[0] = 0;
 
1326
    pq.vals = vals + 1;
 
1327
    pq.idxs = idxs + 1;
 
1328
 
 
1329
        /* Find shortest path of triangles */
 
1330
    sp = triPath(rtr->tg, rtr->tn+2, h_id, t_id, (PQ *) & pq);
 
1331
 
 
1332
    free(vals);
 
1333
    free(idxs);
 
1334
    PQfree(&(pq.pq), 0);
 
1335
 
 
1336
        /* Use path of triangles to generate guiding polygon */
 
1337
    poly = mkPoly(rtr, sp, h_id, t_id, h_p, t_p, &idx);
 
1338
    free(sp);
 
1339
 
 
1340
        /* Generate multiple splines using polygon */
 
1341
    ret = genroute(g, poly, 0, idx, e, doPolyline);
 
1342
    freeTripoly (poly);
 
1343
 
 
1344
    resetGraph(rtr->tg, rtr->tn, ecnt);
 
1345
    return ret;
 
1346
}