How to implement a 3D boolean operation in MATLAB to make an intersection like in Blender (or any other 3D software)?

Question

I am working on creating 3D images from 3 individual binary images, that were taken with 3 cameras. I have a corresponding calibration and know the setup (see below). Since the image preprocessing is mostly done in MATLAB I would like to implement everything there.

The current idea of my code is to extrude the 2D binary image according to the camera calibration. Here's what a typical binary image looks like:

And an extruded image looks like this in MATLAB:

With all 3 cameras extruded and a binning algorithm, I can create my final 3D-shape. This works fine so far, but takes a long time to compute, since I need to create a lot of extrusion steps to get a good surface.

I was thinking now to speed this up by recreating the process I would do in a 3D-modeling software like Blender. There I also extrude the outline of the binary image and create an intersection easily by just creating a spline for the outline, extrude them and use a boolean operator. Here's a Blender example with 2 extruded images:

I have no idea how to implement something like this in MATLAB. I wanted to create two instances of my binary contour at the top and bottom end of the extrusion "tube" and then define faces between the individual points and create an intersection afterwards. The point creation is no problem but the face definition and the intersection (boolean operator) are. Does anyone have an idea how this could be implemented?

Answer 1

This might not be an easy thing to do in MATLAB, but it's possible. I'll outline one set of steps here, using two intersecting cylinders as an example...

Creating a tetrahedral mesh:

The first step is to create a tetrahedral mesh for your extrusion. If your 2D binary image that you're extruding is convex and has no holes, you could do this using the delaunayTriangulation function:

DT = delaunayTriangulation(P);

Here, P contains the coordinate points of the "end caps" of your extrusion (i.e. the faces on each end of your tube). However, when generating tetrahedral meshes, delaunayTriangulation doesn't allow you to specify constrained edges, and as such it can end up filling in holes or concavities in your extrusion. There may be some better mesh-generation alternatives in other toolboxes, such as the Partial Differential Equations Toolbox, but I don't have access to them and can't speak to their applicability.

If automated mesh-generation options don't work, you'll have to build your tetrahedral mesh yourself and pass that data to triangulation. This could be tricky, but I'll show you some steps for how you can do this for a cylinder, which may help you figure it out for more involved shapes. Below, we build a set of coordinate points P1 and an M-by-4 matrix T1 where each row contains indices into the rows of P1 which define one tetrahedron:

% Create circle coordinates for the end caps:
N = 21;
theta = linspace(0, 2*pi, N).';
x = sin(theta(1:(end-1)));
y = cos(theta(1:(end-1)))+0.5;
z = ones(N-1, 1);

% Build tetrahedrons for first cylinder, aligned along the z axis:
P1 = [0 0.5 -1; ...  % Center point of bottom face
      x y -z; ...    % Edge coordinates of bottom face
      0 0.5 1; ...   % Center point of top face
      x y z];        % Edge coordinates of top face
cBottom = ones(N-1, 1);      % Row indices for bottom center coordinate
cEdgeBottom1 = (2:N).';      % Row indices for bottom edge coordinates
cEdgeBottom2 = [3:N 2].';    % Shifted row indices for bottom edge coordinates
cTop = cBottom+N;            % Row indices for top center coordinate
cEdgeTop1 = cEdgeBottom1+N;  % Row indices for top edge coordinates
cEdgeTop2 = cEdgeBottom2+N;  % Shifted row indices for top edge coordinates
% There are 3 tetrahedrons per radial slice of the cylinder: one that includes the
% bottom face and half of the side face (all generated simultaneously by the first row
% below), one that includes the other half of the side face (second row below), and one
% that includes the top face (third row below):
T1 = [cEdgeBottom1 cEdgeBottom2 cEdgeTop1 cBottom; ...
      cEdgeBottom2 cEdgeTop1 cEdgeTop2 cBottom; ...
      cEdgeTop1 cEdgeTop2 cTop cBottom];
TR1 = triangulation(T1, P1);

To better visualize how the cylinder is being divided into tetrahedrons, here's an animation of an exploded view:

Now we can create a second cylinder, offset and rotated so it aligns with the x axis and intersecting the first:

% Build tetrahedrons for second cylinder:
P2 = [P1(:, 3) -P1(:, 2) P1(:, 1)];
T2 = T1;
TR2 = triangulation(T2, P2);

% Plot cylinders:
tetramesh(TR1, 'FaceColor', 'r', 'FaceAlpha', 0.6);
hold on;
tetramesh(TR2, 'FaceColor', 'g', 'FaceAlpha', 0.6);
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');

And here's the plot to visualize them:

Finding the region of intersection:

Once we have tetrahedral representations of the volumes, we can generate a grid of points covering the region of intersection and use the pointLocation function to determine which points are within both cylinders:

nGrid = 101;
[X, Y, Z] = meshgrid(linspace(-1, 1, nGrid));
QP = [X(:) Y(:) Z(:)];
indexIntersect = (~isnan(pointLocation(TR1, QP))) & ...
                 (~isnan(pointLocation(TR2, QP)));
mask = double(reshape(indexIntersect, [nGrid nGrid nGrid]));

We now have volume data mask that contains zeroes and ones, with the ones defining the region of intersection. The finer you make your grid (by adjusting nGrid), the more accurately this will represent the true region of intersection between the cylinders.

Generating a 3D surface:

You may want to create a surface from this data, defining the boundary of the intersection region. There are a couple ways to do this. One is to generate the surface with isosurface, which you could then visualize using featureEdges. For example:

[F, V] = isosurface(mask, 0.5);
TR = triangulation(F, V);
FE = featureEdges(TR, pi/6).';
xV = V(:, 1);
yV = V(:, 2);
zV = V(:, 3);
trisurf(TR, 'FaceColor', 'c', 'FaceAlpha', 0.8, 'EdgeColor', 'none');
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');
hold on;
plot3(xV(FE), yV(FE), zV(FE), 'k');

And the resulting plot:

Another option is to create a "voxelated" Minecraft-like surface, as I illustrate here:

[X, Y, Z, C] = build_voxels(permute(mask, [2 1 3]));
hSurface = patch(X, Y, Z, 'c', ...
                 'AmbientStrength', 0.5, ...
                 'BackFaceLighting', 'unlit', ...
                 'EdgeColor', 'none', ...
                 'FaceLighting', 'flat');
axis equal;
view(-37.5, 30);
set(gca, 'XLim', [0 101], 'YLim', [25 75], 'ZLim', [0 102]);
xlabel('x');
ylabel('y');
zlabel('z');
grid on;
light('Position', get(gca, 'CameraPosition'), 'Style', 'local');

And the resulting plot: