Finding the region in a Voronoi diagram that contains a point

Question

I am trying to build a Voronoi diagram using the code in this link. However, I have a few points and want to know in which region they fall. This code, like the original function in MATLAB (i.e. voronoin) gives two outputs: [vornb,vorvx], one for the vertices and another one for the cells. So, I want to see which region of the Voronoi diagram the point (x, y, z) falls in.

I am actually looking for something like this region masking in 3D.

gnovice gnovice · Accepted Answer · 2017-11-28T17:41:27

Regardless of whether you generate your Voronoi cells using the built-in voronoin (which takes an N-by-D matrix of points X as input) or polybnd_voronoi (the linked File Exchange submission for bounded Voronoi cells, which takes an additional M-by-D matrix of points BX defining a bounding convex polyhedron), you can compute which cell contains a given point [x y z] by using only those same input arguments.

For the bounded Voronoi cell case, you first need to determine if your point is inside the bounding convex hull or not. One way to do this is to create a Delaunay triangulation of the bounding points and use the pointLocation method to determine if the point is within the convex hull:

DT = delaunayTriangulation(BX);
cellIndex = pointLocation(DT, [x y z]);

If cellIndex is NaN, then the point is not within the boundary convex hull. Otherwise, it is in one of the Voronoi cells. To determine which, first consider that a Voronoi cell C{i} represents the set of all points that are closer to point X(i, :) than any other point in X. Therefore, to find out which cell a point [x y z] falls into you just have to find the point in X that it is closest to. You can do this using pdist2 and min like so:

[~, cellIndex] = min(pdist2(X, [x y z]));

If you don't have access to pdist2 (which is in the Statistics Toolbox), you can compute the distances yourself like so:

[~, cellIndex] = min(sqrt(sum(bsxfun(@minus, X, [x y z]).^2, 2)));

Now, cellIndex can be used as an index into the output arguments from voronoin or polybnd_voronoi to get the bounding Voronoi cell.

Generalize to multiple points:

You can generalize the above to more than just one point [x y z], which will allow you to create a 3D region mask:

[PX, PY, PZ] = meshgrid(...);  % Generate regular points in a 3D volume
PXYZ = [PX(:) PY(:) PZ(:)];    % Combine them into one matrix

DT = delaunayTriangulation(BX);      % Create triangulation for boundary
cellMask = pointLocation(DT, PXYZ);  % Find points in boundary

index = ~isnan(cellMask);  % Index of points in boundary
[~, cellMask(index)] = min(pdist2(X, PXYZ(index, :)), [], 1);  % Find cell index
cellMask = reshape(cellMask, size(PX));  % Reshape mask to 3D

And the 3D mask cellMask will contain index values for points within Voronoi cells and NaN for points outside the bounding convex hull.

Finding the region in a Voronoi diagram that contains a point

1 Answers

Generalize to multiple points: