Voronoi diagram

Voronoi diagram
You are encouraged to solve this task according to the task description, using any language you may know.

A Voronoi diagram is a diagram consisting of a number of sites. Each Voronoi site s also has a Voronoi cell consisting of all points closest to s.

The task is to demonstrate how to generate and display a Voroni diagram. See algo K-means++ clustering.

C

C code drawing a color map of a set of Voronoi sites. Image is in PNM P6, written to stdout. Run as a.out > stuff.pnm. <lang c>#include <stdio.h>

1. include <stdlib.h>
2. include <string.h>

1. define N_SITES 150

double site[N_SITES][2]; unsigned char rgb[N_SITES][3];

int size_x = 640, size_y = 480;

inline double sq2(double x, double y) { return x * x + y * y; }

1. define for_k for (k = 0; k < N_SITES; k++)

int nearest_site(double x, double y) { int k, ret = 0; double d, dist = 0; for_k { d = sq2(x - site[k][0], y - site[k][1]); if (!k || d < dist) { dist = d, ret = k; } } return ret; }

/* see if a pixel is different from any neighboring ones */ int at_edge(int *color, int y, int x) { int i, j, c = color[y * size_x + x]; for (i = y - 1; i <= y + 1; i++) { if (i < 0 || i >= size_y) continue;

for (j = x - 1; j <= x + 1; j++) { if (j < 0 || j >= size_x) continue; if (color[i * size_x + j] != c) return 1; } } return 0; }

1. define AA_RES 4 /* average over 4x4 supersampling grid */

void aa_color(unsigned char *pix, int y, int x) { int i, j, n; double r = 0, g = 0, b = 0, xx, yy; for (i = 0; i < AA_RES; i++) { yy = y + 1. / AA_RES * i + .5; for (j = 0; j < AA_RES; j++) { xx = x + 1. / AA_RES * j + .5; n = nearest_site(xx, yy); r += rgb[n][0]; g += rgb[n][1]; b += rgb[n][2]; } } pix[0] = r / (AA_RES * AA_RES); pix[1] = g / (AA_RES * AA_RES); pix[2] = b / (AA_RES * AA_RES); }

1. define for_i for (i = 0; i < size_y; i++)
2. define for_j for (j = 0; j < size_x; j++)

void gen_map() { int i, j, k; int *nearest = malloc(sizeof(int) * size_y * size_x); unsigned char *ptr, *buf, color;

ptr = buf = malloc(3 * size_x * size_y); for_i for_j nearest[i * size_x + j] = nearest_site(j, i);

for_i for_j { if (!at_edge(nearest, i, j)) memcpy(ptr, rgb[nearest[i * size_x + j]], 3); else /* at edge, do anti-alias rastering */ aa_color(ptr, i, j); ptr += 3; }

/* draw sites */ for (k = 0; k < N_SITES; k++) { color = (rgb[k][0]*.25 + rgb[k][1]*.6 + rgb[k][2]*.15 > 80) ? 0 : 255;

for (i = site[k][1] - 1; i <= site[k][1] + 1; i++) { if (i < 0 || i >= size_y) continue;

for (j = site[k][0] - 1; j <= site[k][0] + 1; j++) { if (j < 0 || j >= size_x) continue;

ptr = buf + 3 * (i * size_x + j); ptr[0] = ptr[1] = ptr[2] = color; } } }

printf("P6\n%d %d\n255\n", size_x, size_y); fflush(stdout); fwrite(buf, size_y * size_x * 3, 1, stdout); }

1. define frand(x) (rand() / (1. + RAND_MAX) * x)

int main() { int k; for_k { site[k][0] = frand(size_x); site[k][1] = frand(size_y); rgb [k][0] = frand(256); rgb [k][1] = frand(256); rgb [k][2] = frand(256); }

gen_map(); return 0; }</lang>

C++

<lang cpp>

1. include <windows.h>
2. include <vector>
3. include <string>

using namespace std;

////////////////////////////////////////////////////// struct Point {

 int x, y;

};

////////////////////////////////////////////////////// class MyBitmap {

public:
 MyBitmap() : pen_(nullptr) {}
 ~MyBitmap() {
   DeleteObject(pen_);
   DeleteDC(hdc_);
   DeleteObject(bmp_);
 }

 bool Create(int w, int h) {
   BITMAPINFO	bi;
   ZeroMemory(&bi, sizeof(bi));

   bi.bmiHeader.biSize = sizeof(bi.bmiHeader);
   bi.bmiHeader.biBitCount = sizeof(DWORD) * 8;
   bi.bmiHeader.biCompression = BI_RGB;
   bi.bmiHeader.biPlanes = 1;
   bi.bmiHeader.biWidth = w;
   bi.bmiHeader.biHeight = -h;

   void *bits_ptr = nullptr;
   HDC dc = GetDC(GetConsoleWindow());
   bmp_ = CreateDIBSection(dc, &bi, DIB_RGB_COLORS, &bits_ptr, nullptr, 0);
   if (!bmp_) return false;

   hdc_ = CreateCompatibleDC(dc);
   SelectObject(hdc_, bmp_);
   ReleaseDC(GetConsoleWindow(), dc);

   width_ = w;
   height_ = h;

   return true;
 }

 void SetPenColor(DWORD clr) {
   if (pen_) DeleteObject(pen_);
   pen_ = CreatePen(PS_SOLID, 1, clr);
   SelectObject(hdc_, pen_);
 }

 bool SaveBitmap(const char* path) {
   HANDLE file = CreateFile(path, GENERIC_WRITE, 0, nullptr, CREATE_ALWAYS, FILE_ATTRIBUTE_NORMAL, nullptr);
   if (file == INVALID_HANDLE_VALUE) {
     return false;
   }

   BITMAPFILEHEADER fileheader;
   BITMAPINFO infoheader;
   BITMAP bitmap;    
   GetObject(bmp_, sizeof(bitmap), &bitmap);

   DWORD* dwp_bits = new DWORD[bitmap.bmWidth * bitmap.bmHeight];
   ZeroMemory(dwp_bits, bitmap.bmWidth * bitmap.bmHeight * sizeof(DWORD));
   ZeroMemory(&infoheader, sizeof(BITMAPINFO));
   ZeroMemory(&fileheader, sizeof(BITMAPFILEHEADER));

   infoheader.bmiHeader.biBitCount = sizeof(DWORD) * 8;
   infoheader.bmiHeader.biCompression = BI_RGB;
   infoheader.bmiHeader.biPlanes = 1;
   infoheader.bmiHeader.biSize = sizeof(infoheader.bmiHeader);
   infoheader.bmiHeader.biHeight = bitmap.bmHeight;
   infoheader.bmiHeader.biWidth = bitmap.bmWidth;
   infoheader.bmiHeader.biSizeImage = bitmap.bmWidth * bitmap.bmHeight * sizeof(DWORD);

   fileheader.bfType = 0x4D42;
   fileheader.bfOffBits = sizeof(infoheader.bmiHeader) + sizeof(BITMAPFILEHEADER);
   fileheader.bfSize = fileheader.bfOffBits + infoheader.bmiHeader.biSizeImage;

   GetDIBits(hdc_, bmp_, 0, height_, (LPVOID)dwp_bits, &infoheader, DIB_RGB_COLORS);

   DWORD wb;
   WriteFile(file, &fileheader, sizeof(BITMAPFILEHEADER), &wb, nullptr);
   WriteFile(file, &infoheader.bmiHeader, sizeof(infoheader.bmiHeader), &wb, nullptr);
   WriteFile(file, dwp_bits, bitmap.bmWidth * bitmap.bmHeight * 4, &wb, nullptr);
   CloseHandle(file);

   delete[] dwp_bits;
   return true;
 }

 HDC hdc() { return hdc_; }
 int width() { return width_; }
 int height() { return height_; }

private:
 HBITMAP bmp_;
 HDC hdc_;
 HPEN pen_;
 int width_, height_;

};

static int DistanceSqrd(const Point& point, int x, int y) {

 int xd = x - point.x;
 int yd = y - point.y;
 return (xd * xd) + (yd * yd);

}

////////////////////////////////////////////////////// class Voronoi {

public:
 void Make(MyBitmap* bmp, int count) {
   bmp_ = bmp;
   CreatePoints(count);
   CreateColors();
   CreateSites();
   SetSitesPoints();
 }

private:
 void CreateSites() {
   int w = bmp_->width(), h = bmp_->height(), d;
   for (int hh = 0; hh < h; hh++) {
     for (int ww = 0; ww < w; ww++) {
       int ind = -1, dist = INT_MAX;
       for (size_t it = 0; it < points_.size(); it++) {
         const Point& p = points_[it];
         d = DistanceSqrd(p, ww, hh);
         if (d < dist) {
           dist = d;
           ind = it;
         }
       }

       if (ind > -1)
         SetPixel(bmp_->hdc(), ww, hh, colors_[ind]);
       else
         __asm nop // should never happen!
       }
   }
 }

 void SetSitesPoints() {
   for (const auto& point : points_) {
     int x = point.x, y = point.y;
     for (int i = -1; i < 2; i++)
       for (int j = -1; j < 2; j++)
         SetPixel(bmp_->hdc(), x + i, y + j, 0);
   }
 }

 void CreatePoints(int count) {
   const int w = bmp_->width() - 20, h = bmp_->height() - 20;
   for (int i = 0; i < count; i++) {
     points_.push_back({ rand() % w + 10, rand() % h + 10 });
   }
 }

 void CreateColors() {
   for (size_t i = 0; i < points_.size(); i++) {
     DWORD c = RGB(rand() % 200 + 50, rand() % 200 + 55, rand() % 200 + 50);
     colors_.push_back(c);
   }
 }

 vector<Point> points_;
 vector<DWORD> colors_;
 MyBitmap* bmp_;

};

////////////////////////////////////////////////////// int main(int argc, char* argv[]) {

 ShowWindow(GetConsoleWindow(), SW_MAXIMIZE);
 srand(GetTickCount());

 MyBitmap bmp;
 bmp.Create(512, 512);
 bmp.SetPenColor(0);

 Voronoi v;
 v.Make(&bmp, 50);

 BitBlt(GetDC(GetConsoleWindow()), 20, 20, 512, 512, bmp.hdc(), 0, 0, SRCCOPY);
 bmp.SaveBitmap("v.bmp");

 system("pause");

 return 0;

} </lang>

D

Translation of: Go

<lang d>import std.random, std.algorithm, std.range, bitmap;

struct Point { uint x, y; }

enum randomPoints = (in size_t nPoints, in size_t nx, in size_t ny) =>

   nPoints.iota
   .map!((int) => Point(uniform(0, nx), uniform(0, ny)))
   .array;

Image!RGB generateVoronoi(in Point[] pts,

                         in size_t nx, in size_t ny) /*nothrow*/ {
   // Generate a random color for each centroid.
   immutable rndRBG = (int) => RGB(uniform!"[]"(ubyte.min, ubyte.max),
                                   uniform!"[]"(ubyte.min, ubyte.max),
                                   uniform!"[]"(ubyte.min, ubyte.max));
   const colors = pts.length.iota.map!rndRBG.array;

   // Generate diagram by coloring pixels with color of nearest site.
   auto img = new typeof(return)(nx, ny);
   foreach (immutable x; 0 .. nx)
       foreach (immutable y; 0 .. ny) {
           immutable dCmp = (in Point a, in Point b) pure nothrow =>
               ((a.x - x) ^^ 2 + (a.y - y) ^^ 2) <
               ((b.x - x) ^^ 2 + (b.y - y) ^^ 2);
           // img[x, y] = colors[pts.reduce!(min!dCmp)];
           img[x, y] = colors[pts.length - pts.minPos!dCmp.length];
       }

   // Mark each centroid with a white dot.
   foreach (immutable p; pts)
       img[p.tupleof] = RGB.white;
   return img;

}

void main() {

   enum imageWidth = 640,
        imageHeight = 480;
   randomPoints(150, imageWidth, imageHeight)
   .generateVoronoi(imageWidth, imageHeight)
   .savePPM6("voronoi.ppm");

}</lang>

Go

Output png

<lang go>package main

import (

   "fmt"
   "image"
   "image/color"
   "image/draw"
   "image/png"
   "math/rand"
   "os"
   "time"

)

const (

   imageWidth  = 300
   imageHeight = 200
   nSites      = 10

)

func main() {

   writePngFile(generateVoronoi(randomSites()))

}

func generateVoronoi(sx, sy []int) image.Image {

   // generate a random color for each site
   sc := make([]color.NRGBA, nSites)
   for i := range sx {
       sc[i] = color.NRGBA{uint8(rand.Intn(256)), uint8(rand.Intn(256)),
           uint8(rand.Intn(256)), 255}
   }

   // generate diagram by coloring each pixel with color of nearest site
   img := image.NewNRGBA(image.Rect(0, 0, imageWidth, imageHeight))
   for x := 0; x < imageWidth; x++ {
       for y := 0; y < imageHeight; y++ {
           dMin := dot(imageWidth, imageHeight)
           var sMin int
           for s := 0; s < nSites; s++ {
               if d := dot(sx[s]-x, sy[s]-y); d < dMin {
                   sMin = s
                   dMin = d
               }
           }
           img.SetNRGBA(x, y, sc[sMin])
       }
   }
   // mark each site with a black box
   black := image.NewUniform(color.Black)
   for s := 0; s < nSites; s++ {
       draw.Draw(img, image.Rect(sx[s]-2, sy[s]-2, sx[s]+2, sy[s]+2),
           black, image.ZP, draw.Src)
   }
   return img

}

func dot(x, y int) int {

   return x*x + y*y

}

func randomSites() (sx, sy []int) {

   rand.Seed(time.Now().Unix())
   sx = make([]int, nSites)
   sy = make([]int, nSites)
   for i := range sx {
       sx[i] = rand.Intn(imageWidth)
       sy[i] = rand.Intn(imageHeight)
   }
   return

}

func writePngFile(img image.Image) {

   f, err := os.Create("voronoi.png")
   if err != nil {
       fmt.Println(err)
       return
   }
   if err = png.Encode(f, img); err != nil {
       fmt.Println(err)
   }
   if err = f.Close(); err != nil {
       fmt.Println(err)
   }

}</lang>

Haskell

Uses the repa and repa-io libraries. <lang haskell> -- Compile with: ghc -O2 -fllvm -fforce-recomp -threaded --make {-# LANGUAGE BangPatterns #-} module Main where

import System.Random

import Data.Word import Data.Array.Repa as Repa import Data.Array.Repa.IO.BMP

{-# INLINE sqDistance #-} sqDistance :: Word32 -> Word32 -> Word32 -> Word32 -> Word32 sqDistance !x1 !y1 !x2 !y2 = ((x1-x2)^2) + ((y1-y2)^2)

centers :: Int -> Int -> Array U DIM2 Word32 centers nCenters nCells =

   fromListUnboxed (Z :. nCenters :. 2) $ take (2*nCenters) $ randomRs (0, fromIntegral nCells) (mkStdGen 1)

applyReduce2 arr f =

   traverse arr (\(i :. j) -> i) $ \lookup (Z:.i) ->
       f (lookup (Z:.i:.0)) (lookup (Z:.i:.1))

minimize1D arr = foldS f h t

 where
   indexed arr = traverse arr id (\src idx@(Z :. i) -> (src idx, (fromIntegral i)))        
   (Z :. n) = extent arr
   iarr = indexed arr
   h = iarr ! (Z :. 0)
   t = extract (Z :. 1) (Z :. (n-1)) iarr

   f min@(!valMin, !iMin ) x@(!val, !i) | val < valMin = x
                                        | otherwise = min

voronoi :: Int -> Int -> Array D DIM2 Word32 voronoi nCenters nCells =

   let
     {-# INLINE cellReducer #-}
     cellReducer = applyReduce2 (centers nCenters nCells)
     {-# INLINE nearestCenterIndex #-}
     nearestCenterIndex = snd . (Repa.! Z) . minimize1D
   in        
     Repa.fromFunction (Z :. nCells :. nCells :: DIM2) $ \ (Z:.i:.j) ->
         nearestCenterIndex $ cellReducer (sqDistance (fromIntegral i) (fromIntegral j))

genColorTable :: Int -> Array U DIM1 (Word8, Word8, Word8) genColorTable n = fromListUnboxed (Z :. n) $ zip3 l1 l2 l3

   where
     randoms = randomRs (0,255) (mkStdGen 1)
     (l1, rest1) = splitAt n randoms
     (l2, rest2) = splitAt n rest1
     l3 = take n rest2

colorize :: Array U DIM1 (Word8, Word8, Word8) -> Array D DIM2 Word32 -> Array D DIM2 (Word8, Word8, Word8) colorize ctable = Repa.map $ \x -> ctable Repa.! (Z:. fromIntegral x)

main = do

 let nsites = 150
 let ctable = genColorTable nsites 
 voro <- computeP $ colorize ctable (voronoi nsites 512) :: IO (Array U DIM2 (Word8, Word8, Word8))
 writeImageToBMP "out.bmp" voro

</lang>

Icon and Unicon

The sample images to the right show the screen size, number of sites, and metric used in the title bar.

<lang Icon>link graphics,printf,strings

record site(x,y,colour) # site data position and colour invocable all # needed for string metrics

procedure main(A) # voronoi

&window := open("Voronoi","g","bg=black") | stop("Unable to open window")

WAttrib("canvas=hidden") # figure out maximal size width & height WAttrib(sprintf("size=%d,%d",WAttrib("displaywidth"),WAttrib("displayheight"))) WAttrib("canvas=maximal") height := WAttrib("height") width  := WAttrib("width")

metrics := ["hypot","taxi","taxi3"] # different metrics

while case a := get(A) of { # command line arguments

 "--sites"  | "-s" : sites  := 0 < integer(a := get(A)) | runerr(205,a)
 "--height" | "-h" : height := 0 < (height >= integer(a := get(A))) | runerr(205,a)
 "--width"  | "-w" : width  := 0 < (width  >= integer(a := get(A))) | runerr(205,a)
 "--metric" | "-m" : metric := ((a := get(A)) == !metrics) | runerr(205,a)
 "--help"   | "-?" : write("Usage:\n voronoi [[--sites|-s] n] ",

"[[--height|-h] pixels] [[--width|-w] pixels]", "[[--metric|-m] metric_procedure]", "[--help|-?]\n\n")

 }

/metric := metrics[1] # default to normal /sites := ?(r := integer(.1*width)) + r # sites = random .1 to .2 of width if not given

WAttrib(sprintf("label=Voronoi %dx%d %d %s",width,height,sites,metric)) WAttrib(sprintf("size=%d,%d",width,height))

x := "0123456789abcdef" # hex for random sites (colour) siteL := [] every 1 to sites do # random sites

 put(siteL, site(?width,?height,cat("#",?x,?x,?x,?x,?x,?x))) 

VoronoiDiagram(width,height,siteL,metric) # Voronoi-ize it WDone() end

procedure hypot(x,y,site) # normal metric return sqrt((x-site.x)^2 + (y-site.y)^2) end

procedure taxi(x,y,site) # "taxi" metric return abs(x-site.x)+abs(y-site.y) end

procedure taxi3(x,y,site) # copied from a commented out version (TCL) return (abs(x-site.x)^3+abs(y-site.y)^3)^(.3) end

procedure VoronoiDiagram(width,height,siteL,metric)

  /metric := hypot
  every y := 1 to height & x := 1 to width do {
     dist := width+height         # anything larger than diagonal
     every site := !siteL do {
        if dist < (dt :=  metric(x,y,site)) then next  # skip
        else if dist >:= dt then Fg(site.colour)       # site
        else Fg("#000000")                             # unowned
        DrawPoint(x,y)
        }
     }

  Fg("Black")
  every site := !siteL do                              # mark sites
     DrawCircle(site.x,site.y,1)     

end</lang>

Library: Icon Programming Library

printf.icn provides the printf family graphics.icn provides graphics support strings.icn provides cat

J

Explicit version

A straightforward solution: generate random points and for each pixel find the index of the least distance. Note that the square root is avoided to improve performance. <lang j>NB. (number of points) voronoi (shape) NB. Generates an array of indices of the nearest point voronoi =: 4 :0

 p =. (x,2) ?@$ y
 (i.<./)@:(+/@:*:@:-"1&p)"1 ,"0/&i./ y

)

load'viewmat' viewmat 25 voronoi 500 500</lang>

Another solution generates Voronoi cells from Delaunay triangulation. The page Voronoi diagram/J/Delaunay triangulation also contains a convex hull algorithm.

Tacit version

This a direct reformulation of the explicit version.

<lang j>Voronoi=. ,"0/&i./@:] (i. <./)@:(+/@:*:@:-"1)"1 _ ] ?@$~ 2 ,~ [ 25 viewmat@:([ load bind'viewmat')@:Voronoi 500 500</lang>

Java

Library: Swing

Library: AWT

<lang java>import java.awt.Color; import java.awt.Graphics; import java.awt.Graphics2D; import java.awt.geom.Ellipse2D; import java.awt.image.BufferedImage; import java.io.File; import java.io.IOException; import java.util.Random;

import javax.imageio.ImageIO; import javax.swing.JFrame;

public class Voronoi extends JFrame { static double p = 3; static BufferedImage I; static int px[], py[], color[], cells = 100, size = 1000;

public Voronoi() { super("Voronoi Diagram"); setBounds(0, 0, size, size); setDefaultCloseOperation(EXIT_ON_CLOSE); int n = 0; Random rand = new Random(); I = new BufferedImage(size, size, BufferedImage.TYPE_INT_RGB); px = new int[cells]; py = new int[cells]; color = new int[cells]; for (int i = 0; i < cells; i++) { px[i] = rand.nextInt(size); py[i] = rand.nextInt(size); color[i] = rand.nextInt(16777215);

} for (int x = 0; x < size; x++) { for (int y = 0; y < size; y++) { n = 0; for (byte i = 0; i < cells; i++) { if (distance(px[i], x, py[i], y) < distance(px[n], x, py[n], y)) { n = i;

} } I.setRGB(x, y, color[n]);

} }

Graphics2D g = I.createGraphics(); g.setColor(Color.BLACK); for (int i = 0; i < cells; i++) { g.fill(new Ellipse2D .Double(px[i] - 2.5, py[i] - 2.5, 5, 5)); }

try { ImageIO.write(I, "png", new File("voronoi.png")); } catch (IOException e) {

}

}

public void paint(Graphics g) { g.drawImage(I, 0, 0, this); }

static double distance(int x1, int x2, int y1, int y2) { double d; d = Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); // Euclidian // d = Math.abs(x1 - x2) + Math.abs(y1 - y2); // Manhattan // d = Math.pow(Math.pow(Math.abs(x1 - x2), p) + Math.pow(Math.abs(y1 - y2), p), (1 / p)); // Minkovski return d; }

public static void main(String[] args) { new Voronoi().setVisible(true); } } </lang>

Liberty BASIC

For first site it fills the table with distances to that site. For other sites it looks at vertical lines left and right from its location. If no place on a vertical line is closer to the current site, then there's no point looking further left or right. Don't bother square-rooting to get distances.. <lang lb> WindowWidth =600 WindowHeight =600

sites = 100 xEdge = 400 yEdge = 400 graphicbox #w.gb1, 10, 10, xEdge, yEdge

open "Voronoi neighbourhoods" for window as #w

1. w "trapclose quit"
2. w.gb1 "down ; fill black ; size 4"
3. w.gb1 "font courier_new 12"

dim townX( sites), townY( sites), col$( sites)

for i =1 to sites

   townX( i) =int( xEdge *rnd( 1))
   townY( i) =int( yEdge *rnd( 1))
   col$( i) = int( 256 *rnd( 1)); " "; int( 256 *rnd( 1)); " "; int( 256 *rnd( 1))
   #w.gb1 "color "; col$( i)
   #w.gb1 "set "; townX( i); " "; townY( i)

next i

1. w.gb1 "size 1"

dim nearestIndex(xEdge, yEdge) dim dist(xEdge, yEdge)

start = time$("ms")

'fill distance table with distances from the first site for x = 0 to xEdge - 1

   for y = 0 to yEdge - 1
       dist(x, y) = (townX(1) - x) ^ 2 + (townY(1) - y) ^ 2
       nearestIndex(x, y) = 1
   next y

next x

1. w.gb1 "color darkblue"

'for other towns for i = 2 to sites

   'display some progress
   #w.gb1 "place 0 20"
   #w.gb1 "\computing: "; using("###.#", i / sites * 100); "%"
   'look left
   for x = townX(i) to 0 step -1
       if not(checkRow(i, x,0, yEdge - 1)) then exit for
   next x
   'look right
   for x = townX(i) + 1 to xEdge - 1
       if not(checkRow(i, x, 0, yEdge - 1)) then exit for
   next x
   scan

next i

for x = 0 to xEdge - 1

   for y =0 to yEdge - 1
       #w.gb1 "color "; col$(nearestIndex(x, y))
       startY = y
       nearest = nearestIndex(x, y)
       for y = y + 1 to yEdge
           if nearestIndex(x, y) <> nearest then y = y - 1 : exit for
       next y
       #w.gb1 "line "; x; " "; startY; " "; x; " "; y + 1
   next y

next x

1. w.gb1 "color black; size 4"

for i =1 to sites

   #w.gb1 "set "; townX( i); " "; townY( i)

next i print time$("ms") - start wait

sub quit w$

   close #w$
   end

end sub

function checkRow(site, x, startY, endY)

   dxSquared = (townX(site) - x) ^ 2
   for y = startY to endY
        dSquared = (townY(site) - y) ^ 2 + dxSquared
        if dSquared <= dist(x, y) then
            dist(x, y) = dSquared
            nearestIndex(x, y) = site
            checkRow = 1
        end if
   next y

end function </lang>

Lua

Library: LÖVE

This was adapted from the Python implementation to work with Lua/LÖVE. Credit for the algorithm goes to the author of that section. <lang lua> function love.load()

   love.math.setRandomSeed(os.time()) --set the random seed
   keys = {} --an empty table where we will store key presses
   number_cells = 50 --the number of cells we want in our diagram
   --draw the voronoi diagram to a canvas
   voronoiDiagram = generateVoronoi(love.window.getWidth(), love.window.getHeight(), number_cells)

end

function hypot(x,y)

   return math.sqrt(x*x + y*y)

end

function generateVoronoi(width, height, num_cells)

   canvas = love.graphics.newCanvas(width, height)
   imgx = canvas:getWidth()
   imgy = canvas:getHeight()
   local nx = {}
   local ny = {}
   local nr = {}
   local ng = {}
   local nb = {}
   for a = 1, num_cells do

table.insert(nx, love.math.random(0,imgx)) table.insert(ny, love.math.random(0,imgy)) table.insert(nr, love.math.random(0,255)) table.insert(ng, love.math.random(0,255)) table.insert(nb, love.math.random(0,255))

   end
   love.graphics.setColor({255,255,255})
   love.graphics.setCanvas(canvas)
   for y = 1, imgy do

for x = 1, imgx do

           dmin = hypot(imgx-1, imgy-1)

j = -1 for i = 1, num_cells do d = hypot(nx[i]-x, ny[i]-y) if d < dmin then dmin = d j = i end end love.graphics.setColor({nr[j], ng[j], nb[j]}) love.graphics.point(x, y) end

   end
   --reset color
   love.graphics.setColor({255,255,255})
   --draw points
   for b = 1, num_cells do

love.graphics.point(nx[b], ny[b])

   end
   love.graphics.setCanvas()
   return canvas

end

--RENDER function love.draw()

   --reset color
   love.graphics.setColor({255,255,255})
   --draw diagram
   love.graphics.draw(voronoiDiagram)
   --draw drop shadow text
   love.graphics.setColor({0,0,0})
   love.graphics.print("space: regenerate\nesc: quit",1,1)
   --draw text
   love.graphics.setColor({200,200,0})
   love.graphics.print("space: regenerate\nesc: quit")

end

--CONTROL function love.keyreleased(key)

   if key == ' ' then

voronoiDiagram = generateVoronoi(love.window.getWidth(), love.window.getHeight(), number_cells)

   elseif key == 'escape' then

love.event.quit()

   end

end </lang>

Mathematica

<lang Mathematica>Needs["ComputationalGeometry`"] DiagramPlot[{{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5, 14.9}, {13.2, 11.9}, {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6, 5.1}, {5.3, 2.4}, {8.45, 4.7}, {11.5, 9.6}, {13.8, 7.3}, {12.9, 3.1}, {11, 1.1}}]</lang>

МК-61/52

<lang>0 П4 0 П5 ИП0 1 - x^2 ИП1 1 - x^2 + КвКор П3 9 П6 КИП6 П8 {x} 2 10^x * П9 [x] ИП5 - x^2 ИП9 {x} 2 10^x * ИП4 - x^2 + КвКор П9 ИП3 - x<0 47 ИП9 П3 ИП6 П7 ИП6 ИП2 - 9 - x>=0 17 КИП7 [x] С/П КИП5 ИП5 ИП1 - x>=0 04 КИП4 ИП4 ИП0 - x>=0 02</lang>

Input: Р0 - diagram width; Р1 - diagram height; Р0 - number of the points; РA - РE - coordinates and colors of the points in format C,XXYY (example: 3,0102).

Example of the manually compiled output (graphical output from this class of devices is missing):

·	·	·	·	·	·	·	·	·	·
·	·	·	·	·	·	·	·	·	·
·	·	·	·	·	·	·	·	·	·
·	·	·	•	·	·	·	·	·	·
·	·	·	·	·	·	·	·	·	·
·	·	·	·	·	·	·	·	·	·
·	·	•	·	·	·	·	·	·	·
·	·	·	·	·	·	·	·	·	·
·	·	·	·	·	·	·	•	·	·
·	·	·	·	·	·	·	·	·	·

Prolog

Works with SWI-Prolog and XPCE.
3 Voronoi diagrams are given for the same sites, one with the Manhattan distance, one with the Euclidean distance and the last with the Minkowski distance (order 3).

<lang Prolog>:- dynamic pt/6. voronoi :- V is random(20) + 20, retractall(pt(_,_,_,_)), forall(between(1, V, I), ( X is random(390) + 5, Y is random(390) + 5, R is random(65535), G is random(65535), B is random(65535), assertz(pt(I,X,Y, R, G, B)) )), voronoi(manhattan, V), voronoi(euclide, V), voronoi(minkowski_3, V).

voronoi(Distance, V) :- sformat(A, 'Voronoi 400X400 ~w ~w', [V, Distance]), new(D, window(A)), send(D, size, size(400,400)), new(Img, image(@nil, width := 400, height := 400 , kind := pixmap)),

       % get the list of the sites

bagof((N, X, Y), R^G^B^pt(N, X, Y, R, G, B), L),

forall(between(0,399, I), forall(between(0,399, J), ( get_nearest_site(V, Distance, I, J, L, S), pt(S, _, _, R, G, B), send(Img, pixel(I, J, colour(@default, R, G, B)))))),

new(Bmp, bitmap(Img)), send(D, display, Bmp, point(0,0)), send(D, open).

% define predicatea foldl (functionnal spirit) foldl([], _Pred, R, R).

foldl([H | T], Pred, Acc, R) :- call(Pred, H, Acc, R1), foldl(T, Pred, R1, R).

% predicate for foldl compare(Distance, XP, YP, (N, X, Y), (D, S), R) :- call(Distance, XP, YP, X, Y, DT), ( DT < D -> R = (DT, N) ; R = (D, S)).

% use of a fake site for the init of foldl get_nearest_site(Distance, I, J, L, S) :- foldl(L, compare(Distance, I, J), (65535, nil), (_, S)).

manhattan(X1, Y1, X2, Y2, D) :- D is abs(X2 - X1) + abs(Y2-Y1).

euclide(X1, Y1, X2, Y2, D) :- D is sqrt((X2 - X1)**2 + (Y2-Y1)**2).

minkowski_3(X1, Y1, X2, Y2, D) :- D is (abs(X2 - X1)**3 + abs(Y2-Y1)**3)**0.33. </lang>

PureBasic

Euclidean

Voronoi Diagram in PureBasic

<lang PureBasic>Structure VCoo

 x.i:  y.i
 Colour.i: FillColour.i

EndStructure

Macro RandInt(MAXLIMIT)

 Int(MAXLIMIT*(Random(#MAXLONG)/#MAXLONG))

EndMacro

Macro SQ2(X, Y)

 ((X)*(X) + (Y)*(Y))

EndMacro

Procedure GenRandomPoints(Array a.VCoo(1), xMax, yMax, cnt)

 Protected i, j, k, l
 cnt-1
 Dim a(cnt)
 For i=0 To cnt
   a(i)\x = RandInt(xMax): a(i)\y = RandInt(yMax)
   j = RandInt(255): k = RandInt(255): l = RandInt(255)
   a(i)\Colour = RGBA(j, k, l, 255)
   a(i)\FillColour = RGBA(255-j, 255-k, 255-l, 255)
 Next i
 ProcedureReturn #True

EndProcedure

Procedure MakeVoronoiDiagram(Array a.VCoo(1),xMax, yMax) ; Euclidean

 Protected i, x, y, img, dist.d, dt.d
 img = CreateImage(#PB_Any, xMax+1, yMax+1)
 If StartDrawing(ImageOutput(img))
   For y=0 To yMax
     For x=0 To xMax
       dist = Infinity()
       For i=0 To ArraySize(a())
         dt = SQ2(x-a(i)\x, y-a(i)\y)
         If dt > dist 
           Continue
         ElseIf dt < dist
           dist = dt
           Plot(x,y,a(i)\FillColour)
         Else ; 'Owner ship' is unclear, set pixel to transparent.
           Plot(x,y,RGBA(0, 0, 0, 0))
         EndIf
       Next
     Next
   Next
   For i=0 To ArraySize(a())
     Circle(a(i)\x, a(i)\y, 1, a(i)\Colour)
   Next
   StopDrawing()
 EndIf
 ProcedureReturn img 

EndProcedure

Main code

Define img, x, y, file$ Dim V.VCoo(0) x = 640: y = 480 If Not GenRandomPoints(V(), x, y, 150): End: EndIf img = MakeVoronoiDiagram(V(), x, y) If img And OpenWindow(0, 0, 0, x, y, "Voronoi Diagram in PureBasic", #PB_Window_SystemMenu)

 ImageGadget(0, 0, 0, x, y, ImageID(img))
 Repeat: Until WaitWindowEvent() = #PB_Event_CloseWindow

EndIf

UsePNGImageEncoder() file$ = SaveFileRequester("Save Image?", "Voronoi_Diagram_in_PureBasic.png", "PNG|*.png", 0) If file$ <> ""

 SaveImage(img, file$, #PB_ImagePlugin_PNG)

EndIf</lang>

Taxicab

Voronoi Diagram in PureBasic

<lang PureBasic>Structure VCoo

 x.i:  y.i
 Colour.i: FillColour.i

EndStructure

Macro RandInt(MAXLIMIT)

 Int(MAXLIMIT*(Random(#MAXLONG)/#MAXLONG))

EndMacro

Procedure GenRandomPoints(Array a.VCoo(1), xMax, yMax, cnt)

 Protected i, j, k, l
 cnt-1
 Dim a(cnt)
 For i=0 To cnt
   a(i)\x = RandInt(xMax): a(i)\y = RandInt(yMax)
   j = RandInt(255): k = RandInt(255): l = RandInt(255)
   a(i)\Colour = RGBA(j, k, l, 255)
   a(i)\FillColour = RGBA(255-j, 255-k, 255-l, 255)
 Next i
 ProcedureReturn #True

EndProcedure

Procedure MakeVoronoiDiagram(Array a.VCoo(1),xMax, yMax)

 Protected i, x, y, img, dist, dt, dx, dy
 img = CreateImage(#PB_Any, xMax+1, yMax+1, 32)
 If StartDrawing(ImageOutput(img))
   For y=0 To yMax
     For x=0 To xMax
       dist = #MAXLONG
       For i=0 To ArraySize(a())
         dx = x-a(i)\x
         dy = y-a(i)\y
         dt = Sign(dx)*dx + Sign(dy)*dy
         If dt > dist ; no update
           Continue 
         ElseIf dt < dist  ; an new 'owner' is found
           dist = dt
           Plot(x,y,a(i)\FillColour)
         Else ; dt = dist
           Plot(x,y,RGBA(0,0,0,0)) ; no clear 'owner', make the pixel transparent
         EndIf
       Next
     Next
   Next
   For i=0 To ArraySize(a())
     Circle(a(i)\x, a(i)\y, 1, a(i)\Colour)
   Next
   StopDrawing()
 EndIf
 ProcedureReturn img 

EndProcedure

Main code

Define img, x, y, file$ Dim V.VCoo(0) x = 640: y = 480 If Not GenRandomPoints(V(), x, y, 150): End: EndIf img = MakeVoronoiDiagram(V(), x, y) If img And OpenWindow(0, 0, 0, x, y, "Voronoi Diagram in PureBasic", #PB_Window_SystemMenu)

 ImageGadget(0, 0, 0, x, y, ImageID(img))
 Repeat: Until WaitWindowEvent() = #PB_Event_CloseWindow

EndIf

UsePNGImageEncoder() file$ = SaveFileRequester("Save Image?", "Voronoi_Diagram_in_PureBasic.png", "PNG|*.png", 0) If file$ <> ""

 SaveImage(img, file$, #PB_ImagePlugin_PNG)

EndIf</lang>

Python

This implementation takes in a list of points, each point being a tuple and returns a dictionary consisting of all the points at a given site. <lang python>from PIL import Image import random import math

def generate_voronoi_diagram(width, height, num_cells): image = Image.new("RGB", (width, height)) putpixel = image.putpixel imgx, imgy = image.size nx = [] ny = [] nr = [] ng = [] nb = [] for i in range(num_cells): nx.append(random.randrange(imgx)) ny.append(random.randrange(imgy)) nr.append(random.randrange(256)) ng.append(random.randrange(256)) nb.append(random.randrange(256)) for y in range(imgy): for x in range(imgx): dmin = math.hypot(imgx-1, imgy-1) j = -1 for i in range(num_cells): d = math.hypot(nx[i]-x, ny[i]-y) if d < dmin: dmin = d j = i putpixel((x, y), (nr[j], ng[j], nb[j])) image.save("VoronoiDiagram.png", "PNG")

       image.show()

generate_voronoi_diagram(500, 500, 25)</lang>

Output:

Voronoi Diagram in Python

Racket

Clustering using the nearest neigbour approach.

First approach

<lang racket>

1. lang racket

(require plot)

Performs clustering of points in a grid

using the nearest neigbour approach and shows

clusters in different colors

(define (plot-Voronoi-diagram point-list)

 (define pts
   (for*/list ([x (in-range 0 1 0.005)]
               [y (in-range 0 1 0.005)])
     (vector x y)))
 
 (define clusters (clusterize pts point-list))
 
 (plot
  (append
   (for/list ([r (in-list clusters)] [i (in-naturals)])
     (points (rest r) #:color i #:sym 'fullcircle1))
   (list (points point-list #:sym 'fullcircle5 #:fill-color 'white)))))

Divides the set of points into clusters

using given centroids

(define (clusterize data centroids)

 (for*/fold ([res (map list centroids)]) ([x (in-list data)])
   (define c (argmin (curryr (metric) x) centroids))
   (dict-set res c (cons x (dict-ref res c)))))

</lang>

Different metrics <lang racket> (define (euclidean-distance a b)

 (for/sum ([x (in-vector a)] [y (in-vector b)]) 
   (sqr (- x y))))

(define (manhattan-distance a b)

 (for/sum ([x (in-vector a)] [y (in-vector b)]) 
   (abs (- x y))))

(define metric (make-parameter euclidean-distance)) </lang>

The contour plot of the classification function.

Using Manhattan metric.

Voronoi diagram in 3D space.

Alternative approach

<lang racket>

Plots the Voronoi diagram as a contour plot of

the classification function built for a set of points

(define (plot-Voronoi-diagram2 point-list)

 (define n (length point-list))
 (define F (classification-function point-list))
 (plot 
  (list
   (contour-intervals (compose F vector) 0 1 0 1
                      #:samples 300
                      #:levels n
                      #:colors (range n)
                      #:contour-styles '(solid)
                      #:alphas '(1))
   (points point-list #:sym 'fullcircle3))))

For a set of centroids returns a function

which finds the index of the centroid nearest

to a given point

(define (classification-function centroids)

 (define tbl 
   (for/hash ([p (in-list centroids)] [i (in-naturals)])
     (values p i)))
 (λ (x)
   (hash-ref tbl (argmin (curry (metric) x) centroids))))

</lang>

Output:

<lang racket> (define pts

 (for/list ([i 50]) (vector (random) (random))))

(display (plot-Voronoi-diagram pts))

(display (plot-Voronoi-diagram2 pts))

(parameterize ([metric manhattan-distance])

 (display (plot-Voronoi-diagram2 pts)))

Using the classification function it is possible to plot Voronoi diagram in 3D.

(define pts3d (for/list ([i 7]) (vector (random) (random) (random)))) (plot3d (list

        (isosurfaces3d (compose (classification-function pts3d) vector) 
                       0 1 0 1 0 1
                       #:line-styles '(transparent)
                       #:samples 100
                       #:colors (range 7)
                       #:alphas '(1))
        (points3d pts3d #:sym 'fullcircle3)))

</lang>

Ruby

Uses Raster graphics operations/Ruby

Sample output from Ruby program

<lang ruby>load 'raster_graphics.rb'

class ColourPixel < Pixel

 def initialize(x, y, colour)
   @colour = colour
   super x, y
 end
 attr_accessor :colour

 def distance_to(px, py)
   Math::hypot(px - x, py - y)
 end 

end

width, height = 300, 200 npoints = 20 pixmap = Pixmap.new(width,height)

@bases = npoints.times.collect do |i|

 ColourPixel.new(
     3+rand(width-6), 3+rand(height-6),  # provide a margin to draw a circle
     RGBColour.new(rand(256), rand(256), rand(256))
 )

end

pixmap.each_pixel do |x, y|

 nearest = @bases.min_by {|base| base.distance_to(x, y)}
 pixmap[x, y] = nearest.colour

end

@bases.each do |base|

 pixmap[base.x, base.y] = RGBColour::BLACK
 pixmap.draw_circle(base, 2, RGBColour::BLACK)

end

pixmap.save_as_png("voronoi_rb.png")</lang>

Run BASIC

<lang runbasic>graphic #g, 400,400

1. g flush()

spots = 100 leftSide = 400 rightSide = 400

dim locX(spots) dim locY(spots) dim rgb(spots,3) dim seal(leftSide, rightSide) dim reach(leftSide, rightSide)

for i =1 to spots

   locX(i)	= int(leftSide  * rnd(1))
   locY(i)	= int(rightSide * rnd(1))
   rgb(i,1)	= int(256 * rnd(1))
   rgb(i,2)	= int(256 * rnd(1))
   rgb(i,3)	= int(256 * rnd(1))
   #g color(rgb(i,1),rgb(i,2),rgb(i,3))
   #g set(locX(i),locY(i))

next i

1. g size(1)

' find reach to the first site for x = 0 to leftSide - 1

   for y = 0 to rightSide - 1
       reach(x, y) = (locX(1) - x) ^ 2 + (locY(1) - y) ^ 2
       seal(x, y) = 1
   next y

next x

1. g color("darkblue")

' spots other than 1st spot for i = 2 to spots

   for x = locX(i) to 0 step -1		' looking left
       if not(chkPos(i,x,0, rightSide - 1)) then exit for
   next x
   for x = locX(i) + 1 to leftSide - 1		' looking right
       if not(chkPos(i, x, 0, rightSide - 1)) then exit for
   next x

next i

for x = 0 to leftSide - 1

   for y = 0 to rightSide - 1

c1 = rgb(seal(x, y),1) c2 = rgb(seal(x, y),2) c3 = rgb(seal(x, y),3)

       #g color(c1,c2,c3)
       startY	= y
       nearest	= seal(x, y)
       for y = y + 1 to rightSide
           if seal(x, y) <> nearest then y = y - 1 : exit for
       next y
       #g line(x,startY,x,y + 1)
   next y

next x

1. g color("black")
2. g size(4)

for i =1 to spots

   #g set(locX(i),locY(i))

next i render #g end

function chkPos(site, x, startY, endY) dxSqr = (locX(site) - x) ^ 2 for y = startY to endY dSqr = (locY(site) - y) ^ 2 + dxSqr if dSqr <= reach(x, y) then reach(x,y) = dSqr seal(x,y) = site chkPos = 1 end if next y end function</lang>

Seed7

<lang seed7>$ include "seed7_05.s7i";

 include "draw.s7i";
 include "keybd.s7i";

const type: point is new struct

   var integer: xPos is 0;
   var integer: yPos is 0;
   var color: col is black;
 end struct;

const proc: generateVoronoiDiagram (in integer: width, in integer: height, in integer: numCells) is func

 local
   var array point: points is 0 times point.value;
   var integer: index is 0;
   var integer: x is 0;
   var integer: y is 0;
   var integer: distSquare is 0;
   var integer: minDistSquare is 0;
   var integer: indexOfNearest is 0;
 begin
   screen(width, height);
   points := numCells times point.value;
   for index range 1 to numCells do
     points[index].xPos := rand(0, width);
     points[index].yPos := rand(0, height);
     points[index].col := color(rand(0, 65535), rand(0, 65535), rand(0, 65535));
   end for;
   for y range 0 to height do
     for x range 0 to width do
       minDistSquare := width ** 2 + height ** 2;
       for index range 1 to numCells do
         distSquare := (points[index].xPos - x) ** 2 + (points[index].yPos - y) ** 2;
         if distSquare < minDistSquare then
           minDistSquare := distSquare;
           indexOfNearest := index;
         end if;
       end for;
       point(x, y, points[indexOfNearest].col);
     end for;
   end for;
   for index range 1 to numCells do
     line(points[index].xPos - 2, points[index].yPos, 4, 0, black);
     line(points[index].xPos, points[index].yPos - 2, 0, 4, black);
   end for;
 end func;

const proc: main is func

 begin
   generateVoronoiDiagram(500, 500, 25);
   KEYBOARD := GRAPH_KEYBOARD;
   readln(KEYBOARD);
 end func;</lang>

Original source: [1]

Tcl

Library: Tk

<lang tcl>package require Tk proc r to {expr {int(rand()*$to)}}; # Simple helper

proc voronoi {photo pointCount} {

   for {set i 0} {$i < $pointCount} {incr i} {

lappend points [r [image width $photo]] [r [image height $photo]]

   }
   foreach {x y} $points {

lappend colors [format "#%02x%02x%02x" [r 256] [r 256] [r 256]]

   }
   set initd [expr {[image width $photo] + [image height $photo]}]
   for {set i 0} {$i < [image width $photo]} {incr i} {

for {set j 0} {$j < [image height $photo]} {incr j} { set color black set d $initd foreach {x y} $points c $colors { set h [expr {hypot($x-$i,$y-$j)}] ### Other interesting metrics #set h [expr {abs($x-$i)+abs($y-$j)}] #set h [expr {(abs($x-$i)**3+abs($y-$j)**3)**0.3}] if {$d > $h} {set d $h;set color $c} } $photo put $color -to $i $j } # To display while generating, uncomment this line and the other one so commented #if {$i%4==0} {update idletasks}

   }

}

1. Generate a 600x400 Voronoi diagram with 60 random points

image create photo demo -width 600 -height 400 pack [label .l -image demo]

1. To display while generating, uncomment this line and the other one so commented
2. update

voronoi demo 60</lang>

XPL0

File:VoronoiXPL0.gif

<lang XPL0>include c:\cxpl\codes; \intrinsic 'code' declarations

def N = 15; \number of sites int SiteX(N), SiteY(N), \coordinates of sites

       Dist2, MinDist2, MinI,  \distance squared, and minimums
       X, Y, I;

[SetVid($13); \set 320x200x8 graphics for I:= 0 to N-1 do \create a number of randomly placed sites

       [SiteX(I):= Ran(160);  SiteY(I):= Ran(100)];

for Y:= 0 to 100-1 do \generate Voronoi diagram

   for X:= 0 to 160-1 do       \for all points...
       [MinDist2:= -1>>1;      \find closest site
       for I:= 0 to N-1 do
               [Dist2:= sq(X-SiteX(I)) + sq(Y-SiteY(I));
               if Dist2 < MinDist2 then
                       [MinDist2:= Dist2;  MinI:= I];
               ];
       if MinDist2 then Point(X, Y, MinI+1);   \leave center black
       ];

I:= ChIn(1); \wait for keystroke SetVid($03); \restore normal text screen ]</lang>

zkl

Uses the PPM class from http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm#zkl

Translation of: Python

<lang zkl>fcn generate_voronoi_diagram(width,height,num_cells){

  image,imgx,imgy:=PPM(width,height),width,height;
  nx:=num_cells.pump(List,(0).random.fp(imgx));
  ny:=num_cells.pump(List,(0).random.fp(imgy));
  nr:=num_cells.pump(List,(0).random.fp(256));  // red
  ng:=num_cells.pump(List,(0).random.fp(256));  // blue
  nb:=num_cells.pump(List,(0).random.fp(256));  // green

  foreach y,x in (imgy,imgx){
     dmin:=(imgx-1).toFloat().hypot(imgy-1);
     j:=-1;
     foreach i in (num_cells){
        d:=(nx[i] - x).toFloat().hypot(ny[i] - y);

if(d<dmin) dmin,j = d,i

     }
     image[x,y]=(nr[j]*0xff00 + ng[j])*0xff00 + nb[j];
  }
  image

}</lang> <lang zkl>generate_voronoi_diagram(500,500,25).write(File("VoronoiDiagram.ppm","wb"));</lang>

Output:

http://home.comcast.net/~zenkinetic/Images/VoronoiDiagram.jpg