
Voronoi Tesselation & Delaunay Graph 
Evert
Member #794
November 2000

For a game I'm planning to work on I want to have a playing field (2D) that is build up out of a simple tiling. Don't ask me why but for this game I want to generate the tiling from a set of vertices using Voronoi tesselation to generate the polygons.
After a few hours of thinking and typing I ended up with the following piece of (ugly and somewhat lengthy) code as a first draft: [code]#include <allegro.h> #include <math.h> typedef struct { double x, y; } SITE; typedef struct { double xc, yc; double x1, y1; double x2, y2; double dx, dy; } EDGE; typedef struct { EDGE *edges; int num_edges; int max_edges; } EDGE_LIST; double scale = 50; SITE sites[] = { {1.2, 0.0}, {1.1, 0.0}, {0.0,0.0},{1.0, 2.0}, {1.0, 2.0} }; int red, white, blue, green, black; void draw_sites(SITE *sites, int num_sites) { int n; for (n=0; n<num_sites; n++) putpixel(screen, SCREEN_W/2 + scale * sites[n].x, SCREEN_H/2 + scale * sites[n].y, red); } void add_edge(EDGE_LIST *edge_list, EDGE *edge) { if (edge_list>num_edges >= edge_list>max_edges) { edge_list>max_edges = 2*edge_list>max_edges + 1; edge_list>edges = realloc(edge_list>edges, edge_list>max_edges*sizeof *(edge_list>edges)); } edge_list>edges[edge_list>num_edges] = *edge; edge_list>num_edges++; } void compute_edges(SITE *sites, int num_sites, EDGE_LIST *edge_list, int *num_nn, int **site_nn) { EDGE edge; int n, n2; int c; /* Compute edges for each pair of vertices */ for (n=0; n<num_sites; n++) { for (n2 = 0; n2<num_nn[n]; n2++) { c = site_nn[n][n2]; edge.xc = (sites[n].x + sites[c].x)/2; edge.yc = (sites[n].y + sites[c].y)/2; edge.dx = (sites[n].y  sites[c].y)/2; edge.dy = (sites[n].x  sites[c].x)/2; edge.x1 = edge.xc  50*edge.dx; edge.y1 = edge.yc  50*edge.dy; edge.x2 = edge.xc + 50*edge.dx; edge.y2 = edge.yc + 50*edge.dy; add_edge(edge_list, &edge); } } /* Now constrain edges: for each pair of edges, determine where they * intersect and cutoff the line there */ for (n=0; n<edge_list>num_edges; n++) { for (c=0; c<edge_list>num_edges; c++) if (c!=n) { double t1 = 0.0, t2 = 0.0; double x, y, dx, dy, dx2, dy2; double vx, vy, wx, wy; double wvx, wvy, rx, ry; double det, invdet; /* Find intersection point */ vx = edge_list>edges[n].dx; vy = edge_list>edges[n].dy; wx = edge_list>edges[c].dx; wy = edge_list>edges[c].dy; rx = edge_list>edges[c].xc  edge_list>edges[n].xc; ry = edge_list>edges[c].yc  edge_list>edges[n].yc; /* Make sure that the direction vectors point to a common point */ if (vx*rx + vy*ry < 0) { vx = vx; vy = vy; } if (wx*rx + wy*ry > 0) { wx = wx; wy = wy; } /* Determinant of the matrix for (t t'); nescessary to invert the * matrix */ det = wx*vy  wy*vx; if (fabs(det) > 0.0) { invdet = 1.0/det; t2 = (vy*rx  vx*ry) * invdet; t1 = (wy*rx  wx*ry) * invdet; } else { t1 = t2 = 0; } x = edge_list>edges[n].xc + t1*vx; y = edge_list>edges[n].yc + t1*vy; /* Determine where to intersect */ dx = edge_list>edges[n].x1  edge_list>edges[n].xc; dy = edge_list>edges[n].y1  edge_list>edges[n].yc; dx2 = x  edge_list>edges[n].xc; dy2 = y  edge_list>edges[n].yc; if ( dx2 * dx + dy2 * dy > 0) { if ( dx2*dx2 + dy2*dy2 < dx*dx + dy*dy) { edge_list>edges[n].x1 = x; edge_list>edges[n].y1 = y; } } else { dx = edge_list>edges[n].x2  edge_list>edges[n].xc; dy = edge_list>edges[n].y2  edge_list>edges[n].yc; if ( dx2*dx2 + dy2*dy2 < dx*dx + dy*dy) { edge_list>edges[n].x2 = x; edge_list>edges[n].y2 = y; } } dx = edge_list>edges[c].x1  edge_list>edges[c].xc; dy = edge_list>edges[c].y1  edge_list>edges[c].yc; } } } void draw_edges(EDGE_LIST *edge_list) { EDGE edge; int n; for (n=0; n<edge_list>num_edges; n++) { edge = edge_list>edges[n]; line(screen, SCREEN_W/2 + scale * edge.x1, SCREEN_H/2 + scale * edge.y1, SCREEN_W/2 + scale * edge.x2, SCREEN_H/2 + scale * edge.y2, green); } for (n=0; n<edge_list>num_edges; n++) { edge = edge_list>edges[n]; putpixel(screen, SCREEN_W/2 + scale * edge.xc, SCREEN_H/2 + scale * edge.yc, blue); putpixel(screen, SCREEN_W/2 + scale * edge.x1, SCREEN_H/2 + scale * edge.y1, black); putpixel(screen, SCREEN_W/2 + scale * edge.x2, SCREEN_H/2 + scale * edge.y2, black); } } /* Find nearest neighbours (points in a Dalaunay graph). * Idea: two points are nearest neighbours iff there is a circle through both * of them that does not has any other sites in its interior */ void find_neighbours(SITE *sites, int num_sites, int *num_nn, int **site_nn) { int a, b, n; double r, r1, r2, q1, q2; double x, y; double cx, cy; double nx, ny; double dx, dy; /* Clear lists of nearest neighbours */ memset(num_nn, 0, num_sites * sizeof *num_nn); for (n=0; n<num_sites; n++) memset(site_nn[n], 0, num_sites * sizeof *(site_nn[n])); /* Find circles that intersect two points and don't contain interior points */ for (a=0; a<num_sites; a++) { for (b=a+1; b<num_sites; b++) { int done_plus; int done_neg; /* Initial guess of the radius */ x = sites[b].x  sites[a].x; y = sites[b].y  sites[a].y; r = 0.5*sqrt(x*x + y*y); r1 = r2 = r; done_plus = FALSE; done_neg = FALSE; do { /* Calculate centre of sphere with radius r through points a and b */ /* Relative to point a! */ cx = 0.5*(sites[b].x  sites[a].x); cy = 0.5*(sites[b].y  sites[a].y); /* Normal */ nx = cy / sqrt(cx*cx + cy*cy); ny = cx / sqrt(cx*cx + cy*cy); /* Make sure no other points are inside the circle */ done_plus = TRUE; done_neg = TRUE; for (n=0; n<num_sites; n++) if (n!=a && n!=b) { q1 = sqrt(r1*r1  (cx*cx + cy*cy)); q2 = sqrt(r2*r2  (cx*cx + cy*cy)); x = cx + q1*nx + sites[a].x; y = cy + q1*ny + sites[a].y; dx = x  sites[n].x; dy = y  sites[n].y; if (dx*dx + dy*dy < r1*r1) { done_plus = FALSE; r1 = sqrt(dx*dx + dy*dy); } x = cx  q2*nx + sites[a].x; y = cy  q2*ny + sites[a].y; dx = x  sites[n].x; dy = y  sites[n].y; if (dx*dx + dy*dy < r2*r2) { done_neg = FALSE; r2 = sqrt(dx*dx + dy*dy); } if (!done_neg && !done_plus) break; } if (r1 < r && r2 < r) break; } while (!done_plus && !done_neg); /* Two possibilities: if either of these is ok, then the two points * are nearest neigbours and we can use them for voronoi tesselation */ if (done_plus  done_neg) { site_nn[a][num_nn[a]] = b; site_nn[b][num_nn[b]] = a; num_nn[a]++; num_nn[b]++; } } } } /* Find nearest neighbours (points in a Dalaunay graph). * Idea: two points are nearest neighbours iff there is a circle through both * of them that does not has any other sites in its interior */ void draw_neighbours(SITE *sites, int num_sites, int *num_nn, int **site_nn) { int a, b, n; double r, r1, r2, q1, q2; double x, y; double cx, cy; double nx, ny; double dx, dy; /* Find circles that intersect two points and don't contain interior points */ for (a=0; a<num_sites; a++) { for (n=0; n<num_nn[a]; n++) { b = site_nn[a][n]; line(screen, SCREEN_W/2 + scale * sites[a].x, SCREEN_H/2 + scale * sites[a].y, SCREEN_W/2 + scale * sites[b].x, SCREEN_H/2 + scale * sites[b].y, black); } } } /* Draw empty circles between each pair of vertices */ /* Aid in finding nearest neighbours: two cells are NN iff there is a circle * through both of them that does not contain interior points. */ void draw_circles(SITE *sites, int num_sites) { int a, b, n; double r, r1, r2, q1, q2; double x, y; double cx, cy; double nx, ny; double dx, dy; for (a=0; a<num_sites; a++) { for (b=a+1; b<num_sites; b++) { int done_plus; int done_neg; /* Initial guess of the radius */ /* This won't be enough, but need a better way of getting the radius */ /* with which to try again if the search fails. This works for now. */ x = (sites[b].x  sites[a].x); y = (sites[b].y  sites[a].y); r = 0.5*sqrt(x*x + y*y); r1 = r2 = r; done_plus = FALSE; done_neg = FALSE; do { /* Calculate centre of sphere with radius r through points a and b */ /* Relative to point a! */ cx = 0.5*(sites[b].x  sites[a].x); cy = 0.5*(sites[b].y  sites[a].y); /* Normal */ nx = cy / sqrt(cx*cx + cy*cy); ny = cx / sqrt(cx*cx + cy*cy); /* Make sure no other points are inside the circle */ done_plus = TRUE; done_neg = TRUE; for (n=0; n<num_sites; n++) if (n!=a && n!=b) { q1 = sqrt(r1*r1  (cx*cx + cy*cy)); q2 = sqrt(r2*r2  (cx*cx + cy*cy)); x = cx + q1*nx + sites[a].x; y = cy + q1*ny + sites[a].y; dx = x  sites[n].x; dy = y  sites[n].y; if (dx*dx + dy*dy < r1*r1) { done_plus = FALSE; r1 = sqrt(dx*dx + dy*dy); } x = cx  q2*nx + sites[a].x; y = cy  q2*ny + sites[a].y; dx = x  sites[n].x; dy = y  sites[n].y; if (dx*dx + dy*dy < r2*r2) { done_neg = FALSE; r2 = sqrt(dx*dx + dy*dy); } if (!done_neg && !done_plus) break; } if (r1 < r && r2 < r) break; } while (!done_plus && !done_neg); /* Two possibilities: if either of these is ok, then the two points * are nearest neigbours and we can use them for voronoi tesselation */ if (done_plus) { x = cx + q1*nx + sites[a].x; y = cy + q1*ny + sites[a].y; circle(screen, SCREEN_W/2 + scale * x, SCREEN_H/2 + scale * y, scale * r1, blue); } if (done_neg) { x = cx  q2*nx + sites[a].x; y = cy  q2*ny + sites[a].y; circle(screen, SCREEN_W/2 + scale * x, SCREEN_H/2 + scale * y, scale * r2, green); } } } } int main(void) { EDGE_LIST edge_list = { NULL, 0, 0 }; int *site_num_neighbours; int **site_neighbours; int n, c; n = sizeof sites / sizeof *sites; site_num_neighbours = malloc(n * sizeof *site_num_neighbours); site_neighbours = malloc(n * sizeof *site_neighbours); for (c=0; c<n; c++) { site_neighbours[c] = malloc(n * sizeof *site_neighbours); } allegro_init(); install_keyboard(); set_color_depth(16); set_gfx_mode(GFX_AUTODETECT_WINDOWED, 640, 480, 0,0); red = makecol(255, 0, 0); white = makecol(255, 255, 255); blue = makecol(0, 0, 255); green = makecol(0, 255, 0); black = makecol(0, 0, 0); clear_to_color(screen, white); /* Find nearest neighbours for all sites */ find_neighbours(sites, n, site_num_neighbours, site_neighbours); /* Compute voronoi edges in the set */ compute_edges(sites, n, &edge_list, site_num_neighbours, site_neighbours); /* Graphics output */ /* Draw empty circles */ draw_circles(sites, n); /* Draw lines between neighbours (Delaunay graph) */ draw_neighbours(sites, n, site_num_neighbours, site_neighbours); /* Draw edges */ draw_edges(&edge_list); /* Draw all voronoi sites */ draw_sites(sites, n); readkey(); return 0; } END_OF_MAIN()[/code] What I do is construct the lines that separate two points (called sites in the code) and cut them off in places where they intersect eachother. For this, I need to find the nearest neighbours of each point, which I take from the Delaunay graph. This works fine in principle but I've hit a few snags. First, my way of calculating the Delaunay graph is a bit whacky, probably slower than it has to be and I haven't coded it in the most flexible way, but it's the best I could think of for now. It works for the grids I've tried it with anyway. It works by looking for the smallest circle (well, here I was sloppy and just did something I know won't always work, but it works for now) that connects two points and if the circle contains no points in its interior, then the two points are marked as nearest neighbours. Second, for each particle I calculate the dividing line between it and all it's neighbours. The lines are then clipped when they intersect eachother. The intention was that I'd end up with a set of edges that describe polygons with one point inside. Here's the snag: when an edge is clipped, it can happen that the edge it's clipped against is itself clipped, so that the first edge now ends in the middle of nowhere instead of joining with other edges in a vertex. I can of course look for such occurances afterward and extend the line again to the nearest vertex along the extension of the edge, but that feels a bit klunky and I'm sure it can be done better. So my question(s) to people who know more about this subject than I do: what is a better way of finding the Voronoi tesselation from the Delaunay graph? Is there a more efficient way to get the proper edges? For that matter, is there a more reliable and straightforward algorithm to find the Delaunay graph? I googled around a bit but things I could find were either scientific papers or lowlevel stuff that didn't contain much new information. I found a nice Java applet that let me draw and play with these things (which is how I got convinced that I want to use this for the game in question) but I can't make heads or tails of the code. Hmm... I just remembered that I have a book on graph theory... I'll have a look in that as well. In the meantime all help is appreciated! 
Rash
Member #2,374
May 2002

I'm convinced that if Allegro.cc ever gets an official FAQ, it should make prominent mention of the comp.graphics.algorithms FAQ, as it often holds the solution to many of the complex questions asked here. In your case there is the following entry: 
Inphernic
Member #1,111
March 2001

