Voronoi diagram: Difference between revisions
Line 644: | Line 644: | ||
This a direct reformulation of the explicit version. |
This a direct reformulation of the explicit version. |
||
<lang j>Voronoi=. ,"0/&i./@:] (i. <./)@:(+/@:*:@:-"1)"1 _ ] ?.@$~ 2 ,~ [ |
<lang j>Voronoi=. ,"0/&i./@:] (i. <./)@:(+/@:*:@:-"1)"1 _ ] ?.@$~ 2 ,~ [ |
||
25 viewmat@:([ load bind'viewmat')@:Voronoi 500 500</lang> |
|||
=={{header|Liberty BASIC}}== |
=={{header|Liberty BASIC}}== |
Revision as of 22:32, 25 June 2014
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>
- include <stdlib.h>
- include <string.h>
- 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; }
- 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; }
- 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); }
- define for_i for (i = 0; i < size_y; i++)
- 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); }
- 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++
- include <windows.h>
- include <vector>
- include <string>
//-------------------------------------------------------------------------------------------------- using namespace std;
//-------------------------------------------------------------------------------------------------- class point { public:
int x, y;
point() { x = y = 0; } point( int a, int b ) { x = a; y = b; } int distanceSqrd( const point& p ) {
int xd = p.x - x, yd = p.y - y;
return xd * xd + yd * yd;
}
}; //-------------------------------------------------------------------------------------------------- class myBitmap { public:
myBitmap() : pen( NULL ) {} ~myBitmap() {
DeleteObject( pen ); DeleteDC( hdc ); DeleteObject( bmp );
}
bool create( int w, int h ) {
BITMAPINFO bi; void *pBits; 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;
HDC dc = GetDC( GetConsoleWindow() ); bmp = CreateDIBSection( dc, &bi, DIB_RGB_COLORS, &pBits, NULL, 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 );
}
void saveBitmap( string path ) {
BITMAPFILEHEADER fileheader; BITMAPINFO infoheader; BITMAP bitmap; DWORD* dwpBits; DWORD wb; HANDLE file;
GetObject( bmp, sizeof( bitmap ), &bitmap );
dwpBits = new DWORD[bitmap.bmWidth * bitmap.bmHeight]; ZeroMemory( dwpBits, 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 )dwpBits, &infoheader, DIB_RGB_COLORS );
file = CreateFile( path.c_str(), GENERIC_WRITE, 0, NULL, CREATE_ALWAYS, FILE_ATTRIBUTE_NORMAL, NULL ); WriteFile( file, &fileheader, sizeof( BITMAPFILEHEADER ), &wb, NULL ); WriteFile( file, &infoheader.bmiHeader, sizeof( infoheader.bmiHeader ), &wb, NULL ); WriteFile( file, dwpBits, bitmap.bmWidth * bitmap.bmHeight * 4, &wb, NULL ); CloseHandle( file );
delete [] dwpBits;
}
HDC getDC() { return hdc; } int getWidth() { return width; } int getHeight() { return height; }
private:
HBITMAP bmp; HDC hdc; HPEN pen; int width, height;
}; //-------------------------------------------------------------------------------------------------- class Voronoi { public:
void make( myBitmap* bmp, int count ) {
_bmp = bmp; createPoints( count ); createColors(); createSites(); setSitesPoints();
}
private:
void createSites() {
int w = _bmp->getWidth(), h = _bmp->getHeight(), d; for( int hh = 0; hh < h; hh++ ) { for( int ww = 0; ww < w; ww++ ) { point bpt( ww, hh ); int ind = -1, dist = INT_MAX; for( int it = 0; it < points.size(); it++ ) { d = ( points[it] ).distanceSqrd( bpt ); if( d < dist ) { dist = d; ind = it; } }
if( ind > -1 ) SetPixel( _bmp->getDC(), ww, hh, colors[ind] ); else __asm nop // should never happen! }
} }
void setSitesPoints() {
for( vector<point>::iterator it = points.begin(); it < points.end(); it++ ) { int x = ( *it ).x, y = ( *it ).y; for( int i = -1; i < 2; i++ ) for( int j = -1; j < 2; j++ ) SetPixel( _bmp->getDC(), x + i, y + j, 0 ); }
}
void createPoints( int count ) {
int w = _bmp->getWidth() - 20, h = _bmp->getHeight() - 20; for( int i = 0; i < count; i++ ) { point p( rand() % w + 10, rand() % h + 10 ); points.push_back( p ); }
}
void createColors() {
for( int 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.getDC(), 0, 0, SRCCOPY ); bmp.saveBitmap( "f://rc//v.bmp" );
system( "pause" );
return 0;
} //-------------------------------------------------------------------------------------------------- </lang>
D
<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
<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>
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>
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
- w "trapclose quit"
- w.gb1 "down ; fill black ; size 4"
- 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
- 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
- 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
- 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>
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
<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
<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>
Sample Output:
Racket
First approach
<lang racket>
- 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>
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
<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>
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
<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}
}
}
- Generate a 600x400 Voronoi diagram with 60 random points
image create photo demo -width 600 -height 400 pack [label .l -image demo]
- To display while generating, uncomment this line and the other one so commented
- update
voronoi demo 60</lang>
XPL0
<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>