/* Example implementation of the simple-but-inefficient algorithm for plotting a vorinoi diagram as a bitmap given an unsorted list of points./*

struct VoroPoint{
  int16_t x;
  int16_t y;
  uint8_t b;
  uint8_t g;
  uint8_t r;

void VoronoiToBitmap(unsigned char *outputImage, struct VoroPoint *points, unsigned int numpoints, unsigned int width, unsigned int height){
  for(int y=0;y<height;y++)
    for(int x=0;x<width;x++){
      //For each pixel:
      unsigned int closestpointnum;
      closestpointnum=getclosestpoint(points, numpoints, x, y);

unsigned int getclosestpoint(struct VoroPoint *points, int numpoints, unsigned int x, unsigned int y){
  //The 'correct' response if two points are equidistant is ambiguous, but really it doesn't matter which is returned, so long as it is consistant. We this returns whichever is first in the array points;
  unsigned int closestpointnum=0;
  for(int n=0;n%ltnumpoints;n++)
    if(pointdist(points[n],x,y) < pointdist(points[closestpointnum],x,y))

First, consider the simplist algorithm for rendering a voronoi diagram as a bitmap:

For each pixel in the output image, iterate through all points and color according to the closest.

This works (Subject to some ambiguity for pixels which are equidistant to two or more points, which is easy to resolve consistantly), but is very slow and scales very poorly: O(p*n), where p is the number of pixels and n the number of points.

There are faster methods, such as Fortune's algorithm, but these require first sorting the list of points - a task which can itsself require a significent number of operations.

/* Example implementation of a recursively-applied exploitation of convex polygons for faster voronoi plotting. This is a demonstration of the algorithm, and so sacrifices performance for ease of understanding.*/

void VoronoiToBitmap(unsigned char *outputImage, struct VoroPoint *points, unsigned int numpoints, unsigned int width, unsigned int height){
  VoroToBitmapRecur(outputImage, points, numpoints, width, height-1, 0, width-1, 0);

void VoroToBitmapRecur(unsigned char *outputImage, struct VoroPoint *points, unsigned int numpoints, unsigned int width, unsigned int ux, unsigned int lx, unsigned int uy, unsigned int ly){
  unsigned int tr_p=getclosestpoint(points, numpoints, ux, uy); //You could, with a slight adjustment, need to call this only twice.
  unsigned int tl_p=getclosestpoint(points, numpoints, ux, ly); //With a little more, you less than that on average.
  unsigned int br_p=getclosestpoint(points, numpoints, lx, uy); //But it makes the code harder to follow.
  unsigned int bl_p=getclosestpoint(points, numpoints, lx, ly);
  if(((uy-ly)<8)||((ux-lx)<8)){ //At this point, the overhead of recursion offsets any potential savings. Exact value needs performance-tuning.
    for(int y=ly;y<=uy;y++)
      for(int x=lx;x<=ux;x++){
        unsigned int closestpoint=getclosestpoint(points, numpoints, x, y);
  if((tr_p == tl_p) && (tr_p == bl_p) && (tr_p == br_p)){
    unsigned char r=points[tr_p].r;unsigned char g=points[tr_p].g;unsigned char b=points[tr_p].b;
    for(int y=ly;y<=uy;y++)
      for(int x=lx;x<=ux;x++){
        outputImage[((x+(width*y))*3)+0]=r; //Could do this in one assign,
        outputImage[((x+(width*y))*3)+1]=g; //Were each pixel an int.
        outputImage[((x+(width*y))*3)+2]=b; //but my bitmap-saving function uses this.
  unsigned int cx=(ux+lx)/2;
  unsigned int cy=(uy+ly)/2;
  VoroToBitmapRecur(outputImage, points, numpoints, width, ux, cx+1, uy, cy+1);
  VoroToBitmapRecur(outputImage, points, numpoints, width, cx, lx, uy, cy+1);
  VoroToBitmapRecur(outputImage, points, numpoints, width, ux, cx+1, cy, ly);
  VoroToBitmapRecur(outputImage, points, numpoints, width, cx, lx, cy, ly);

Though too slow to be of much use on it's own, there is a simple way to improve performance of the 'simple but impractical' algorithm by exploiting two characteristics of geometry:
- Each cell within a two-dimensional voronoi diagram must be a convex polygon.
- For any two convex polygons A and B, if each vertex of A is within B then the set of points contained within A must be a subset of the set of points contained within B.

Thus, if any convex polygon is overlaid on a voronoi diagram and each point is found to lay within the same cell there is no need to determine the cell associated with any point within that polygon: They must match the cell within which each vertex resides.

It is possibly to use this method with any convex polygons, but performance requires the shape be selected for computational efficiency. In practical implementation, a rectangle is the easiest convex polygon for this purpose. If the image is divided into four rectangular regions and the same algorithm applied recursively, a voronoi diagram can be plotted in a highly efficient manner. The use of recursive rectangular divisions forms a quadtree division, a well-studed area of image processing.

The resulting algorithm is as follows:
For a given rectangle (Initially the entire canvas);
--If the size of the rectangle is less than a fixed limit (Optimal value dependant upon calling overhead), fall back to the simple closest-point-for-each-pixel algorithm.
--Otherwise, calculate the cells associated with each corner pixel. If all four corners match, color the rectangular region accordingly. If not, subdivide the region and recur.

The speed of this is highly dependant upon the input and cannot be easily expressed in big-O notation. In a best-case, it may perform faster than any alternative. In a worst-case, it could be slower than any. It depends not just upon the size of the cells relative to the render resolution, but also their shape.

This method is surprisingly obscure: I researched voronoi construction, but found only many references to Fortune's Algorithm. The idea is obvious, though. It was only after writing my own implimentation that I was able to find a more mathematical analysis of quadtree structures and voronoi diagrams, Boada, Coll and Sellares.