I was recently asked: “How can you find the intersecting volume of two polyhedra?”
R2017b introduced polyshape to MATLAB to make working with two dimensional polygons easy. However, it does not support 3D.
One way to get an approximation to the volume would be to convert the polyhedra to a three dimensional image and find the overlapping voxels. Dirk-Jan’s polygon2voxel helps with the first step of this.
First, let’s create two 3D shapes, one convex, and one not.
xyz = gallery('uniformdata', [30 3], 0)*100; DT = delaunayTriangulation(xyz); figure tetramesh(DT)
xyz = gallery('uniformdata', [30 3], 1)*100; AS = alphaShape(xyz, 40); figure plot(AS)
Now let’s convert these to voxels. I will use just the freeBoundary, i.e. the surface, of the delaunay triamgulation.
[f, t] = freeBoundary(DT); ConvexImage = polygon2voxel(struct('faces', f, 'vertices', t), [100 100 100], 'auto');
For an alphaShape, boundaryFacets will return the surface.
[f, t] = boundaryFacets(AS); ConcaveImage = polygon2voxel(struct('faces', f, 'vertices', t), [100 100 100], 'auto');
We can watch both volumes as a movie.
The holes in the volumes have not been filled; imfill can do that for us.
ConvexImage = imfill(ConvexImage, 'holes'); ConcaveImage = imfill(ConcaveImage, 'holes'); implay([ConvexImage ConcaveImage])
Now logical indexing can be used to find the overlapping volume:
VoxelIntersection = ConvexImage & ConcaveImage; figure isosurface(VoxelIntersection)
The volume could be computed as just the sum of the voxels:
volumevoxels = sum(VoxelIntersection, 'all'); disp("The intersecting volume is: " + volumevoxels + "units")
The intersecting volume is: 272741units
Published with MATLAB® R2019b
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.