# Voronoi diagrams for generators other than points

Any suggestions on how to determine a Voronoi diagram for sites other than points, as e.g. in the picture below? The input is a raster image.

• Answered at gis.stackexchange.com/a/53435 (another site, but the answer uses Mathematica. – whuber Mar 5 '13 at 20:15
• thanks for the link and for the answer below. – DeeDee Mar 5 '13 at 20:39
• – Szabolcs Mar 5 '13 at 21:15

Obtain the image:

i = Import["http://i.stack.imgur.com/iab6u.png"];


Compute the distance transform:

k = DistanceTransform[ColorNegate[i]] // ImageAdjust;
ReliefPlot[Reverse@ImageData[k]] (* To illustrate *) Identify the "peaks," which must bound the Voronoi cells:

l = ColorNegate[Binarize[ColorNegate[LaplacianGaussianFilter[k, 2] // ImageAdjust]]];


Clean the result and identify its connected components (the cells):

m = Erosion[Dilation[MorphologicalComponents[l] // Colorize, 2], 1];


Show this with the original features:

ImageMultiply[m, ColorNegate[i]]


## ### Edit

A cleaner solution--albeit one that takes substantially more processing time--exploits WatershedComponents (new in Version 8):

l = WatershedComponents[k];
m = Dilation[MorphologicalComponents[l] // Colorize, 1] (* Needs little or no cleaning *)
ImageMultiply[m, ColorNegate[i]] (* As before *) I like this one better, but fear it might take too much processing for large complex images.

• thank you very much for your answer. – DeeDee Mar 5 '13 at 20:38
• (+1) very nice & simple ;) – Vitaliy Kaurov Mar 5 '13 at 21:07

Here is a Nearest-based method. This is quite similar to what @Mr. Wizard did for approximating 3D (ordinary) Voronoi.

comps     = MorphologicalComponents[img];
cmap      = Flatten[MapIndexed[#2 -> #1 &, comps, {2}]];
comparray = DeleteCases[cmap, _ -> 0];
nf        = Nearest[comparray];


Now we build the table giving Voronoi components.

Timing[
voronoi2 =
Array[
First @ nf[{##}] &,
Length /@ {comps, comps[]}
];
]

(* Out= {0.640000, Null} *)


The picture:

MatrixPlot[voronoi2 + comps, ColorFunction -> "BrightBands", ImageSize -> 500] • Good one. Thanks, Daniel Lichtblau. – DeeDee Mar 6 '13 at 4:37
• +1 Nearest is flexible and squarely within the spirit of Voronoi tessellations. It is nice to see that it can be reasonably efficient with this problem. – whuber Mar 6 '13 at 18:12
• +1 as well, and thanks for the mention, but why can't you just use comparray = DeleteCases[cmap, _ -> 0];? – Mr.Wizard Mar 6 '13 at 18:16
• @Mr. Wizard Yes, DeleteCases is the better way to do this. I was floundering a bit there and just ran with the first thing that worked for me. – Daniel Lichtblau Mar 6 '13 at 18:50

OK, you can get the Voronoi diagram using raster-graphics tricks, but what if you want to do it the old-fashioned way?

img = Import["http://i.stack.imgur.com/iab6u.png"];
morph = MorphologicalComponents[img];
boundary =
N[{{0, 0}, {#1, 0}, {#1, #2}, {0, #2}} & @@ (ImageDimensions[img] +
1)];
comps = (boundary[[-1]] + {1, -1} Reverse[#] & /@
Position[morph, #]) & /@ Range[Max[morph]];


Take the hulls of the solid shapes. You don't have to do this, but it should make things faster.

hulls = Parallelize[
Function[set,
Cases[set,
x_ /; Count[
MemberQ[set, x + #] & /@
DeleteCases[Tuples[{-1, 0, 1}, {2}], {0, 0}], True] <= 5]] /@
comps];


Deploy ComputationalGeometry! Start with the Delaunay triangulation, and then get the bounded diagram. (I had to cheat here. Due to some badly conditioned matrices, BoundedDiagram blows up if the boundary is too close.)

Needs["ComputationalGeometry"]
del = DelaunayTriangulation[pts = N[Join @@ hulls]];
vor = BoundedDiagram[
boundary /. {0. -> -400., x_?Positive :> x + 400}, pts, del];


Now we just need to paste these back together. This function will fuse adjacent polygons.

fuse[p_, q_] :=
Module[{l = Intersection[p, q], al, nal, qrot, prot, qdrop},
al = Alternatives @@ l;
nal = Except[Alternatives @@ l]; {qrot, prot} =
RotateLeft[#, (-Length[l] +
Position[
Differences@Position[#, Alternatives @@ l] -
1, {_Integer?
Positive}] /. {} -> {{Min[
Position[#, Alternatives @@ l]] - 1}})[[1,
1]]] & /@ {DeleteDuplicates[q], p};
If[Length[DeleteDuplicates[q]] == Length[l],
qrot = Take[prot, Length[l]]];
qdrop = Drop[RotateLeft[qrot], Length[l] - 2];
prot /. {al .., A : nal ...} :>
Join[If[Take[prot, Length[l]] =!= Take[qrot, Length[l]], Identity,
Reverse][qdrop], {A}]]

mergevor = {vor[],
MapIndexed[
Function[{span, i}, {i[],
Fold[fuse, First[#], Rest[#]] &@
SortBy[vor[[2, Span @@ (span + {1, 0})]],
N@Norm[pts[[#[]]] - pts[[vor[[2, span[] + 1, 1]]]]] &][[
All, 2]]}],
Partition[Prepend[Accumulate[Length /@ hulls], 0], 2, 1]]}


And we're done:

Graphics[MapThread[{Hue[#1/Length[hulls]], Opacity[0.5],
Polygon[mergevor[[1, #3]]], Opacity, PointSize[0.01],
Point[#2]} &, {Range[Length[hulls]], hulls,
mergevor[[2, All, 2]]}],
PlotRange -> ({0, #} & /@ ImageDimensions[img])] • I like the approach. Did you notice a slight difference of your result and that of whuber? – DeeDee Mar 6 '13 at 0:33
• I hope not; they should be identical. If I flip back and forth between our final images, I don't see any difference. – Xerxes Mar 6 '13 at 0:36
• One difference is that whuber's method will handle Voronoi cells with holes; mine does not. My method will give you a list of Polygons when you're done; his will need further processing. – Xerxes Mar 6 '13 at 0:39
• It is not the same as mine. Any ideas why? – DeeDee Mar 6 '13 at 0:50
• I suppose because the Meyer algorithm for finding watersheds is not exactly the same as a Voronoi tessellation. I'm not familiar with the details, but you might be able to find them in this paper. – Xerxes Mar 6 '13 at 0:57

While I cannot match @whuber's simple elegance, I will show a bit of brutishness by using Fast Marching from scratch. This finds distances from a specified boundary. I'll modify it so that, for each pixel, it returns the value of the nearest boundary components.

The code is a bit long but mostly cribbed, from this blog. The only modification is the extra bookkeeping and alteration in the returned result noted above. I include it below the example.

For the example itself, there's not much work for me there either because I cribbed that from @whuber.

img = Import["http://i.stack.imgur.com/iab6u.png"];

comps = MorphologicalComponents[img];

negcomps = -comps;
Timing[voronoi = findNearestIndexC[negcomps];]

(* Out= {2.200000, Null} *)


(Not as fast as @whuber's, but not bad either *)

So let's have a look.

MatrixPlot[voronoi + comps, ColorFunction -> "BrightBands",
ImageSize -> 500] Not too bad. There is a bit of jaggedness which might be from resolving seeming ties in the "wrong" way. or maybe they really should be there, I'm not sure.

--- edit ---

Or, more likely, the jaggedness is due to the nature of the algorithm. This is a way to handle certain diffusion-like problems; that is, it is in effect solving a PDE of some sort. Since all steps are essentially "local" (that is, based directly on what territory we have recently traversed but not on past regions), my guess is we get some jaggedness due to global accumulation of error.

--- end edit ---

Code used:

frozen = -1.;
frozenQ[aa_] := aa < 0.
unseen = 0.;
far = 30000.;
outofbounds = 100000.;
bigstate = 10000;
band = 0.;

Clear[FastCompile];
SetAttributes[FastCompile, HoldAll];
FastCompile[stuff__] :=
Compile[stuff,
CompilationOptions -> {"InlineCompiledFunctions" -> False,
"InlineExternalDefinitions" -> True}, RuntimeOptions -> "Speed",
CompilationTarget -> "C"];

state = FastCompile[{{states, _Real,
2}, {x, _Integer}, {y, _Integer}},
If[x > Length[states] || x < 1 || y > Length[states[[x]]] || y < 1,
bigstate, states[[x, y]]]];

distance =
FastCompile[{{distances, _Real, 2}, {x, _Integer}, {y, _Integer}},
If[x > Length[distances] || x < 1 || y > Length[distances[[x]]] ||
y < 1, outofbounds, distances[[x, y]]]];

neighborValue =
FastCompile[{{l1, _Integer, 1}, {l2, _Integer, 1}, {states, _Real,
2}, {distances, _Real, 2}},
Module[{s1, s2, d1, l11, l12, l21, l22}, {l11, l12} = l1;
{l21, l22} = l2;
s1 = state[states, l1[], l1[]];
s2 = state[states, l2[], l2[]];
d1 = distance[distances, l1[], l1[]];
Which[s1 >= 0. && s2 >= 0., outofbounds, s1 <= -1. && s2 <= -1.,
Min[distance[distances, l1[], l1[]],
distance[distances, l2[], l2[]]], s1 <= -1.,
distance[distances, l1[], l1[]], True,
distance[distances, l2[], l2[]]]]];

distanceToBoundary2 =
FastCompile[{v1, v2,
f}, (Sqrt[(-f^2)*(-2 + f^2*(v1 - v2)^2)] + f^2*(v1 + v2))/(2*f^2)];

distanceToBoundary1 = FastCompile[{v1, f}, v1 + 1/f];

newDistance =
FastCompile[{{x, _Integer}, {y, _Integer}, {states, _Real,
2}, {distances, _Real, 2}},
Module[{up, down, left, right, f = 1., res, xvalue, yvalue},
up = {x, y + 1}; down = {x, y - 1}; left = {x - 1, y};
right = {x + 1, y};
xvalue = neighborValue[right, left, states, distances];
yvalue = neighborValue[up, down, states, distances];
res = Which[xvalue == yvalue == outofbounds, outofbounds,
xvalue != outofbounds && yvalue != outofbounds,
distanceToBoundary2[xvalue, yvalue, f], xvalue != outofbounds,
distanceToBoundary1[xvalue, f],
True, distanceToBoundary1[yvalue, f]];
res]];

findNearestIndexC =
FastCompile[{{ll, _Real, 2}},
Module[{hindex = 0, dist, j1, j2, nbrs, pt, x, y, x1, y1, x2, y2,
next, prev, done, cond = False, len, wid, hsize, distancetable,
statetable, statetable2, heaptable, bandheap},
len = Length[ll];
wid = Length[ll[]];
hsize = len*wid;
distancetable = ll;
statetable = Map[If[TrueQ[# == unseen], far, #] &, ll, {2}];
statetable2 = ll;
heaptable = Table[0, {len}, {wid}];
bandheap = Table[{0., 0., 0.}, {hsize}];
Do[If[statetable[[ii, jj]] >= 0., Continue[]];
nbrs = {{ii, jj + 1}, {ii, jj - 1}, {ii - 1, jj}, {ii + 1, jj}};
Do[{x, y} = nbrs[[kk]];
If[! (0 < x <= len && 0 < y <= wid && statetable[[x, y]] == far),
Continue[]];
hindex++;
statetable[[x, y]] = band;
statetable2[[x, y]] = statetable2[[ii, jj]];
dist = newDistance[x, y, statetable, distancetable];
distancetable[[x, y]] = dist;
bandheap[[hindex]] = {dist, N[x], N[y]};
j1 = hindex;
While[(j2 = Floor[j1/2]) >= 1 &&
bandheap[[j2, 1]] > bandheap[[j1, 1]],
bandheap[[{j1, j2}]] = bandheap[[{j2, j1}]];
{x1, y1} = Round[Rest[bandheap[[j1]]]];
heaptable[[x1, y1]] = j1;
j1 = j2;];
heaptable[[x, y]] = j1, {kk, Length[nbrs]}], {ii, len}, {jj,
wid}];
While[hindex > 0, pt = bandheap[];
{x, y} = Round[Rest[pt]];
statetable[[x, y]] = frozen;
bandheap[] = bandheap[[hindex]];
done = False;
prev = 1; next = 1;
{j1, j2} = 2*prev + {0, 1};
While[j1 < hindex && ! done,
If[j2 < hindex,
If[TrueQ[bandheap[[j1, 1]] <= bandheap[[j2, 1]]], next = j1,
next = j2], next = j1];
cond = bandheap[[prev, 1]] > bandheap[[next, 1]];
If[TrueQ[cond],
bandheap[[{prev, next}]] = bandheap[[{next, prev}]];
{x1, y1} = Round[Rest[bandheap[[prev]]]];
heaptable[[x1, y1]] = prev;
prev = next;
{j1, j2} = 2*prev + {0, 1};
, done = True];
];
{x1, y1} = Round[Rest[bandheap[[prev]]]];
heaptable[[x1, y1]] = prev;
nbrs = {{x, y + 1}, {x, y - 1}, {x - 1, y}, {x + 1, y}};
Do[{x2, y2} = nbrs[[kk]];
If[! (0 < x2 <= len && 0 < y2 <= wid &&
statetable[[x2, y2]] == band),
Continue[]];
dist = newDistance[x2, y2, statetable, distancetable];
distancetable[[x2, y2]] = dist;
statetable2[[x2, y2]] = statetable2[[x, y]];
j1 = heaptable[[x2, y2]];
bandheap[[j1]] = {dist, N[x2], N[y2]};
While[(j2 = Floor[j1/2]) >= 1 &&
bandheap[[j2, 1]] > bandheap[[j1, 1]],
bandheap[[{j1, j2}]] = bandheap[[{j2, j1}]];
{x1, y1} = Round[Rest[bandheap[[j1]]]];
heaptable[[x1, y1]] = j1;
j1 = j2;];
heaptable[[x2, y2]] = j1, {kk, Length[nbrs]}];
hindex--;
Do[{x2, y2} = nbrs[[kk]];
If[! (0 < x2 <= len && 0 < y2 <= wid &&
statetable[[x2, y2]] == far),
Continue[]];
hindex++;
statetable[[x2, y2]] = band;
dist = newDistance[x2, y2, statetable, distancetable];
distancetable[[x2, y2]] = dist;
statetable2[[x2, y2]] = statetable2[[x, y]];
bandheap[[hindex]] = {dist, N[x2], N[y2]};
j1 = hindex;
While[(j2 = Floor[j1/2]) >= 1 &&
bandheap[[j2, 1]] > bandheap[[j1, 1]],
bandheap[[{j1, j2}]] = bandheap[[{j2, j1}]];
{x1, y1} = Round[Rest[bandheap[[j1]]]];
heaptable[[x1, y1]] = j1;
j1 = j2;];
heaptable[[x2, y2]] = j1, {kk, Length[nbrs]}];];
statetable2(*distancetable*)]];


I think this works too.

i = Import["http://i.stack.imgur.com/iab6u.png"];
cn = ColorNegate[i];
iws = Image[WatershedComponents[cn]];
ImageMultiply[cn, iws] • I believe this method was already covered in whuber's answer before you posted, was it not? – Mr.Wizard Mar 6 '13 at 0:33
• I see now. I didn't notice it before. – DeeDee Mar 6 '13 at 0:37
• You can delete it if you wish but no one is forcing your hand. Usually we try not to post answers that duplicate earlier ones but it does happen; further, sometimes the simply-stated answer is best, and obviously five people liked this answer. Mostly I wondered if you felt this was a different approach to the one whuber showed. – Mr.Wizard Mar 6 '13 at 0:41
• By the way the result is slightly different. E.g. the border between "S" element in the upper left corner and the element above the pentagon is slightly different than in other solutions. And so is border around the ellipse. Which result is more accurate? – DeeDee Mar 6 '13 at 0:47
• Your method appears to be in error. The key difference is whuber's use of the DistanceTransform prior to using the watershed method. – Xerxes Mar 6 '13 at 1:30

Using ImageMesh[] to convert the image into polygons, along with a strategy similar to this answer, here is how to generate a fake Voronoi diagram for the OP's image:

img = Import["http://i.stack.imgur.com/iab6u.png"];
polys = Flatten[MeshPrimitives[#, 2] & /@ ConnectedMeshComponents[ImageMesh[img]]];
rdf = RegionDistance /@ polys;

(* Schlick's "bias" function,
bias[a_, t_] := t/((1/a - 2) (1 - t) + 1)

DensityPlot[bias[0.6, (HarmonicMean[#] - First[#]) &[
MinimalBy[Through[rdf[{x, y}]], Identity, 2]]],
{x, 5, 451}, {y, 13, 287}, AspectRatio -> Automatic,
ColorFunction -> GrayLevel, Exclusions -> None,
PlotPoints -> 105, PlotRange -> All] As it turns out, there is an easy way to get a MeshRegion corresponding to the Voronoi diagram derived from the OP's image. The method, as described in Okabe et al.'s "Spatial Tessellations", is very easily implemented in Mathematica (altho I had to use a bunch of undocumented routines):

img = Import["http://i.stack.imgur.com/iab6u.png"];
polys = Flatten[MeshPrimitives[#, 2] & /@
ConnectedMeshComponents[ImageMesh[img, Method -> "LinearSeparable"]]];

ptset = polys /. Polygon -> Sequence;
vm = VoronoiMesh[Flatten[ptset, 1]];

facs = MeshPrimitives[vm, 2];
cif = RegionMeshMeshMemberCellIndex[vm];

vmp = DiscretizeGraphics[Graphics[GraphicsPolygonUtilsPolygonCombine[facs[[#]]] & /@
InternalPartitionRagged[cif[Flatten[ptset, 1]][[All, -1]], Length /@ ptset]]] Show Voronoi diagram with the generators:

Show[vmp, Graphics[{Red, polys}], PlotRange -> {{5, 451}, {13, 287}}] 