Conference rusure::math

Depends on what you mean by "random".

As I recall, someone about 40-50 years ago found a way to tile the plane 
with convex 7-gons, except that the further you got from the origin, the 
bigger they get. So the random "centers" would have to be distributed a bit
wildly, perhaps. Certainly not uniformly.

I think it's also not possible to tile a finite plane convex region with
7-gons (or even a combination of 7-gons and hexagons), in which case the
average value of 6 can probably be proven. 

Lynn 

What you've constructed is the Vornoi diagram for the set of points.

A Vornoi diagram for a set of points divides the plane into regions
such that each region contains one of the original n points and all
the points that are closer to that point than any other original point.

Assuming that no four points are co-circular, the number of edges for
a Vornoi diagram for n points is (perhaps) a simple function of n,
and the number of points on the convex hull.

    The following program may be of some interest.  It computes the
    dual of the Voronoi via a divide and conquer algorithm that was
    published a while back by Leo Guibas in the ACM TOG, and shows the
    progress of the triangulation.  It would be very easy to show the
    Voronoi too by clipping the bisector of each triangular edge against
    neighboring bisectors (and the graphics window bounds.)

    It uses UIS, but the graphics calls can easily be changed to output
    REGIS, GKS, XLIB or other graphics protocols.

    - Jim

/*
 *  Experiment with Voronoi diagrams and Delaunay triangulations in the plane
 *
 *  J. Roth	CADM Advanced Development	    May 1986
 *
 *  This program is based on the paper:
 *
 *	"Primitives for the Manipulation of General Subdivisions and the
 *	 Computation of Voronoi Diagrams"
 *
 *	L. Guibas, J. Stolfi, ACM TOG, April 1985
 *
 *  Unfinished work:
 *
 *  o	Compute Voronoi from the Delaunay after triangualation is done
 *  o	Compute EMST
 *  o	Modify for GPX color displays
 */

#include stdio
#include math
#include descrip

#include "vsi$library:uisentry.h"
#include "vsi$library:uisusrdef.h"

#define DEFAULT_ATTR 0
#define OVER_ATTR 1
#define ERAS_ATTR 2
#define COMP_ATTR 3
#define TEXT_ATTR 4

#define PI 3.14159265358979323846264338

int vd_id, wd_id;
float wc_xmin, wc_xmax, wc_ymin, wc_ymax;
float vd_width, vd_height;

int seed = 1234567;

float mth$random();

$DESCRIPTOR(wsname, "SYS$WORKSTATION");
$DESCRIPTOR(title, "Delaunay Triangulation");

float x, y;
float x0, y0, r;
float x1, y1, x2, y2, x3, y3;

float xcur, ycur;

int box_flag = 0;

plot_box(x, y, s)
float x, y, s;
{
  float xl = x-s, yl = y-s;
  float xh = x+s, yh = y+s;

  uis$plot(&vd_id, &COMP_ATTR,
	   &xl, &yl, &xh, &yl,
	   &xh, &yl, &xh, &yh,
	   &xh, &yh, &xl, &yh,
	   &xl, &yh, &xl, &yl);

  box_flag = 1-box_flag;
}

move_ast()
{
  if (uis$get_pointer_position(&vd_id, &wd_id, &x, &y)) {
    }
}

exit_ast()
{
}

setup_window(range, size)
float range, size;
{
  wc_xmin = wc_ymin = -range;
  wc_xmax = wc_ymax =  range;
  vd_width = size;
  vd_height = size;
  vd_id = uis$create_display(&wc_xmin, &wc_ymin, &wc_xmax, &wc_ymax,
			     &vd_width, &vd_height);
  wd_id = uis$create_window(&vd_id, &wsname, &title);

  uis$set_writing_mode(&vd_id, &DEFAULT_ATTR, &OVER_ATTR, &UIS$C_MODE_OVER);
  uis$set_writing_mode(&vd_id, &DEFAULT_ATTR, &ERAS_ATTR, &UIS$C_MODE_ERAS);
  uis$set_writing_mode(&vd_id, &DEFAULT_ATTR, &COMP_ATTR, &UIS$C_MODE_COMP);

  uis$set_pointer_ast(&vd_id, &wd_id,
		      &move_ast, 0,
		      0, 0, 0, 0,
		      &exit_ast, 0);
  box_flag = 0;
}

clnup_window()
{
  uis$delete_display(&vd_id);
}

typedef struct site_struct {
  float x, y;
  } site_struct;

typedef struct edge_struct {
  int next[4];
  site_struct *data[4];
  } edge_struct;

#define NSITES 1000

site_struct sites[NSITES];
int nsites = 100;

/*
 *  Methods to access edge records
 */

#define ROT(e) (((e)&0xfffffffc)+(((e)+1)&3))
#define SYM(e) (((e)&0xfffffffc)+(((e)+2)&3))
#define ROT3(e) (((e)&0xfffffffc)+(((e)+3)&3))
#define ONEXT(e) ((edge_struct *)((e)&0xfffffffc))->next[(e)&3]
#define RNEXT(e) ((edge_struct *)((e)&0xfffffffc))->next[((e)+1)&3]
#define DNEXT(e) ((edge_struct *)((e)&0xfffffffc))->next[((e)+2)&3]
#define LNEXT(e) ((edge_struct *)((e)&0xfffffffc))->next[((e)+3)&3]
#define ORG(e) ((edge_struct *)((e)&0xfffffffc))->data[(e)&3]
#define RIGHT(e) ((edge_struct *)((e)&0xfffffffc))->data[((e)+1)&3]
#define DEST(e) ((edge_struct *)((e)&0xfffffffc))->data[((e)+2)&3]
#define LEFT(e) ((edge_struct *)((e)&0xfffffffc))->data[((e)+3)&3]
#define OPREV(e) ROT(RNEXT(e))
#define RPREV(e) ROT(DNEXT(e))

int norm_flag = 0;

float ca, sa;

void delaunay();

int makedge();
void deletedge();
void splice();

int incircle();
float ccw();
int rightof();
int leftof();

void makesites();
void sortsites();

int zone_id;

int lib$get_vm();
int lib$free_vm();
int lib$create_vm_zone();
int lib$delete_vm_zone();

/*
 *  Make the initial set of random sites
 */

void makesites(nsites)
int nsites;
{
  int i, j;
  float s, t;

  if (norm_flag) {	/* central limit approx to normal distribuition */
    for (i=0; i<nsites; i++) {
      for (s=t=0.0, j=0; j<9; j++,
	   s += mth$random(&seed)-0.5, t += mth$random(&seed)-0.5) ;
      sites[i].x = s/5.0;
      sites[i].y = t/5.0;
      }
    }
  else {		/* uniform distribution */
    for (i=0; i<nsites; i++) {
      sites[i].x = mth$random(&seed)-0.5;
      sites[i].y = mth$random(&seed)-0.5;
      }
    }
}

/*
 *  Sort the sites into (x, y) order along a given axis; triangulation
 *  should be the same regardless of sort axis.  Shell sort.
 */

void sortsites(nsites)
int nsites;
{
  int gap, i, j;
  site_struct tmp;

  for (gap = nsites/2; gap > 0; gap /= 2)
    for (i = gap; i < nsites; i++)
      for (j= i-gap; j >= 0 && 
	(ca*sites[j].x-sa*sites[j].y != ca*sites[j+gap].x-sa*sites[j+gap].y ?
	 ca*sites[j].x-sa*sites[j].y > ca*sites[j+gap].x-sa*sites[j+gap].y :
	 sa*sites[j].x+ca*sites[j].y > sa*sites[j+gap].x+ca*sites[j+gap].y) ; j -= gap) {
	tmp = sites[j]; sites[j] = sites[j+gap]; sites[j+gap] = tmp;
	}
}

/*
 *  Show an edge of the Delaunay in the window
 */

showedge(e)
int e;
{
  site_struct *o, *d;

  o = ORG(e); d = DEST(e);
  uis$plot(&vd_id, &COMP_ATTR, &o->x, &o->y, &d->x, &d->y);
}

/*
 *  Erase an edge of the Delaunay from the window
 */

erasedge(e)
int e;
{
  site_struct *o, *d;

  o = ORG(e); d = DEST(e);
  uis$plot(&vd_id, &COMP_ATTR, &o->x, &o->y, &d->x, &d->y);
}

/*
 *  Recursively create the Delaunay triangulation of a set of sites
 *  sorted in (x,y) order; divide and conquer algorithm.
 */

void delaunay(sl, sh, le, re)
int sl, sh;
int *le, *re;
{
  int a, b, c;
  int t;
  float ct;
  int ldo, ldi, rdi, rdo;
  int basel;
  int lcand, rcand;
  int sm;

  if (sh == sl+2) {
    a = makedge();
    ORG(a) = &sites[sl]; DEST(a) = &sites[sl+1];
    showedge(a);
    *le = a;  *re = SYM(a);
    return;
    }

  if (sh == sl+3) {
    a = makedge();
    b = makedge();
    splice(SYM(a), b);
    ORG(a) = &sites[sl]; DEST(a) = &sites[sl+1];
    ORG(b) = &sites[sl+1];  DEST(b) = &sites[sl+2];
    ct = ccw(&sites[sl], &sites[sl+1], &sites[sl+2]);
    if (ct == 0.0) {
      showedge(a); showedge(b);
      *le = a; *re = SYM(b);
      return;
      }
    else {
      c = connect(b, a);
      if (ct > 0.0) { *le = a; *re = SYM(b); }
      else { *le = SYM(c); *re = c; }
      showedge(a); showedge(b);
      return;
      }
    }

  sm = (sl+sh)/2;
  delaunay(sl, sm, &ldo, &ldi);
  delaunay(sm, sh, &rdi, &rdo);

  while (1) {
    if (leftof(ORG(rdi), ldi)) ldi = ROT(LNEXT(ldi));
    else if (rightof(ORG(ldi), rdi)) rdi = ONEXT(SYM(rdi));
    else break;
    }

  basel = connect(SYM(rdi), ldi);
  if (ORG(ldi) == ORG(ldo)) ldo = SYM(basel);
  if (ORG(rdi) == ORG(rdo)) rdo = basel;

  while (1) {

    lcand = ONEXT(SYM(basel));
    if (rightof(DEST(lcand), basel))
      while (incircle(DEST(basel), ORG(basel), DEST(lcand), DEST(ONEXT(lcand)))) {
	t = ONEXT(lcand); deletedge(lcand); lcand = t;
	}

    rcand = OPREV(basel);
    if (rightof(DEST(rcand), basel))
      while (incircle(DEST(basel), ORG(basel), DEST(rcand), DEST(OPREV(rcand)))) {
	t = OPREV(rcand); deletedge(rcand); rcand = t;
	}

    if (!rightof(DEST(lcand), basel) && !rightof(DEST(rcand), basel)) break;

    if (!rightof(DEST(lcand), basel) ||
	(rightof(DEST(rcand), basel) && incircle(DEST(lcand),
						 ORG(lcand),
						 ORG(rcand),
						 DEST(rcand))))
      basel = connect(rcand, SYM(basel));
    else
      basel = connect(SYM(basel), SYM(lcand));
    }
  *le = ldo; *re = rdo;
  return;
}

/*
 *  Splice primitive
 */

void splice(a, b)
int a, b;
{
  int ta, tb;
  int alpha = ROT(ONEXT(a));
  int beta = ROT(ONEXT(b));

  ta = ONEXT(a);
  tb = ONEXT(b);
  ONEXT(a) = tb;
  ONEXT(b) = ta;
  ta = ONEXT(alpha);
  tb = ONEXT(beta);
  ONEXT(alpha) = tb;
  ONEXT(beta) = ta;    
}

/*
 *  Make a new edge
 */

int makedge()
{
  int e;

  lib$get_vm(&sizeof(edge_struct), &e, &zone_id);
  ONEXT(e) = e;
  DNEXT(e) = SYM(e);
  RNEXT(e) = ROT3(e);
  LNEXT(e) = ROT(e);
  return e;
}

/*
 *  Delete an edge
 */

void deletedge(e)
int e;
{
  erasedge(e);
  splice(e, OPREV(e));
  splice(SYM(e), OPREV(SYM(e)));  
}

/*
 *  Connect two vertices with a new edge
 */

int connect(a, b)
int a, b;
{
  int e;

  e = makedge();
  ORG(e) = DEST(a);
  DEST(e) = ORG(b);
  splice(e, ROT(LNEXT(a)));
  splice(SYM(e), b);
  showedge(e);
  return e;
}

/*
 *  Test if point to right of given edge
 */

int rightof(s, e)
site_struct *s;
int e;
{
  return ccw(s, DEST(e), ORG(e)) > 0.0;
}

/*
 *  Test if point to left of given edge
 */

int leftof(s, e)
site_struct *s;
int e;
{
  return ccw(s, ORG(e), DEST(e)) > 0.0;
}

/*
 *  Counterclockwise triangle predicate
 */

float ccw(a, b, c)
site_struct *a, *b, *c;
{
  double x1 = a->x, y1 = a->y;
  double x2 = b->x, y2 = b->y;
  double x3 = c->x, y3 = c->y;

  return (x2*y3-y2*x3) - (x1*y3-y1*x3) + (x1*y2-y1*x2);
}

/*
 *  InCircle predicate
 */

int incircle(a, b, c, d)
site_struct *a, *b, *c, *d;
{
  double x1 = a->x, y1 = a->y;
  double x2 = b->x, y2 = b->y;
  double x3 = c->x, y3 = c->y;
  double x4 = d->x, y4 = d->y;

  return ((y4-y1)*(x2-x3)+(x4-x1)*(y2-y3))*((x4-x3)*(x2-x1)-(y4-y3)*(y2-y1)) >
	 ((y4-y3)*(x2-x1)+(x4-x3)*(y2-y1))*((x4-x1)*(x2-x3)-(y4-y1)*(y2-y3));
}

float ask(query)
char query[];
{
  char reply[132];
  float a;

  do {
    printf(query);
    gets(reply);
    } while(sscanf(reply, "%f", &a) != 1);
  return a;
}

main()
{
  char dummy[132];
  int i;
  int le, re;
  float a;

  lib$create_vm_zone(&zone_id);

  seed = ask("enter seed: ");
  for (i=0; i<250; i++) mth$random(&seed);

  nsites = ask("enter number of sites: ");
  if (nsites < 10) nsites = 10;
  if (nsites > NSITES) nsites = NSITES;

  norm_flag = ask("enter 1 for normal distribution: ");

  a = ask("enter angle of sort axis: ");
  ca = cos(PI/180.0*a); sa = sin(PI/180.0*a);

  makesites(nsites);
  sortsites(nsites);

  setup_window(0.6, 30.0);

  uis$erase(&vd_id);

  for (i=0; i<nsites; i++) {
    plot_box(sites[i].x, sites[i].y, .005);
    }

  delaunay(0, nsites, &le, &re);

  printf("delaunay returned %x %x\n", le, re);

  gets(dummy);

  lib$delete_vm_zone(&zone_id);

  clnup_window();
}

    
    (I used REPLY/NOEXTRACT so notes wouldn't copy all 500 lines of
    your program over the net merely to let me reply!)
    
    ANyway, I've got DECwindows, so if you convert to xlib, I'll try it.
    
    Hey, has anyone got a general uis-to-decwindows converter ?  Some
    sort of transfer vector perhaps with jacket routines ?
    
    I'd like to address the "depends on what you mean by random" question.
    
    Perhaps one usable definition might be:
    
    	Pick a square on a plane.  Define (x,y) pairs by
    	picking n uniformly distributed random numbers from 0 through the
    	the length of the side of the square.  Call these the x's.  Do the
    	same for the y's.
    
    	Now let the circles expand into the polygons.  Ignore the polygons
    	whose edges share the edge of the square.
    
    	Average the number of sides of the complete polygons.
    
    	What does the average approach as we let n approach infinity ?
    
    /Eric

    There's an amazingly elegant algorithm for computing the Voronoi
    diagram in O(n*log(n)) time that can be thought of as inspired by
    the observation in the base note about expanding circles about each
    site.

    Consider identical opaque cones growing out above each site in the plane.
    They will intersect in hyperbolic curves.

    The edges of the Voronoi diagram are simply the projection of the
    pieces of the hyperbolas which are visible from below!

    This can be used to get a beautiful algorithm for computing the
    Voronoi.  Consider sweeping a plane inclinded at the same angle as
    the sides of the cones.  It will intersect the cones in parabolas.
    Projecting these parabolas on the plane containing the sites will
    result in a parabolic wavefront advancing across the plane, if only
    the visible pieces of each parabola are kept.  If we stop each time
    a new site encounters the sweep plane, and each time a visible vertex
    of the parabolic wave front changes, we can trace out the voronoi
    diagram.  It turns out to be fairly easy to queue up these possible
    events in the order of their appearance - since such queue manipulation
    can be done in O(log(n)) time, and such manipulation has to be done
    for each site or Voronoi vertex (there are O(n) of them) and O(n*log(n))
    algorithm results.

    The divide and conquer algorithm posted earler also uses 3-D geometry -
    in that case the points are lifted up to a paraboloid of revolution
    above the plane.  Then making the Delaunay diagram is equivalent to
    calculating the convex hull of the projected points in 3-D.

    I'd like to try out the plane sweep/cone algorithm (which was written
    up by S. Fortune in an ACM conference on Computational Geometry, 1986)
    and compare it to the divide and conquer method...

    - Jim