From: erg Date: Fri, 31 Oct 2008 20:29:18 +0000 (+0000) Subject: Commit initial version of orthogonal edge code X-Git-Tag: LAST_LIBGRAPH~32^2~2914 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=27f8c7f49e2413fab629dd089725599f679a0922;p=graphviz Commit initial version of orthogonal edge code --- diff --git a/configure.ac b/configure.ac index eb91f5433..2e3c4a66b 100644 --- a/configure.ac +++ b/configure.ac @@ -2928,6 +2928,7 @@ AC_CONFIG_FILES(Makefile lib/twopigen/Makefile lib/patchwork/Makefile lib/pack/Makefile + lib/ortho/Makefile lib/expr/Makefile lib/expr/libexpr.pc lib/common/Makefile diff --git a/debian/rules b/debian/rules index 6d6d30ccd..ed9521e6d 100644 --- a/debian/rules +++ b/debian/rules @@ -64,6 +64,7 @@ configure-stamp: --with-rsvg \ --with-devil \ --with-gts \ + --without-ortho \ --without-sfdp \ --without-smyrna \ --disable-sharp \ diff --git a/graphviz.spec.in b/graphviz.spec.in index 67958df29..b92dfd95b 100644 --- a/graphviz.spec.in +++ b/graphviz.spec.in @@ -696,6 +696,7 @@ CFLAGS="$RPM_OPT_FLAGS" \ --with%{!?MING:out}-ming \ --with%{!?PANGOCAIRO:out}-pangocairo \ --with%{!?RSVG:out}-rsvg \ + --with%{!?ORTHO:out}-ortho \ --with%{!?SFDP:out}-sfdp \ --with%{!?SMYRNA:out}-smyrna make %{?_smp_mflags} diff --git a/lib/Makefile.am b/lib/Makefile.am index 7d7ff8411..45f1e4425 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -2,7 +2,7 @@ ## Process this file with automake to produce Makefile.in SUBDIRS = cdt graph cgraph gd pathplan sfio vmalloc ast \ - vpsc rbtree sparse patchwork expr common \ + vpsc rbtree ortho sparse patchwork expr common \ pack gvc xdot ingraphs topfish glcomp \ circogen dotgen fdpgen neatogen twopigen sfdpgen diff --git a/lib/dotgen/Makefile.am b/lib/dotgen/Makefile.am index 33c11a50c..c243e80a3 100644 --- a/lib/dotgen/Makefile.am +++ b/lib/dotgen/Makefile.am @@ -11,6 +11,7 @@ AM_CPPFLAGS = \ -I$(top_srcdir) \ -I$(top_srcdir)/lib/common \ -I$(top_srcdir)/lib/gvc \ + -I$(top_srcdir)/lib/ortho \ -I$(top_srcdir)/lib/$(GRAPH) \ -I$(top_srcdir)/lib/cdt \ -I$(top_srcdir)/lib/pathplan diff --git a/lib/gvc/Makefile.am b/lib/gvc/Makefile.am index 89f504254..0188d9da3 100644 --- a/lib/gvc/Makefile.am +++ b/lib/gvc/Makefile.am @@ -55,6 +55,10 @@ libgvc_C_la_LIBADD = \ libgvc_C_la_DEPENDENCIES = \ $(top_builddir)/lib/pack/libpack_C.la \ $(top_builddir)/lib/common/libcommon_C.la +if WITH_ORTHO +libgvc_C_la_LIBADD += $(top_builddir)/lib/ortho/libortho_C.la +libgvc_C_la_DEPENDENCIES += $(top_builddir)/lib/ortho/libortho_C.la +endif #For use with plugins. # so it is linked with an empty table of builtins. libgvc_la_LDFLAGS = -version-info $(VERSION_INFO) -no-undefined diff --git a/lib/neatogen/Makefile.am b/lib/neatogen/Makefile.am index e989ff876..b07c68ba0 100644 --- a/lib/neatogen/Makefile.am +++ b/lib/neatogen/Makefile.am @@ -12,6 +12,7 @@ AM_CPPFLAGS = \ -I$(top_srcdir)/lib/common \ -I$(top_srcdir)/lib/gvc \ -I$(top_srcdir)/lib/pack \ + -I$(top_srcdir)/lib/ortho \ -I$(top_srcdir)/lib/pathplan \ -I$(top_srcdir)/lib/$(GRAPH) \ -I$(top_srcdir)/lib/sfdpgen \ diff --git a/lib/ortho/Makefile.am b/lib/ortho/Makefile.am new file mode 100644 index 000000000..ba15bab0e --- /dev/null +++ b/lib/ortho/Makefile.am @@ -0,0 +1,27 @@ +# $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/gvc \ + -I$(top_srcdir)/lib/neatogen \ + -I$(top_srcdir)/lib/pack \ + -I$(top_srcdir)/lib/pathplan \ + -I$(top_srcdir)/lib/graph \ + -I$(top_srcdir)/lib/cdt + +if WITH_ORTHO +AM_CFLAGS = -D_BLD_common=1 +endif + +noinst_HEADERS = fPQ.h intset.h maze.h partition.h rawgraph.h sgraph.h structures.h trap.h + +if WITH_ORTHO +noinst_LTLIBRARIES = libortho_C.la +endif + +libortho_C_la_SOURCES = fPQ.c intset.c maze.c ortho.c partition.c rawgraph.c sgraph.c trapezoid.c + +EXTRA_DIST = Makefile.old + diff --git a/lib/ortho/Makefile.old b/lib/ortho/Makefile.old new file mode 100644 index 000000000..be8064beb --- /dev/null +++ b/lib/ortho/Makefile.old @@ -0,0 +1,29 @@ +all: libortho.a +ROOT=../.. +include $(ROOT)/Config.mk +include $(ROOT)/makearch/$(ARCH) + +INCS = -I. -I.. + +OBJS = partition.o trapezoid.o main.o ortho.o maze.o sgraph.o fPQ.o \ + rawgraph.o intset.o + +libortho.a: $(OBJS) + $(RM) libortho.a + $(AR) cr libortho.a $(OBJS) + $(RANLIB) libortho.a + +install: libortho.a + $(MKPATH) $(LIBDIR) + $(INSTALL) libortho.a $(LIBDIR) +# $(SHLIB_LD) -o $(LIBDIR)/libortho.so.$(VERSION) $(OBJS) + $(MKPATH) $(INCDIR) + $(INSTALL) pathgeom.h $(INCDIR) + $(MKPATH) $(LIBMANDIR) + $(INSTALL) ortho.3 $(LIBMANDIR) + +clean: + $(RM) *.o core + +distclean: clean + $(RM) *.a lib*.so.* diff --git a/lib/ortho/fPQ.c b/lib/ortho/fPQ.c new file mode 100644 index 000000000..dced61bb4 --- /dev/null +++ b/lib/ortho/fPQ.c @@ -0,0 +1,149 @@ +/* vim:set shiftwidth=4 ts=8: */ +/* Priority Queue Code for shortest path in graph */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +static snode** pq; +static int PQcnt; +static snode guard; +static int PQsize; + +void +PQgen(int sz) +{ + if (!pq) { + pq = N_NEW(sz+1,snode*); + pq[0] = &guard; + PQsize = sz; + } + PQcnt = 0; +} + +void +PQfree() +{ + free (pq); + pq = NULL; + PQcnt = 0; +} + +void +PQinit() +{ + PQcnt = 0; +} + +void +PQcheck () +{ + int i; + + for (i = 1; i <= PQcnt; i++) { + if (N_IDX(pq[i]) != i) { + assert (0); + } + } +} + +void +PQupheap(int k) +{ + snode* x = pq[k]; + int v = x->n_val; + int next = k/2; + snode* n; + + while (N_VAL(n = pq[next]) < v) { + pq[k] = n; + N_IDX(n) = k; + k = next; + next /= 2; + } + pq[k] = x; + N_IDX(x) = k; +} + +void +PQ_insert(snode* np) +{ + if (PQcnt == PQsize) { + fprintf (stderr, "Heap overflow\n"); + exit (1); + } + PQcnt++; + pq[PQcnt] = np; + PQupheap (PQcnt); + PQcheck(); +} + +void +PQdownheap (int k) +{ + snode* x = pq[k]; + int v = N_VAL(x); + int lim = PQcnt/2; + snode* n; + int j; + + while (k <= lim) { + j = k+k; + n = pq[j]; + if (j < PQcnt) { + if (N_VAL(n) < N_VAL(pq[j+1])) { + j++; + n = pq[j]; + } + } + if (v >= N_VAL(n)) break; + pq[k] = n; + N_IDX(n) = k; + k = j; + } + pq[k] = x; + N_IDX(x) = k; +} + +snode* +PQremove () +{ + snode* n; + + if (PQcnt) { + n = pq[1]; + pq[1] = pq[PQcnt]; + PQcnt--; + if (PQcnt) PQdownheap (1); + PQcheck(); + return n; + } + else return 0; +} + +void +PQupdate (snode* n, int d) +{ + N_VAL(n) = d; + PQupheap (n->n_idx); + PQcheck(); +} + +void +PQprint () +{ + int i; + snode* n; + + fprintf (stderr, "Q: "); + for (i = 1; i <= PQcnt; i++) { + n = pq[i]; + fprintf (stderr, "%d(%d:%d) ", + n->index, N_IDX(n), N_VAL(n)); + } + fprintf (stderr, "\n"); +} diff --git a/lib/ortho/fPQ.h b/lib/ortho/fPQ.h new file mode 100644 index 000000000..31c478779 --- /dev/null +++ b/lib/ortho/fPQ.h @@ -0,0 +1,172 @@ +/* vim:set shiftwidth=4 ts=8: */ +/* Priority Queue Code for shortest path in graph */ + +#include +/* typedef snode** PQ; */ + +#define N_VAL(n) (n)->n_val +#define N_IDX(n) (n)->n_idx +#define N_DAD(n) (n)->n_dad +#define N_EDGE(n) (n)->n_edge +#define E_WT(e) (e->weight) +#define E_INCR(e) (e->incr) + +#ifdef INLINEPQ + +#include "assert.h" +static snode** pq; +static int PQcnt; +static snode guard; +static int PQsize; + +static void +PQgen(int sz) +{ + if (!pq) { + pq = N_NEW(sz+1,snode*); + pq[0] = &guard; + PQsize = sz; + } + PQcnt = 0; +} + +static void +PQfree() +{ + free (pq); + pq = NULL; + PQcnt = 0; +} + +static void +PQinit() +{ + PQcnt = 0; +} + +static void +PQcheck () +{ + int i; + + for (i = 1; i <= PQcnt; i++) { + if (N_IDX(pq[i]) != i) { + assert (0); + } + } +} + +static void +PQupheap(int k) +{ + snode* x; + x = pq[k]; + int v = x->n_val; + int next = k/2; + snode* n; + + while (N_VAL(n = pq[next]) < v) { + pq[k] = n; + N_IDX(n) = k; + k = next; + next /= 2; + } + pq[k] = x; + N_IDX(x) = k; +} + +static void +PQ_insert(snode* np) +{ + if (PQcnt == PQsize) { + fprintf (stderr, "Heap overflow\n"); + exit (1); + } + PQcnt++; + pq[PQcnt] = np; + PQupheap (PQcnt); + PQcheck(); +} + +static void +PQdownheap (int k) +{ + snode* x = pq[k]; + int v = N_VAL(x); + int lim = PQcnt/2; + snode* n; + int j; + + while (k <= lim) { + j = k+k; + n = pq[j]; + if (j < PQcnt) { + if (N_VAL(n) < N_VAL(pq[j+1])) { + j++; + n = pq[j]; + } + } + if (v >= N_VAL(n)) break; + pq[k] = n; + N_IDX(n) = k; + k = j; + } + pq[k] = x; + N_IDX(x) = k; +} + +static snode* +PQremove () +{ + snode* n; + + if (PQcnt) { + n = pq[1]; + pq[1] = pq[PQcnt]; + PQcnt--; + if (PQcnt) PQdownheap (1); + PQcheck(); + return n; + } + else return 0; +} + +static void +PQupdate (snode* n, int d) +{ + N_VAL(n) = d; + PQupheap (n->n_idx); + PQcheck(); +} + +static void +PQprint () +{ + int i; + snode* n; + + fprintf (stderr, "Q: "); + for (i = 1; i <= PQcnt; i++) { + n = pq[i]; + fprintf (stderr, "%s(%d:%d) ", + n->index, N_IDX(n), N_VAL(n)); + } + fprintf (stderr, "\n"); +} +#else + +#ifndef FPQ_H +#define FPQ_H + +void PQgen(int sz); +void PQfree(); +void PQinit(); +void PQcheck (); +void PQupheap(int); +void PQ_insert(snode* np); +void PQdownheap (int k); +snode* PQremove (); +void PQupdate (snode* n, int d); +void PQprint (); +#endif +#endif diff --git a/lib/ortho/intset.c b/lib/ortho/intset.c new file mode 100644 index 000000000..14bfa94b9 --- /dev/null +++ b/lib/ortho/intset.c @@ -0,0 +1,49 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +static Void_t* +mkIntItem(Dt_t* d,intitem* obj,Dtdisc_t* disc) +{ + intitem* np = NEW(intitem); + np->id = obj->id; + return (Void_t*)np; +} + +static void +freeIntItem(Dt_t* d,intitem* obj,Dtdisc_t* disc) +{ + free (obj); +} + +static int +cmpid(Dt_t* d, int* key1, int* key2, Dtdisc_t* disc) +{ + if (*key1 > *key2) return 1; + else if (*key1 < *key2) return -1; + else return 0; +} + +static Dtdisc_t intSetDisc = { + offsetof(intitem,id), + sizeof(int), + offsetof(intitem,link), + (Dtmake_f)mkIntItem, + (Dtfree_f)freeIntItem, + (Dtcompar_f)cmpid, + 0, + 0, + 0 +}; + +Dt_t* +openIntSet () +{ + return dtopen(&intSetDisc,Dtoset); +} + diff --git a/lib/ortho/intset.h b/lib/ortho/intset.h new file mode 100644 index 000000000..b9f240309 --- /dev/null +++ b/lib/ortho/intset.h @@ -0,0 +1,13 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifndef INTSET_H +#define INTSET_H + +#include + +typedef struct { + int id; + Dtlink_t link; +} intitem; + +extern Dt_t* openIntSet (); +#endif diff --git a/lib/ortho/maze.c b/lib/ortho/maze.c new file mode 100644 index 000000000..2a65bc9b9 --- /dev/null +++ b/lib/ortho/maze.c @@ -0,0 +1,398 @@ +/* vim:set shiftwidth=4 ts=8: */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include +#include +/* #include */ + +#define MARGIN 36; + +#ifdef DEBUG +char* pre = "%!PS-Adobe-2.0\n\ +/node {\n\ + /Y exch def\n\ + /X exch def\n\ + /y exch def\n\ + /x exch def\n\ + newpath\n\ + x y moveto\n\ + x Y lineto\n\ + X Y lineto\n\ + X y lineto\n\ + closepath fill\n\ +} def\n\ +/cell {\n\ + /Y exch def\n\ + /X exch def\n\ + /y exch def\n\ + /x exch def\n\ + newpath\n\ + x y moveto\n\ + x Y lineto\n\ + X Y lineto\n\ + X y lineto\n\ + closepath stroke\n\ +} def\n"; + +char* post = "showpage\n"; + +static void +psdump (cell* gcells, int n_gcells, boxf BB, boxf* rects, int nrect) +{ + int i; + boxf bb; + box absbb; + + absbb.LL.y = absbb.LL.x = 10; + absbb.UR.x = absbb.LL.x + BB.UR.x - BB.LL.x; + absbb.UR.y = absbb.LL.y + BB.UR.y - BB.LL.y; + fputs (pre, stderr); + fprintf (stderr, "%%%%Page: 1 1\n%%%%PageBoundingBox: %d %d %d %d\n", + absbb.LL.x, absbb.LL.y, absbb.UR.x, absbb.UR.y); + + + fprintf (stderr, "%f %f translate\n", 10-BB.LL.x, 10-BB.LL.y); + fputs ("0 0 1 setrgbcolor\n", stderr); + for (i = 0; i < n_gcells; i++) { + bb = gcells[i].bb; + fprintf (stderr, "%f %f %f %f node\n", bb.LL.x, bb.LL.y, bb.UR.x, bb.UR.y); + } + fputs ("0 0 0 setrgbcolor\n", stderr); + for (i = 0; i < nrect; i++) { + bb = rects[i]; + fprintf (stderr, "%f %f %f %f cell\n", bb.LL.x, bb.LL.y, bb.UR.x, bb.UR.y); + } + fputs ("1 0 0 setrgbcolor\n", stderr); + fprintf (stderr, "%f %f %f %f cell\n", BB.LL.x, BB.LL.y, BB.UR.x, BB.UR.y); + fputs (post, stderr); +} +#endif + +static int +vcmpid(Dt_t* d, pointf* key1, pointf* key2, Dtdisc_t* disc) +{ + if (key1->x > key2->x) return 1; + else if (key1->x < key2->x) return -1; + else if (key1->y > key2->y) return 1; + else if (key1->y < key2->y) return -1; + else return 0; +} + +static int +hcmpid(Dt_t* d, pointf* key1, pointf* key2, Dtdisc_t* disc) +{ + if (key1->y > key2->y) return 1; + else if (key1->y < key2->y) return -1; + else if (key1->x > key2->x) return 1; + else if (key1->x < key2->x) return -1; + else return 0; +} + +typedef struct { + snode* np; + pointf p; + Dtlink_t link; +} snodeitem; + +static Dtdisc_t vdictDisc = { + offsetof(snodeitem,p), + sizeof(pointf), + offsetof(snodeitem,link), + 0, + 0, + (Dtcompar_f)vcmpid, + 0, + 0, + 0 +}; +static Dtdisc_t hdictDisc = { + offsetof(snodeitem,p), + sizeof(pointf), + offsetof(snodeitem,link), + 0, + 0, + (Dtcompar_f)hcmpid, + 0, + 0, + 0 +}; + +#define delta 1 /* weight of length */ +#define mu 500 /* weight of bends */ + +#define BEND(g,e) ((g->nodes + e->v1)->isVert != (g->nodes + e->v2)->isVert) +#define HORZ(g,e) ((g->nodes + e->v1)->isVert) +#define BIG 16384 +#define CHANSZ(w) (((w)-3)/2) + +static void +updateWt (cell* cp, sedge* ep, int sz) +{ + ep->cnt++; + if (ep->cnt > sz) { + ep->cnt = 0; + ep->weight += BIG; + } +} + +/* updateWts: + * Iterate over edges in a cell, adjust weights as necessary. + * It always updates the bent edges belonging to a cell. + * A horizontal/vertical edge is updated only if the edge traversed + * is bent, or if it is the traversed edge. + */ +void +updateWts (sgraph* g, cell* cp, sedge* ep) +{ + int i; + sedge* e; + int isBend = BEND(g,ep); + int hsz = CHANSZ (cp->bb.UR.y - cp->bb.LL.y); + int vsz = CHANSZ (cp->bb.UR.x - cp->bb.LL.x); + int minsz = MIN(hsz, vsz); + + /* Bend edges are added first */ + for (i = 0; i < cp->nedges; i++) { + e = cp->edges[i]; + if (!BEND(g,e)) break; + updateWt (cp, e, minsz); + } + + for (; i < cp->nedges; i++) { + e = cp->edges[i]; + if (isBend || (e == ep)) updateWt (cp, e, (HORZ(g,e)?hsz:vsz)); + } +} + +static void +createSEdges (cell* cp, sgraph* g) +{ + boxf bb = cp->bb; + double hwt = delta*(bb.UR.x-bb.LL.x); + double vwt = delta*(bb.UR.y-bb.LL.y); + double wt = (hwt + vwt)/2.0 + mu; + + if (CHANSZ(bb.UR.y-bb.LL.y) < 2) { + hwt = BIG; + wt = BIG; + } + if (CHANSZ(bb.UR.x-bb.LL.x) < 2) { + vwt = BIG; + wt = BIG; + } + + if (cp->sides[M_LEFT] && cp->sides[M_TOP]) + cp->edges[cp->nedges++] = createSEdge (g, cp->sides[M_LEFT], cp->sides[M_TOP], wt); + if (cp->sides[M_TOP] && cp->sides[M_RIGHT]) + cp->edges[cp->nedges++] = createSEdge (g, cp->sides[M_TOP], cp->sides[M_RIGHT], wt); + if (cp->sides[M_LEFT] && cp->sides[M_BOTTOM]) + cp->edges[cp->nedges++] = createSEdge (g, cp->sides[M_LEFT], cp->sides[M_BOTTOM], wt); + if (cp->sides[M_BOTTOM] && cp->sides[M_RIGHT]) + cp->edges[cp->nedges++] = createSEdge (g, cp->sides[M_BOTTOM], cp->sides[M_RIGHT], wt); + if (cp->sides[M_TOP] && cp->sides[M_BOTTOM]) + cp->edges[cp->nedges++] = createSEdge (g, cp->sides[M_TOP], cp->sides[M_BOTTOM], vwt); + if (cp->sides[M_LEFT] && cp->sides[M_RIGHT]) + cp->edges[cp->nedges++] = createSEdge (g, cp->sides[M_LEFT], cp->sides[M_RIGHT], hwt); +} + +static snode* +findSVert (sgraph* g, Dt_t* cdt, pointf p, snodeitem* ditems, boolean isVert) +{ + snodeitem* n = dtmatch (cdt, &p); + + if (!n) { + snode* np = createSNode (g); + assert(ditems); + n = ditems + np->index; + n->p = p; + n->np = np; + np->isVert = isVert; + dtinsert (cdt, n); + } + + return n->np; +} + +static sgraph* +mkMazeGraph (maze* mp, boxf bb) +{ + int nsides, i, ncnt, maxdeg; + int bound = 4*mp->ncells; + sgraph* g = createSGraph (bound + 2); + Dt_t* vdict = dtopen(&vdictDisc,Dtoset); + Dt_t* hdict = dtopen(&hdictDisc,Dtoset); + snodeitem* ditems = N_NEW(bound, snodeitem); + snode** sides; + + /* For each cell, create if necessary and attach a node in search + * corresponding to each internal face. The node also gets + * a pointer to the cell. + */ + sides = N_NEW(4*mp->ncells, snode*); + ncnt = 0; + for (i = 0; i < mp->ncells; i++) { + cell* cp = mp->cells+i; + snode* np; + pointf pt; + + cp->nsides = 4; + cp->sides = sides + 4*i; + if (cp->bb.UR.x < bb.UR.x) { + pt.x = cp->bb.UR.x; + pt.y = cp->bb.LL.y; + np = findSVert (g, vdict, pt, ditems, TRUE); + np->cells[0] = cp; + cp->sides[M_RIGHT] = np; + } + if (cp->bb.UR.y < bb.UR.y) { + pt.x = cp->bb.LL.x; + pt.y = cp->bb.UR.y; + np = findSVert (g, hdict, pt, ditems, FALSE); + np->cells[0] = cp; + cp->sides[M_TOP] = np; + } + if (cp->bb.LL.x > bb.LL.x) { + np = findSVert (g, vdict, cp->bb.LL, ditems, TRUE); + np->cells[1] = cp; + cp->sides[M_LEFT] = np; + } + if (cp->bb.LL.y > bb.LL.y) { + np = findSVert (g, hdict, cp->bb.LL, ditems, FALSE); + np->cells[1] = cp; + cp->sides[M_BOTTOM] = np; + } + } + + /* For each gcell, corresponding to a node in the input graph, + * connect it to its corresponding search nodes. + */ + maxdeg = 0; + sides = N_NEW(g->nnodes, snode*); + nsides = 0; + for (i = 0; i < mp->ngcells; i++) { + cell* cp = mp->gcells+i; + pointf pt; + cp->sides = sides+nsides; + pt = cp->bb.LL; + snodeitem* np = dtmatch (hdict, &pt); + for (; np && np->p.x < cp->bb.UR.x; np = dtnext (hdict, np)) { + cp->sides[cp->nsides++] = np->np; + np->np->cells[1] = cp; + } + np = dtmatch (vdict, &pt); + for (; np && np->p.y < cp->bb.UR.y; np = dtnext (vdict, np)) { + cp->sides[cp->nsides++] = np->np; + np->np->cells[1] = cp; + } + pt.y = cp->bb.UR.y; + np = dtmatch (hdict, &pt); + for (; np && np->p.x < cp->bb.UR.x; np = dtnext (hdict, np)) { + cp->sides[cp->nsides++] = np->np; + np->np->cells[0] = cp; + } + pt.x = cp->bb.UR.x; + pt.y = cp->bb.LL.y; + np = dtmatch (vdict, &pt); + for (; np && np->p.y < cp->bb.UR.y; np = dtnext (vdict, np)) { + cp->sides[cp->nsides++] = np->np; + np->np->cells[0] = cp; + } + nsides += cp->nsides; + if (cp->nsides > maxdeg) maxdeg = cp->nsides; + } + RALLOC (nsides, sides, snode*); + + /* Set index of two dummy nodes used for real nodes */ + g->nodes[g->nnodes].index = g->nnodes; + g->nodes[g->nnodes+1].index = g->nnodes+1; + + /* create edges + * For each ordinary cell, there can be at most 6 edges. + * At most 2 gcells will be used at a time, and each of these + * can have at most degree maxdeg. + */ + initSEdges (g, maxdeg); + for (i = 0; i < mp->ncells; i++) { + cell* cp = mp->cells+i; + createSEdges (cp, g); + } + + /* tidy up memory */ + g->nodes = RALLOC (g->nnodes+2, g->nodes, snode); + g->edges = RALLOC (g->nedges+2*maxdeg, g->edges, sedge); + dtclose (vdict); + dtclose (hdict); + free (ditems); + + /* save core graph state */ + gsave(g); + return g; +} + +maze* +mkMaze (graph_t* g) +{ + node_t* n; + maze* mp = NEW(maze); + boxf* rects; + int i, nrect; + cell* cp; + double w2, h2; + boxf bb, BB; + + mp->ngcells = agnnodes(g); + cp = mp->gcells = N_NEW(mp->ngcells, cell); + + BB.LL.x = BB.LL.y = MAXDOUBLE; + BB.UR.x = BB.UR.y = -MAXDOUBLE; + for (n = agfstnode (g); n; n = agnxtnode(g,n)) { + w2 = 36.0*ND_width(n); + h2 = 36.0*ND_height(n); + bb.LL.x = ND_coord(n).x - w2; + bb.UR.x = ND_coord(n).x + w2; + bb.LL.y = ND_coord(n).y - h2; + bb.UR.y = ND_coord(n).y + h2; + BB.LL.x = MIN(BB.LL.x, bb.LL.x); + BB.LL.y = MIN(BB.LL.y, bb.LL.y); + BB.UR.x = MAX(BB.UR.x, bb.UR.x); + BB.UR.y = MAX(BB.UR.y, bb.UR.y); + cp->bb = bb; + cp->flags |= MZ_ISNODE; + ND_alg(n) = cp; + cp++; + } + + BB.LL.x -= MARGIN; + BB.LL.y -= MARGIN; + BB.UR.x += MARGIN; + BB.UR.y += MARGIN; + rects = partition (mp->gcells, mp->ngcells, &nrect, BB); + +#ifdef DEBUG + psdump (mp->gcells, mp->ngcells, BB, rects, nrect); +#endif + mp->cells = N_NEW(nrect, cell); + mp->ncells = nrect; + for (i = 0; i < nrect; i++) { + mp->cells[i].bb = rects[i]; + } + free (rects); + + mp->sg = mkMazeGraph (mp, BB); + return mp; +} + +void freeMaze (maze* mp) +{ + free (mp->cells[0].sides); + free (mp->gcells[0].sides); + free (mp->cells); + free (mp->gcells); +} + diff --git a/lib/ortho/maze.h b/lib/ortho/maze.h new file mode 100644 index 000000000..3415f705c --- /dev/null +++ b/lib/ortho/maze.h @@ -0,0 +1,40 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifndef MAZE_H +#define MAZE_H + +#include + +enum {M_RIGHT=0, M_TOP, M_LEFT, M_BOTTOM}; + +#define MZ_ISNODE 1 +#define MZ_VSCAN 2 +#define MZ_HSCAN 4 + +#define IsNode(cp) (cp->flags & MZ_ISNODE) /* cell corresponds to node */ + /* cell already inserted in vertical channel */ +#define IsVScan(cp) (cp->flags & MZ_VSCAN) + /* cell already inserted in horizontal channel */ +#define IsHScan(cp) (cp->flags & MZ_HSCAN) + +typedef struct cell { + int flags; + int nedges; + sedge* edges[6]; + int nsides; + snode** sides; + boxf bb; +} cell; + +typedef struct { + int ncells, ngcells; + cell* cells; /* cells not corresponding to graph nodes */ + cell* gcells; /* cells corresponding to graph nodes */ + sgraph* sg; + Dt_t* hchans; + Dt_t* vchans; +} maze; + +extern maze* mkMaze (graph_t*); +extern void freeMaze (maze*); +void updateWts (sgraph* g, cell* cp, sedge* ep); +#endif diff --git a/lib/ortho/ortho.c b/lib/ortho/ortho.c new file mode 100644 index 000000000..9202b44e1 --- /dev/null +++ b/lib/ortho/ortho.c @@ -0,0 +1,1422 @@ +/* vim:set shiftwidth=4 ts=8: */ + +/* TODO: + * Increase weight as edges in sgraph are used, so that alternate + * paths will be found. + * In dot, prefer bottom or top routing + * Ports/compass points + * Edge labels + * Loops + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include "fPQ.h" +#include +#include +#include +#include + +#define DEBUG + +#define CELL(n) ((cell*)ND_alg(n)) +#define MID(a,b) (((a)+(b))/2.0) +#define SC 1 + +/* cellOf: + * Given 2 snodes sharing a cell, return the cell. + */ +static cell* +cellOf (snode* p, snode* q) +{ + cell* cp = p->cells[0]; + if ((cp == q->cells[0]) || (cp == q->cells[1])) return cp; + else return p->cells[1]; +} + +static pointf +midPt (cell* cp) +{ + pointf p; + p.x = MID(cp->bb.LL.x,cp->bb.UR.x); + p.y = MID(cp->bb.LL.y,cp->bb.UR.y); + return p; +} + +/* sidePt: + * Given a cell and an snode on one of its sides, return the + * midpoint of the side. + */ +static pointf +sidePt (snode* ptr, cell* cp) +{ + pointf pt; + if (cp == ptr->cells[1]) { + if (ptr->isVert) { + pt.x = cp->bb.LL.x; + pt.y = MID(cp->bb.LL.y,cp->bb.UR.y); + } + else { + pt.x = MID(cp->bb.LL.x,cp->bb.UR.x); + pt.y = cp->bb.LL.y; + } + } + else { + if (ptr->isVert) { + pt.x = cp->bb.UR.x; + pt.y = MID(cp->bb.LL.y,cp->bb.UR.y); + } + else { + pt.x = MID(cp->bb.LL.x,cp->bb.UR.x); + pt.y = cp->bb.UR.y; + } + } + return pt; +} + +/* setSet: + * Initialize and normalize segments. + * p1 stores smaller value + * Assume b1 != b2 + */ +static void +setSeg (segment* sp, int dir, double fix, double b1, double b2, int l1, int l2) +{ + sp->isVert = dir; + sp->comm_coord = fix; + if (b1 < b2) { + sp->p.p1 = b1; + sp->p.p2 = b2; + sp->l1 = l1; + sp->l2 = l2; + sp->flipped = 0; + } + else { + sp->p.p2 = b1; + sp->p.p1 = b2; + sp->l2 = l1; + sp->l1 = l2; + sp->flipped = 1; + } +} + +/* Convert route in shortest path graph to route + * of segments. This records the first and last cells, + * plus cells where the path bends. + * Note that the shortest path will always have at least 4 nodes: + * the two dummy nodes representing the center of the two real nodes, + * and the two nodes on the boundary of the two real nodes. + */ +#define PUSH(rte,P) (rte.p[rte.n++] = P) + +static route +convertSPtoRoute (sgraph* g, snode* fst, snode* lst) +{ + route rte; + snode* ptr; + snode* next; + snode* prev; /* node in shortest path just previous to next */ + int i, sz = 0; + cell* cp; + cell* ncp; + segment seg; + double fix, b1, b2; + int l1, l2; + pointf bp1, bp2, prevbp; /* bend points */ + + /* count no. of nodes in shortest path */ + for (ptr = fst; ptr; ptr = N_DAD(ptr)) sz++; + rte.n = 0; + rte.segs = N_NEW(sz-2, segment); /* at most sz-2 segments */ + + seg.prev = seg.next = 0; + ptr = prev = N_DAD(fst); + next = N_DAD(ptr); + if (IsNode(ptr->cells[0])) + cp = ptr->cells[1]; + else + cp = ptr->cells[0]; + bp1 = sidePt (ptr, cp); + while (N_DAD(next)!=NULL) { + ncp = cellOf (prev, next); + updateWts (g, ncp, N_EDGE(ptr)); + + /* add seg if path bends or at end */ + if ((ptr->isVert != next->isVert) || (N_DAD(next) == lst)) { + if (ptr->isVert != next->isVert) + bp2 = midPt (ncp); + else + bp2 = sidePt(next, ncp); + if (ptr->isVert) { /* horizontal segment */ + if (ptr == N_DAD(fst)) l1 = B_NODE; + else if (prevbp.y > bp1.y) l1 = B_UP; + else l1 = B_DOWN; + if (ptr->isVert != next->isVert) { + if (next->cells[0] == ncp) l2 = B_UP; + else l2 = B_DOWN; + } + else l2 = B_NODE; + fix = cp->bb.LL.y; + b1 = cp->bb.LL.x; + b2 = ncp->bb.LL.x; + } + else { /* vertical segment */ + if (ptr == N_DAD(fst)) l1 = B_NODE; + else if (prevbp.x > bp1.x) l1 = B_RIGHT; + else l1 = B_LEFT; + if (ptr->isVert != next->isVert) { + if (next->cells[0] == ncp) l2 = B_RIGHT; + else l2 = B_LEFT; + } + else l2 = B_NODE; + fix = cp->bb.LL.x; + b1 = cp->bb.LL.y; + b2 = ncp->bb.LL.y; + } + setSeg (&seg, !ptr->isVert, fix, b1, b2, l1, l2); + rte.segs[rte.n++] = seg; + cp = ncp; + prevbp = bp1; + bp1 = bp2; + if ((ptr->isVert != next->isVert) && (N_DAD(next) == lst)) { + bp2 = sidePt(next, ncp); + l2 = B_NODE; + if (next->isVert) { /* horizontal segment */ + if (prevbp.y > bp1.y) l1 = B_UP; + else l1 = B_DOWN; + fix = cp->bb.LL.y; + b1 = cp->bb.LL.x; + b2 = ncp->bb.LL.x; + } + else { + if (prevbp.x > bp1.x) l1 = B_RIGHT; + else l1 = B_LEFT; + fix = cp->bb.LL.x; + b1 = cp->bb.LL.y; + b2 = ncp->bb.LL.y; + } + setSeg (&seg, !next->isVert, fix, b1, b2, l1, l2); + rte.segs[rte.n++] = seg; + } + ptr = next; + } + prev = next; + next = N_DAD(next); + } + + rte.segs = realloc (rte.segs, rte.n*sizeof(segment)); + for (i=0; i 0) + rte.segs[i].prev = rte.segs + (i-1); + if (i < rte.n-1) + rte.segs[i].next = rte.segs + (i+1); + } + + return rte; +} + +typedef struct { + Dtlink_t link; + double v; + Dt_t* chans; +} chanItem; + +static void +freeChannel (Dt_t* d, channel* cp, Dtdisc_t* disc) +{ + free_graph (cp->G); + free (cp); +} + +static void +freeChanItem (Dt_t* d, chanItem* cp, Dtdisc_t* disc) +{ + dtclose (cp->chans); + free (cp); +} + +/* chancmpid: + * Compare intervals. Two intervals are equal if one contains + * the other. Otherwise, the one with the smaller p1 value is + * less. + * This combines two separate functions into one. Channels are + * disjoint, so we really only need to key on p1. + * When searching for a channel containing a segment, we rely on + * interval containment to return the correct channel. + */ +static int +chancmpid(Dt_t* d, paird* key1, paird* key2, Dtdisc_t* disc) +{ + if (key1->p1 > key2->p1) { + if (key1->p2 <= key2->p2) return 0; + else return 1; + } + else if (key1->p1 < key2->p1) { + if (key1->p2 >= key2->p2) return 0; + else return -1; + } + else return 0; +} + +static int +dcmpid(Dt_t* d, double* key1, double* key2, Dtdisc_t* disc) +{ + if (*key1 > *key2) return 1; + else if (*key1 < *key2) return -1; + else return 0; +} + +static Dtdisc_t chanDisc = { + offsetof(channel,p), + sizeof(paird), + offsetof(channel,link), + 0, + (Dtfree_f)freeChannel, + (Dtcompar_f)chancmpid, + 0, + 0, + 0 +}; + +static Dtdisc_t chanItemDisc = { + offsetof(chanItem,v), + sizeof(double), + offsetof(chanItem,link), + 0, + (Dtfree_f)freeChanItem, + (Dtcompar_f)dcmpid, + 0, + 0, + 0 +}; + +static void +addChan (Dt_t* chdict, channel* cp, double j) +{ + chanItem* subd = dtmatch (chdict, &j); + + if (!subd) { + subd = NEW (chanItem); + subd->v = j; + subd->chans = dtopen (&chanDisc, Dtoset); + dtinsert (chdict, subd); + } + dtinsert (subd->chans, cp); +} + +static Dt_t* +extractHChans (maze* mp) +{ + int i; + snode* np; + Dt_t* hchans = dtopen (&chanItemDisc, Dtoset); + + for (i = 0; i < mp->ncells; i++) { + channel* chp; + cell* cp = mp->cells+i; + cell* nextcp; + if (IsHScan(cp)) continue; + + /* move left */ + while ((np = cp->sides[M_LEFT]) && (nextcp = np->cells[0]) && + !IsNode(nextcp)) { + cp = nextcp; + } + + chp = NEW(channel); + chp->cp = cp; + chp->p.p1 = cp->bb.LL.x; + + /* move right */ + cp->flags |= MZ_HSCAN; + while ((np = cp->sides[M_RIGHT]) && (nextcp = np->cells[1]) && + !IsNode(nextcp)) { + cp = nextcp; + cp->flags |= MZ_HSCAN; + } + + chp->p.p2 = cp->bb.UR.x; + addChan (hchans, chp, chp->cp->bb.LL.y); + } + return hchans; +} + +static Dt_t* +extractVChans (maze* mp) +{ + int i; + snode* np; + Dt_t* vchans = dtopen (&chanItemDisc, Dtoset); + + for (i = 0; i < mp->ncells; i++) { + channel* chp; + cell* cp = mp->cells+i; + cell* nextcp; + if (IsVScan(cp)) continue; + + /* move down */ + while ((np = cp->sides[M_BOTTOM]) && (nextcp = np->cells[0]) && + !IsNode(nextcp)) { + cp = nextcp; + } + + chp = NEW(channel); + chp->cp = cp; + chp->p.p1 = cp->bb.LL.y; + + /* move up */ + cp->flags |= MZ_VSCAN; + while ((np = cp->sides[M_TOP]) && (nextcp = np->cells[1]) && + !IsNode(nextcp)) { + cp = nextcp; + cp->flags |= MZ_VSCAN; + } + + chp->p.p2 = cp->bb.UR.y; + addChan (vchans, chp, chp->cp->bb.LL.x); + } + return vchans; +} + +static void +insertChan (channel* chan, segment* seg) +{ + seg->ind_no = chan->cnt++; + chan->seg_list = ALLOC(chan->cnt, chan->seg_list, segment*); + chan->seg_list[chan->cnt-1] = seg; +} + +static channel* +chanSearch (Dt_t* chans, segment* seg) +{ + channel* cp; + chanItem* chani = dtmatch (chans, &seg->comm_coord); + assert (chani); + cp = dtmatch (chani->chans, &seg->p); + assert (cp); + return cp; +} + +static void +assignSegs (int nrtes, route* route_list, maze* mp) +{ + channel* chan; + int i, j; + + for (i=0;iisVert) + chan = chanSearch(mp->vchans, seg); + else + chan = chanSearch(mp->hchans, seg); + insertChan (chan, seg); + } + } +} + +/* addNodeEdges: + * Add temporary node to sgraph corresponding to cell cp, represented + * by np. + */ +static void +addNodeEdges (sgraph* sg, cell* cp, snode* np) +{ + int i; + for (i = 0; i < cp->nsides; i++) { + snode* onp = cp->sides[i]; + cell* ocp; + pointf midp, p; + double wt; + if (onp->cells[0] == cp) + ocp = onp->cells[1]; + else + ocp = onp->cells[0]; + p = sidePt (onp, ocp); + midp = midPt (cp); + wt = abs(p.x - midp.x) + abs(p.y - midp.y); + createSEdge (sg, np, onp, 0); /* FIX weight */ + } + sg->nnodes++; +#ifdef DEBUG + np->cells[0] = np->cells[1] = cp; +#endif +} + +#ifdef DEBUG + +#include +static char* bendToStr (bend b) +{ + char* s; + switch (b) { + case B_NODE : + s = "B_NODE"; + break; + case B_UP : + s = "B_UP"; + break; + case B_LEFT : + s = "B_LEFT"; + break; + case B_DOWN : + s = "B_DOWN"; + break; + case B_RIGHT : + s = "B_RIGHT"; + break; + } + return s; +} + +static void putSeg (FILE* fp, segment* seg) +{ + if (seg->isVert) + fprintf (fp, "((%f,%f),(%f,%f)) %s %s", seg->comm_coord, seg->p.p1, + seg->comm_coord, seg->p.p2, bendToStr (seg->l1), bendToStr (seg->l2)); + else + fprintf (fp, "((%f,%f),(%f,%f)) %s %s", seg->p.p1,seg->comm_coord, + seg->p.p2, seg->comm_coord, bendToStr (seg->l1), bendToStr (seg->l2)); +} + +static void +dumpChanG (channel* cp, int v) +{ + int k; + intitem* ip; + Dt_t* adj; + + if (cp->cnt < 2) return; + fprintf (stderr, "channel %d (%f,%f)\n", v, cp->p.p1, cp->p.p2); + for (k=0;kcnt;k++) { + adj = cp->G->vertices[k].adj_list; + if (dtsize(adj) == 0) continue; + putSeg (stderr, cp->seg_list[k]); + fputs (" ->\n", stderr); + for (ip = (intitem*)dtfirst(adj); ip; ip = (intitem*)dtnext(adj, ip)) { + fputs (" ", stderr); + putSeg (stderr, cp->seg_list[ip->id]); + fputs ("\n", stderr); + } + } +} +#endif + +static void +assignTrackNo (Dt_t* chans) +{ + Dt_t* lp; + Dtlink_t* l1; + Dtlink_t* l2; + channel* cp; + int k; + + for (l1 = dtflatten (chans); l1; l1 = dtlink(chans,l1)) { + lp = ((chanItem*)l1)->chans; + for (l2 = dtflatten (lp); l2; l2 = dtlink(lp,l2)) { + cp = (channel*)l2; + if (cp->cnt) { +#ifdef DEBUG +/* dumpChanG (cp, ((chanItem*)l1)->v); */ +#endif + top_sort (cp->G); + for (k=0;kcnt;k++) + cp->seg_list[k]->track_no = cp->G->vertices[k].topsort_order+1; + } + } + } +} + +static void +create_graphs(Dt_t* chans) +{ + Dt_t* lp; + Dtlink_t* l1; + Dtlink_t* l2; + channel* cp; + + for (l1 = dtflatten (chans); l1; l1 = dtlink(chans,l1)) { + lp = ((chanItem*)l1)->chans; + for (l2 = dtflatten (lp); l2; l2 = dtlink(lp,l2)) { + cp = (channel*)l2; + cp->G = make_graph (cp->cnt); + } + } +} + +static int +eqEndSeg (bend S1l2, bend S2l2, bend T1, bend T2) +{ + if (((S1l2==T2)&&(S2l2=!T2)) + || ((S1l2==B_NODE)&&(S2l2==T1))) + return(0); + else + return(-1); +} + +static int +overlapSeg (segment* S1, segment* S2, bend T1, bend T2) +{ + if(S1->p.p2p.p2) { + if(S1->l2==T1&&S2->l1==T2) return(-1); + else if(S1->l2==T2&&S2->l1==T1) return(1); + else return(0); + } + else if(S1->p.p2==S2->p.p2) { + if(S2->l1==T2) return eqEndSeg (S1->l2, S2->l2, T1, T2); + else return -1*eqEndSeg (S2->l2, S1->l2, T1, T2); + } + else { /* S1->p.p2>S2->p.p2 */ + if(S2->l1==T2&&S2->l2==T2) return(-1); + else if (S2->l1==T1&&S2->l2==T1) return(1); + else return(0); + } +} + +static int +ellSeg (bend S1l1, bend S1l2, bend T) +{ + if (S1l1 == T) { + if (S1l2== T) return -1; + else return 0; + } + else return 1; +} + +static int +segCmp (segment* S1, segment* S2, bend T1, bend T2) +{ + /* no overlap */ + if((S1->p.p2p.p1)||(S1->p.p1>S2->p.p2)) return(0); + /* left endpoint of S2 inside S1 */ + if(S1->p.p1p.p1&&S2->p.p1p.p2) + return overlapSeg (S1, S2, T1, T2); + /* left endpoint of S1 inside S2 */ + else if(S2->p.p1p.p1&&S1->p.p1p.p2) + return -1*overlapSeg (S2, S1, T1, T2); + else if(S1->p.p1==S2->p.p1) { + if(S1->p.p2==S2->p.p2) { + if((S1->l1==S2->l1)&&(S1->l2==S2->l2)) + return(0); + else if (S2->l1==S2->l2) { + if(S2->l1==T1) return(1); + else if(S2->l1==T2) return(-1); + else if ((S1->l1!=T1)&&(S1->l2!=T1)) return (1); + else if ((S1->l1!=T2)&&(S1->l2!=T2)) return (-1); + else return 0; + } + else if ((S2->l1==T1)&&(S2->l2==T2)) { + if ((S1->l1!=T1)&&(S1->l2==T2)) return 1; + else if ((S1->l1==T1)&&(S1->l2!=T2)) return -1; + else return 0; + } + else if ((S2->l2==T1)&&(S2->l1==T2)) { + if ((S1->l2!=T1)&&(S1->l1==T2)) return 1; + else if ((S1->l2==T1)&&(S1->l1!=T2)) return -1; + else return 0; + } + else if ((S2->l1==B_NODE)&&(S2->l2==T1)) { + return ellSeg (S1->l1, S1->l2, T1); + } + else if ((S2->l1==B_NODE)&&(S2->l2==T2)) { + return -1*ellSeg (S1->l1, S1->l2, T2); + } + else if ((S2->l1==T1)&&(S2->l2==B_NODE)) { + return ellSeg (S1->l2, S1->l1, T1); + } + else { /* ((S2->l1==T2)&&(S2->l2==B_NODE)) */ + return -1*ellSeg (S1->l2, S1->l1, T2); + } + } + else if(S1->p.p2p.p2) { + if(S1->l2==T1) + return eqEndSeg (S2->l1, S1->l1, T1, T2); + else + return -1*eqEndSeg (S2->l1, S1->l1, T1, T2); + } + else { /* S1->p.p2>S2->p.p2 */ + if(S2->l2==T2) + return eqEndSeg (S1->l1, S2->l1, T1, T2); + else + return -1*eqEndSeg (S1->l1, S2->l1, T1, T2); + } + } + else if(S1->p.p2==S2->p.p1) { + if(S1->l2==S2->l1) return(0); + else if(S1->l2==T2) return(1); + else return(-1); + } + else { /* S1->p.p1==S2->p.p2 */ + if(S1->l1==S2->l2) return(0); + else if(S1->l1==T2) return(1); + else return(-1); + } + assert(0); + return 0; +} + +/* Function seg_cmp returns + * -1 if S1 HAS TO BE to the right/below S2 to avoid a crossing, + * 0 if a crossing is unavoidable or there is no crossing at all or + * the segments are parallel, + * 1 if S1 HAS TO BE to the left/above S2 to avoid a crossing + * + * Note: This definition means horizontal segments have track numbers + * increasing as y decreases, while vertical segments have track numbers + * increasing as x increases. It would be good to make this consistent, + * with horizontal track numbers increasing with y. This can be done by + * switching B_DOWN and B_UP in the first call to segCmp. At present, + * though, I'm not sure what assumptions are made in handling parallel + * segments, so we leave the code alone for the time being. + */ +static int +seg_cmp(segment* S1, segment* S2) +{ + if(S1->isVert!=S2->isVert||S1->comm_coord!=S2->comm_coord) { + fprintf (stderr, "incomparable segments !! -- Aborting\n"); + exit(1); + } + if(S1->isVert) + return segCmp (S1, S2, B_RIGHT, B_LEFT); + else + return segCmp (S1, S2, B_DOWN, B_UP); +} + +static void +add_edges_in_G(channel* cp) +{ + int x,y; + segment** seg_list = cp->seg_list; + int size = cp->cnt; + rawgraph* G = cp->G; + + for(x=0;x+1chans; + for (l2 = dtflatten (lp); l2; l2 = dtlink(lp,l2)) { + cp = (channel*)l2; + if (cp->cnt) + add_edges_in_G(cp); + } + } +} + +static segment* +next_seg(segment* seg, int dir) +{ + assert(seg); + if (!dir) + return(seg->prev); + else + return(seg->next); +} + +/* propagate_prec propagates the precedence relationship along + * a series of parallel segments on 2 edges + */ +static int +propagate_prec(segment* seg, int prec, int hops, int dir) +{ + int x; + int ans=prec; + segment* next; + segment* current; + + current = seg; + for(x=1;x<=hops;x++) { + next = next_seg(current, dir); + if(!current->isVert) { + if(next->comm_coord==current->p.p1) { + if(current->l1==B_UP) ans *= -1; + } + else { + if(current->l2==B_DOWN) ans *= -1; + } + } + else { + if(next->comm_coord==current->p.p1) { + if(current->l1==B_RIGHT) ans *= -1; + } + else { + if(current->l2==B_LEFT) ans *= -1; + } + } + current = next; + } + return(ans); +} + +static int +is_parallel(segment* s1, segment* s2) +{ + assert (s1->comm_coord==s2->comm_coord); + return ((s1->p.p1==s2->p.p1)&& + (s1->p.p2==s2->p.p2)&& + (s1->l1==s2->l1)&& + (s1->l2==s2->l2)); +} + +/* decide_point returns the number of hops needed in the given directions + * along the 2 edges to get to a deciding point (or NODES) and also puts + * into prec the appropriate dependency (follows same convention as seg_cmp) + */ +static pair +decide_point(segment* si, segment* sj, int dir1, int dir2) +{ + int prec, ans = 0, temp; + pair ret; + segment* np1; + segment* np2; + + while ((np1 = next_seg(si,dir1)) && (np2 = next_seg(sj,dir2)) && + is_parallel(np1, np2)) { + ans++; + si = np1; + sj = np2; + } + if (!np1) + prec = 0; + else if (!np2) + assert(0); /* FIXME */ + else { + temp = seg_cmp(np1, np2); + prec = propagate_prec(np1, temp, ans+1, 1-dir1); + } + + ret.a = ans; + ret.b = prec; + return(ret); +} + +/* sets the edges for a series of parallel segments along two edges starting + * from segment i, segment j. It is assumed that the edge should be from + * segment i to segment j - the dependency is appropriately propagated + */ +static void +set_parallel_edges (segment* seg1, segment* seg2, int dir1, int dir2, int hops, + maze* mp) +{ + int x; + channel* chan; + channel* nchan; + segment* prev1; + segment* prev2; + + if (seg1->isVert) + chan = chanSearch(mp->vchans, seg1); + else + chan = chanSearch(mp->hchans, seg1); + insert_edge(chan->G, seg1->ind_no, seg2->ind_no); + + for (x=1;x<=hops;x++) { + prev1 = next_seg(seg1, dir1); + prev2 = next_seg(seg2, dir2); + if(!seg1->isVert) { + nchan = chanSearch(mp->vchans, prev1); + if(prev1->comm_coord==seg1->p.p1) { + if(seg1->l1==B_UP) { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + else + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + } + else { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + else + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + } + } + else { + if(seg1->l2==B_UP) { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G,prev1->ind_no, prev2->ind_no); + else + insert_edge(nchan->G,prev2->ind_no, prev1->ind_no); + } + else { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + else + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + } + } + } + else { + nchan = chanSearch(mp->hchans, prev1); + if(prev1->comm_coord==seg1->p.p1) { + if(seg1->l1==B_LEFT) { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + else + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + } + else { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + else + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + } + } + else { + if(seg1->l2==B_LEFT) { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + else + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + } + else { + if(edge_exists(chan->G, seg1->ind_no, seg2->ind_no)) + insert_edge(nchan->G, prev1->ind_no, prev2->ind_no); + else + insert_edge(nchan->G, prev2->ind_no, prev1->ind_no); + } + } + } + chan = nchan; + seg1 = prev1; + seg2 = prev2; + } +} + +/* removes the edge between segments after the resolution of a conflict + */ +static void +removeEdge(segment* seg1, segment* seg2, int dir, maze* mp) +{ + segment* ptr1; + segment* ptr2; + channel* chan; + + ptr1 = seg1; + ptr2 = seg2; + while(is_parallel(ptr1, ptr2)) { + ptr1 = next_seg(ptr1, 1); + ptr2 = next_seg(ptr2, dir); + } + if(ptr1->isVert) + chan = chanSearch(mp->vchans, ptr1); + else + chan = chanSearch(mp->hchans, ptr1); + remove_redge (chan->G, ptr1->ind_no, ptr2->ind_no); +} + +static void +addPEdges (channel* cp, maze* mp) +{ + int i,j; + /* dir[1,2] are used to figure out whether we should use prev + * pointers or next pointers -- 0 : decrease, 1 : increase + */ + int dir; + /* number of hops along the route to get to the deciding points */ + pair hops; + /* precedences of the deciding points : same convention as + * seg_cmp function + */ + int prec1, prec2; + pair p; + rawgraph* G = cp->G; + segment** segs = cp->seg_list; + + for(i=0;i+1cnt;i++) { + for(j=i+1;jcnt;j++) { + if (!edge_exists(G,i,j) && !edge_exists(G,j,i)) { + if (is_parallel(segs[i], segs[j])) { + /* get_directions */ + if(segs[i]->prev==0) { + if(segs[j]->prev==0) + dir = 0; + else + dir = 1; + } + else if(segs[j]->prev==0) { + assert (0); /* FIXME */ + } + else { + if(segs[i]->prev->comm_coord==segs[j]->prev->comm_coord) + dir = 0; + else + dir = 1; + } + + p = decide_point(segs[i], segs[j], 0, dir); + hops.a = p.a; + prec1 = p.b; + p = decide_point(segs[i], segs[j], 1, 1-dir); + hops.b = p.a; + prec2 = p.b; + + switch(prec1) { + case -1 : + set_parallel_edges (segs[j], segs[i], dir, 0, hops.a, mp); + set_parallel_edges (segs[j], segs[i], 1-dir, 1, hops.b, mp); + if(prec2==1) + removeEdge (segs[i], segs[j], 1-dir, mp); + break; + case 0 : + switch(prec2) { + case -1: + set_parallel_edges (segs[j], segs[i], dir, 0, hops.a, mp); + set_parallel_edges (segs[j], segs[i], 1-dir, 1, hops.b, mp); + break; + case 0 : + set_parallel_edges (segs[i], segs[j], 0, dir, hops.a, mp); + set_parallel_edges (segs[i], segs[j], 1, 1-dir, hops.b, mp); + break; + case 1: + set_parallel_edges (segs[i], segs[j], 0, dir, hops.a, mp); + set_parallel_edges (segs[i], segs[j], 1, 1-dir, hops.b, mp); + break; + } + break; + case 1 : + set_parallel_edges (segs[i], segs[j], 0, dir, hops.a, mp); + set_parallel_edges (segs[i], segs[j], 1, 1-dir, hops.b, mp); + if(prec2==-1) + removeEdge (segs[i], segs[j], 1-dir, mp); + break; + } + } + } + } + } +} + +static void +add_p_edges (Dt_t* chans, maze* mp) +{ + Dt_t* lp; + Dtlink_t* l1; + Dtlink_t* l2; + + for (l1 = dtflatten (chans); l1; l1 = dtlink(chans,l1)) { + lp = ((chanItem*)l1)->chans; + for (l2 = dtflatten (lp); l2; l2 = dtlink(lp,l2)) { + addPEdges ((channel*)l2, mp); + } + } +} + +static void +assignTracks (int nrtes, route* route_list, maze* mp) +{ + /* Create the graphs for each channel */ + create_graphs(mp->hchans); + create_graphs(mp->vchans); + + /* add edges between non-parallel segments */ + add_np_edges(mp->hchans); + add_np_edges(mp->vchans); + + /* add edges between parallel segments + remove appropriate edges */ + add_p_edges(mp->hchans, mp); + add_p_edges(mp->vchans, mp); + + /* Assign the tracks after a top sort */ + assignTrackNo (mp->hchans); + assignTrackNo (mp->vchans); +} + +#ifdef DEBUG +static void emitGraph (FILE* fp, maze* mp, Agraph_t* g, route* route_list); +#endif + +static double +vtrack (segment* seg, maze* m) +{ + channel* chp = chanSearch(m->vchans, seg); + double f = ((double)seg->track_no)/(chp->cnt+1); + double left = chp->cp->bb.LL.x; + double right = chp->cp->bb.UR.x; + return left + f*(right-left); +} + +static int +htrack (segment* seg, maze* m) +{ + channel* chp = chanSearch(m->hchans, seg); + double f = 1.0 - ((double)seg->track_no)/(chp->cnt+1); + double lo = chp->cp->bb.LL.y; + double hi = chp->cp->bb.UR.y; + return lo + f*(hi-lo); +} + +typedef struct { + int d; + Agedge_t* e; +} epair_t; + +static pointf +addPoints(pointf p0, pointf p1) +{ + p0.x += p1.x; + p0.y += p1.y; + return p0; +} + +static void +attachOrthoEdges (maze* mp, int n_edges, route* route_list, splineInfo *sinfo, epair_t es[]) +{ + int irte = 0; + int i, ipt, npts;; + pointf* ispline = 0; + int splsz = 0; + pointf p, p1, q1; + route rte; + segment* seg; + Agedge_t* e; + + for (; irte < n_edges; irte++) { + e = es[irte].e; + p1 = addPoints(ND_coord(e->tail), ED_tail_port(e).p); + q1 = addPoints(ND_coord(e->head), ED_head_port(e).p); + + rte = route_list[irte]; + npts = 1 + 3*rte.n; + if (npts > splsz) { + if (ispline) free (ispline); + ispline = N_GNEW(npts, pointf); + splsz = npts; + } + + seg = rte.segs; + if (seg->isVert) { + p.x = vtrack(seg, mp); + p.y = p1.y; + } + else { + p.y = htrack(seg, mp); + p.x = p1.x; + } + ispline[0] = ispline[1] = p; + ipt = 2; + + for (i = 1;iisVert) + p.x = vtrack(seg, mp); + else + p.y = htrack(seg, mp); + ispline[ipt+2] = ispline[ipt+1] = ispline[ipt] = p; + ipt += 3; + } + + if (seg->isVert) { + p.x = vtrack(seg, mp); + p.y = q1.y; + } + else { + p.y = htrack(seg, mp); + p.x = q1.x; + } + ispline[ipt] = ispline[ipt+1] = p; + + if (Verbose > 1) + fprintf(stderr, "ortho %s %s\n", e->tail->name, e->head->name); + clip_and_install(e, e->head, ispline, npts, sinfo); + } + free(ispline); +} + +static int +edgeLen (Agedge_t* e) +{ + pointf p = ND_coord(e->tail); + pointf q = ND_coord(e->head); + return (int)DIST2(p,q); +} + +static int edgecmp(epair_t* e0, epair_t* e1) +{ + return (e0->d - e1->d); +} + +void +orthoEdges (Agraph_t* g, int useLbls, splineInfo* sinfo) +{ + sgraph* sg; + maze* mp; + int n_edges = agnedges(g); + route* route_list = N_NEW (n_edges, route); + int i, gstart; + Agnode_t* n; + Agedge_t* e; + snode* sn; + snode* dn; + epair_t* es = N_GNEW(n_edges, epair_t); + cell* start; + cell* dest; + + mp = mkMaze (g); + sg = mp->sg; + + i = 0; + for (n = agfstnode (g); n; n = agnxtnode(g, n)) { + for (e = agfstout(g, n); e; e = agnxtout(g,e)) { + es[i].e = e; + es[i].d = edgeLen (e); + i++; + } + } + + qsort((char *)es, n_edges, sizeof(epair_t), (qsort_cmpf) edgecmp); + + fflush (stdout); + gstart = sg->nnodes; + PQgen (sg->nnodes+2); + sn = &sg->nodes[gstart]; + dn = &sg->nodes[gstart+1]; + for (i = 0; i < n_edges; i++) { + e = es[i].e; + start = CELL(e->tail); + dest = CELL(e->head); + + addNodeEdges (sg, dest, dn); + addNodeEdges (sg, start, sn); + shortPath (sg, dn, sn); + + route_list[i] = convertSPtoRoute(sg, sn, dn); + reset (sg); + } + PQfree (); + + mp->hchans = extractHChans (mp); + mp->vchans = extractVChans (mp); + assignSegs (agnedges(g), route_list, mp); + assignTracks (agnedges(g), route_list, mp); + /* emitGraph (stdout, mp, g, route_list); */ + attachOrthoEdges (mp, n_edges, route_list, sinfo, es); + + freeSGraph (sg); + dtclose (mp->hchans); + dtclose (mp->vchans); + for (i=0; i < agnedges(g); i++) + free (route_list[i].segs); + free (route_list); + freeMaze (mp); +} + +#ifdef DEBUG +#include +/* #include */ +#define TRANS 10 + +static char* prolog2 = +"%%!PS-Adobe-2.0\n\ +%%%%BoundingBox: (atend)\n\ +/point {\n\ + /Y exch def\n\ + /X exch def\n\ + newpath\n\ + X Y 3 0 360 arc fill\n\ +} def\n\ +/cell {\n\ + /Y exch def\n\ + /X exch def\n\ + /y exch def\n\ + /x exch def\n\ + newpath\n\ + x y moveto\n\ + x Y lineto\n\ + X Y lineto\n\ + X y lineto\n\ + closepath stroke\n\ +} def\n\ +/node {\n\ + /u exch def\n\ + /r exch def\n\ + /d exch def\n\ + /l exch def\n\ + newpath l d moveto\n\ + r d lineto r u lineto l u lineto\n\ + closepath fill\n\ +} def\n\ +\n"; + +static char* epilog2 = +"showpage\n\ +%%%%Trailer\n\ +%%%%BoundingBox: %d %d %d %d\n"; + +static point +coordOf (cell* cp, snode* np) +{ + point p; + if (cp->sides[M_TOP] == np) { + p.x = (cp->bb.LL.x + cp->bb.UR.x)/2; + p.y = cp->bb.UR.y; + } + else if (cp->sides[M_BOTTOM] == np) { + p.x = (cp->bb.LL.x + cp->bb.UR.x)/2; + p.y = cp->bb.LL.y; + } + else if (cp->sides[M_LEFT] == np) { + p.y = (cp->bb.LL.y + cp->bb.UR.y)/2; + p.x = cp->bb.LL.x; + } + else if (cp->sides[M_RIGHT] == np) { + p.y = (cp->bb.LL.y + cp->bb.UR.y)/2; + p.x = cp->bb.UR.x; + } + return p; +} + +static boxf +emitEdge (FILE* fp, Agedge_t* e, route rte, maze* m, int ix, boxf bb) +{ + int i, x, y; + boxf n = CELL(e->head)->bb; + segment* seg = rte.segs; + if (seg->isVert) { + x = vtrack(seg, m); + y = (n.UR.y + n.LL.y)/2; + } + else { + y = htrack(seg, m); + x = (n.UR.x + n.LL.x)/2; + } + bb.LL.x = MIN(bb.LL.x, SC*x); + bb.LL.y = MIN(bb.LL.y, SC*y); + bb.UR.x = MAX(bb.UR.x, SC*x); + bb.UR.y = MAX(bb.UR.y, SC*y); + fprintf (fp, "newpath %d %d moveto\n", SC*x, SC*y); + + for (i = 1;iisVert) { + x = vtrack(seg, m); + } + else { + y = htrack(seg, m); + } + bb.LL.x = MIN(bb.LL.x, SC*x); + bb.LL.y = MIN(bb.LL.y, SC*y); + bb.UR.x = MAX(bb.UR.x, SC*x); + bb.UR.y = MAX(bb.UR.y, SC*y); + fprintf (fp, "%d %d lineto\n", SC*x, SC*y); + } + + n = CELL(e->tail)->bb; + if (seg->isVert) { + x = vtrack(seg, m); + y = (n.UR.y + n.LL.y)/2; + } + else { + y = htrack(seg, m); + x = (n.LL.x + n.UR.x)/2; + } + bb.LL.x = MIN(bb.LL.x, SC*x); + bb.LL.y = MIN(bb.LL.y, SC*y); + bb.UR.x = MAX(bb.UR.x, SC*x); + bb.UR.y = MAX(bb.UR.y, SC*y); + fprintf (fp, "%d %d lineto stroke\n", SC*x, SC*y); + + return bb; +} + +static void +emitSearchGraph (FILE* fp, sgraph* sg) +{ + cell* cp; + snode* np; + sedge* ep; + point p; + int i; + fputs ("graph G {\n", fp); + fputs (" node[shape=point]\n", fp); + for (i = 0; i < sg->nnodes; i++) { + np = sg->nodes+i; + cp = np->cells[0]; + if (cp == np->cells[1]) { + pointf pf = midPt (cp); + p.x = pf.x; + p.y = pf.y; + } + else { + if (IsNode(cp)) cp = np->cells[1]; + p = coordOf (cp, np); + } + fprintf (fp, " %d [pos=\"%d,%d\"]\n", i, p.x, p.y); + } + for (i = 0; i < sg->nedges; i++) { + ep = sg->edges+i; + fprintf (fp, " %d -- %d\n", ep->v1, ep->v2); + } + fputs ("}\n", fp); +} + +static void +emitGraph (FILE* fp, maze* mp, Agraph_t* g, route* route_list) +{ + int i; + boxf bb, absbb; + box bbox; + Agnode_t* n; + Agedge_t* e; + + absbb.LL.x = absbb.LL.y = MAXDOUBLE; + absbb.UR.x = absbb.UR.y = -MAXDOUBLE; + + fprintf (fp, prolog2); + fprintf (fp, "%d %d translate\n", TRANS, TRANS); + + fputs ("0 0 1 setrgbcolor\n", fp); + for (i = 0; i < mp->ngcells; i++) { + bb = mp->gcells[i].bb; + printf ("%f %f %f %f node\n", bb.LL.x, bb.LL.y, bb.UR.x, bb.UR.y); + } + + i = 0; + for (n = agfstnode (g); n; n = agnxtnode(g, n)) { + for (e = agfstout(g, n); e; e = agnxtout(g,e)) { + absbb = emitEdge (fp, e, route_list[i], mp, i, absbb); + i++; + } + } + + fputs ("0.8 0.8 0.8 setrgbcolor\n", fp); + for (i = 0; i < mp->ncells; i++) { + bb = mp->cells[i].bb; + printf ("%f %f %f %f cell\n", bb.LL.x, bb.LL.y, bb.UR.x, bb.UR.y); + absbb.LL.x = MIN(absbb.LL.x, bb.LL.x); + absbb.LL.y = MIN(absbb.LL.y, bb.LL.y); + absbb.UR.x = MAX(absbb.UR.x, bb.UR.x); + absbb.UR.y = MAX(absbb.UR.y, bb.UR.y); + } + + bbox.LL.x = absbb.LL.x + TRANS; + bbox.LL.y = absbb.LL.y + TRANS; + bbox.UR.x = absbb.UR.x + TRANS; + bbox.UR.y = absbb.UR.y + TRANS; + fprintf (fp, epilog2, bbox.LL.x, bbox.LL.y, bbox.UR.x, bbox.UR.y); +} +#endif diff --git a/lib/ortho/ortho.h b/lib/ortho/ortho.h new file mode 100644 index 000000000..c25983dcf --- /dev/null +++ b/lib/ortho/ortho.h @@ -0,0 +1,8 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifndef ORTHO_H +#define ORTHO_H +#include + +void orthoEdges (Agraph_t* g, int useLbls, splineInfo* sinfo); + +#endif diff --git a/lib/ortho/partition.c b/lib/ortho/partition.c new file mode 100644 index 000000000..09ce6741f --- /dev/null +++ b/lib/ortho/partition.c @@ -0,0 +1,750 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include +#include + +#define NPOINTS 4 /* only rectangles */ +#define TRSIZE(ss) (5*(ss)+1) + +#define TR_FROM_UP 1 /* for traverse-direction */ +#define TR_FROM_DN 2 + +#define SP_SIMPLE_LRUP 1 /* for splitting trapezoids */ +#define SP_SIMPLE_LRDN 2 +#define SP_2UP_2DN 3 +#define SP_2UP_LEFT 4 +#define SP_2UP_RIGHT 5 +#define SP_2DN_LEFT 6 +#define SP_2DN_RIGHT 7 +#define SP_NOSPLIT -1 + +#define DOT(v0, v1) ((v0).x * (v1).x + (v0).y * (v1).y) +#define CROSS_SINE(v0, v1) ((v0).x * (v1).y - (v1).x * (v0).y) +#define LENGTH(v0) (sqrt((v0).x * (v0).x + (v0).y * (v0).y)) + +typedef struct { + int vnum; + int next; /* Circularly linked list */ + int prev; /* describing the monotone */ + int marked; /* polygon */ +} monchain_t; + +typedef struct { + pointf pt; + int vnext[4]; /* next vertices for the 4 chains */ + int vpos[4]; /* position of v in the 4 chains */ + int nextfree; +} vertexchain_t; + +static int chain_idx, mon_idx; + /* Table to hold all the monotone */ + /* polygons . Each monotone polygon */ + /* is a circularly linked list */ +static monchain_t* mchain; + /* chain init. information. This */ + /* is used to decide which */ + /* monotone polygon to split if */ + /* there are several other */ + /* polygons touching at the same */ + /* vertex */ +static vertexchain_t* vert; + /* contains position of any vertex in */ + /* the monotone chain for the polygon */ +static int* mon; + +/* return a new mon structure from the table */ +#define newmon() (++mon_idx) +/* return a new chain element from the table */ +#define new_chain_element() (++chain_idx) + +static void +convert (boxf bb, int flip, int ccw, pointf* pts) +{ + pts[0] = bb.LL; + pts[2] = bb.UR; + if (ccw) { + pts[1].x = bb.UR.x; + pts[1].y = bb.LL.y; + pts[3].x = bb.LL.x; + pts[3].y = bb.UR.y; + } + else { + pts[1].x = bb.LL.x; + pts[1].y = bb.UR.y; + pts[3].x = bb.UR.x; + pts[3].y = bb.LL.y; + } + if (flip) { + int i; + for (i = 0; i < NPOINTS; i++) { + double tmp = pts[i].y; + pts[i].y = pts[i].x; + pts[i].x = -tmp; + } + } +} + +static int +store (segment_t* seg, int first, pointf* pts) +{ + int i, last = first + NPOINTS - 1; + int j = 0; + + for (i = first; i <= last; i++, j++) { + if (i == first) { + seg[i].next = first+1; + seg[i].prev = last; + } + else if (i == last) { + seg[i].next = first; + seg[i].prev = last-1; + } + else { + seg[i].next = i+1; + seg[i].prev = i-1; + } + seg[i].is_inserted = FALSE; + seg[seg[i].prev].v1 = seg[i].v0 = pts[j]; + } + return (last+1); +} + +static void +genSegments (cell* cells, int ncells, boxf bb, segment_t* seg, int flip) +{ + int j = 0, i = 1; + pointf pts[4]; + + convert (bb, flip, 1, pts); + i = store (seg, i, pts); + for (j = 0; j < ncells; j++) { + convert (cells[j].bb, flip, 0, pts); + i = store (seg, i, pts); + } +} + +/* Generate a random permutation of the segments 1..n */ +static void +generateRandomOrdering(int n, int* permute) +{ + int i, j, tmp; + for (i = 0; i <= n; i++) permute[i] = i; + + for (i = 1; i <= n; i++) { + j = i + drand48() * (n + 1 - i); + if (j != i) { + tmp = permute[i]; + permute [i] = permute[j]; + permute [j] = tmp; + } + } +} + +/* Function returns TRUE if the trapezoid lies inside the polygon */ +static int +inside_polygon (trap_t *t, segment_t* seg) +{ + int rseg = t->rseg; + + if (t->state == ST_INVALID) + return 0; + + if ((t->lseg <= 0) || (t->rseg <= 0)) + return 0; + + if (((t->u0 <= 0) && (t->u1 <= 0)) || + ((t->d0 <= 0) && (t->d1 <= 0))) /* triangle */ + return (_greater_than(&seg[rseg].v1, &seg[rseg].v0)); + + return 0; +} + +static double +get_angle (pointf *vp0, pointf *vpnext, pointf *vp1) +{ + pointf v0, v1; + + v0.x = vpnext->x - vp0->x; + v0.y = vpnext->y - vp0->y; + + v1.x = vp1->x - vp0->x; + v1.y = vp1->y - vp0->y; + + if (CROSS_SINE(v0, v1) >= 0) /* sine is positive */ + return DOT(v0, v1)/LENGTH(v0)/LENGTH(v1); + else + return (-1.0 * DOT(v0, v1)/LENGTH(v0)/LENGTH(v1) - 2); +} + +/* (v0, v1) is the new diagonal to be added to the polygon. Find which */ +/* chain to use and return the positions of v0 and v1 in p and q */ +static int +get_vertex_positions (int v0, int v1, int *ip, int *iq) +{ + vertexchain_t *vp0, *vp1; + register int i; + double angle, temp; + int tp, tq; + + vp0 = &vert[v0]; + vp1 = &vert[v1]; + + /* p is identified as follows. Scan from (v0, v1) rightwards till */ + /* you hit the first segment starting from v0. That chain is the */ + /* chain of our interest */ + + angle = -4.0; + for (i = 0; i < 4; i++) + { + if (vp0->vnext[i] <= 0) + continue; + if ((temp = get_angle(&vp0->pt, &(vert[vp0->vnext[i]].pt), + &vp1->pt)) > angle) + { + angle = temp; + tp = i; + } + } + + *ip = tp; + + /* Do similar actions for q */ + + angle = -4.0; + for (i = 0; i < 4; i++) + { + if (vp1->vnext[i] <= 0) + continue; + if ((temp = get_angle(&vp1->pt, &(vert[vp1->vnext[i]].pt), + &vp0->pt)) > angle) + { + angle = temp; + tq = i; + } + } + + *iq = tq; + + return 0; +} + +/* v0 and v1 are specified in anti-clockwise order with respect to + * the current monotone polygon mcur. Split the current polygon into + * two polygons using the diagonal (v0, v1) + */ +static int +make_new_monotone_poly (int mcur, int v0, int v1) +{ + int p, q, ip, iq; + int mnew = newmon(); + int i, j, nf0, nf1; + vertexchain_t *vp0, *vp1; + + vp0 = &vert[v0]; + vp1 = &vert[v1]; + + get_vertex_positions(v0, v1, &ip, &iq); + + p = vp0->vpos[ip]; + q = vp1->vpos[iq]; + + /* At this stage, we have got the positions of v0 and v1 in the */ + /* desired chain. Now modify the linked lists */ + + i = new_chain_element(); /* for the new list */ + j = new_chain_element(); + + mchain[i].vnum = v0; + mchain[j].vnum = v1; + + mchain[i].next = mchain[p].next; + mchain[mchain[p].next].prev = i; + mchain[i].prev = j; + mchain[j].next = i; + mchain[j].prev = mchain[q].prev; + mchain[mchain[q].prev].next = j; + + mchain[p].next = q; + mchain[q].prev = p; + + nf0 = vp0->nextfree; + nf1 = vp1->nextfree; + + vp0->vnext[ip] = v1; + + vp0->vpos[nf0] = i; + vp0->vnext[nf0] = mchain[mchain[i].next].vnum; + vp1->vpos[nf1] = j; + vp1->vnext[nf1] = v0; + + vp0->nextfree++; + vp1->nextfree++; + +#ifdef DEBUG + fprintf(stderr, "make_poly: mcur = %d, (v0, v1) = (%d, %d)\n", + mcur, v0, v1); + fprintf(stderr, "next posns = (p, q) = (%d, %d)\n", p, q); +#endif + + mon[mcur] = p; + mon[mnew] = i; + return mnew; +} + +/* recursively visit all the trapezoids */ +static int +traverse_polygon (int* visited, boxf* decomp, int size, segment_t* seg, trap_t* tr, + int mcur, int trnum, int from, int flip, int dir) +{ + trap_t *t = &tr[trnum]; + int mnew; + int v0, v1; + int retval; + int do_switch = FALSE; + + if ((trnum <= 0) || visited[trnum]) + return size; + + visited[trnum] = TRUE; + + if ((t->hi.y > t->lo.y) && + (seg[t->lseg].v0.x == seg[t->lseg].v1.x) && + (seg[t->rseg].v0.x == seg[t->rseg].v1.x)) { + if (flip) { + decomp[size].LL.x = t->lo.y; + decomp[size].LL.y = -seg[t->rseg].v0.x; + decomp[size].UR.x = t->hi.y; + decomp[size].UR.y = -seg[t->lseg].v0.x; + } else { + decomp[size].LL.x = seg[t->lseg].v0.x; + decomp[size].LL.y = t->lo.y; + decomp[size].UR.x = seg[t->rseg].v0.x; + decomp[size].UR.y = t->hi.y; + } + size++; + } + + /* We have much more information available here. */ + /* rseg: goes upwards */ + /* lseg: goes downwards */ + + /* Initially assume that dir = TR_FROM_DN (from the left) */ + /* Switch v0 and v1 if necessary afterwards */ + + + /* special cases for triangles with cusps at the opposite ends. */ + /* take care of this first */ + if ((t->u0 <= 0) && (t->u1 <= 0)) + { + if ((t->d0 > 0) && (t->d1 > 0)) /* downward opening triangle */ + { + v0 = tr[t->d1].lseg; + v1 = t->lseg; + if (from == t->d1) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + } + } + else + { + retval = SP_NOSPLIT; /* Just traverse all neighbours */ + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + } + } + + else if ((t->d0 <= 0) && (t->d1 <= 0)) + { + if ((t->u0 > 0) && (t->u1 > 0)) /* upward opening triangle */ + { + v0 = t->rseg; + v1 = tr[t->u0].rseg; + if (from == t->u1) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + } + } + else + { + retval = SP_NOSPLIT; /* Just traverse all neighbours */ + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + } + } + + else if ((t->u0 > 0) && (t->u1 > 0)) + { + if ((t->d0 > 0) && (t->d1 > 0)) /* downward + upward cusps */ + { + v0 = tr[t->d1].lseg; + v1 = tr[t->u0].rseg; + retval = SP_2UP_2DN; + if (((dir == TR_FROM_DN) && (t->d1 == from)) || + ((dir == TR_FROM_UP) && (t->u1 == from))) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + } + } + else /* only downward cusp */ + { + if (_equal_to(&t->lo, &seg[t->lseg].v1)) + { + v0 = tr[t->u0].rseg; + v1 = seg[t->lseg].next; + + retval = SP_2UP_LEFT; + if ((dir == TR_FROM_UP) && (t->u0 == from)) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + } + } + else + { + v0 = t->rseg; + v1 = tr[t->u0].rseg; + retval = SP_2UP_RIGHT; + if ((dir == TR_FROM_UP) && (t->u1 == from)) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + } + } + } + } + else if ((t->u0 > 0) || (t->u1 > 0)) /* no downward cusp */ + { + if ((t->d0 > 0) && (t->d1 > 0)) /* only upward cusp */ + { + if (_equal_to(&t->hi, &seg[t->lseg].v0)) + { + v0 = tr[t->d1].lseg; + v1 = t->lseg; + retval = SP_2DN_LEFT; + if (!((dir == TR_FROM_DN) && (t->d0 == from))) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + } + } + else + { + v0 = tr[t->d1].lseg; + v1 = seg[t->rseg].next; + + retval = SP_2DN_RIGHT; + if ((dir == TR_FROM_DN) && (t->d1 == from)) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + } + } + } + else /* no cusp */ + { + if (_equal_to(&t->hi, &seg[t->lseg].v0) && + _equal_to(&t->lo, &seg[t->rseg].v0)) + { + v0 = t->rseg; + v1 = t->lseg; + retval = SP_SIMPLE_LRDN; + if (dir == TR_FROM_UP) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + } + } + else if (_equal_to(&t->hi, &seg[t->rseg].v1) && + _equal_to(&t->lo, &seg[t->lseg].v1)) + { + v0 = seg[t->rseg].next; + v1 = seg[t->lseg].next; + + retval = SP_SIMPLE_LRUP; + if (dir == TR_FROM_UP) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->d0, trnum, flip, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mnew, t->u1, trnum, flip, TR_FROM_DN); + } + } + else /* no split possible */ + { + retval = SP_NOSPLIT; + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u0, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d0, trnum, flip, TR_FROM_UP); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->u1, trnum, flip, TR_FROM_DN); + size = traverse_polygon (visited, decomp, size, seg, tr, mcur, t->d1, trnum, flip, TR_FROM_UP); + } + } + } + + return size; +} + +static int +monotonate_trapezoids(int nsegs, segment_t*seg, trap_t* tr, + int flip, boxf* decomp) +{ + int i, size; + int tr_start; + int tr_size = TRSIZE(nsegs); + int* visited = N_NEW(tr_size,int); + + mchain = N_NEW(tr_size, monchain_t); + vert = N_NEW(nsegs+1,vertexchain_t); + mon = N_NEW(nsegs, int); + + /* First locate a trapezoid which lies inside the polygon */ + /* and which is triangular */ + for (i = 0; i < TRSIZE(nsegs); i++) + if (inside_polygon(&tr[i], seg)) break; + tr_start = i; + + /* Initialise the mon data-structure and start spanning all the */ + /* trapezoids within the polygon */ + + for (i = 1; i <= nsegs; i++) { + mchain[i].prev = seg[i].prev; + mchain[i].next = seg[i].next; + mchain[i].vnum = i; + vert[i].pt = seg[i].v0; + vert[i].vnext[0] = seg[i].next; /* next vertex */ + vert[i].vpos[0] = i; /* locn. of next vertex */ + vert[i].nextfree = 1; + } + + chain_idx = nsegs; + mon_idx = 0; + mon[0] = 1; /* position of any vertex in the first */ + /* chain */ + + /* traverse the polygon */ + if (tr[tr_start].u0 > 0) + size = traverse_polygon (visited, decomp, 0, seg, tr, 0, tr_start, tr[tr_start].u0, flip, TR_FROM_UP); + else if (tr[tr_start].d0 > 0) + size = traverse_polygon (visited, decomp, 0, seg, tr, 0, tr_start, tr[tr_start].d0, flip, TR_FROM_DN); + + free (visited); + free (mchain); + free (vert); + free (mon); + + /* return the number of rects created */ + return size; +} + +static int +rectIntersect (boxf *d, const boxf *r0, const boxf *r1) +{ + double t; + + t = (r0->LL.x > r1->LL.x) ? r0->LL.x : r1->LL.x; + d->UR.x = (r0->UR.x < r1->UR.x) ? r0->UR.x : r1->UR.x; + d->LL.x = t; + + t = (r0->LL.y > r1->LL.y) ? r0->LL.y : r1->LL.y; + d->UR.y = (r0->UR.y < r1->UR.y) ? r0->UR.y : r1->UR.y; + d->LL.y = t; + + if ((d->LL.x >= d->UR.x) || + (d->LL.y >= d->UR.y)) + return 0; + + return 1; +} + +#ifdef DEBUG +static void +dumpTrap (trap_t* tr, int n) +{ + int i; + for (i = 1; i <= n; i++) { + tr++; + fprintf (stderr, "%d : %d %d (%f,%f) (%f,%f) %d %d %d %d\n", i, + tr->lseg, tr->rseg, tr->hi.x, tr->hi.y, tr->lo.x, tr->lo.y, + tr->u0, tr->u1, tr->d0, tr->d1); + fprintf (stderr, " %d %d %d %d\n", tr->sink, tr->usave, + tr->uside, tr->state); + } + fprintf (stderr, "====\n"); +} + +static void +dumpSegs (segment_t* sg, int n) +{ + int i; + for (i = 1; i <= n; i++) { + sg++; + fprintf (stderr, "%d : (%f,%f) (%f,%f) %d %d %d %d %d\n", i, + sg->v0.x, sg->v0.y, sg->v1.x, sg->v1.y, + sg->is_inserted, sg->root0, sg->root1, sg->next, sg->prev); + } + fprintf (stderr, "====\n"); +} +#endif + +boxf* +partition (cell* cells, int ncells, int* nrects, boxf bb) +{ + int nsegs = 4*(ncells+1); + segment_t* segs = N_GNEW(nsegs+1, segment_t); + int* permute = N_NEW(nsegs+1, int); + int hd_size, vd_size; + int i, j, cnt = 0; + boxf* rs; + int ntraps = TRSIZE(nsegs); + trap_t* trs = N_GNEW(ntraps, trap_t); + boxf* hor_decomp = N_NEW(ntraps, boxf); + boxf* vert_decomp = N_NEW(ntraps, boxf); + int nt; + + /* fprintf (stderr, "cells = %d segs = %d traps = %d\n", ncells, nsegs, ntraps); */ + genSegments (cells, ncells, bb, segs, 0); +#if 0 +fprintf (stderr, "%d\n\n", ncells+1); +for (i = 1; i<= nsegs; i++) { + if (i%4 == 1) fprintf(stderr, "4\n"); + fprintf (stderr, "%f %f\n", segs[i].v0.x, segs[i].v0.y); + if (i%4 == 0) fprintf(stderr, "\n"); +} +#endif + srand48(173); + generateRandomOrdering (nsegs, permute); + nt = construct_trapezoids(nsegs, segs, permute, ntraps, trs); + /* fprintf (stderr, "hor traps = %d\n", nt); */ + hd_size = monotonate_trapezoids (nsegs, segs, trs, 0, hor_decomp); + + genSegments (cells, ncells, bb, segs, 1); + generateRandomOrdering (nsegs, permute); + nt = construct_trapezoids(nsegs, segs, permute, ntraps, trs); + /* fprintf (stderr, "ver traps = %d\n", nt); */ + vd_size = monotonate_trapezoids (nsegs, segs, trs, 1, vert_decomp); + + rs = N_NEW (hd_size*vd_size, boxf); + for (i=0; i + +extern boxf* partition (cell*, int, int*, boxf); + +#endif diff --git a/lib/ortho/rawgraph.c b/lib/ortho/rawgraph.c new file mode 100644 index 000000000..a3e9fa7f9 --- /dev/null +++ b/lib/ortho/rawgraph.c @@ -0,0 +1,153 @@ +/* vim:set shiftwidth=4 ts=8: */ + /* Implements graph.h */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "rawgraph.h" +#include "memory.h" +#include "intset.h" + +#define UNSCANNED 0 +#define SCANNING 1 +#define SCANNED 2 + +rawgraph* +make_graph(int n) +{ + int i; + rawgraph* g = NEW(rawgraph); + g->nvs = n; + g->vertices = N_NEW(n, vertex); + for(i=0;ivertices[i].adj_list = openIntSet (); + g->vertices[i].color = UNSCANNED; + } + return g; +} + +void +free_graph(rawgraph* g) +{ + int i; + for(i=0;invs;i++) + dtclose(g->vertices[i].adj_list); + free (g->vertices); + free (g); +} + +void +insert_edge(rawgraph* g, int v1, int v2) +{ + intitem obj; + + obj.id = v2; + dtinsert(g->vertices[v1].adj_list,&obj); +} + +void +remove_redge(rawgraph* g, int v1, int v2) +{ + intitem obj; + obj.id = v2; + dtdelete (g->vertices[v1].adj_list, &obj); + obj.id = v1; + dtdelete (g->vertices[v2].adj_list, &obj); +} + +int +edge_exists(rawgraph* g, int v1, int v2) +{ + return (dtmatch (g->vertices[v1].adj_list, &v2) != 0); +} + +typedef struct { + int top; + int* vals; +} stack; + +static stack* +mkStack (int i) +{ + stack* sp = NEW(stack); + sp->vals = N_NEW(i,int); + sp->top = -1; + return sp; +} + +static void +freeStack (stack* s) +{ + free (s->vals); + free (s); +} + +static void +pushStack (stack* s, int i) +{ + s->top++; + s->vals[s->top] = i; +} + +static int +popStack (stack* s) +{ + int v; + + if (s->top == -1) return -1; + v = s->vals[s->top]; + s->top--; + return v; +} + +static int +DFS_visit(rawgraph* g, int v, int time, stack* sp) +{ + Dt_t* adj; + Dtlink_t* link; + int id; + vertex* vp; + + vp = g->vertices + v; + vp->color = SCANNING; + /* g->vertices[v].d = time; */ + adj = vp->adj_list; + time = time + 1; + + for(link = dtflatten (adj); link; link = dtlink(adj,link)) { + id = ((intitem*)dtobj(adj,link))->id; + if(g->vertices[id].color == UNSCANNED) + time = DFS_visit(g, id, time, sp); + } + vp->color = SCANNED; + /* g->vertices[v].f = time; */ + pushStack (sp, v); + return (time + 1); +} + +void +top_sort(rawgraph* g) +{ + int i, v; + int time = 0; + int count = 0; + stack* sp; + + if (g->nvs == 0) return; + if (g->nvs == 1) { + g->vertices[0].topsort_order = count; + return; + } + + sp = mkStack (g->nvs); + for(i=0;invs;i++) { + if(g->vertices[i].color == UNSCANNED) + time = DFS_visit(g, i, time, sp); + } + while((v = popStack(sp)) >= 0) { + g->vertices[v].topsort_order = count; + count++; + } + freeStack (sp); +} diff --git a/lib/ortho/rawgraph.h b/lib/ortho/rawgraph.h new file mode 100644 index 000000000..fc5daae34 --- /dev/null +++ b/lib/ortho/rawgraph.h @@ -0,0 +1,29 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifndef RAWGRAPH_H +#define RAWGRAPH_H + +#include + +typedef struct { + int color; + int topsort_order; + Dt_t* adj_list; /* adj_list */ +} vertex; + +typedef struct { + int nvs; + vertex* vertices; +} rawgraph; + +extern rawgraph* make_graph(int n); /* makes a graph wih n vertices, 0 edges */ +extern void free_graph(rawgraph*); + /* inserts edge FROM v1 to v2 */ +extern void insert_edge(rawgraph*, int v1, int v2); + /* removes any edge between v1 to v2 -- irrespective of direction */ +extern void remove_redge(rawgraph*, int v1, int v2); + /* tests if there is an edge FROM v1 TO v2 */ +extern int edge_exists(rawgraph*, int v1, int v2); + /* topologically sorts the directed graph */ +extern void top_sort(rawgraph*); + +#endif diff --git a/lib/ortho/sgraph.c b/lib/ortho/sgraph.c new file mode 100644 index 000000000..d350d84d6 --- /dev/null +++ b/lib/ortho/sgraph.c @@ -0,0 +1,252 @@ +/* vim:set shiftwidth=4 ts=8: */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include "memory.h" +#include "sgraph.h" +#include "fPQ.h" + +#if 0 +/* Max. number of maze segments around a node */ +static int MaxNodeBoundary = 100; + +typedef struct { + int left, right, up, down; +} irect; + +/* nodes in the search graph correspond to line segments in the + * grid formed by n_hlines horizontal lines and n_vlines vertical lines. + * The vertical segments are enumerated first, top to bottom, left to right. + * Then the horizontal segments left to right, top to bottom. For example, + * with an array of 4 vertical and 3 horizontal lines, we have + * + * |--14--|--15--|--16--| + * 1 3 5 7 + * |--11--|--12--|--13--| + * 0 2 4 6 + * |-- 8--|-- 9--|--10--| + */ +static irect +get_indices(orthograph* OG,int i, int j) +{ + irect r; + int hl = OG->n_hlines-1; + int vl = OG->n_vlines-1; + r.left = i*hl + j; + r.right = r.left + hl; + r.down = (vl+1)*hl + j*vl + i; + r.up = r.down + vl; + return r; +} + +static irect +find_boundary(orthograph* G, int n) +{ + rect R = G->Nodes[n]; + irect r; + int i; + + for (i = 0; i < G->n_vlines; i++) { + if (R.left == G->v_lines[i]) { + r.left = i; + break; + } + } + for (; i < G->n_vlines; i++) { + if (R.right == G->v_lines[i]) { + r.right = i; + break; + } + } + for (i = 0; i < G->n_hlines; i++) { + if (R.down == G->h_lines[i]) { + r.down = i; + break; + } + } + for (; i < G->n_hlines; i++) { + if (R.up == G->h_lines[i]) { + r.up = i; + break; + } + } + return r; +} +#endif + +void +gsave (sgraph* G) +{ + int i; + G->save_nnodes = G->nnodes; + G->save_nedges = G->nedges; + for (i = 0; i < G->nnodes; i++) + G->nodes[i].save_n_adj = G->nodes[i].n_adj; +} + +void +reset(sgraph* G) +{ + int i; + G->nnodes = G->save_nnodes; + G->nedges = G->save_nedges; + for (i = 0; i < G->nnodes; i++) + G->nodes[i].n_adj = G->nodes[i].save_n_adj; + for (; i < G->nnodes+2; i++) + G->nodes[i].n_adj = 0; +} + +void +initSEdges (sgraph* g, int maxdeg) +{ + int i; + int* adj = N_NEW (6*g->nnodes + 2*maxdeg, int); + g->edges = N_NEW (3*g->nnodes + maxdeg, sedge); + for (i = 0; i < g->nnodes; i++) { + g->nodes[i].adj_edge_list = adj; + adj += 6; + } + for (; i < g->nnodes+2; i++) { + g->nodes[i].adj_edge_list = adj; + adj += maxdeg; + } +} + +sgraph* +createSGraph (int nnodes) +{ + sgraph* g = NEW(sgraph); + + /* create the nodes vector in the search graph */ + g->nnodes = 0; + g->nodes = N_NEW(nnodes, snode); + return g; +} + +snode* +createSNode (sgraph* g) +{ + snode* np = g->nodes+g->nnodes; + np->index = g->nnodes; + g->nnodes++; + return np; +} + +static void +addEdgeToNode (snode* np, sedge* e, int idx) +{ + np->adj_edge_list[np->n_adj] = idx; + np->n_adj++; +} + +sedge* +createSEdge (sgraph* g, snode* v1, snode* v2, double wt) +{ + sedge* e; + int idx = g->nedges++; + + e = g->edges + idx; + e->v1 = v1->index; + e->v2 = v2->index; + e->weight = wt; + e->cnt = 0; + + addEdgeToNode (v1, e, idx); + addEdgeToNode (v2, e, idx); + + return e; +} + +void +freeSGraph (sgraph* g) +{ + free (g->nodes[0].adj_edge_list); + free (g->nodes); + free (g->edges); + free (g); +} + +#include "fPQ.h" + +/* shortest path: + * Constructs the path of least weight between from and to. + * + * Assumes graph, node and edge type, and that nodes + * have associated values N_VAL, N_IDX, and N_DAD, the first two + * being ints, the last being a node*. Edges have a E_WT function + * to specify the edge length or weight. + * + * Assumes there are functions: + * agnnodes: graph -> int number of nodes in the graph + * agfstnode, agnxtnode : iterators over the nodes in the graph + * agfstedge, agnxtedge : iterators over the edges attached to a node + * adjacentNode : given an edge e and an endpoint n of e, returns the + * other endpoint. + * + * The path is given by + * to, N_DAD(to), N_DAD(N_DAD(to)), ..., from + */ + +#define UNSEEN INT_MIN + +static snode* +adjacentNode(sgraph* g, sedge* e, snode* n) +{ + if (e->v1==n->index) + return (&(g->nodes[e->v2])); + else + return (&(g->nodes[e->v1])); +} + +void +shortPath (sgraph* g, snode* from, snode* to) +{ + snode* n; + sedge* e; + snode* adjn; + int d; + int x, y; + + for (x = 0; xnnodes; x++) { + snode* temp; + temp = &(g->nodes[x]); + N_VAL(temp) = UNSEEN; + } + + PQinit(); + PQ_insert (from); + N_DAD(from) = NULL; + N_VAL(from) = 0; + + while ((n = PQremove())) { + N_VAL(n) *= -1; + if (n == to) break; + for (y=0; yn_adj; y++) { + e = &(g->edges[n->adj_edge_list[y]]); + adjn = adjacentNode(g, e, n); + if (N_VAL(adjn) < 0) { + d = -(N_VAL(n) + E_WT(e)); + if (N_VAL(adjn) == UNSEEN) { + N_VAL(adjn) = d; + PQ_insert(adjn); + N_DAD(adjn) = n; + N_EDGE(adjn) = e; + } + else { + if (N_VAL(adjn) < d) { + PQupdate(adjn, d); + N_DAD(adjn) = n; + N_EDGE(adjn) = e; + } + } + } + } + } + + /* PQfree(); */ + return ; +} + diff --git a/lib/ortho/sgraph.h b/lib/ortho/sgraph.h new file mode 100644 index 000000000..6c7d787ba --- /dev/null +++ b/lib/ortho/sgraph.h @@ -0,0 +1,51 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifndef SEARCH_G_H +#define SEARCH_G_H + +#include "structures.h" + +typedef struct snode snode; +typedef struct sedge sedge; + +struct snode { + int n_val, n_idx; + snode* n_dad; + sedge* n_edge; + short n_adj; + short save_n_adj; + struct cell* cells[2]; + + /* edges incident on this node + * -- stored as indices of the nodes vector in the graph + */ + int* adj_edge_list; + int index; + boolean isVert; /* true if node corresponds to vertical segment */ +}; + +struct sedge { + double weight; /* weight of edge */ + int cnt; /* paths using edge */ + /* end-points of the edge + * -- stored as indices of the nodes vector in the graph + */ + int v1, v2; +}; + +typedef struct { + int nnodes, nedges; + int save_nnodes, save_nedges; + snode* nodes; + sedge* edges; +} sgraph; + +extern void reset(sgraph*); +extern void gsave(sgraph*); +extern sgraph* createSGraph(int); +extern void freeSGraph (sgraph*); +extern void initSEdges (sgraph* g, int maxdeg); +extern void shortPath (sgraph* g, snode* from, snode* to); +extern snode* createSNode (sgraph*); +extern sedge* createSEdge (sgraph* g, snode* v0, snode* v1, double wt); + +#endif diff --git a/lib/ortho/structures.h b/lib/ortho/structures.h new file mode 100644 index 000000000..55cc2efe2 --- /dev/null +++ b/lib/ortho/structures.h @@ -0,0 +1,82 @@ +/* vim:set shiftwidth=4 ts=8: */ +#ifndef STRUCTURES_H +#define STRUCTURES_H + +#include "types.h" +#include "graph.h" +#include "rawgraph.h" + +typedef struct { + double p1, p2; +} paird; + +typedef struct { + int a,b; +} pair; + +typedef struct { + pair t1, t2; +} pair2; + +typedef enum {B_NODE, B_UP, B_LEFT, B_DOWN, B_RIGHT} bend; + +/* Example : segment connecting maze point (3,2) + * and (3,8) has isVert = 1, common coordinate = 3, p1 = 2, p2 = 8 + */ +typedef struct segment { + boolean isVert; + boolean flipped; + double comm_coord; /* the common coordinate */ + paird p; /* end points */ + bend l1, l2; + int ind_no; /* index number of this segment in its channel */ + int track_no; /* track number assigned in the channel */ + struct segment* prev; + struct segment* next; +} segment; + +typedef struct { + int n; + segment* segs; +} route; + +typedef struct { + Dtlink_t link; + paird p; /* extrema of channel */ + int cnt; /* number of segments */ + segment** seg_list; /* array of segment pointers */ + rawgraph* G; + struct cell* cp; +} channel; + +#if 0 +typedef struct { + int i1, i2, j; + int cnt; + int* seg_list; /* list of indices of the segment list */ + + rawgraph* G; +} hor_channel; + +typedef struct { + hor_channel* hs; + int cnt; +} vhor_channel; + +typedef struct { + int i, j1, j2; + int cnt; + int* seg_list; /* list of indices of the segment list */ + + rawgraph* G; +} vert_channel; + +typedef struct { + vert_channel* vs; + int cnt; +} vvert_channel; +#endif + +#define N_DAD(n) (n)->n_dad + +#endif diff --git a/lib/ortho/trap.h b/lib/ortho/trap.h new file mode 100644 index 000000000..7d664a0d0 --- /dev/null +++ b/lib/ortho/trap.h @@ -0,0 +1,45 @@ +#ifndef TRAP_H +#define TRAP_H + +/* Segment attributes */ + +typedef struct { + pointf v0, v1; /* two endpoints */ + int is_inserted; /* inserted in trapezoidation yet ? */ + int root0, root1; /* root nodes in Q */ + int next; /* Next logical segment */ + int prev; /* Previous segment */ +} segment_t; + + +/* Trapezoid attributes */ + +typedef struct { + int lseg, rseg; /* two adjoining segments */ + pointf hi, lo; /* max/min y-values */ + int u0, u1; + int d0, d1; + int sink; /* pointer to corresponding in Q */ + int usave, uside; /* I forgot what this means */ + int state; +} trap_t; + +#define ST_VALID 1 /* for trapezium state */ +#define ST_INVALID 2 + +#define C_EPS 1.0e-7 /* tolerance value: Used for making */ + /* all decisions about collinearity or */ + /* left/right of segment. Decrease */ + /* this value if the input points are */ + /* spaced very close together */ +#define FP_EQUAL(s, t) (fabs(s - t) <= C_EPS) + +#define _equal_to(v0,v1) \ + (FP_EQUAL((v0)->y, (v1)->y) && FP_EQUAL((v0)->x, (v1)->x)) + +#define _greater_than(v0, v1) \ + (((v0)->y > (v1)->y + C_EPS) ? TRUE : (((v0)->y < (v1)->y - C_EPS) ? FALSE : ((v0)->x > (v1)->x))) + +extern int construct_trapezoids(int, segment_t*, int*, int, trap_t*); + +#endif diff --git a/lib/ortho/trapezoid.c b/lib/ortho/trapezoid.c new file mode 100644 index 000000000..a518d91c9 --- /dev/null +++ b/lib/ortho/trapezoid.c @@ -0,0 +1,1066 @@ +/* vim:set shiftwidth=4 ts=8: */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#define __USE_ISOC99 +#include +#include +#include +#include +#include + +/* Node types */ + +#define T_X 1 +#define T_Y 2 +#define T_SINK 3 + +#define FIRSTPT 1 /* checking whether pt. is inserted */ +#define LASTPT 2 + +#define S_LEFT 1 /* for merge-direction */ +#define S_RIGHT 2 + +#define INF 1<<30 + +#define CROSS(v0, v1, v2) (((v1).x - (v0).x)*((v2).y - (v0).y) - \ + ((v1).y - (v0).y)*((v2).x - (v0).x)) + +typedef struct { + int nodetype; /* Y-node or S-node */ + int segnum; + pointf yval; + int trnum; + int parent; /* doubly linked DAG */ + int left, right; /* children */ +} qnode_t; + +/* static int chain_idx, op_idx, mon_idx; */ +static int q_idx; +static int tr_idx; +static int QSIZE; +static int TRSIZE; + +/* Return a new node to be added into the query tree */ +static int newnode() +{ + if (q_idx < QSIZE) + return q_idx++; + else { + fprintf(stderr, "newnode: Query-table overflow\n"); + assert(0); + return -1; + } +} + +/* Return a free trapezoid */ +static int newtrap(trap_t* tr) +{ + if (tr_idx < TRSIZE) { + tr[tr_idx].lseg = -1; + tr[tr_idx].rseg = -1; + tr[tr_idx].state = ST_VALID; + return tr_idx++; + } + else { + fprintf(stderr, "newtrap: Trapezoid-table overflow %d\n", tr_idx); + assert(0); + return -1; + } +} + +/* Return the maximum of the two points into the yval structure */ +static int _max (pointf *yval, pointf *v0, pointf *v1) +{ + if (v0->y > v1->y + C_EPS) + *yval = *v0; + else if (FP_EQUAL(v0->y, v1->y)) + { + if (v0->x > v1->x + C_EPS) + *yval = *v0; + else + *yval = *v1; + } + else + *yval = *v1; + + return 0; +} + +/* Return the minimum of the two points into the yval structure */ +static int _min (pointf *yval, pointf *v0, pointf *v1) +{ + if (v0->y < v1->y - C_EPS) + *yval = *v0; + else if (FP_EQUAL(v0->y, v1->y)) + { + if (v0->x < v1->x) + *yval = *v0; + else + *yval = *v1; + } + else + *yval = *v1; + + return 0; +} + +static int _greater_than_equal_to (pointf *v0, pointf *v1) +{ + if (v0->y > v1->y + C_EPS) + return TRUE; + else if (v0->y < v1->y - C_EPS) + return FALSE; + else + return (v0->x >= v1->x); +} + +static int _less_than (pointf *v0, pointf *v1) +{ + if (v0->y < v1->y - C_EPS) + return TRUE; + else if (v0->y > v1->y + C_EPS) + return FALSE; + else + return (v0->x < v1->x); +} + +/* Initilialise the query structure (Q) and the trapezoid table (T) + * when the first segment is added to start the trapezoidation. The + * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes + * + * 4 + * ----------------------------------- + * \ + * 1 \ 2 + * \ + * ----------------------------------- + * 3 + */ + +static int +init_query_structure(int segnum, segment_t* seg, trap_t* tr, qnode_t* qs) +{ + int i1, i2, i3, i4, i5, i6, i7, root; + int t1, t2, t3, t4; + segment_t *s = &seg[segnum]; + + i1 = newnode(); + qs[i1].nodetype = T_Y; + _max(&qs[i1].yval, &s->v0, &s->v1); /* root */ + root = i1; + + qs[i1].right = i2 = newnode(); + qs[i2].nodetype = T_SINK; + qs[i2].parent = i1; + + qs[i1].left = i3 = newnode(); + qs[i3].nodetype = T_Y; + _min(&qs[i3].yval, &s->v0, &s->v1); /* root */ + qs[i3].parent = i1; + + qs[i3].left = i4 = newnode(); + qs[i4].nodetype = T_SINK; + qs[i4].parent = i3; + + qs[i3].right = i5 = newnode(); + qs[i5].nodetype = T_X; + qs[i5].segnum = segnum; + qs[i5].parent = i3; + + qs[i5].left = i6 = newnode(); + qs[i6].nodetype = T_SINK; + qs[i6].parent = i5; + + qs[i5].right = i7 = newnode(); + qs[i7].nodetype = T_SINK; + qs[i7].parent = i5; + + t1 = newtrap(tr); /* middle left */ + t2 = newtrap(tr); /* middle right */ + t3 = newtrap(tr); /* bottom-most */ + t4 = newtrap(tr); /* topmost */ + + tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval; + tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval; + tr[t4].hi.y = (double) (INF); + tr[t4].hi.x = (double) (INF); + tr[t3].lo.y = (double) -1* (INF); + tr[t3].lo.x = (double) -1* (INF); + tr[t1].rseg = tr[t2].lseg = segnum; + tr[t1].u0 = tr[t2].u0 = t4; + tr[t1].d0 = tr[t2].d0 = t3; + tr[t4].d0 = tr[t3].u0 = t1; + tr[t4].d1 = tr[t3].u1 = t2; + + tr[t1].sink = i6; + tr[t2].sink = i7; + tr[t3].sink = i4; + tr[t4].sink = i2; + + tr[t1].state = tr[t2].state = ST_VALID; + tr[t3].state = tr[t4].state = ST_VALID; + + qs[i2].trnum = t4; + qs[i4].trnum = t3; + qs[i6].trnum = t1; + qs[i7].trnum = t2; + + s->is_inserted = TRUE; + return root; +} + +/* Retun TRUE if the vertex v is to the left of line segment no. + * segnum. Takes care of the degenerate cases when both the vertices + * have the same y--cood, etc. + */ +static int +is_left_of (int segnum, segment_t* seg, pointf *v) +{ + segment_t *s = &seg[segnum]; + double area; + + if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */ + { + if (FP_EQUAL(s->v1.y, v->y)) + { + if (v->x < s->v1.x) + area = 1.0; + else + area = -1.0; + } + else if (FP_EQUAL(s->v0.y, v->y)) + { + if (v->x < s->v0.x) + area = 1.0; + else + area = -1.0; + } + else + area = CROSS(s->v0, s->v1, (*v)); + } + else /* v0 > v1 */ + { + if (FP_EQUAL(s->v1.y, v->y)) + { + if (v->x < s->v1.x) + area = 1.0; + else + area = -1.0; + } + else if (FP_EQUAL(s->v0.y, v->y)) + { + if (v->x < s->v0.x) + area = 1.0; + else + area = -1.0; + } + else + area = CROSS(s->v1, s->v0, (*v)); + } + + if (area > 0.0) + return TRUE; + else + return FALSE; +} + +/* Returns true if the corresponding endpoint of the given segment is */ +/* already inserted into the segment tree. Use the simple test of */ +/* whether the segment which shares this endpoint is already inserted */ +static int inserted (int segnum, segment_t* seg, int whichpt) +{ + if (whichpt == FIRSTPT) + return seg[seg[segnum].prev].is_inserted; + else + return seg[seg[segnum].next].is_inserted; +} + +/* This is query routine which determines which trapezoid does the + * point v lie in. The return value is the trapezoid number. + */ +static int +locate_endpoint (pointf *v, pointf *vo, int r, segment_t* seg, qnode_t* qs) +{ + qnode_t *rptr = &qs[r]; + + switch (rptr->nodetype) { + case T_SINK: + return rptr->trnum; + + case T_Y: + if (_greater_than(v, &rptr->yval)) /* above */ + return locate_endpoint(v, vo, rptr->right, seg, qs); + else if (_equal_to(v, &rptr->yval)) /* the point is already */ + { /* inserted. */ + if (_greater_than(vo, &rptr->yval)) /* above */ + return locate_endpoint(v, vo, rptr->right, seg, qs); + else + return locate_endpoint(v, vo, rptr->left, seg, qs); /* below */ + } + else + return locate_endpoint(v, vo, rptr->left, seg, qs); /* below */ + + case T_X: + if (_equal_to(v, &seg[rptr->segnum].v0) || + _equal_to(v, &seg[rptr->segnum].v1)) + { + if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */ + { + if (vo->x < v->x) + return locate_endpoint(v, vo, rptr->left, seg, qs); /* left */ + else + return locate_endpoint(v, vo, rptr->right, seg, qs); /* right */ + } + + else if (is_left_of(rptr->segnum, seg, vo)) + return locate_endpoint(v, vo, rptr->left, seg, qs); /* left */ + else + return locate_endpoint(v, vo, rptr->right, seg, qs); /* right */ + } + else if (is_left_of(rptr->segnum, seg, v)) + return locate_endpoint(v, vo, rptr->left, seg, qs); /* left */ + else + return locate_endpoint(v, vo, rptr->right, seg, qs); /* right */ + + default: + fprintf(stderr, "unexpected case in locate_endpoint\n"); + assert (0); + break; + } + return 1; /* stop warning */ +} + +/* Thread in the segment into the existing trapezoidation. The + * limiting trapezoids are given by tfirst and tlast (which are the + * trapezoids containing the two endpoints of the segment. Merges all + * possible trapezoids which flank this segment and have been recently + * divided because of its insertion + */ +static void +merge_trapezoids (int segnum, int tfirst, int tlast, int side, trap_t* tr, + qnode_t* qs) +{ + int t, tnext, cond; + int ptnext; + + /* First merge polys on the LHS */ + t = tfirst; + while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo)) + { + if (side == S_LEFT) + cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) || + (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum))); + else + cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) || + (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum))); + + if (cond) + { + if ((tr[t].lseg == tr[tnext].lseg) && + (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */ + { /* merge them */ + /* Use the upper node as the new node i.e. t */ + + ptnext = qs[tr[tnext].sink].parent; + + if (qs[ptnext].left == tr[tnext].sink) + qs[ptnext].left = tr[t].sink; + else + qs[ptnext].right = tr[t].sink; /* redirect parent */ + + + /* Change the upper neighbours of the lower trapezoids */ + + if ((tr[t].d0 = tr[tnext].d0) > 0) { + if (tr[tr[t].d0].u0 == tnext) + tr[tr[t].d0].u0 = t; + else if (tr[tr[t].d0].u1 == tnext) + tr[tr[t].d0].u1 = t; + } + + if ((tr[t].d1 = tr[tnext].d1) > 0) { + if (tr[tr[t].d1].u0 == tnext) + tr[tr[t].d1].u0 = t; + else if (tr[tr[t].d1].u1 == tnext) + tr[tr[t].d1].u1 = t; + } + + tr[t].lo = tr[tnext].lo; + tr[tnext].state = ST_INVALID; /* invalidate the lower */ + /* trapezium */ + } + else /* not good neighbours */ + t = tnext; + } + else /* do not satisfy the outer if */ + t = tnext; + + } /* end-while */ + +} + +/* Add in the new segment into the trapezoidation and update Q and T + * structures. First locate the two endpoints of the segment in the + * Q-structure. Then start from the topmost trapezoid and go down to + * the lower trapezoid dividing all the trapezoids in between . + */ +static int +add_segment (int segnum, segment_t* seg, trap_t* tr, qnode_t* qs) +{ + segment_t s; + int tu, tl, sk, tfirst, tlast; + int tfirstr, tlastr, tfirstl, tlastl; + int i1, i2, t, tn; + pointf tpt; + int tritop = 0, tribot = 0, is_swapped; + int tmptriseg; + + s = seg[segnum]; + if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */ + { + int tmp; + tpt = s.v0; + s.v0 = s.v1; + s.v1 = tpt; + tmp = s.root0; + s.root0 = s.root1; + s.root1 = tmp; + is_swapped = TRUE; + } + else is_swapped = FALSE; + + if ((is_swapped) ? !inserted(segnum, seg, LASTPT) : + !inserted(segnum, seg, FIRSTPT)) /* insert v0 in the tree */ + { + int tmp_d; + + tu = locate_endpoint(&s.v0, &s.v1, s.root0, seg, qs); + tl = newtrap(tr); /* tl is the new lower trapezoid */ + tr[tl].state = ST_VALID; + tr[tl] = tr[tu]; + tr[tu].lo.y = tr[tl].hi.y = s.v0.y; + tr[tu].lo.x = tr[tl].hi.x = s.v0.x; + tr[tu].d0 = tl; + tr[tu].d1 = 0; + tr[tl].u0 = tu; + tr[tl].u1 = 0; + + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + /* Now update the query structure and obtain the sinks for the */ + /* two trapezoids */ + + i1 = newnode(); /* Upper trapezoid sink */ + i2 = newnode(); /* Lower trapezoid sink */ + sk = tr[tu].sink; + + qs[sk].nodetype = T_Y; + qs[sk].yval = s.v0; + qs[sk].segnum = segnum; /* not really reqd ... maybe later */ + qs[sk].left = i2; + qs[sk].right = i1; + + qs[i1].nodetype = T_SINK; + qs[i1].trnum = tu; + qs[i1].parent = sk; + + qs[i2].nodetype = T_SINK; + qs[i2].trnum = tl; + qs[i2].parent = sk; + + tr[tu].sink = i1; + tr[tl].sink = i2; + tfirst = tl; + } + else /* v0 already present */ + { /* Get the topmost intersecting trapezoid */ + tfirst = locate_endpoint(&s.v0, &s.v1, s.root0, seg, qs); + tritop = 1; + } + + + if ((is_swapped) ? !inserted(segnum, seg, FIRSTPT) : + !inserted(segnum, seg, LASTPT)) /* insert v1 in the tree */ + { + int tmp_d; + + tu = locate_endpoint(&s.v1, &s.v0, s.root1, seg, qs); + + tl = newtrap(tr); /* tl is the new lower trapezoid */ + tr[tl].state = ST_VALID; + tr[tl] = tr[tu]; + tr[tu].lo.y = tr[tl].hi.y = s.v1.y; + tr[tu].lo.x = tr[tl].hi.x = s.v1.x; + tr[tu].d0 = tl; + tr[tu].d1 = 0; + tr[tl].u0 = tu; + tr[tl].u1 = 0; + + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + /* Now update the query structure and obtain the sinks for the */ + /* two trapezoids */ + + i1 = newnode(); /* Upper trapezoid sink */ + i2 = newnode(); /* Lower trapezoid sink */ + sk = tr[tu].sink; + + qs[sk].nodetype = T_Y; + qs[sk].yval = s.v1; + qs[sk].segnum = segnum; /* not really reqd ... maybe later */ + qs[sk].left = i2; + qs[sk].right = i1; + + qs[i1].nodetype = T_SINK; + qs[i1].trnum = tu; + qs[i1].parent = sk; + + qs[i2].nodetype = T_SINK; + qs[i2].trnum = tl; + qs[i2].parent = sk; + + tr[tu].sink = i1; + tr[tl].sink = i2; + tlast = tu; + } + else /* v1 already present */ + { /* Get the lowermost intersecting trapezoid */ + tlast = locate_endpoint(&s.v1, &s.v0, s.root1, seg, qs); + tribot = 1; + } + + /* Thread the segment into the query tree creating a new X-node */ + /* First, split all the trapezoids which are intersected by s into */ + /* two */ + + t = tfirst; /* topmost trapezoid */ + + while ((t > 0) && + _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo)) + /* traverse from top to bot */ + { + int t_sav, tn_sav; + sk = tr[t].sink; + i1 = newnode(); /* left trapezoid sink */ + i2 = newnode(); /* right trapezoid sink */ + + qs[sk].nodetype = T_X; + qs[sk].segnum = segnum; + qs[sk].left = i1; + qs[sk].right = i2; + + qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */ + qs[i1].trnum = t; + qs[i1].parent = sk; + + qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */ + qs[i2].trnum = tn = newtrap(tr); + tr[tn].state = ST_VALID; + qs[i2].parent = sk; + + if (t == tfirst) + tfirstr = tn; + if (_equal_to(&tr[t].lo, &tr[tlast].lo)) + tlastr = tn; + + tr[tn] = tr[t]; + tr[t].sink = i1; + tr[tn].sink = i2; + t_sav = t; + tn_sav = tn; + + /* error */ + + if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */ + { + fprintf(stderr, "add_segment: error\n"); + break; + } + + /* only one trapezoid below. partition t into two and make the */ + /* two resulting trapezoids t and tn as the upper neighbours of */ + /* the sole lower trapezoid */ + + else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0)) + { /* Only one trapezoid below */ + if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) + { /* continuation of a chain from abv. */ + if (tr[t].usave > 0) /* three upper neighbours */ + { + if (tr[t].uside == S_LEFT) + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = -1; + tr[tn].u1 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[tn].u0].d0 = tn; + tr[tr[tn].u1].d0 = tn; + } + else /* intersects in the right */ + { + tr[tn].u1 = -1; + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[t].u0; + tr[t].u0 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } + + tr[t].usave = tr[tn].usave = 0; + } + else /* No usave.... simple case */ + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d0 = tn; + } + } + else + { /* fresh seg. or upward cusp */ + int tmp_u = tr[t].u0; + int td0, td1; + if (((td0 = tr[tmp_u].d0) > 0) && + ((td1 = tr[tmp_u].d1) > 0)) + { /* upward cusp */ + if ((tr[td0].rseg > 0) && + !is_left_of(tr[td0].rseg, seg, &s.v1)) + { + tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d1 = tn; + } + else /* cusp going leftwards */ + { + tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1; + tr[tr[t].u0].d0 = t; + } + } + else /* fresh segment */ + { + tr[tr[t].u0].d0 = t; + tr[tr[t].u0].d1 = tn; + } + } + + if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && + FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) + { /* bottom forms a triangle */ + + if (is_swapped) + tmptriseg = seg[segnum].prev; + else + tmptriseg = seg[segnum].next; + + if ((tmptriseg > 0) && is_left_of(tmptriseg, seg, &s.v0)) + { + /* L-R downward cusp */ + tr[tr[t].d0].u0 = t; + tr[tn].d0 = tr[tn].d1 = -1; + } + else + { + /* R-L downward cusp */ + tr[tr[tn].d0].u1 = tn; + tr[t].d0 = tr[t].d1 = -1; + } + } + else + { + if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0)) + { + if (tr[tr[t].d0].u0 == t) /* passes thru LHS */ + { + tr[tr[t].d0].usave = tr[tr[t].d0].u1; + tr[tr[t].d0].uside = S_LEFT; + } + else + { + tr[tr[t].d0].usave = tr[tr[t].d0].u0; + tr[tr[t].d0].uside = S_RIGHT; + } + } + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = tn; + } + + t = tr[t].d0; + } + + + else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0)) + { /* Only one trapezoid below */ + if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) + { /* continuation of a chain from abv. */ + if (tr[t].usave > 0) /* three upper neighbours */ + { + if (tr[t].uside == S_LEFT) + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = -1; + tr[tn].u1 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[tn].u0].d0 = tn; + tr[tr[tn].u1].d0 = tn; + } + else /* intersects in the right */ + { + tr[tn].u1 = -1; + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[t].u0; + tr[t].u0 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } + + tr[t].usave = tr[tn].usave = 0; + } + else /* No usave.... simple case */ + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d0 = tn; + } + } + else + { /* fresh seg. or upward cusp */ + int tmp_u = tr[t].u0; + int td0, td1; + if (((td0 = tr[tmp_u].d0) > 0) && + ((td1 = tr[tmp_u].d1) > 0)) + { /* upward cusp */ + if ((tr[td0].rseg > 0) && + !is_left_of(tr[td0].rseg, seg, &s.v1)) + { + tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d1 = tn; + } + else + { + tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1; + tr[tr[t].u0].d0 = t; + } + } + else /* fresh segment */ + { + tr[tr[t].u0].d0 = t; + tr[tr[t].u0].d1 = tn; + } + } + + if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && + FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) + { /* bottom forms a triangle */ + /* int tmpseg; */ + + if (is_swapped) + tmptriseg = seg[segnum].prev; + else + tmptriseg = seg[segnum].next; + + /* if ((tmpseg > 0) && is_left_of(tmpseg, seg, &s.v0)) */ + if ((tmptriseg > 0) && is_left_of(tmptriseg, seg, &s.v0)) + { + /* L-R downward cusp */ + tr[tr[t].d1].u0 = t; + tr[tn].d0 = tr[tn].d1 = -1; + } + else + { + /* R-L downward cusp */ + tr[tr[tn].d1].u1 = tn; + tr[t].d0 = tr[t].d1 = -1; + } + } + else + { + if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0)) + { + if (tr[tr[t].d1].u0 == t) /* passes thru LHS */ + { + tr[tr[t].d1].usave = tr[tr[t].d1].u1; + tr[tr[t].d1].uside = S_LEFT; + } + else + { + tr[tr[t].d1].usave = tr[tr[t].d1].u0; + tr[tr[t].d1].uside = S_RIGHT; + } + } + tr[tr[t].d1].u0 = t; + tr[tr[t].d1].u1 = tn; + } + + t = tr[t].d1; + } + + /* two trapezoids below. Find out which one is intersected by */ + /* this segment and proceed down that one */ + + else + { + /* int tmpseg = tr[tr[t].d0].rseg; */ + double y0, yt; + pointf tmppt; + int tnext, i_d0, i_d1; + + i_d0 = i_d1 = FALSE; + if (FP_EQUAL(tr[t].lo.y, s.v0.y)) + { + if (tr[t].lo.x > s.v0.x) + i_d0 = TRUE; + else + i_d1 = TRUE; + } + else + { + tmppt.y = y0 = tr[t].lo.y; + yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y); + tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x); + + if (_less_than(&tmppt, &tr[t].lo)) + i_d0 = TRUE; + else + i_d1 = TRUE; + } + + /* check continuity from the top so that the lower-neighbour */ + /* values are properly filled for the upper trapezoid */ + + if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) + { /* continuation of a chain from abv. */ + if (tr[t].usave > 0) /* three upper neighbours */ + { + if (tr[t].uside == S_LEFT) + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = -1; + tr[tn].u1 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[tn].u0].d0 = tn; + tr[tr[tn].u1].d0 = tn; + } + else /* intersects in the right */ + { + tr[tn].u1 = -1; + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[t].u0; + tr[t].u0 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } + + tr[t].usave = tr[tn].usave = 0; + } + else /* No usave.... simple case */ + { + tr[tn].u0 = tr[t].u1; + tr[tn].u1 = -1; + tr[t].u1 = -1; + tr[tr[tn].u0].d0 = tn; + } + } + else + { /* fresh seg. or upward cusp */ + int tmp_u = tr[t].u0; + int td0, td1; + if (((td0 = tr[tmp_u].d0) > 0) && + ((td1 = tr[tmp_u].d1) > 0)) + { /* upward cusp */ + if ((tr[td0].rseg > 0) && + !is_left_of(tr[td0].rseg, seg, &s.v1)) + { + tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d1 = tn; + } + else + { + tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1; + tr[tr[t].u0].d0 = t; + } + } + else /* fresh segment */ + { + tr[tr[t].u0].d0 = t; + tr[tr[t].u0].d1 = tn; + } + } + + if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && + FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) + { + /* this case arises only at the lowest trapezoid.. i.e. + tlast, if the lower endpoint of the segment is + already inserted in the structure */ + + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = -1; + tr[tr[t].d1].u0 = tn; + tr[tr[t].d1].u1 = -1; + + tr[tn].d0 = tr[t].d1; + tr[t].d1 = tr[tn].d1 = -1; + + tnext = tr[t].d1; + } + else if (i_d0) + /* intersecting d0 */ + { + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = tn; + tr[tr[t].d1].u0 = tn; + tr[tr[t].d1].u1 = -1; + + /* new code to determine the bottom neighbours of the */ + /* newly partitioned trapezoid */ + + tr[t].d1 = -1; + + tnext = tr[t].d0; + } + else /* intersecting d1 */ + { + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = -1; + tr[tr[t].d1].u0 = t; + tr[tr[t].d1].u1 = tn; + + /* new code to determine the bottom neighbours of the */ + /* newly partitioned trapezoid */ + + tr[tn].d0 = tr[t].d1; + tr[tn].d1 = -1; + + tnext = tr[t].d1; + } + + t = tnext; + } + + tr[t_sav].rseg = tr[tn_sav].lseg = segnum; + } /* end-while */ + + /* Now combine those trapezoids which share common segments. We can */ + /* use the pointers to the parent to connect these together. This */ + /* works only because all these new trapezoids have been formed */ + /* due to splitting by the segment, and hence have only one parent */ + + tfirstl = tfirst; + tlastl = tlast; + merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT, tr, qs); + merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT, tr, qs); + + seg[segnum].is_inserted = TRUE; + return 0; +} + +/* Update the roots stored for each of the endpoints of the segment. + * This is done to speed up the location-query for the endpoint when + * the segment is inserted into the trapezoidation subsequently + */ +static void +find_new_roots(int segnum, segment_t* seg, trap_t* tr, qnode_t* qs) +{ + segment_t *s = &seg[segnum]; + + if (s->is_inserted) return; + + s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0, seg, qs); + s->root0 = tr[s->root0].sink; + + s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1, seg, qs); + s->root1 = tr[s->root1].sink; +} + +/* Get log*n for given n */ +static int math_logstar_n(int n) +{ + register int i; + double v; + + for (i = 0, v = (double) n; v >= 1; i++) + v = log2(v); + + return (i - 1); +} + +static int math_N(int n, int h) +{ + register int i; + double v; + + for (i = 0, v = (int) n; i < h; i++) + v = log2(v); + + return (int) ceil((double) 1.0*n/v); +} + +/* Main routine to perform trapezoidation */ +int +construct_trapezoids(int nseg, segment_t* seg, int* permute, int ntraps, + trap_t* tr) +{ + int i; + int root, h; + int segi = 1; + qnode_t* qs; + + QSIZE = 2*ntraps; + TRSIZE = ntraps; + qs = N_NEW (2*ntraps, qnode_t); + q_idx = tr_idx = 1; + memset((void *)tr, 0, ntraps*sizeof(trap_t)); + + /* Add the first segment and get the query structure and trapezoid */ + /* list initialised */ + + root = init_query_structure(permute[segi++], seg, tr, qs); + + for (i = 1; i <= nseg; i++) + seg[i].root0 = seg[i].root1 = root; + + for (h = 1; h <= math_logstar_n(nseg); h++) { + for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++) + add_segment(permute[segi++], seg, tr, qs); + + /* Find a new root for each of the segment endpoints */ + for (i = 1; i <= nseg; i++) + find_new_roots(i, seg, tr, qs); + } + + for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++) + add_segment(permute[segi++], seg, tr, qs); + + free (qs); + return tr_idx; +} +