如何在 MATLAB 中实现 3D 布尔运算以像在 Blender（或任何其他 3D 软件）中那样创建交集？
How to implement a 3D boolean operation in MATLAB to make an intersection like in Blender (or any other 3D software)?
我正致力于从 3 个单独的二进制图像创建 3D 图像，这些图像是用 3 个相机拍摄的。我有相应的校准并知道设置（见下文）。由于图像预处理主要是在 MATLAB 中完成的，所以我想在那里实现所有内容。
我的代码目前的想法是根据相机标定来拉伸二维二值图像。这是典型的二进制图像的样子：
在 MATLAB 中，拉伸图像如下所示：
在所有 3 个相机都被挤压并采用合并算法后，我可以创建最终的 3D 形状。到目前为止效果很好，但需要很长时间来计算，因为我需要创建很多挤压步骤才能获得良好的表面。
我现在想通过重新创建我将在 Blender 等 3D 建模软件中执行的过程来加快速度。在那里我还挤出了二值图像的轮廓并通过为轮廓创建样条线轻松创建交集，挤出它们并使用布尔运算符。这是一个带有 2 个挤压图像的 Blender 示例：
我不知道如何在 MATLAB 中实现类似的东西。我想在挤压 "tube" 的顶端和底端创建我的二进制轮廓的两个实例，然后定义各个点之间的面，然后创建一个交点。点创建没问题，但面定义和交集（布尔运算符）是。有人知道如何实施吗？
这在 MATLAB 中可能不是一件容易的事，但它是可能的。我将在这里概述一组步骤，以两个相交的圆柱体为例...
创建四面体网格：
第一步是为挤压创建一个四面体网格。如果您正在挤压的 2D 二值图像是凸面的并且没有孔，您可以使用 delaunayTriangulation 函数执行此操作：
DT = delaunayTriangulation(P);
此处，P 包含挤压的 "end caps" 的坐标点（即管子两端的面）。但是，在生成四面体网格时，delaunayTriangulation 不允许您指定 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，但我无法访问它们，也无法说出它们的适用性。
如果自动网格生成选项不起作用，您必须自己构建四面体网格并将该数据传递给 triangulation。这可能很棘手，但我会向您展示一些步骤，说明如何对圆柱体执行此操作，这可能会帮助您找出更多复杂的形状。下面，我们构建了一组坐标点 P1 和一个 M×4 矩阵 T1，其中每一行都包含 P1 行的索引，这些行定义了一个四面体：
% 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);
为了更好地形象化圆柱体是如何被分成四面体的，下面是分解图的动画：
现在我们可以创建第二个圆柱体，偏移并旋转使其与 x 轴对齐并与第一个圆柱体相交：
% 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');
下面是 visualize them 的情节：
寻找相交区域：
一旦我们有了体积的四面体表示，我们就可以 generate a grid of points covering the region of intersection and use the pointLocation 函数来确定哪些点在两个圆柱体内：
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]));
我们现在有体积数据 mask，其中包含零和一，这些数据定义了相交区域。您制作的网格越精细（通过调整 nGrid），这将更准确地代表圆柱体之间的真实交叉区域。
正在生成 3D 表面：
您可能希望根据该数据创建一个曲面，定义交叉区域的边界。有几种方法可以做到这一点。一种是用 isosurface, which you could then visualize using featureEdges 生成表面。例如：
[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');
结果图：
另一种选择是创建一个 "voxelated" 类似 Minecraft 的表面，如我所示 :
[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');
结果图：
我正致力于从 3 个单独的二进制图像创建 3D 图像，这些图像是用 3 个相机拍摄的。我有相应的校准并知道设置（见下文）。由于图像预处理主要是在 MATLAB 中完成的，所以我想在那里实现所有内容。
我的代码目前的想法是根据相机标定来拉伸二维二值图像。这是典型的二进制图像的样子：
在 MATLAB 中，拉伸图像如下所示：
在所有 3 个相机都被挤压并采用合并算法后，我可以创建最终的 3D 形状。到目前为止效果很好，但需要很长时间来计算，因为我需要创建很多挤压步骤才能获得良好的表面。
我现在想通过重新创建我将在 Blender 等 3D 建模软件中执行的过程来加快速度。在那里我还挤出了二值图像的轮廓并通过为轮廓创建样条线轻松创建交集，挤出它们并使用布尔运算符。这是一个带有 2 个挤压图像的 Blender 示例：
我不知道如何在 MATLAB 中实现类似的东西。我想在挤压 "tube" 的顶端和底端创建我的二进制轮廓的两个实例，然后定义各个点之间的面，然后创建一个交点。点创建没问题，但面定义和交集（布尔运算符）是。有人知道如何实施吗？
这在 MATLAB 中可能不是一件容易的事，但它是可能的。我将在这里概述一组步骤，以两个相交的圆柱体为例...
创建四面体网格：
第一步是为挤压创建一个四面体网格。如果您正在挤压的 2D 二值图像是凸面的并且没有孔，您可以使用 delaunayTriangulation 函数执行此操作：
DT = delaunayTriangulation(P);
此处，P 包含挤压的 "end caps" 的坐标点（即管子两端的面）。但是，在生成四面体网格时，delaunayTriangulation 不允许您指定 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，但我无法访问它们，也无法说出它们的适用性。
如果自动网格生成选项不起作用，您必须自己构建四面体网格并将该数据传递给 triangulation。这可能很棘手，但我会向您展示一些步骤，说明如何对圆柱体执行此操作，这可能会帮助您找出更多复杂的形状。下面，我们构建了一组坐标点 P1 和一个 M×4 矩阵 T1，其中每一行都包含 P1 行的索引，这些行定义了一个四面体：
% 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);
为了更好地形象化圆柱体是如何被分成四面体的，下面是分解图的动画：
现在我们可以创建第二个圆柱体，偏移并旋转使其与 x 轴对齐并与第一个圆柱体相交：
% 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');
下面是 visualize them 的情节：
寻找相交区域：
一旦我们有了体积的四面体表示，我们就可以 generate a grid of points covering the region of intersection and use the pointLocation 函数来确定哪些点在两个圆柱体内：
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]));
我们现在有体积数据 mask，其中包含零和一，这些数据定义了相交区域。您制作的网格越精细（通过调整 nGrid），这将更准确地代表圆柱体之间的真实交叉区域。
正在生成 3D 表面：
您可能希望根据该数据创建一个曲面，定义交叉区域的边界。有几种方法可以做到这一点。一种是用 isosurface, which you could then visualize using featureEdges 生成表面。例如：
[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');
结果图：
另一种选择是创建一个 "voxelated" 类似 Minecraft 的表面，如我所示
[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');
结果图：