% This file is part of the Stanford GraphBase (c) Stanford University 1990 \nocon @* Introduction. This is a hastily written implementation of treehull, using treaps to guarantee good average access time. @f Graph int /* |gb_graph| defines the |Graph| type and a few others */ @f Vertex int @f Arc int @f Area int @p #include "gb_graph.h" #include "gb_miles.h" #include "gb_rand.h" @@; int n=128; @@; @@; @# main() { @@; Graph *g=miles(128,0,0,0,0,0,0); @# mems=ccs=0; @; printf("Total of %d mems and %d calls on ccw.\n",mems,ccs); } @ I'm instrumenting this in a simple way. @d o mems++ @d oo mems+=2 @d ooo mems+=3 @= int mems; /* memory accesses */ int ccs; /* calls on |ccw| */ @*Data structures. For now, each vertex is represented by two coordinates stored in the utility fields |x.I| and |y.I|. I'm also putting a serial number into |z.I|, so that I can check whether different algorithms generate identical hulls. We use separate nodes for the current convex hull. These nodes have a bunch of fields: |p->vert| points to the vertex; |p->succ| and |p->pred| point to next and previous nodes in a circular list; |p->left| and |p->right| point to left and right children in a tree that's superimposed on the list; |p->parent| is present too, it points to the parent node; |p->prio| is the priority if we are implementing the tree as a treap. The |head| node has the root of the tree in its |right| field, and it represents the special vertex that isn't in the tree. @= typedef struct node_struct { struct vertex_struct *vert; struct node_struct *succ,*pred,*left,*right,*parent; long prio; } node; @ @= head=(node*)gb_alloc((g->n)*sizeof(node),working_storage); if (head==NULL) return(1); /* fixthis */ next_node=head; @ @= node *head; /* beginning of the hull data structure */ node *next_node; /* first unused slot in that array */ Area working_storage; int serial_no=1; /* used to disambiguate entries with equal coordinates */ @ We assume that the vertices have been given to us in a GraphBase-type graph. The algorithm begins with a trivial hull that contains only the first two vertices. @= init_area(working_storage); @; o,u=g->vertices; v=u+1; u->z.I=0; v->z.I=1; p=++next_node; ooo,head->succ=head->pred=head->right=p; oo,p->succ=p->pred=head; o,p->parent=head; oo,p->left=p->right=NULL; gb_init_rand(110); o,p->prio=gb_next_rand(); o,head->vert=u; o,p->vert=v; next_node++; if (n<150) printf("Beginning with (%s; %s)\n",u->name,v->name); @ We'll probably need a bunch of local variables to do elementary operations on data structures. @= Vertex *u,*v,*vv,*w; node *p,*pp,*q, *qq, *qqq, *r, *rr, *s, *ss, *tt, **par, **ppar, *prepar, *preppar; int replaced; /* will be nonzero if we've just replaced a hull element */ @ Here's a routine I used when debugging (in fact I should have written it sooner than I did). @= p=head; count=0; do { count++; p->prio=(p->prio&0xffff0000)+count; if (p->succ->pred!=p) printf("succ/pred failure at %s!\n",p->vert->name); if (p->left!=NULL && p->left->parent!=p) printf("parent/lchild failure at %s!\n",p->vert->name); if (p->right!=NULL && p->right->parent!=p) printf("parent/rchild failure at %s!\n",p->vert->name); p=p->succ; } while (p!=head); count=1; inorder(head->right); @ @= int count; inorder(p) node *p; { if (p) { inorder(p->left); if ((p->prio&0xffff)!=++count) { printf("tree node %d is missing at %d: %s!\n",count,p->prio&0xffff, p->vert->name); count=p->prio&0xffff; } inorder(p->right); } } @*Hull updating. The main loop of the algorithm updates the data structure incrementally by adding one new vertex at a time. If the new vertex lies outside the current convex hull, we put it into the cycle and possibly delete some vertices that were previously part of the hull. @= @; for (oo,vv=g->vertices+2;vvvertices+g->n;vv++) { vv->z.I=++serial_no; o,q=head->pred; replaced=0; o,u=head->vert; if (o,ccw(vv,u,q->vert)) @@; else @; @; } @; @ Let me do the easy part first, since it's bedtime and I can worry about the rest tomorrow. @= p=head; printf("The convex hull is:\n"); do { printf(" %s\n",p->vert->name); p=p->succ; } while (p!=head); @ In Case 1 we don't need the tree structure since we've already found that the new vertex is outside the hull at the tree root position. @= {@+qqq=head; while (1) { o,r=qqq->succ; if (r==q) break; /* can't eliminate any more */ if (oo,ccw(vv,qqq->vert,r->vert)) break; @; qqq=r; } qq=qqq; qqq=q; while (1) { o,r=qqq->pred; if (r==qq) break; if (oo,ccw(vv,r->vert,qqq->vert)) break; @; qqq=r; } q=qqq; if (!replaced) @; if (n<150) printf("New hull sequence (%s; %s; %s)\n", q->vert->name,vv->name,qq->vert->name); } @ At this point |q==head->pred| is the tree's rightmost node. @= { tt=next_node++; o,tt->vert=vv; o,tt->succ=head; o,tt->pred=q; o,head->pred=tt; o,q->succ=tt; oo,tt->left=tt->right=NULL; o,tt->prio=gb_next_rand(); if (n<150) printf("(Inserting %s at right of tree, prio=%d)\n",vv->name,tt->prio); if (o,tt->prioprio) rotup(q,&(q->right),tt,tt->prio); else { /* easy case, no rotation necessary */ o,tt->parent=q; o,q->right=tt; } } @ The link from parent to child hasn't been set when the priorities indicate necessary rotation. @= rotup(p,pp,q,qp) node *p; /* parent of inserted node */ node **pp; /* link field in parent */ node *q; /* inserted node */ long qp; /* its priority */ {@+node *pr, **ppr; /* grandparent */ node *qq; /* child who is reparented */ while (1) { o,pr=p->parent; if (o,pr->right==p) ppr=&(pr->right); else ppr=&(pr->left); if (pp==&(p->right)) { /* we should rotate left */ if (n<150) printf("...(rotating left)\n"); o,qq=q->left; o,q->left=p; o,p->parent=q; o,p->right=qq; if (qq!=NULL) o,qq->parent=p; } else { /* we should rotate right */ if (n<150) printf("...(rotating right)\n"); o,qq=q->right; o,q->right=p; o,p->parent=q; o,p->left=qq; if (qq!=NULL) o,qq->parent=p; } if (o,qp>=pr->prio) break; p=pr; pp=ppr; } o,q->parent=pr; o,*ppr=q; } @ Nodes don't need to be recycled. @= if (replaced) { o,pp=qqq->pred; o,tt=qqq->succ; o,pp->succ=tt; o,tt->pred=pp; o,prepar=qqq->parent; if (o,prepar->right==qqq) par=&(prepar->right); else par=&(prepar->left); o,pp=qqq->left; if (o,(ss=qqq->right)==NULL) { if (n<150) printf("(Deleting %s from tree, case 1)\n",qqq->vert->name); o,*par=pp; if (pp!=NULL) o,pp->parent=prepar; } else if (pp==NULL) { if (n<150) printf("(Deleting %s from tree, case 2)\n",qqq->vert->name); o,*par=ss; o,ss->parent=prepar; } else { if (n<150) printf("(Deleting %s from tree, hard case)\n",qqq->vert->name); oo,deldown(prepar,par,pp,ss,pp->prio,ss->prio); } } else { o,qqq->vert=vv; replaced=1; } @ @= deldown(p,pp,ql,qr,qlp,qrp) node *p; /* parent of deleted node */ node **pp; /* link field in that parent */ node *ql, *qr; /* children of deleted node */ int qlp, qrp; /* their priorities */ {@+node *qq; /* grandchild of deleted node */ if (qlpparent=p; o,*pp=ql; o,qq=ql->right; if (qq!=NULL) o,deldown(ql,&(ql->right),qq,qr,qq->prio,qrp); /* tail recursion */ else { o,ql->right=qr; o,qr->parent=ql; } } else { if (n<150) printf("...(moving right child up)\n"); o,qr->parent=p; o,*pp=qr; o,qq=qr->left; if (qq!=NULL) o,deldown(qr,&(qr->left),ql,qq,qlp,qq->prio); /* tail recursion */ else { o,qr->left=ql; o,ql->parent=qr; } } } @ @= {@+o,qq=head->right; while (1) { if (qq==q || (o,ccw(u,vv,qq->vert))) { o,r=qq->left; if (r==NULL) { preppar=qq; o,ppar=&(qq->left); break; } } else { o,r=qq->right; if (r==NULL) { preppar=qq; o,ppar=&(qq->right); o,qq=qq->succ; break; } } qq=r; } if (o,(r=qq->pred)==head || (oo,ccw(vv,qq->vert,r->vert))) { if (r!=head) { while (1) { qqq=r; o,r=qqq->pred; if (r==head) break; if (oo,ccw(vv,r->vert,qqq->vert)) break; @; } r=qqq; } qqq=qq; while (1) { if (qqq==q) break; oo,rr=qqq->succ; if (oo,ccw(vv,qqq->vert,rr->vert)) break; @; qqq=rr; } if (!replaced) @; if (n<150) printf("New hull sequence (%s; %s; %s)\n", r->vert->name,vv->name,qqq->vert->name); } } @ @= { tt=next_node++; o,tt->vert=vv; o,tt->succ=qq; o,tt->pred=r; o,qq->pred=tt; o,r->succ=tt; oo,tt->left=tt->right=NULL; o,tt->prio=gb_next_rand(); if (n<150) printf("(Inserting %s at bottom of tree, prio=%d)\n",vv->name,tt->prio); if (o,tt->prioprio) rotup(preppar,ppar,tt,tt->prio); else { /* easy case, no rotation needed */ o,tt->parent=preppar; o,*ppar=tt; } } @*Determinants. I need code for the primitive function |ccw|. Floating-point arithmetic suffices for my purposes. We want to evaluate the determinant $$ccw(u,v,w)=\left\vert\matrix{u(x)&u(y)&1\cr v(x)&v(y)&1\cr w(x)&w(y)&1\cr} \right\vert=\left\vert\matrix{u(x)-w(x)&u(y)-w(y)\cr v(x)-w(x)&v(y)-w(y)\cr} \right\vert\,.$$ @= int ccw(u,v,w) Vertex *u,*v,*w; {@+register double wx=(double)w->x.I, wy=(double)w->y.I; register double det=((double)u->x.I-wx)*((double)v->y.I-wy) -((double)u->y.I-wy)*((double)v->x.I-wx); Vertex *uu=u,*vv=v,*ww=w,*t; if (det==0) { det=1; if (u->x.I>v->x.I || (u->x.I==v->x.I && (u->y.I>v->y.I || (u->y.I==v->y.I && u->z.I>v->z.I)))) { t=u;@+u=v;@+v=t;@+det=-det; } if (v->x.I>w->x.I || (v->x.I==w->x.I && (v->y.I>w->y.I || (v->y.I==w->y.I && v->z.I>w->z.I)))) { t=v;@+v=w;@+w=t;@+det=-det; } if (u->x.I>v->x.I || (u->x.I==v->x.I && (u->y.I>v->y.I || (u->y.I==v->y.I && u->z.Iz.I)))) { det=-det; } } if (n<150) printf("cc(%s; %s; %s) is %s\n",uu->name,vv->name,ww->name, det>0? "true": "false"); ccs++; return (det>0); }