File Exchange

## Polytope bounded Voronoi diagram in 2D and 3D

version 1.15 (12.6 KB) by
The function cacluates arbitrary polytope bounded Voronoi diagram in 2D/3D

Updated 03 Jun 2020

View Version History

Polytope-bounded-Voronoi-diagram

This is a MATLAB script

What is this for?

The function calculates Voronoi diagram with the finite set of points that are bounded by an arbitrary polytope. The Voronoi diagram is obtained using linear ineqaulities formed with perpendicular bisecters between any two connected points in the Deluanay triangulation.

Description

File name Description
demo.m an example script
polybnd_voronoi.m main function that obtains polytope bounded Voronoi diagram
pbisec.m a function computes perpendicular bisectors of two points
MY_con2vert.m inequality constraints to set of vertices (written by Michael Keder)
vert2lcon.m a function is used to find linear inequalities from a polyhedron
(written by Matt Jacobson and Michael Keder)
inhull.m a test function to see if a set of points are inside some convex hull
(written by John D'Errico)
MY_setdiff.m, MY_intersect.m fuctions which are much faster than MATLAB built-in functions
(written by Nick, see http://www.mathworks.com/matlabcentral/profile
authors/1739467-nick)
Note: This is still for experimental use ONLY!

### Cite As

Hyongju Park (2021). Polytope bounded Voronoi diagram in 2D and 3D (https://github.com/hyongju/Polytope-bounded-Voronoi-diagram/releases/tag/1.15), GitHub. Retrieved .

Emma De Cocker

@Hyongju Park, thank you for your reply. Unfortunately, the second output of the convhull command seems to return the volume 1 convex hull when I want the volume of each polyhedron. I am new to Matlab and I am not sure I understand fully. Does the function create polyhedron of the same volume? Thank you!

Hyongju Park

@Emma De Cocker: Hi. It looks like you can use Matlab's convhull function to obtain the volume of the given convex hull from the function's 2nd output. For this see, e.g., https://www.mathworks.com/help/matlab/ref/convhull.html

Emma De Cocker

Hello @Hyongju Park ! Great code and very useful output! I have adapted your code for mine and I need the volumes of each polyhedrons. Is there a way to obtain these? Thanks a lot

Alessio Trebbi

@Hyongju Park. Thanks for your interest. I figured out a solution. I first check if two faces are parallel by looking if their orthogonal vectors are parallel. Then I check if they have two points in common. And finally I check if the two respectively third points are on the same half plane. If this is the case then the face is removed.

Hyongju Park

@Alessio Trebbi: Sorry about the late reply. Yes, I am aware of the issue, and on top of my head, I cannot think of any elegant method of locating the repeated faces either. One dirty way to do is to iterate over all faces to generate list of unique faces as represented as points in CC/CCW order.

Alessio Trebbi

Hi Hyongju! Amazing work! I ask you a suggestion:
This code generates polyhedrons and where they have a face in common, this face is repeated twice (One for each polyhedron). In particular, if there is a quadrangular face of vertices 1 2 3 4, one polyhedron will have his faces 1 2 3 and 3 4 2, and the other one 1 3 4 and 1 2 4. Since the two polyhedrons are using different combinations of points, its hard to find when there is a repeated face. I would like to remove these repeated faces. How would you do it?
Thank you a lot

Hansong Ji

Thanks for the code. There are some irreguar cases that some patches violates the boundary constrainswe when I use it.
This bug comes from the function MY_con2vert, the code "c=A\B" cannot get the proper solution sometimes.
I fix this proble by make the seedpoint as the initial value of c.

Trung Hoang Dinh

Dear Mr.Park
I would like to create voronoi diagram for a ellipsoid with axes: a=b=1.2, longer axis c=2.
Could you please give me some clues how to do it? I am just started studying Matlab few days ago. I will be really helpful for me if you can give me suggestions.
Thank you very much.

Hyongju Park

@Binghui Tian: Hi. Thanks for the comment and the suggestion. I will be testing this more later to correct the issue; but you can also try changing the "tol" to be a smaller value then the current 1e-07. Would you mind posting one of your example codes here which has resulted in the irregular case?

Binghui Tian

Thanks for the code, But when i use it, there are some irregular cases that some patches violates the boundary constrainswe, I think we should check the value of favl not ef in "MY_con2vert.m" ([c,fval,ef] = fminsearch( @(x)obj(x, {A,b}),c,options)); and if favl is not 0, do FMINSEARCH again, two or three times will be good.

Hyongju Park

@ankit agrawal: (cont') one thing I would do to manually identify the *true* voronoi neighbor; is to obtain all the non-boundary edges (that is not the part of the arbitrary polytope) from the given voronoi cell (e.g., 20), and check if each a line connecting from the node 20 to neighbour (16,12,30,26, 27,22) is perpendicular to the any given non-boundary edge. If it does then the given neighbor is *true* voronoi neighbor (in this case it would be 16,30,12,26) and *invalid* voronoi neighbor (22,27).

Hyongju Park

@ankit agrawal: Hi. Actually both 22 and 27 do share their edges with cell 20; but this is not obvious from within polytope (two edges 22-20, 27-20 should probably meet far away from the polytope). I think the notion of Voronoi neighbor only holds for those *unbounded* Voronoi diagrams (where some of the cells can be unbounded). Since by its definition Voronoi diagram has a dual relationship with Delaunay Triangulation, I believe that if there is an edge between any two generators (nodes) from the Triangulation, the related cells should be neighbors from its Voronoi diagam.

ankit agrawal

Hi Hyongju, Thank you very much for uploading this. Very nice script you wrote. I think there is a problem in identifying voronoi neighbours. The neighbor calculation is extrapolated. I am attaching the image in this link (https://imgur.com/ohEO6M1). According to script cell 20 has 6 neighbour (16,12,30,26, 27,22). But we can see visually that 22 and 27 does share any edge so cannot be neighbour. Blue lines are delaunary triangulation. Black lines are voronoi partition.

albert gu

@Hyongju Park：Thank you very much!

Hyongju Park

albert gu: Hi. The inequality is to tell if a given test point x' in R^d is inside the convex hull (polytope) whose boundary is defined by Ax = b, where x in R^d. In this case, it might be used to test if x' is in the lower (closed) half-space, again defined by Ax <= b.

albert gu

Thanks for your code!! However, I still get confused with the “inequality Ax <= b” in pbisec.m. Could you please give me a hint on it? Thanks a lot!

Lukos

Unfortunately not.

Besides a volume exceeding the boundary, sometimes a volume within the boundary is not plotted/displayed.

This is most likely related to each other. Thanks for helping me out!

Hyongju Park

@ Luuk Altenburg: Hi. I have just updated the script "MY_con2vert.m" based on previous comment from William Warriner (see his comments below), and this seems to fix the issue you have mentioned. Please confirm. Thanks!

Lukos

The vertices sometimes exceed the boundaries, creating a faulty trisurf plot. I used a cube defined by 8 points as as the boundary polytope:

bnd0 = [0 0 0;...
0 0 1;...
0 1 0;...
1 0 0;...
1 1 0;...
1 1 1;...
0 1 1;...
1 0 1];

Does anyone know what went wrong? or could it be a bug.

Hyongju Park

@ Luuk Yes, you are correct. I am glad that you were able to figure it out yourself!

Lukos

@Hyongju Thanks for the quick reply! I was just fooling around, but I think the code is easily adjusted to add the value of the volume. In the demo it is in the following section:

case 3
figure('position',[0 0 600 600],'Color',[1 1 1]);
for i = 1:size(pos,1)
[K VVV(i)]= convhulln(vorvx{i}); <<<<<<<<< THIS LINE
trisurf(K,vorvx{i}(:,1),vorvx{i}(:,2),vorvx{i}(:,3),'FaceColor',col(i,:),'FaceAlpha',0.5,'EdgeAlpha',1)
hold on;
end
scatter3(pos(:,1),pos(:,2),pos(:,3),'Marker','o','MarkerFaceColor',[0 .75 .75], 'MarkerEdgeColor','k');
axis('equal')
axis([0 1 0 1 0 1]);
set(gca,'xtick',[0 1]);
set(gca,'ytick',[0 1]);
set(gca,'ztick',[0 1]);
set(gca,'FontSize',16);
xlabel('X');ylabel('Y');zlabel('Z');

I added VVV(i) in the convhull statement. I am not sure if I am correct. What do you think?

Hyongju Park

@Luuk: In fact, you do not need the previous script. The matlab function "convhull" returns the volume for any convex polyhydren as well, see https://www.mathworks.com/help/matlab/ref/convhull.html

Hyongju Park

@Luuk: It is not possible with the current code; but you can try an existing script that can compute the volume of convex bodies from File Exchange at the following URL: https://www.mathworks.com/matlabcentral/fileexchange/43596-volume-computation-of-convex-bodies

Lukos

Is it possible to extract the volume of each Voronion cell with this code?

Thanks for publishing these files!

Hyongju Park

@William: Thanks for your feedback. I am glad that you were able to resolve the issue. I will try to update my code based on your comments!

William Warriner

Hi Hyongju, I sent an email about an error on line 153 of "MY_con2vert.m". Changing the line to: if ~all(abs(A*c - b) < tol) resolved the issue for me. I posted on that submission as well, but I don't hold out hope it will get fixed there. See: https://www.mathworks.com/matlabcentral/fileexchange/7894-con2vert-constraints-to-vertices

Hyongju Park

zhengya si: Yes, you can send me an email. Thanks.

zhengya si

Thanks for your response. Maybe your code is so adaptable that it can also be compatible for lower versions.
In the MPT2/MPT3, it has a function mpt_voronoi which is used for generating voronoi tessellations. If you don't mind, may I send an email to you for discussing about voronoi? I want to show my figures to you.

Hyongju Park

@zhengya si: I don't think the result will differ in R2015b. I have only used MPT2 before; not MPT3.

zhengya si

Thanks for your code. When I run your code in R2015b(Linux), I haven’t got any errors. So I wonder to know whether your codes are also compatible with R2015b or lower version. Do you think the result which obtain with R2015b will different from R2016~2018?
By the way, I am researching about voronoi tessellation. In the beginning, I created the voronoi tessellation by toolbox MPT3. And I can obtain Voronoi diagram in 2D/3D by using MPT3 in Windows. But when I run the same program in Linux, the Voronoi diagram is asunder. Have you used MPT3? And If you have, could you give me some suggestions about my problem?

Zhang Shi

@Hyongju Park, Thanks for your quick response. I used a box boundary in my case such as [1 1 1;1 0 0;0 1 0;1 0 1; 1 1 0;0 1 1;0 0 1;0 0 0; I will try to increase the boundary points and see how it goes.

Hyongju Park

@ Zhang Shi - Try downloading the code&run the demo.m again. The issue is merely due to the smaller number of boundary points (m) relative to the test points (n). I have increased the value m to 50. After the fix, I have ran the demo 100 times for both 2D and 3D and encountered 0 error. Hope this addresses your issue.

Zhang Shi

Thanks for your code. But I am always having this error : on-bounding constraints detected. (Consider box constraints on variables.) Error in polybnd_voronoi (line 75)
V{i}= MY_con2vert(Aaug{i},baug{i});

Could you help me on how to fix it?
Thanks

Hyongju Park

Geeta Monpara-You can do so by replacing 'bnd0 = rand(m,d)' with, e.g., 'bnd0 = [0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1]' in the 'demo.m'

Geeta Monpara

Thank you for the code. I was wondering if there is any way to fit voronoi tessellation produced bounded in a cube? I am interested in using generated vertices to create a finite element mesh.

mengke yuan

The package is perfect. But when i use it to visualize the 3d voronoi diagram, there are some irregular cases that some patches violates the boundary constrains, Some times it gives "Error in My_con2vert(line 170) Non-bounding constraints detected."

Hyongju Park

Michael: thanks for the feedback! In fact, I have corrected the issue you have mentioned and others (e.g., Voronoi cells near the boundary not displaying correctly); however failed to merge to the master. Things are all working now as it should. Also, I will really appreciate if you can test with the updated files.

Michael

I think the functionality of fminserch has changed, line 153 of MY_con2vert needs changing in 2016b(at least) from:
[c,~,ef] = fminsearch(@obj,c,'params',{A,b});

to:

objf=@(c)obj(c,{A,b});
[c,~,ef] = fminsearch(objf,c);

Other wise you get this error:
Error using fminsearch (line 105)
Argument 3 must be an options structure.

shasha

There're bugs in these programs. Some tiles are covering other tiles, and one tile or maybe more are missing in the vorvx.

This is exactly what I was looking for. Thank you very much.
I'd like to know if this works even if my polygon isn't convex (I'm only using the 2D part of the code currently.)

##### MATLAB Release Compatibility
Created with R2016a
Compatible with R2016a and later releases
##### Platform Compatibility
Windows macOS Linux