File Exchange

image thumbnail

Polytope bounded Voronoi diagram in 2D and 3D

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

23 Downloads

Updated 03 Nov 2018

GitHub view license on GitHub

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 persendicular bisecters between any two connected points in the Deluanay triangulation.
Here are the description of the uploads.
"demo.m" an example
"polybnd_voronoi.m" main function that obtains polytope bounded Voronoi diagram
"pbisec.m" obtains half space created with perpendicular bisector of two points in the form Ax <= b

"MY_con2vert.m" convert a convex set of constraint inequalities into the set of vertices at the intersections of those inequalities (written by Michael Keder)

"vert2lcon.m" used for finding the %linear constraints defining a polyhedron in R^n given its vertices (written by Matt Jacobson and Michael Keder)

"inhull.m" tests if a set of points are inside a convex hull (written by John D'Errico)

"MY_setdiff.m", "MY_intersect.m" are much fasten than MATLAB built-in "setdiff.m", "intersect.m". Two functions are written by Nick (http://www.mathworks.com/matlabcentral/profile/authors/1739467-nick)

Cite As

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

Comments and Ratings (33)

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?

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.

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!

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!

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.
Again, thanks for your help.

@ 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'

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

Updates

1.14.0.0

Increase size of m to 50 in the 'demo.m'

1.13.0.0

Minor bug fixes

1.12.0.0

corrected a few errors found in the DEMO file (see GitHub page for details). You will no longer see cell overlapping issues, etc.

1.11.0.0

Title is updated.

1.11.0.0

Image is replaced.

1.1.0.0

Image uploaded.

MATLAB Release Compatibility
Created with R2016a
Compatible with R2016a to R2018b
Platform Compatibility
Windows macOS Linux