Starts with a point set, repeatedly moves each point to centroid of Voronoi cell.
LLOYDSALGORITHM runs Lloyd's algorithm on the particles at xy positions
(Px,Py) within the boundary polygon crs for numIterations iterations
showPlot = true will display the results graphically.
Lloyd's algorithm starts with an initial distribution of samples or
points and consists of repeatedly executing one relaxation step:
1. The Voronoi diagram of all the points is computed.
2. Each cell of the Voronoi diagram is integrated and the centroid is computed.
3. Each point is then moved to the centroid of its Voronoi cell.
Inspired by http://www.mathworks.com/matlabcentral/fileexchange/34428voronoilimit
Requires the Polybool function of the mapping toolbox to run.
Run with no input to see example. To initialize a square with 50 robots
in left middle, run:
lloydsAlgorithm(0.01*rand(50,1),zeros(50,1)+1/2, [0,0;0,1;1,1;1,0], 200, true)
1.6  corrected error about some vertices not being in clockwise error by sorting the vertices 

1.5  Corrected error when called as a function by adding a calculation for n, removed warning about order of vertices by sorting them in CW order, added example call as a function. Thanks for the feedback Kirill Smirnov! 

1.3  When run with no arguments, starts the robots in a small grid pattern. 

1.1  removed dependency on www.mathworks.com/matlabcentral/fileexchange/34428 (which sometimes miscalculates the Voronoi cell boundaries), robots never leave the boundary, reduces to a single file. 
Phuc Ho (view profile)
Ivaylo Pantev (view profile)
Thank you Aaron Becker  the script works a treat! As for the ring domain...I could not see a straightforward way of doing it, so approached in the following way *beware ugly hack to follow*:
1. Generate vertices of regular polygon on the external ring edge (x=ext_radius*cos(t)+x_centre; y=ext_radius*sin(t)+y_centre, where t=linspace(2pi,0,n) and n is the ngon vertices you want to approximate the circle shape).
2. Repeat for vertices of internal ring edge (this time t=linspace(0,2pi,n)!!)
3. For the initial placement of robots  generate points on an (n1)gon on a circle with r=(r_ext+r_int)/2
If you use n big enough (180+) and are willing to live with some errors  this can work. There is surely a better way to do this of course!
Vinh Phu Nguyen (view profile)
Dear Becker,
Thanks for sharing your code. I am wondering what I should do to have a Centroidal Voronoi cells for a ring domain? I do not know how to specify the boundary polygon in this case.
Thanks.
Phu
Aaron T. Becker's Robot Swarm Lab (view profile)
Duc Van Le,
You are almost there. You just need to change the face color of the patches.
After the line:
verCellHandle(i) = patch(Px(i),Py(i),cellColors(i,:));
Add the line:
set(verCellHandle(i), 'FaceAlpha',0.2);
This will make semitransparent patches. You can also make the patches with no face color:
verCellHandle(i) = patch(Px(i),Py(i),cellColors(i,:)); % use color i  no robot assigned yet
set(verCellHandle(i), 'FaceColor','none');
Duc Van Le (view profile)
Dear Becker,
Thank you for your suggestions.
I have done to calculate the centroids of polygons with nonconstant density by applying the trapezoidal method to calculate the integral of density function over polygon area.
Now, I have a problem in Visualization. I want to show the color of each polygon based on the density distribution of its position on the area (similar to the case we use "contour" or "surfc" to plot the density distribution of the area).
Actually, I have added the following code before line 65 (Visualization step in your code) to plot the density function and hope to keep it as background of the figure.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Area of interest
[x,y] = meshgrid(0:1:xrange, 0:1:yrange);
%% Interest density function
z = exp((x800).^2/(2*200^2)(y800).^2/(2*200^2));
contourf(x,y,z);
colorbar
colormap hsv
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
However, when the Lloyd algorithm is running to plot Voronoi cells and their centroids, the background with the density distribution is removed.
Could you give me any advices to solve this problem?
Thank you so much!
Aaron T. Becker's Robot Swarm Lab (view profile)
Duc Van Le, what you would need to do is calculate at "weighted Voronoi digram", and find the center of mass of this partition. This would require changing lines 89 and 102 of this code. I couldn't find any readymade code for this, but you can start with
https://en.wikipedia.org/wiki/Weighted_Voronoi_diagram
for a discussion.
Another place to look is 'power diagrams': http://www.mathworks.com/matlabcentral/fileexchange/44385powerdiagrams
Duc Van Le (view profile)
Duc Van Le (view profile)
Dear Becker,
Fist of all, I would like to thank you and really appreciate your contribution on this work.
Now I am also studying to implement Lloyd algorithm. In my work, it is assumed that the area is given with a distribution density function f(x,y) that is function of position (x,y). So, the centroids of Voronoi cells will depend on function f(x,y).
In my understanding, your work simply assumed that f(x,y) = 1 (correct?) for every position in the area.
I would appreciate if you could give me any advices on how to implement the calculation of centroids with the given density function f(x,y).
Thank you so much!
Aaron T. Becker's Robot Swarm Lab (view profile)
VeronoiLimi, http://www.mathworks.com/matlabcentral/fileexchange/34428voronoilimit
has been updated to handle both internal and external boundaries. If this is needed for your implementation of Lloyd's algorithm, please comment to let me know.
liu (view profile)
Aaron T. Becker's Robot Swarm Lab (view profile)
Corrected error when called as a function by adding a calculation for n, removed warning about order of vertices by sorting them in CW order, added example call as a function. Thanks for the feedback Kirill Smirnov!
Kirill Smirnov (view profile)
I just downloaded this file but when I specify the arguments, I have an error message:
'Undefined function or variable "n".
Error in lloydsAlgorithm (line 63)'.
Because 'n' is initialized only in the 'demo' part of the code. If I put at the beginning something like n = length(Px), everything works fine. Anyway, thank you for your effort!
Aaron T. Becker's Robot Swarm Lab (view profile)
Certainly Lindsey  the code was designed so you can supply the boundary of your environment by supplying the variable crs. Alternatively, to get rid of the white rectangle in the top middle just comment out lines 42,43,44, and 45:
crs = [ 0, 0;
0, yrange;
% 1/3*xrange, yrange; % a world with a narrow passage
% 1/3*xrange, 1/4*yrange;
% 2/3*xrange, 1/4*yrange;
% 2/3*xrange, yrange;
xrange, yrange;
xrange, 0];
Aaron T. Becker
Lindsey (view profile)
Dear All,
First of all, thank you for sharing the Lloyd Algorithm with everyone. I have a question regarding the algorithm. Is it possible to remove the white space in the code?
Thank you and I appreciate if you can help me on this regard,
Lindsey