Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

share|improve this question
2  
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
    
Related: stackoverflow.com/a/7941325/695132 –  Szabolcs Mar 5 '13 at 21:15
add comment

5 Answers 5

up vote 28 down vote accepted

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 *)

Relief plot

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]]

Image

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 *)

Solution 2

I like this one better, but fear it might take too much processing for large complex images.

share|improve this answer
    
thank you very much for your answer. –  DeeDee Mar 5 '13 at 20:38
3  
(+1) very nice & simple ;) –  Vitaliy Kaurov Mar 5 '13 at 21:07
1  
+1 Enlightened badge –  Mr.Wizard Mar 5 '13 at 21:12
add comment

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[[1]]}
  ];
]

(* Out[138]= {0.640000, Null} *)

The picture:

MatrixPlot[voronoi2 + comps, ColorFunction -> "BrightBands", ImageSize -> 500]

enter image description here

share|improve this answer
    
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
add comment

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 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 returned result noted above. I include it below the example.

For the example itself, 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[84]= {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]

enter image description here

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[[1]], l1[[2]]];
    s2 = state[states, l2[[1]], l2[[2]]];
    d1 = distance[distances, l1[[1]], l1[[2]]];
    Which[s1 >= 0. && s2 >= 0., outofbounds, s1 <= -1. && s2 <= -1., 
     Min[distance[distances, l1[[1]], l1[[2]]], 
      distance[distances, l2[[1]], l2[[2]]]], s1 <= -1., 
     distance[distances, l1[[1]], l1[[2]]], True, 
     distance[distances, l2[[1]], l2[[2]]]]]];

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[[1]]];
    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[[1]];
     {x, y} = Round[Rest[pt]];
     statetable[[x, y]] = frozen;
     bandheap[[1]] = 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*)]];
share|improve this answer
add comment

I think this works too.

i = Import["http://i.stack.imgur.com/iab6u.png"];
cn = ColorNegate[i];
iws = Image[WatershedComponents[cn]];
ImageMultiply[cn, iws]

enter image description here

share|improve this answer
1  
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
1  
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
4  
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
show 1 more comment

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

Download the graphic from the net and parse its components:

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[[1]], 
  MapIndexed[
   Function[{span, i}, {i[[1]], 
     Fold[fuse, First[#], Rest[#]] &@
      SortBy[vor[[2, Span @@ (span + {1, 0})]], 
        N@Norm[pts[[#[[1]]]] - pts[[vor[[2, span[[1]] + 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[1], PointSize[0.01], 
    Point[#2]} &, {Range[Length[hulls]], hulls, 
   mergevor[[2, All, 2]]}], 
 PlotRange -> ({0, #} & /@ ImageDimensions[img])]

Voronoi diagram of squiggles

share|improve this answer
    
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
1  
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
1  
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
show 4 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.