# Efficient Lloyd Sampling of Images

I want to sample a grayscale image so that every voronoi cell contains approximately the same total intensity using lloyd sampling.

My current code is kind of slow and I was hoping for some advice to improve my current implementation.

The setup:

img = ColorConvert[
DensityPlot[
Exp[-0.2 (x^2 + y^2)], {x, -5, 5}, {y, -5, 5},
PlotTheme -> "Monochrome", Frame -> False, PlotRangePadding -> 0,
ImageSize -> 200, PlotPoints -> 100
], "Grayscale"];

SeedRandom["Lloyd"]
seeds = RandomReal[{1, 200}, {20, 2}];
Show[img,Graphics@MeshPrimitives[VoronoiMesh[seeds, {{1, 200}, {1, 200}}], 1]]


My implementation:

1. Calculate the voronoi diagram for the initial seeds.
2. Segment the image according to the voronoi cells.
3. Calculate the intensity centroid for every cell and move the seed to the centroid.

.

InPolyQ[poly_, pt_] := GraphicsPolygonUtilsPointWindingNumber[poly, pt] =!= 0

LloydIteration[img_, seeds_] := Module[{
imgD = First@ImageDimensions@img,
polygones,
polygones = Cases[MeshPrimitives[VoronoiMesh[seeds, {{1, imgD}, {1, imgD}}], 2],
t : Polygon[x_] :> x];
masks = Table[Boole[InPolyQ[#, {i, j}]], {i, 1, imgD}, {j, 1, imgD}] & /@ polygones;
ImageMeasurements[ImageMultiply[Image@#, img], "IntensityCentroid"] & /@ masks]


30 Iteration of LloydIteration yield the following result after ~120 seconds. The result looks promising but the execution time is pretty slow especially with a higher number of seed points.

• Do you need "Voronoi" cells or just "cells". A quick method for partitioning an image in "cells" approximating a Voronoi diagram is shown at the end of the help on Thinning[] Feb 6 '15 at 20:27
• Your algorithm converges to a Voronoi diagram in which every site is the intensity centroid of its cell, which is not the same as a Voronoi diagram in which every cell contains approximately the same total intensity. Observe that even for a constant intensity, the result is a centroidal Voronoi tessellation whose cells can have different areas.
– user484
Feb 6 '15 at 20:44
• @belisarius sorry, but I don't get how this example relates to my problem at all.
– paw
Feb 6 '15 at 21:01
• @Rahul ups, you are right! Will revise!
– paw
Feb 6 '15 at 21:02
• Apparently the weighted CVT works pretty well nevertheless. See the paper "Weighted Voronoi stippling" by Adrian Secord, which I think you will find quite interesting and relevant.
– user484
Feb 6 '15 at 23:07

Edit: At the end of this post, you will find an implementation with Nearest which is as fast as Nikies solution. The unfortunate thing is, that my first idea was to use Nearest but I somehow did not time it correctly. Since I wanted to answer to this comment and show that my implementation is faster, I timed it again - this time correctly - and I have to admit, that Nearest is faster and about a Zillion lines shorter than my parallel implementation.

As nikie pointed out, it is often easily possible that you run into problems when you use things like Rasterize or high-level visualization functions for image processing tasks that require pixel-exact solutions. Theses things are rather made for visually pleasing results.

Let me give an alternative for creating the Voronoi-Mask from a list of seeds which uses non of these functions at all. In fact, it uses only very elementary functions so that we can compile it down and apply it in parallel. It will create results that a visually very similar to Nikie's implementation:

The timing here is 1.9s compared to Nikies 0.7s.

## Idea

The idea is simple and to quote nikie

Building a mask by checking every single pixel for every single polygon is terribly inefficient.

Correct, and still it is exactly what I'm doing here. The complexity is $O(n^2)$ but you'll see that we will reach similar runtime for n=200 as in this example.

In your image, every pixel has a position {x,y}. To apply the voronoi-creation in parallel, it is crucial that we have the matrix of positions up-front. The good thing is that we have to create this exactly once:

makePositions[{nx_Integer, ny_Integer}] :=
DeveloperToPackedArray[Reverse@Table[{x, y}, {y, ny}, {x, nx}]]


The reason for the Reverse is the same as what nikie said in the first paragraph.

When you have a list of seed-points and a position from the above matrix, we go through the seeds-points and calculate its distance to the position. We remember the index of the seed-point with the smallest distance and return it. The following function does exactly this:

parallelVoronoiMask =
Compile[{{point, _Integer, 1}, {seeds, _Real, 2}},
Module[{min = 10.^9, minPos = 0, id = 0, p = {0, 0}},
Do[
id++;
With[{d = #.# &[point - c]},
If[d < min,
min = d;
minPos = id;]
],
{c, seeds}];
minPos
], Parallelization -> True, RuntimeAttributes -> {Listable},
RuntimeOptions -> "Speed"]


Note that this function only takes one position but it is parallelized. Therefore, when we apply it to a matrix of positions, it will automatically be executed for each position in the matrix.

## Main function

The main wrapper function does nothing more than initially creating the matrix of positions and then it uses Nest to iterate as often as we want. The usage of ComponentMeasurements is similar to the implementation of Nikie

Lloyd[img_?ImageQ, seeds_, iter_] := With[{
pos = makePositions[ImageDimensions[img]],
r = Range[Length[seeds]]},
"IntensityCentroid"]) &, seeds, iter]
];

Lloyd[img, seeds, 30]


## Using Nearest

If you want to replace my parallelVoronoiMask with a solution using Nearest then you just have to define

voronoiMaskNearest[ptMat_, seeds_] :=
With[{nf = Nearest[seeds -> Range[Length[seeds]]]},
Map[First, nf[ptMat], {2}]
]


and then replace parallelVoronoiMask with voronoiMaskNearest in the implementation of Lloyd. The overall timing for the example given by the OP is then 0.7s which is exactly as fast as Nikies implementation.

• Calculating the distance to the nearest seed for each pixel, good idea, I didn't think of that! It might get slower if nx and ny get large (10^4? 10^5?), but for 200, it's pretty fast. I don't have time to update my answer right now, but I think Nearest[seeds -> Range[Length[seeds]], pos][[All, All, 1]] is even faster than your Compiled code. Feb 7 '15 at 16:09
• @nikie Hmm, Nearest was my first try but I skipped it because the overall algorithm needed 10s, but I had my own implementation of the ComponentMeasurement part. I guess I forgot to time it correctly because Nearest really is faster (as I would have expected it too). Feb 7 '15 at 16:35

That's because Mathematica array indices are row, column, 1-based, and data[[1,1]] is at the top-left corner of the image, while coordinates are x,y, 0-based and {0,0} is at the bottom-left corner. So the right way to build these masks would be:

masks = Table[
Boole[InPolyQ[#, {j - 1, imgD - i}]], {i, 1, imgD}, {j, 1, imgD}] & /@ polygones


Where {j - 1, imgD - i} is the coordinate for the pixel at i,j. Yes, it's confusing.

Second: Building a mask by checking every single pixel for every single polygon is terribly inefficient. It's much faster to simply draw the polygons:

polygonsToLabels = Function[polygons,
Module[{graphic, n = 255},
graphic =
Rasterize[
Style[Graphics[
MapIndexed[{GrayLevel[#2[[1]]/n], #1} &, polygons],
PlotRange -> {{1, 200}, {1, 200}-1}, ImageSize -> 200],
Antialiasing -> False]];
Round[ImageData[graphic][[All, All, 1]]*n]]];


this generates a "label matrix", where the pixels inside the nth polygon are set to n.

With this, we can use ComponentMeasurements instead of ImageMeasurements to calculate the measurements of all the components at once:

lloydIteration[img_, seeds_] :=
Module[{imgD = First@ImageDimensions@img, polygones, labels},
polygones =
MeshPrimitives[VoronoiMesh[seeds, {{1, imgD}, {1, imgD}}], 2];
labels = polygonsToLabels[polygones];
Range[Length[seeds]] /.
ComponentMeasurements[{labels, img}, "IntensityCentroid"]]


On my PC, this takes about .015s, instead of 5s per iteration for the original version. 30 iterations take less than a second:

(
res = Nest[lloydIteration[img, #] &, seeds, 30];
Show[img,
Graphics@MeshPrimitives[VoronoiMesh[res, {{1, 200}, {1, 200}}], 1]]
) // Timing


It's even fast enough to animate it:

Monitor[Nest[Function[seeds,
temp =
Show[img,
Graphics@
MeshPrimitives[VoronoiMesh[seeds, {{1, 200}, {1, 200}}], 1]];
lloydIteration[img, seeds]], seeds, 250], temp]


To make results comparable, these are the seeds and the results of the first iteration:

SeedRandom["Lloyd"]
seeds = RandomReal[{1, 200}, {20, 2}]


{{33.836, 73.7841}, {90.5691, 79.3203}, {131.345, 191.074}, {109.001, 33.9445}, {108.521, 157.942}, {34.5898, 146.338}, {39.8998, 163.298}, {179.774, 137.209}, {177.902, 169.573}, {64.9409, 56.6131}, {116.162, 102.832}, {44.0173, 110.004}, {42.1332, 174.476}, {2.76295, 88.0033}, {139.134, 145.829}, {16.7267, 130.617}, {194.23, 63.942}, {156.263, 73.6834}, {59.3779, 54.5255}, {50.8815, 36.7287}}

lloydIteration[img, seeds]


{{53.9621, 55.3154}, {42.4279, 142.645}, {9.05126, 84.5329}, {90.519, 81.8095}, {178.08, 126.159}, {114.133, 105.155}, {185.286, 44.604}, {43.3518, 159.704}, {43.5777, 183.988}, {13.6021, 130.025}, {129.24, 186.91}, {177.955, 175.638}, {71.7956, 58.5852}, {139.295, 138.999}, {152.249, 71.7255}, {31.572, 69.1432}, {116.826, 29.0632}, {39.814, 22.3326}, {97.6739, 152.166}, {56.9203, 111.011}}

• Is it possible that your include you own seed with a fixed SeedRandom so that one can compare results? It's strange, but I get different results than the OP when I use "Lloyd". Feb 7 '15 at 9:30
• Thanks for the corrections! Your implementation is much more efficient. There might still be something wrong with the algorithm though. If I increase the seed number to 200 and the iterations to 100 I get a very asymmetric distribution of cells and some cell are clearly moving away from their intensity centroid. The Algorithm also fails with 200 seeds and 200 iterations.
– paw
Feb 7 '15 at 11:56
• @paw: It seems as if Rasterize[Graphics[...]] creates an 8-bit image, and using GrayLevel[index/200] leads to rounding errors. GrayLevel[index/255] seems to work better; If you have more than 255 seeds, you'll have to use RGBColor to encode the seed index. Also, I've noticed that the seed centers seem to be "moving downwards" - I'm guessing Polygon and ComponentMeasurements are rounding in different directions - that's why I've added the (admittedly, ugly) -1` in the PlotRange Feb 7 '15 at 12:22