gb_plane.w

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

% This file is part of the Stanford GraphBase (c) Stanford University 1993
@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
@i gb_types.w

\def\title{GB\_\,PLANE}

\prerequisite{GB\_\,MILES}
@* Introduction. This GraphBase module contains the |plane| subroutine,
which constructs undirected planar graphs from vertices located randomly
in a rectangle,
as well as the |plane_miles| routine, which constructs planar graphs
based on the mileage and coordinate data in \.{miles.dat}. Both
routines use a general-purpose |delaunay| subroutine,
which computes the Delaunay triangulation of a given set of points.

@d plane_miles p_miles /* abbreviation for Procrustean external linkage */

@(gb_plane.h@>=
#define plane_miles p_miles
extern Graph *plane();
extern Graph *plane_miles();
extern void delaunay();

@ The subroutine call |plane(n,x_range,y_range,extend,prob,seed)| constructs
a planar graph whose vertices have integer coordinates
uniformly distributed in the rectangle
$$\{\,(x,y)\;\mid\;0\le x<|x_range|, \;0\le y<|y_range|\,\}\,.$$
The values of |x_range| and |y_range| must be at most $2^{14}=16384$; the
latter value is the default, which is substituted if |x_range| or |y_range|
is given as zero. If |extend==0|, the graph will have |n| vertices; otherwise
it will have |n+1| vertices, where the |(n+1)|st is assigned the coordinates
$(-1,-1)$ and may be regarded as a point at~$\infty$.
Some of the |n|~finite vertices might have identical coordinates, particularly
if the point density |n/(x_range*y_range)| is not very small.

The subroutine works by first constructing the Delaunay triangulation
of the points, then discarding
each edge of the resulting graph with probability |prob/65536|. Thus,
for example, if |prob| is zero the full Delaunay triangulation will be
returned; if |prob==32768|, about half of the Delaunay edges will remain.
Each finite edge is assigned a length equal to the Euclidean distance between
points, multiplied by $2^{10}$ and
rounded to the nearest integer. If |extend!=0|, the
Delaunay triangulation will also contain edges between $\infty$ and
all points of the convex hull; such edges, if not discarded, are
assigned length $2^{28}$, otherwise known as |INFTY|.

If |extend!=0| and |prob==0|, the graph will have $n+1$ vertices and
$3(n-1)$ edges; this is the maximum number of edges that a planar graph
on $n+1$ vertices can have. In such a case the average degree of a vertex will
be $6(n-1)/(n+1)$, slightly less than~6; hence, if |prob==32768|,
the average degree of a vertex will usually be near~3.

As with all other GraphBase routines that rely on random numbers,
different values of |seed| will produce different graphs, in a
machine-independent fashion that is reproducible on many different
computers. Any |seed| value between 0 and $2^{31}-1$ is permissible.

@d INFTY 0x10000000L /* ``infinite'' length */

@(gb_plane.h@>=
#define INFTY @t\quad@> 0x10000000L

@ If the |plane| routine encounters a problem, it returns |NULL|
(\.{NULL}), after putting a code number into the external variable
|panic_code|. This code number identifies the type of failure.
Otherwise |plane| returns a pointer to the newly created graph, which
will be represented with the data structures explained in {\sc GB\_\,GRAPH}.
(The external variable |panic_code| is itself defined in {\sc GB\_\,GRAPH}.)

@d panic(c) @+{@+panic_code=c;@+gb_trouble_code=0;@+return NULL;@+}

@ Here is the overall shape of the \CEE/ file \.{gb\_plane.c}\kern.2em:

@p
#include "gb_flip.h"
 /* we will use the {\sc GB\_\,FLIP} routines for random numbers */
#include "gb_graph.h" /* we will use the {\sc GB\_\,GRAPH} data structures */
#include "gb_miles.h" /* and we might use {\sc GB\_\,MILES} for mileage data */
#include "gb_io.h"
 /* and {\sc GB\_\,MILES} uses {\sc GB\_\,IO}, which has |str_buf| */
@h@#
@<Type declarations@>@;
@<Global variables@>@;
@<Subroutines for arithmetic@>@;
@<Other subroutines@>@;
@<The |delaunay| routine@>@;
@<The |plane| routine@>@;
@<The |plane_miles| routine@>@;

@ @<The |plane| routine@>=
Graph *plane(n,x_range,y_range,extend,prob,seed)
  unsigned long n; /* number of vertices desired */
  unsigned long x_range,y_range; /* upper bounds on rectangular coordinates */
  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| */
  register Vertex *v; /* the current vertex of interest */
  register long k; /* the canonical all-purpose index */
  gb_init_rand(seed);
  if (x_range>16384 || y_range>16384) panic(bad_specs); /* range too large */
  if (n<2) panic(very_bad_specs); /* don't make |n| so small, you fool */
  if (x_range==0) x_range=16384; /* default */
  if (y_range==0) y_range=16384; /* default */
  @<Set up a graph with |n| uniformly distributed vertices@>;
  @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with their Euclidean length@>;
  if (gb_trouble_code) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  }
  if (extend) new_graph->n++; /* make the ``infinite'' vertex legitimate */
  return new_graph;
}

@ The coordinates are placed into utility fields |x_coord| and |y_coord|.
A random ID number is also stored in utility field~|z_coord|; this number is
used by the |delaunay| subroutine to break ties when points are equal or
collinear or cocircular. No two vertices have the same ID number.
(The header file \.{gb\_miles.h} defines |x_coord|, |y_coord|, and
|index_no| to be |x.I|, |y.I|, and |z.I| respectively.)

@d z_coord z.I

@<Set up a graph with |n| uniform...@>=
if (extend) extra_n++; /* allocate one more vertex than usual */
new_graph=gb_new_graph(n);
if (new_graph==NULL)
  panic(no_room); /* out of memory before we're even started */
sprintf(new_graph->id,"plane(%lu,%lu,%lu,%lu,%lu,%ld)",
  n,x_range,y_range,extend,prob,seed);
strcpy(new_graph->util_types,"ZZZIIIZZZZZZZZ");
for (k=0,v=new_graph->vertices; k<n; k++,v++) {
  v->x_coord=gb_unif_rand(x_range);
  v->y_coord=gb_unif_rand(y_range);
  v->z_coord=((long)(gb_next_rand()/n))*n+k;
  sprintf(str_buf,"%ld",k);@+v->name=gb_save_string(str_buf);
}
if (extend) {
  v->name=gb_save_string("INF");
  v->x_coord=v->y_coord=v->z_coord=-1;
  extra_n--;
}

@ @(gb_plane.h@>=
#define x_coord @t\quad@> x.I
#define y_coord @t\quad@> y.I
#define z_coord @t\quad@> z.I

@* Delaunay triangulation. The Delaunay triangulation of a set of
vertices in the plane consists of all line segments $uv$ such that
there exists a circle passing through $u$ and~$v$ containing no other
vertices. Equivalently, $uv$ is a Delaunay edge if and only if the
Voronoi regions for $u$ and $v$ are adjacent; the Voronoi region of a
vertex~$u$ is the polygon with the property that all points inside it
are closer to $u$ than to any other vertex. In this sense, we can say
that Delaunay edges connect vertices with their ``neighbors.''
@^Delaunay [Delone], Boris Nikolaevich@>

The definitions in the previous paragraph assume that no two vertices are
equal, that no three vertices lie on a straight line, and that no four vertices
lie on a circle. If those nondegeneracy conditions aren't satisfied, we can
perturb the points very slightly so that the assumptions do hold.

Another way to characterize the Delaunay triangulation is to consider
what happens when we map a given set of
points onto the unit sphere via stereographic projection: Point $(x,y)$ is
mapped to
$$\bigl(2x/(r^2+1),2y/(r^2+1),(r^2-1)/(r^2+1)\bigr)\,,$$ where $r^2=x^2+y^2$.
If we now extend the configuration by adding $(0,0,1)$,
which is the limiting point on the sphere when $r$ approaches infinity,
the Delaunay edges of the original points
turn out to be edges of the polytope defined by the mapped
points. This polytope, which is the 3-dimensional convex hull of $n+1$ points
on the sphere, also has edges from $(0,0,1)$ to the mapped points
that correspond to the 2-dimensional convex hull of the original points. Under
our assumption of nondegeneracy, the faces of this polytope are all
triangles; hence its edges are said to form a triangulation.

A self-contained presentation of all the relevant theory, together with
an exposition and proof of correctness of the algorithm below, can be found
in the author's monograph {\sl Axioms and Hulls}, Lecture Notes in
Computer Science {\bf606} (Springer-Verlag, 1992).
@^Axioms and Hulls@>
For further references, see Franz Aurenhammer, {\sl ACM Computing Surveys\/
@^Aurenhammer, Franz@>
\bf23} (1991), 345--405.

@ The |delaunay| procedure, which finds the Delaunay triangulation of
a given set of vertices, is the key ingredient in \\{gb\_plane}'s
algorithms for generating planar graphs. The given vertices should
appear in a GraphBase graph~|g| whose edges, if any, are ignored by
|delaunay|. The coordinates of each vertex appear in utility fields
|x_coord| and~|y_coord|, which must be nonnegative and less than
$2^{14}=16384$. The utility field~|z_coord| must contain a unique ID
number, distinct for every vertex, so that the algorithm can break
ties in cases of degeneracy. (Note: These assumptions about the input
data are the responsibility of the calling procedure; |delaunay| does not
double-check them. If they are violated, catastrophic failure is possible.)

Instead of returning the Delaunay triangulation as a graph, |delaunay|
communicates its answer implicitly by performing the procedure call
|f(u,v)| on every pair of vertices |u| and~|v| joined by a Delaunay edge.
Here |f|~is a procedure supplied as a parameter; |u| and~|v| are either
pointers to vertices or |NULL| (i.e., \.{NULL}), where |NULL| denotes the
vertex ``$\infty$.'' As remarked above, edges run between $\infty$ and all
vertices on the convex hull of the given points. The graph of all edges,
including the infinite edges, is planar.

For example, if the vertex at infinity is being ignored, the user can
declare
$$\vcenter{\halign{#\hfil\cr
|void ins_finite(u,v)|\cr
\qquad|Vertex *u,*v;|\cr
|{@+if (u&&v)@+gb_new_edge(u,v,1L);@+}|\cr}}$$
Then the procedure call |delaunay(g,ins_finite)| will add all the finite
Delaunay edges to the current graph~|g|, giving them all length~1.

If |delaunay| is unable to allocate enough storage to do its work, it
will set |gb_trouble_code| nonzero and there will be no edges in
the triangulation.

@<The |delaunay| routine@>=
void delaunay(g,f)
  Graph *g; /* vertices in the plane */
  void @[@] (*f)(); /* procedure that absorbs the triangulated edges */
{@+@<Local variables for |delaunay|@>;@#
  @<Find the Delaunay triangulation of |g|, or return with |gb_trouble_code|
    nonzero if out of memory@>;
  @<Call |f(u,v)| for each Delaunay edge \\{uv}@>;
  gb_free(working_storage);
}

@ The procedure passed to |delaunay| will communicate with |plane| via
global variables called |gprob| and |inf_vertex|.

@<Glob...@>=
static unsigned long gprob; /* copy of the |prob| parameter */
static Vertex *inf_vertex; /* pointer to the vertex $\infty$, or |NULL| */

@ @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with their Euclidean length@>=
gprob=prob;
if (extend) inf_vertex=new_graph->vertices+n;
else inf_vertex=NULL;
delaunay(new_graph,new_euclid_edge);

@ @<Other...@>=
static void new_euclid_edge(u,v)
  Vertex *u,*v;
{@+register long dx,dy;
  if ((gb_next_rand()>>15)>=gprob) {
    if (u) {
      if (v) {
        dx=u->x_coord-v->x_coord;
        dy=u->y_coord-v->y_coord;
        gb_new_edge(u,v,int_sqrt(dx*dx+dy*dy));
      }@+else if (inf_vertex) gb_new_edge(u,inf_vertex,INFTY);
    }@+else if (inf_vertex) gb_new_edge(inf_vertex,v,INFTY);
  }
}

@* Arithmetic. Before we lunge into the world of geometric algorithms,
let's build up some confidence by polishing off some subroutines that
will be needed to ensure correct results. We assume that |long| integers
are less than $2^{31}$.

First is a routine to calculate $s=\lfloor2^{10}\sqrt x+{1\over2}\rfloor$,
the nearest integer to $2^{10}$ times the square root of a given nonnegative
integer~|x|. If |x>0|, this
is the unique integer such that $2^{20}x-s\le s^2<2^{20}x+s$.

The following routine appears to work by magic, but the mystery goes
away when one considers the invariant relations
$$ m=\lfloor 2^{2k-21}\rfloor,\qquad
   0<y=\lfloor 2^{20-2k}x\rfloor-s^2+s\le q=2s.$$
(Exception: We might actually have $y=0$ for a short time when |q=2|.)

@<Subroutines for arith...@>=
static long int_sqrt(x)
  long x;
{@+register long y, m, q=2; long k;
  if (x<=0) return 0;
  for (k=25,m=0x20000000;x<m;k--,m>>=2) ; /* find the range */
  if (x>=m+m) y=1;
  else y=0;
  do @<Decrease |k| by 1, maintaining the invariant relations
       between |x|, |y|, |m|, and |q|@>@;
  while (k);
  return q>>1;
}

@ @<Decrease |k| by 1, maintaining the invariant relations...@>=
{
  if (x&m) y+=y+1;
  else y+=y;
  m>>=1;
  if (x&m) y+=y-q+1;
  else y+=y-q;
  q+=q;
  if (y>q)
    y-=q,q+=2;
  else if (y<=0)
    q-=2,y+=q;
  m>>=1;
  k--;
}

@ We are going to need multiple-precision arithmetic in order to
calculate certain geometric predicates properly, but it turns out
that we do not need to implement general-purpose subroutines for bignums.
It suffices to have a single special-purpose routine called
$|sign_test|(|x1|,|x2|,|x3|,\allowbreak|y1|,|y2|,|y3|)$, which
computes a single-precision integer having the same sign as the dot product
$$\hbox{|x1*y1+x2*y2+x3*y3|}$$
when we have $-2^{29}<|x1|,|x2|,|x3|<2^{29}$ and $0\le|y1|,|y2|,|y3|<2^{29}$.

@<Subroutines for arith...@>=
static long sign_test(x1,x2,x3,y1,y2,y3)
  long x1,x2,x3,y1,y2,y3;
{@+long s1,s2,s3; /* signs of individual terms */
  long a,b,c; /* components of a redundant representation of the dot product */
  register long t; /* temporary register for swapping */
  @<Determine the signs of the terms@>;
  @<If the answer is obvious, return it without further ado; otherwise,
    arrange things so that |x3*y3| has the opposite sign to |x1*y1+x2*y2|@>;
  @<Compute a redundant representation of |x1*y1+x2*y2+x3*y3|@>;
  @<Return the sign of the redundant representation@>;
}