3D Voxel mask from geometry/ does plane intersect with volume

Hi Does anyone know if there is a way to efficently convert from geometry (for example a quadrilateral defined by its corners) to a binary voxel mask. I want somthing like the function poly2mask but which will work in 3D and return true voxel values if the input plane or line touches the voxel.

My ultimate aim is to use this with the output from voronoi() to generate x-ray CT like data for finite element modelling.

I've seen this: https://uk.mathworks.com/matlabcentral/fileexchange/37863-blended-3d-poly2mask But i want to be able to cope with polys that arent aligned with the voxel planes.

Answer 1

[Ahmet Cecen]

Well, the first problem you define is slightly more problematic, but the actual task you mention is easier.

A Voronoi tessellation is simply an L2 distance boundary. This means if you generate your random centers, you can voxelize it by assigning each voxel to the closest center.

It would look something like:

[X,Y,Z] = meshgrid(xgv,ygv,zgv)
Space = [X(:),Y(:),Z(:)];
IDX = knnsearch(CENTERS,Space);
Voxelized = reshape(IDX,sizeofVolume);

[Michael]

Thanks a lot for your answer. I need just the boundaries between the Voronoi cells so you end up with something analogous to a closed cell foam (or an open cell foam if you use just the edges of the boundaries. I suppose for the closed cell foam I could make a labelled image and off set it from its self or use diff() to find transitions. But I'm still stuck on the open cell version.

[Ahmet Cecen]

This sounds like something you might need to explain to me in detail. Provide detailed description of this problem with figures, clearly explaining what open cell / closed cell foam is.

Only thing I can say for now is, you can find the boundaries of such a matrix where each region is labeled differently, by using a diff operator. Basically find the locations where the gradient is positive.

You can do this by using either 'imfilter', or 'circshift', or 'gradient'.

[M.S. Khan]

[X,Y,Z] = meshgrid(xgv,ygv,zgv)
Space = [X(:),Y(:),Z(:)];
IDX = knnsearch(CENTERS,Space);
Voxelized = reshape(IDX,sizeofVolume);

where are values of xgv, ygv and zgv?

could you plz explain the terminology of this coding to me. will be very thankful

[George Abrahams]

Hi @M.S. Khan. What @Ahmet Cecen provided is a simple method of generating a Voronoi diagram in a 3D volume. Hopefully my example below is clearer.

% Options.
numberOfVoronoiRegions = 10;
volumeSize = [ 100, 150, 200 ];
% Create an Nx3 matrix of all voxel coordinates in the volume. Each row
% gives the subscripts for a voxel.
volumeXvalues = 1 : volumeSize(1);
volumeYvalues = 1 : volumeSize(2);
volumeZvalues = 1 : volumeSize(3);
[ X, Y, Z ] = ndgrid( volumeXvalues, volumeYvalues, volumeZvalues );
voxelCoordinates = [ X(:), Y(:), Z(:) ];
% Randomly select some voxels as Voronoi seeds, i.e., the center of the 
% regions or cells.
rng( 1 )
voronoiSeedCoordinates = ...
    datasample( voxelCoordinates, numberOfVoronoiRegions );
% Create the Voronoi diagram by find the closest Voronoi seed to each 
% voxel, and then reshaping this to again form a volume.
indexOfClosestSeed = knnsearch( voronoiSeedCoordinates, voxelCoordinates );
voronoiDiagram = reshape( indexOfClosestSeed, volumeSize );
% Find the boundaries where the label value changes, i.e., the edges.
% Uncomment the line below if you want the "outer walls" of the volume.
% voronoiDiagram = padarray( voronoiDiagram, [1 1 1], 0 );
[ gX, gY, gZ ] = gradient( voronoiDiagram );
voronoiBoundaries = ( gX .^2 + gY .^2 + gZ .^2 ) > 0;
% Show the Voronoi boundaries.
isosurface( voronoiBoundaries )
axis tight equal