gb_plane.w

模拟器提供了一个简单易用的平台

in the worst case. So our program will allocate blocks of nodes dynamically
instead of assuming a maximum size.

@d nodes_per_block 127 /* on most computers we want it $\equiv 15$ (mod 16) */
@d new_node(x)
  if (next_node==max_node) {
    x=gb_typed_alloc(nodes_per_block,node,working_storage);
    if (x==NULL) {
      gb_free(working_storage); /* release |delaunay|'s auxiliary memory */
      return; /* |gb_trouble_code| is nonzero */
    }
    next_node=x+1; max_node=x+nodes_per_block;
  }@+else x=next_node++;
@#
@d terminal_node(x,p) {@+new_node(x); /* allocate a new node */
   x->v=(Vertex*)(p); /* make it point to a given arc from the triangle */
}  /* note that |x->u==NULL|, representing a terminal node */

@<Local variables for |delaunay|@>=
node *next_node; /* the first yet-unused node slot
   in the current block of nodes */
node *max_node; /* address of nonexistent node following the current
   block of nodes */
node root_node; /* start here to locate a vertex in its triangle */
Area working_storage; /* where |delaunay| builds its triangulation */

@ The algorithm begins with a trivial triangulation that contains
only the first two vertices, together with two ``triangles'' extending
to infinity at their left and right.

@<Initialize the data structures@>=
next_node=max_node=NULL;
init_area(working_storage);
@<Initialize the array of arcs@>;
u=g->vertices;
v=u+1;
@<Make two ``triangles'' for |u|, |v|, and $\infty$@>;

@ We'll need a bunch of local variables to do elementary operations on
data structures.

@<Local variables for |delaunay|@>=
Vertex *p, *q, *r, *s, *t, *tp, *tpp, *u, *v;
arc *a,*aa,*b,*c,*d, *e;
node *x,*y,*yp,*ypp;

@ @<Make two ``triangles'' for |u|, |v|, and $\infty$@>=
root_node.u=u; root_node.v=v;
a=next_arc;
terminal_node(x,a+1);
root_node.l=x;
a->vert=v;@+a->next=a+1;@+a->inst=x;
(a+1)->next=a+2;@+(a+1)->inst=x;
 /* |(a+1)->vert=NULL|, representing $\infty$ */
(a+2)->vert=u;@+(a+2)->next=a;@+(a+2)->inst=x;
mate(a,b);
terminal_node(x,b-2);
root_node.r=x;
b->vert=u;@+b->next=b-2;@+b->inst=x;
(b-2)->next=b-1;@+(b-2)->inst=x;
 /* |(b-2)->vert=NULL|, representing $\infty$ */
(b-1)->vert=v;@+(b-1)->next=b;@+(b-1)->inst=x;
next_arc+=3;

@*Delaunay updating.
The main loop of the algorithm updates the data structure incrementally
by adding one new vertex at a time. The new vertex will always be connected
by an edge (i.e., by two arcs) to each of the vertices of the triangle that
previously enclosed it. It might also deserve to be connected to other
nearby vertices.

@<Find the Delaunay triangulation...@>=
if (g->n<2) return; /* no edges unless there are at least 2 vertices */
@<Initialize the data structures@>;
for (p=g->vertices+2;p<g->vertices+g->n;p++) {
  @<Find an arc |a| on the boundary of the triangle containing |p|@>;
  @<Divide the triangle left of |a| into three triangles surrounding |p|@>;
  @<Explore the triangles surrounding |p|, ``flipping'' their neighbors
    until all triangles that should touch |p| are found@>;
}

@ We have set up the branch nodes so that they solve the triangle location
problem.

@<Find an arc |a| on the boundary of the triangle containing |p|@>=
x=&root_node;
do@+{
  if (ccw(x->u,x->v,p))
    x = x->l;
  else x = x->r;
}@+while (x->u);
a = (arc*) x->v; /* terminal node points to the arc we want */

@ Subdividing a triangle is an easy exercise in data structure manipulation,
except that we must do something special when one of the vertices is
infinite. Let's look carefully at what needs to be done.

Suppose the triangle containing |p| has the vertices |q|, |r|, and |s|
in counterclockwise order. Let |x| be the terminal node that points to
the triangle~$\Delta qrs$. We want to change |x| so that we will be
able to locate a future point of $\Delta qrs$ within either $\Delta pqr$,
$\Delta prs$, or $\Delta psq$.

If |q|, |r|, and |s| are finite, we will change |x| and add five new nodes
as follows:
$$\vbox{\halign{\hfil#:\enspace&#\hfil\cr
$x$&if left of $rp$, go to $x''$, else go to $x'$;\cr
$x'$&if left of $sp$, go to $y$, else go to $y'$;\cr
$x''$&if left of $qp$, go to $y'$, else go to $y''$;\cr
$y$&you're in $\Delta prs$;\cr
$y'$&you're in $\Delta psq$;\cr
$y''$&you're in $\Delta pqr$.\cr}}$$

But if, say, $q=\infty$, such instructions make no sense,
because there are lines in all directions that run from $\infty$ to any point.
In such a case we use ``wedges'' instead of triangles, as explained below.

At the beginning of the following code, we have |x==a->inst|.

@<Divide the triangle left of |a| into three triangles surrounding |p|@>=
b=a->next;@+c=b->next;
q=a->vert;@+r=b->vert;@+s=c->vert;
@<Create new terminal nodes |y|, |yp|, |ypp|, and new arcs pointing to them@>;
if (q==NULL) @<Compile instructions to update convex hull@>
else {@+register node *xp;
  x->u=r;@+x->v=p;
  new_node(xp);
  xp->u=q;@+xp->v=p;@+xp->l=yp;@+xp->r=ypp; /* instruction $x''$ above */
  x->l=xp;
  new_node(xp);
  xp->u=s;@+xp->v=p;@+xp->l=y;@+xp->r=yp; /* instruction $x'$ above */
  x->r=xp;
}

@ The only subtle point here is that |q=a->vert| might be |NULL|. A terminal
node must point to the proper arc of an infinite triangle.

@<Create new terminal nodes |y|, |yp|, |ypp|, and new arcs pointing to them@>=
terminal_node(yp,a);@+terminal_node(ypp,next_arc);@+terminal_node(y,c);
c->inst=y;@+a->inst=yp;@+b->inst=ypp;
mate(next_arc,e);
a->next=e;@+b->next=e-1;@+c->next=e-2;
next_arc->vert=q;@+next_arc->next=b;@+next_arc->inst=ypp;
(next_arc+1)->vert=r;@+(next_arc+1)->next=c;@+(next_arc+1)->inst=y;
(next_arc+2)->vert=s;@+(next_arc+2)->next=a;@+(next_arc+2)->inst=yp;
e->vert=(e-1)->vert=(e-2)->vert=p;
e->next=next_arc+2;@+(e-1)->next=next_arc;@+(e-2)->next=next_arc+1;
e->inst=yp;@+(e-1)->inst=ypp;@+(e-2)->inst=y;
next_arc += 3;

@ Outside of the current convex hull, we have ``wedges'' instead of
triangles. Wedges are exterior angles whose points lie outside an
edge $rs$ of the convex hull, but not outside the next edge on the other
side of point |r|. When a new point lies in such a wedge, we have to
see if it also lies outside the edges $st$, $tu$, etc., in the
clockwise direction, in which case the convex hull loses points
$s$, $t$, etc., and we must update the new wedges accordingly.

This was the hardest part of the program to prove correct; a complete
proof can be found in {\sl Axioms and Hulls}.

@<Compile...@>=
{@+register node *xp;
  x->u=r;@+x->v=p;@+x->l=ypp;
  new_node(xp);
  xp->u=s;@+xp->v=p;@+xp->l=y;@+xp->r=yp;
  x->r=xp;
  mate(a,aa);@+d=aa->next;@+t=d->vert;
  while (t!=r && (ccw(p,s,t))) {@+register node *xpp;
    terminal_node(xpp,d);
    xp->r=d->inst;
    xp=d->inst;
    xp->u=t;@+xp->v=p;@+xp->l=xpp;@+xp->r=yp;
    flip(a,aa,d,s,NULL,t,p,xpp,yp);
    a=aa->next;@+mate(a,aa);@+d=aa->next;
    s=t;@+t=d->vert;
    yp->v=(Vertex*)a;
  }
  terminal_node(xp,d->next);
  x=d->inst;@+x->u=s;@+x->v=p;@+x->l=xp;@+x->r=yp;
  d->inst=xp;@+d->next->inst=xp;@+d->next->next->inst=xp;
  r=s; /* this value of |r| shortens the exploration step that follows */
}

@ The updating process finishes by walking around the triangles
that surround |p|, making sure that none of them are adjacent to
triangles containing |p| in their circumcircle. (Such triangles are
no longer in the Delaunay triangulation, by definition.)

@<Explore...@>=
while(1) {
  mate(c,d);@+e=d->next;
  t=d->vert;@+tp=c->vert;@+tpp=e->vert;
  if (tpp && incircle(tpp,tp,t,p)) { /* triangle $tt''t'$ no longer Delaunay */
    register node *xp, *xpp;
    terminal_node(xp,e);
    terminal_node(xpp,d);
    x=c->inst;@+x->u=tpp;@+x->v=p;@+x->l=xp;@+x->r=xpp;
    x=d->inst;@+x->u=tpp;@+x->v=p;@+x->l=xp;@+x->r=xpp;
    flip(c,d,e,t,tp,tpp,p,xp,xpp);
    c=e;
  }
  else if (tp==r) break;
  else {
    mate(c->next,aa);
    c=aa->next;
  }
}

@ Here |d| is the mate of |c|, |e=d->next|, |t=d->vert|, |tp=c->vert|,
and |tpp=e->vert|. The triangles $\Delta tt'p$ and $\Delta t'tt''$ to the
left and right of arc~|c| are being replaced in the current triangulation
by $\Delta ptt''$ and $\Delta t''t'p$, corresponding to terminal nodes
|xp| and |xpp|. (The values of |t| and |tp| are not actually used, so
some optimization is possible.)

@<Other...@>=
static void flip(c,d,e,t,tp,tpp,p,xp,xpp)
  arc *c,*d,*e;
  Vertex *t,*tp,*tpp,*p;
  node *xp,*xpp;
{@+register arc *ep=e->next, *cp=c->next, *cpp=cp->next;
  e->next=c;@+c->next=cpp;@+cpp->next=e;
  e->inst=c->inst=cpp->inst=xp;
  c->vert=p;
  d->next=ep;@+ep->next=cp;@+cp->next=d;
  d->inst=ep->inst=cp->inst=xpp;
  d->vert=tpp;
}

@*Use of mileage data. The |delaunay| routine is now complete, and the
only missing piece of code is the promised routine that generates
planar graphs based on data from the real world.

The subroutine call {\advance\thinmuskip 0mu plus 2mu
|plane_miles(n,north_weight,west_weight,pop_weight, extend,prob,seed)|}
will construct a planar graph with min$(128,n)$ vertices, where the
vertices are exactly the same as the cities produced by the subroutine call
|miles(n,north_weight,west_weight, pop_weight,0,0,seed)|. (As
explained in module {\sc GB\_\,MILES}, the weight parameters |north_weight|,
|west_weight|, and |pop_weight| are used to rank the cities by
location and/or population.)  The edges of the new graph are obtained
by first constructing the Delaunay triangulation of those cities,
based on a simple projection onto the plane using their latitude and
longitude, then discarding each Delaunay edge with probability
|prob/65536|. The length of each surviving edge is the same as the
mileage between cities that would appear in the complete graph
produced by |miles|.

If |extend!=0|, an additional vertex representing $\infty$ is also
included. The Delaunay triangulation includes edges of length |INFTY|
connecting this vertex with all cities on the convex hull; these edges,
like the others, are subject to being discarded with probability |prob/65536|.
(See the description of |plane| for further comments about using
|prob| to control the sparseness of the graph.)

The weight parameters must satisfy
$$ \vert|north_weight|\vert\le100{,}000,\quad
   \vert|west_weight|\vert\le100{,}000,\quad
   \vert|pop_weight|\vert\le100.$$
Vertices of the graph will appear in order of decreasing weight.
The |seed| parameter defines the pseudo-random numbers used wherever
a ``random'' choice between equal-weight vertices needs to be made,
or when deciding whether to discard a Delaunay edge.

@<The |plane_miles| routine@>=
Graph *plane_miles(n,north_weight,west_weight,pop_weight,extend,prob,seed)
  unsigned long n; /* number of vertices desired */
  long north_weight; /* coefficient of latitude in the weight function */
  long west_weight; /* coefficient of longitude in the weight function */
  long pop_weight; /* coefficient of population in the weight function */
  unsigned long extend; /* should a point at infinity be included? */
  unsigned long prob; /* probability of rejecting a Delaunay edge */
  long seed; /* random number seed */
{@+Graph *new_graph; /* the graph constructed by |plane_miles| */
  @<Use |miles| to set up the vertices of a graph@>;
  @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with the road length in miles@>;
  if (gb_trouble_code) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  }
  gb_free(new_graph->aux_data); /* recycle special memory used by |miles| */
  if (extend) new_graph->n++; /* make the ``infinite'' vertex legitimate */
  return new_graph;
}

@ By setting the |max_distance| parameter to~1, we cause |miles|
to produce a graph having the desired vertices but no edges.
The vertices of this graph will have appropriate coordinate fields
|x_coord|, |y_coord|, and~|z_coord|.

@<Use |miles|...@>=
if (extend) extra_n++; /* allocate one more vertex than usual */
if (n==0 || n>MAX_N) n=MAX_N; /* compute true number of vertices */
new_graph=miles(n,north_weight,west_weight,pop_weight,1L,0L,seed);
if (new_graph==NULL) return NULL; /* |panic_code| has been set by |miles| */
sprintf(new_graph->id,"plane_miles(%lu,%ld,%ld,%ld,%lu,%lu,%ld)",
  n,north_weight,west_weight,pop_weight,extend,prob,seed);
if (extend) extra_n--; /* restore |extra_n| to its previous value */

@ @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with the road length in miles@>=
gprob=prob;
if (extend) {
  inf_vertex=new_graph->vertices+new_graph->n;
  inf_vertex->name=gb_save_string("INF");
  inf_vertex->x_coord=inf_vertex->y_coord=inf_vertex->z_coord= -1;
}@+else inf_vertex=NULL;
delaunay(new_graph,new_mile_edge);

@ The mileages will all have been negated by |miles|, so we make them
positive again.

@<Other...@>=
static void new_mile_edge(u,v)
  Vertex *u,*v;
{
  if ((gb_next_rand()>>15)>=gprob) {
    if (u) {
      if (v) gb_new_edge(u,v,-miles_distance(u,v));
      else if (inf_vertex) gb_new_edge(u,inf_vertex,INFTY);
    }@+else if (inf_vertex) gb_new_edge(inf_vertex,v,INFTY);
  }
}

@* Index. As usual, we close with an index that
shows where the identifiers of \\{gb\_plane} are defined and used.