File Exchange

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

version 1.14.0.0 (12.6 KB) by Hyongju Park

### Hyongju Park (view profile)

The function cacluates arbitrary polytope bounded Voronoi diagram in 2D/3D

Updated 03 Nov 2018

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 .

Trung Hoang Dinh

### Trung Hoang Dinh (view profile)

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

### Hyongju Park (view profile)

@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

### Binghui Tian (view profile)

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

### Hyongju Park (view profile)

@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

### Hyongju Park (view profile)

@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

### ankit agrawal (view profile)

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

### albert gu (view profile)

@Hyongju Park：Thank you very much!

Hyongju Park

### Hyongju Park (view profile)

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

### albert gu (view profile)

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

### Lukos (view profile)

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

### Hyongju Park (view profile)

@ 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

### Lukos (view profile)

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

### Hyongju Park (view profile)

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

Lukos

### Lukos (view profile)

@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

### Hyongju Park (view profile)

@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

### Hyongju Park (view profile)

@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

### Lukos (view profile)

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

Thanks for publishing these files!

Hyongju Park

### Hyongju Park (view profile)

@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

### William Warriner (view profile)

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

### Hyongju Park (view profile)

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

zhengya si

### zhengya si (view profile)

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

### Hyongju Park (view profile)

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

zhengya si

### zhengya si (view profile)

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

### Zhang Shi (view profile)

@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

### Hyongju Park (view profile)

@ 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

### Zhang Shi (view profile)

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

### Hyongju Park (view profile)

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

### Geeta Monpara (view profile)

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

### mengke yuan (view profile)

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

### Hyongju Park (view profile)

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

### Michael (view profile)

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

### shasha (view profile)

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

Jean-Philippe Morissette

### Jean-Philippe Morissette (view profile)

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