File Exchange

image thumbnail


version (10.5 KB) by Jakob Sievers
Constrain the vertices of a Voronoi decomposition to the domain of the input data.


Updated 08 Jun 2020

View Version History

View License

The routine performs a Voronoi decomposition of an input dataset and constrains the vertices to the domain of the data themselves, such that even unbounded Voronoi cells become useful polygons (See attached figure). The function also supports internal/external user-defined boundaries.

Cite As

Jakob Sievers (2021). VoronoiLimit(varargin) (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (88)

sahil bansal

Very robust function. Can we get an output of neighborhood cells. List of cells that share boundary with a particular cell. Thanks

Noam Chen

Excellent function!
Is there a way to receive the cells in the same order as the input points?
Noam Chen

Etienne Jessen

Jakob Sievers

@Ehud Haimov As indicated in the "help" of the script, the vertice-indices within V for each individual cell can be found within the cell structure of the output C. That is: V(C{ij},:) gives you the xy coordinates of all vertices of the ij'th cell. As for the other question, I have tried to recreate the problem and fixed a small related problem (a new updated version has been uploaded) but I wasn't able to recreate the conditions you mention. I.e.: that the four vertices of the boundary are somehow included in the output V. Could you please check (1) if the problem persists with the new version and (2) if so: send me more specific description or input values used so I can recreate it and try to understand it.

Ehud Haimov

And another small question following the former - The output "V" of this code includes the four vertices of the rectangular boundary together with the "real" Voronoi vertices. Can I change something in the code to automatically dispose of these four vertices in "V"? Thank you.

Ehud Haimov

Thank you for this wonderful code. I was wondering, using this code, how can I extract the 'bonds' matrix? i.e. which vertex is connected to which.
Because "V" gives only the coordinates of all vertices but doesn't give any information as to who is connected with whom. Thank you for any answer.

Jakob Sievers

@Milan Sztilkovics: You're welcome. And thank you for pointing this problem out so I had a chance to improve the function.

Milan Sztilkovics

@Jakob Sievers: Works like a charm! Thank you, highly appreciated! :)

Jakob Sievers

@Milan Sztilkovics: I've spent a few hours looking through the code and discovered several minor errors which grew in significance as fewer points were used. I have corrected the errors and have tested it on a number of cases and they all seem to run now. Try downloading the new version ( that I just uploaded.

Milan Sztilkovics

@Jakob Sievers: My problem is that the Voronoi cells are not expanded to the external limiting square. The right edges between the points are there (however sometimes there are invalid additional lines), but the my problem is that the final vertices did not expand to the limiting region. See the attached figures below:
( , )

The example code for these images:
coord = rand(3,2)
bs_ext=[0 0 1 1 0; 0 1 1 0 0]';
for i=1:size(coord,1)

Yes, I use the latest version.

Thanks for your help!

Jakob Sievers

@Milan Sztilkovics: Hm... it comes out right when I run it. Are you sure you're using the most recent version of the code? (version

Milan Sztilkovics

@Jakob Sievers:
Here is the sample code that generates faulty samples. (3 points in this case):
coord = rand(3,2)
bs_ext=[0 0 1 1 0; 0 1 1 0 0]';

Thanks for your quick response!

Jakob Sievers

@Milan Sztilkovics: Could you write or send the xy coordinates of a few such cases so I can investigate?

Milan Sztilkovics

Thanks for the upload, very useful!

However I found that for low number of input points (3,4) the function tends to fail.
(fgures: , , , )
Does anyone has a hint on how to solve this behavior?


Thank you so much.


Hello everyone

Jakob Sievers

@AnindyaG : Thank you. No it does not and I don't envision such extension to the code anytime soon. Cheers

Anindya G

Thank you. This is very useful. Does this also work for the 3-D Voronoi diagrams?

Dimosthenis Floros

Excellent function. It works perfectly on my case of a square mesh with a hole.

Jakob Sievers

Examples of use added to help section for easier application

Jakob Sievers

@Tsinuel Geleta: Thank you for bringing these cases to my attention. I have used them to improve on the edge handling of the script and now it handles all the cases perfectly. Let me know if you experience any more problems. /Cheers

Tsinuel Geleta

I was using your function many times and it worked fine. However, I got the same problem as some other users in which it fails to create some corner and edge polygons. I have put a google drive link to a matlab data file that includes all the cases that failed. If the link fails, I can put some samples here.

Jiaqi Yang

@Jakob Sievers Thanks! I think my problem got solved by the new version.

Jakob Sievers

Does this still happen with the new version that I recently uploaded? If so, please send me the points in questions so I can investigate the problem.

Jiaqi Yang

I wonder could someone solve my problem as follows: When I generate five or less points, and then call VoronoiLimit function on them, it sometimes returns the error: Index exceeds the number of array elements (0).
Error in VoronoiLimit (line 239)
And this problem gets solved when I generate more points, could someone explain to me why this happened and how I can get away from this error only with five points.

Jakob Sievers

@Vanessa Sanchez: Your points and defined boundaries do not overlap. That's your problem.
@Mahmut Okatan: Solved as per version
@Reiner Vanselow: Solved as per version
@Bharath Attoti: No idea. Not with this routine for sure

Jonas Allgeier

Great function, works very well! If it is important to keep the points in the same order (e.q. for plotting with "patch" to predefined colors), it is necessary to add 'stable' to the "unique" command in line 97.


For those who like to use this function in the future in combination with patching of polygons, please note that in the documentation ( the example is stated with with the line `if all(c{i}~=1)`. If you want to make sure all polygons are patched arter this function, remove that line.


After carefully checking results, this indeed changes the sides/corners correctly, but it created another problem. Patching the polygons lead to fine results in sides/corners and most of the midpoints, but for some reason there were a few polygons that weren't patched.


Just want to say this is really useful! I was looking for a way to close the polygons at sides/corners and I came upon your question/answer on the mathworks forum. This gave me exactly what I was looking for.

Bharath Attoti

How can I do it in N dimensions (not just 2d)?

Reiner Vanselow

A useful program, though I found a case where it failed to do Voronoi tessellation properly.

X = P(:,1);
Y = P(:,2);

Error in VoronoiLimit (line 175)
mvx=vx(2,mix); %missing vx

Who can help?

void faceless

Mahmut Okatan

Thanks a lot for this useful program, though I found a case where it failed to do Voronoi tessellation properly.



[X,Y] = ndgrid(XX,YY);

bbox=[0 -1; 29 -1; 29 16.5; 0 16.5; 0 -1];


Two corner cells turns out be missing.

Vanessa Sanchez

I am a new user of matlab, and have a problem with the function. I used the next lines:
x = [1177.3918,1199.0418, 1212.0418, 1189.2218, 1202.2218, 1179.4118, 1192.4118];
y = [-17.25, -17.25, -17.25, -34.25, -34.25, -51.25, -51.25];

[V,C,XY]=VoronoiLimit(x,y,'bs_ext',[0 0; 0 -40; 13 -40; 13 0], 'figura', 'on');

The result is only the figure, but the contents of V, C and XY are empty. There are in my workspace but without data.
Please help me.

Kaijun Ren

yo lecat

This is very usefull thank you for the work.
However, I am using your code to demodulate a noisy signal. For low modulations as QPSK, if there is a lot of noise, the constellation points are received beyond the boundaries of the polygon. I can't figure out how to put the boundaries far enough, it seems to be have a limit. Can someone help me on this?
Thank you


yo lecat

Martin Hofmann

Really useful, thanks a lot!!

Jakob Sievers

Hi all, sorry for the silence. Replies:

@Laura Veldenz:
The help text has been corrected. The right way to initiate the function is:

@Ignacio Bordeu and @Michal Kvasnicka:
It seems that the old algorithm linked the wrong points for a very small number of the cases. I have added some extra basic requirements. Namely that (1) each proposed cell contains the original centroid (duh) and (2) that no side of a cell can intersect with any other part of the cell (again.. duh). These extra requirements seem to have corrected the problem for your specific data points. It also might have slowed the function a bit but I'll look at that if anyone are experiencing a serious slowdown. Please try it out and get back to me if it fails for some reason.

Simple. Try running e.g.:
VoronoiLimit(x,y,'bs_ext',[0 0;0 1;1 1;1 0],'figure','on');




Hi, Can anyone help me to understand how to get the vertices for unit square (for ex:defined by x=0,y=0,x=1 and y=1) using this function? How can I define the boundaries as x=0,y=0,x=1 and y=1?

Michal Kvasnicka

I tried to contact Jacob Sievers directly with this problem, but he is not ready to make any bug fixing now.

Ignacio Bordeu

@Michal Kvasnicka: I am trying to find a solution, the problem exists only for certain sets of (xc,yc) points and region sizes r.

Michal Kvasnicka

@Ignacio Bordeu: I can confirm this problem on my data sets. Any help?

Ignacio Bordeu

Hey guys, I have found an issue when considering certain sets of points. When trying the following example:

xc = [25.0362, 619.2890, 137.8061, 387.3650, 2.7294, 85.5450]

yc = [643.6328, 8.6653, 621.3243, 55.4894, 545.6431, 193.8842]

r = 500;

You can see that the limits are not drawn correctly, this issue disappears for smaller values of r.. any ideas on how to fix it?

Laura Veldenz

there is a difference between the function and the comment about the function (first line in the comments), as seen below, the input is in different orders:

function [V,C,XY]=VoronoiLimi(varagin)


so which one is correct? Or does it not make any difference>


Michal Kvasnicka

good piece of code ...
well done!!!


Thank you for this great function!
Unfortunately, I sometimes get this error:

Subscript indices must either be real positive integers or logicals.

Error in VoronoiLimit (line 152)
mvx=vx(2,mix); %missing vx

I used the unit square as an external boundary. Any idea how I can avoid this error? Thanks in advance!

Morten Rishøj

Works as advertised, simply put.

Matt J

Just a question. Are the polygons obtained by VoronoiLimit guaranteed to be the same as what you get if you take the unbounded Voronoi cells (as computed by VORONOIN, for example) and intersect them with a large rectangle?

Paul Martin

Excellent, thanks Jakob,

it worked just fine for me!

I used the convex hull of my data as a boundary, as follows:

% Get convex hull of data
K=convhull(X, Y);

% Turn the hull into X,Y boundary values

BX = X(K); BY = Y(K);

% Do the special Voronoi limited by the data
[V , C, XY] = VoronoiLimit(X, Y, 'bs_ext', [BX; BY], 'plot', 'on');


This fails for me when I use large bounds, for example in the given examples using the following bounds:

world_bounds_x = [-8,10];
world_bounds_y = [-4,14];
bs_ext = [world_bounds_x([1 1 2 2 1]);world_bounds_y([1 2 2 1 1])]';

Also, when shrinking this bound, it sometimes makes slightly innaccurate bounds (i.e. one of the bounded-rectangle corners is slightly shifted).


Jakob Sievers

Hi Preeti

Please send me an email with a description of the problem along with the specific xy coordinates you use (.mat file) so that I may have a look at it and figure out what the problem is.


Preeti Goel

I wanted to get the points as per the external boundary, I got bugs with some examples. It was fixed when I used VoronoiBounded from the Lloyds Algorithm(Px,Py, Crs, Num Iterations, Show Plot).

Preeti Goel

I believe there is a bug in this section of the code

if length(C{ij})<3

for one of my examples:
with 3 points
one voronoi vertex in ext boundary
3 end points computed
the voronoi cells are not computed correctly and all are the same. Please help. I wanted to add files, but I do not see an option to do that here, I can send them an email id if needed. Thanks

Jakob Sievers

Hi Zu
I have added a triangle example to the current version of the file. Please download and test to see how it works. You can use any form of external boundary as long as it contains all of the x,y input points.

Zu Mailand

could you give me an example? the external boundary is belonging to triangle.

Jakob Sievers

The input variable "bs_ext" which is responsible for the external boundary, takes any number of vertices, meaning that you could just as easily describe a circle or whatever else you need :)


Hi Jakob,

The current limited domain is only a rectangular or a square. If the limited domain is a circle or a triangle, can we modify more for your Matlab code? Thanks

Jakob Sievers

Hi francesco
I'm not sure I understand your question. The function supports both internal an external bounds such as illustrated (e.g. run with no input to see how it works). Is the bound you are thinking of not covered by any of these options?


Thanks Jakob. I figured out :)

We should use dt=DelaunayTri(x(:),y(:)).

A question, if I want a bound domain which is rectangular, triangle or any more complicated geometry. How could I do for your function?

Jakob Sievers

Hi Francesco. The DelaunayTri function will be removed from Matlab in favor of the DelaunayTriangulation function ( I don't know when the latter was implemented but it might be after the 2012 edition of Matlab, which is why it does not show up in your version of Matlab. Unless you upgrade to current matlab I suggest you use DelaunayTri if that works for you.



Jakob Sievers

Hi Francesco

I forgot to add the link to the information in question:


It seems that it should be the function 'dt=DelaunayTri(x(:),y(:))', right?


I have an error when compiling

Undefined function 'delaunayTriangulation' for input arguments of type 'double'.

Error in VoronoiLimit (line 96)

As checked, I cannot find the function delaunayTriangulation in my Matlab R2012b.

Please help.


Jakob Sievers

One of my comments was removed so I'll reiterate briefly: three points is at the lower limit of when delaunayTriangulation can actually provide voronoi cells. if length(x)>3, however, cells begin to form.

Jakob Sievers

Follow-up comment: You may be able to solve this by: (1) deriving the centroid of the three points, (2) extrapolating each point based on this centroid far beyond the external bounds you have specified to add three extra points, (3) using all 6 pointsto create the voronoi cells. The program should remove the outer three points automatically. This is just off the top of my head so get back to me if you try this and it works or doesn't.


Hi Jakob
That problem has been sovled. However, I find another bug. When you just choose three points as the generator, for example, (-1,0), (1,0) and (0,1) are generators, and choose the boundary as 'boundary_x=[-2;2;2;-2;-2];boundary_y=[-2;-2;2;2;-2];'. The program cannot give the right result. At this stage, I have not found the reason and I am trying to solve it. Do you know why?

Yuhang Wang

Jakob Sievers

New version uploaded. Please let me know if this solves the problem for you!

Jakob Sievers

AH. I just noticed I had made a "range" function back in the day and its output was a vector of size=[1 2] where the values where min/max of the input vector respectively. This explains everything. Matlab must have used my "range" function instead of the one native to the program. I will fix this immediately and upload a new version ASAP. Stay tuned for the update.


Hi Jakob
Another question. Could one vertex connect with three or more infinite vertices? If so, it seems that the program does not consider this situation. Is that right? Thank you.


Hi Jakob
I have check the result of the code. It seems that "range function gives the difference between the minimum and maximum of a vector". So the code "rx=range(x);" just gets rx of 1 times 1 dimension. But if I change this code to "rx=[min(x),max(x)];", the program could run. Is that right?

Jakob Sievers

Dear Michele
This is all a bit puzzling to me. When run without input x and y are both vectors. Consequently rx=range(x) and ry=range(y) gives two-element min/max ranges (size=[1,2]) and bnd=[rx ry]; becomes a four-element vector (size=[1,4]) on which crs is based. I cant figure out what might be going wrong based on your (and Yuanzhe's) description. Could you elaborate a bit? Are you running with no input or with user-defined input, and if so; what kind?

Michele Guidi

I downloaded the last verion today since this function would really help me a lot and I found the same error of Yuanzhe.
It is in the creation of the crs matrix (lines 77-82), where using the function RANGE(X), with X being a vector, simply return the range between the max and min value and not a vector, so it attempt to access an index out of bound.
I tried to fix the problem myself but made some avalanche consecutive errors.

Jakob Sievers

Dear Yuanzhe
I do not get that error, but I have been doing a number of bugfixes these past few days so maybe if you try downloading the current version you won't experience that error? Please get back to me in case that didn't work for you.


Dear Jakob
When run the program without input, there is an error. As x and y are vectors, bnd only has two elements and bnd(4) does not exist. So how to solve it? Thanks.


Excellent function. I agree with Michele that an edge constraint functionality would be great to handle internal polygonal edges in the domain, e.g. plate with a hole cutout. The default DT = DelaunayTri(points,EC) function has the EC input variable to handle this. From this Delaunay Triangulation, the Voronoi tesselation can be calculated. This could be a good starting point.


Dear Jakob,
There is only a small bug in the code you may resolve that. for some Voronoi cells (specially those near edges) vertices are not correct. Would you mind to check that.

Good luck,


Jakob Sievers

Hi Michele
That is an interesting case. I am very limited timewise at the moment but when I do return to this function I will definitely take it into consideration. Could be interesting to add options for empty spaces of arbitrary polygon shapes within the bounding box.


Great function!
I have a question though. Is it possible to add new internal bounds, for example a "hole" inside a square region?


Great function !
I upgraded it for using it with any bounding shape. Just replacing "crs" variable by any polygon.
Works fine !



I 2nd Richard. I'm trying to get extended boundary, but still have no luck. Let me know if you have chance to upgrade the code. However, I don't have the mapping toolbox either :"(

Jakob Sievers

I have indeed considered adding an option to describe extended boundaries yourself. I will see if I get the time to make those changes.
I am not aware of any file-exchange alternatives to what polybool does, so I'm afraid the mapping toolbox remains a requirement.

Richard Garner

Thanks for this, as I have run into the same problem, with a slight variation. I want the Voronoi diagram limited by a larger boundary outside of the original set of vertices from which the Voronoi diagram is calculated. HOWEVER, unfortunatey I cannot use your function because I do not have the mapping toolbox.

MATLAB Release Compatibility
Created with R2019a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!