用Matlab从表面网格创建3D二进制矩阵/图像

Question

如何在Matlab中从表面网格创建3D二进制矩阵/图像？

例如，当我使用以下代码创建椭球时：

[x, y, z] = ellipsoid(0,0,0,5.9,3.25,3.25,30);

X、Y 和 Z 都是大小为 31 x 31 的二维矩阵。

function Create_Mask_Basedon_Ellapsoid3()
         close all
        SurroundingVol  = [50, 50, 20];
        %DATA
        [MatX,MatY,MatZ] = meshgrid(-24:1:25, -24:1:25, -9:1:10);
        [mask1, x, y, z] = DrawEllipsoid([0, -10, 0], [6, 3, 3], MatX,MatY,MatZ);
        [mask2, x2, y2, z2] = DrawEllipsoid([15, 14, 6], [6, 3, 3], MatX,MatY,MatZ);
        mask =  mask1 + mask2;

        %Surface PLOT
        figure('Color', 'w');
        subplot(1,2,1);

        %help: Ideally I would like to generate surf plot directly from combined mask= mask1 + mask2;
        s = surf(x,y,z); hold on;
        s2 = surf(x2,y2,z2); hold off;     
        title('SURFACE', 'FontSize', 16);
        view(-78,22)

        subplot(1,2,2);
        xslice = median(MatX(:));
        yslice = median(MatY(:));
        zslice = median(MatZ(:));

        %help: Also how do I decide correct "slice" and angles to 3D visualization.
        h = slice(MatX, MatY, MatZ, double(mask), xslice, yslice, zslice)
        title('BINARY MASK - SLICE VOLUME', 'FontSize', 16);
        set(h, 'EdgeColor','none');   
        view(-78, 22)
        %az = 0; el = 90;
        %view(az, el);

    end

    function [mask, Ellipsoid_x, Ellipsoid_y, Ellipsoid_z] = DrawEllipsoid(CenterEllipsoid, SizeEllipsoid, MatX, MatY, MatZ)
        [Ellipsoid_x, Ellipsoid_y, Ellipsoid_z] =  ellipsoid(CenterEllipsoid(1), CenterEllipsoid(2), CenterEllipsoid(3), SizeEllipsoid(1)/2 , SizeEllipsoid(2)/2 , SizeEllipsoid(3)/2 ,30);
        v = [Ellipsoid_x(:), Ellipsoid_y(:), Ellipsoid_z(:)];                       %3D points
        %v = [x(:), y(:), z(:)];                                                    %3D points
        tri = DelaunayTri(v);                                                       %triangulation
        SI = pointLocation(tri,MatX(:),MatY(:),MatZ(:));                            %index of simplex (returns NaN for all points outside the convex hull)
        mask = ~isnan(SI);                                                          %binary
        mask = reshape(mask,size(MatX));                                            %reshape the mask 
    end

- Garima Singh

---

你是否只是想要绘制plot3(x,y,z)或其他东西？您可以使用cat(3,x,y,z)，其中3表示3D连接。 - Vuwox

我无法理解你的解决方案。我想知道如何从椭球体制作一个3D面罩。 - Garima Singh

你所说的“mask”是什么意思？过滤器？我不确定理解了。 - Vuwox

所谓“掩码”，是指对于所有位于椭球体内部的坐标（x，y，z），掩码（x，y，z）= 0；否则，掩码（x，y，z）= 1。 - Garima Singh

@Alexandre Bizeau，您可能会对下面的凸包解决方案感兴趣，它适用于任何类型的体积——这是来自Gnovice的一种相当聪明的解决方案。 - marsei

Answer 1

这是您需要的:

%// Points you want to test. Define as you need. This example uses a grid of 1e6
%// points on a cube of sides [-10,10]:
[x y z] = meshgrid(linspace(-10,10,100)); 
x = x(:);
y = y(:);
z = z(:); %// linearize

%// Ellipsoid data
center = [0 0 0]; %// center
semiaxes = [5 4 3]; %// semiaxes

%// Actual computation:
inner = (x-center(1)).^2/semiaxes(1).^2 ...
      + (y-center(2)).^2/semiaxes(2).^2 ...
      + (z-center(3)).^2/semiaxes(3).^2 <= 1;

对于网格的第n个点，其坐标为x(n)，y(n)，z(n)，如果该点位于椭球体内部，则inner(n)为1，否则为0。

例如：绘制内部点：

plot3(x(inner), y(inner), z(inner), '.' , 'markersize', .5)

- Luis Mendo

---

我喜欢使用椭球的参数形式 :) 好主意。 - Vuwox

@AlexandreBizeau 实际上，我在 help ellipsoid 中看到了它 :-) - Luis Mendo

感谢@ADonda提出“线性化”这个术语（http://stackoverflow.com/questions/19408612/insert-new-values-into-an-array/19408989#comment28770801_19408989） - Luis Mendo

Answer 2

这是一种从椭球创建二进制掩模的方法。它创建一个相应的体积，并将椭球外的点设置为NaN（内部为ones）。

它不考虑椭球的公式，而是使用凸包。实际上，对于任何可以正确描述为3D凸包的体积，它都可以工作。在这里，convexhulln步骤被绕过，因为椭球已经是凸包。

所有功劳归功于将凸包转换为二进制掩模。

以下是绘图：

是由...生产的

%DATA
[x, y, z] = ellipsoid(0,0,0,5.9,3.25,3.25,30);

%METHOD
v = [x(:), y(:), z(:)];                       %3D points
[X,Y,Z] = meshgrid(min(v(:)):0.1:max(v(:)));  %volume mesh
tri = DelaunayTri(v);                         %triangulation
SI = pointLocation(tri,X(:),Y(:),Z(:));       %index of simplex (returns NaN for all points outside the convex hull)
mask = ~isnan(SI);                            %binary
mask = reshape(mask,size(X));                 %reshape the mask 

%PLOT
figure('Color', 'w');
subplot(1,2,1);
s = surf(x,y,z);
title('SURFACE', 'FontSize', 16);
view(-78,22)
subplot(1,2,2);
xslice = median(X(:));
yslice = median(Y(:));
zslice = median(Z(:));
h = slice(X, Y, Z, double(mask), xslice, yslice, zslice)
title('BINARY MASK - SLICE VOLUME', 'FontSize', 16);
set(h, 'EdgeColor','none');
view(-78,22)

---

多个椭球

如果您有多个椭球，可以对每个椭球使用此遮罩方法，然后使用&组合生成的遮罩。

切片和角度选择

"正确"是个人选择的问题。您可以：

• 创建未旋转的遮罩，然后在旋转后进行旋转（在matlab中旋转3D数组）。
• 在已旋转的椭球上创建遮罩。
• 在稍微旋转的椭球上创建遮罩（这样可以选择"正确"的切片），然后将其进一步旋转到最终位置。

- marsei

---

谢谢。这正是我在寻找的。不过还有一个问题。我怎样才能让每个体素只有一个单位？也就是说，如果我想要我的椭球体是6 x 3 x 3体素^3，并且位于100 x 100 x 20体素^3的周围体积中心为（i，j，k）的位置，该怎么办？ - Garima Singh

一种方法是使用 ellipsoid(0,0,0, 6/2 , 3/2 , 3/2 ,30) 定义一个椭球体，然后使用 meshgrid(-49:50, -49:50, -9:10)（-49:50 共100个元素）来生成网格。由于网格的每个体素都是单位体素，因此您的椭球体（x_length=6）将在x轴上具有6个正元素。 - marsei

非常感谢您愿意提供帮助。 我尝试根据您的建议编写了代码。请查看我编辑过的原始帖子中的代码。 然而，我仍有两个问题： 1）如果我有两个椭球体，根据您的建议将其转换为3D数据点。如何使用它来显示正确大小的周围环境？ 2）还有，我该如何确定正确的“切片”和角度以进行3D可视化？ 我的代码中这两个问题都标记为%help:。如果您能回答这些问题，我将不胜感激。 - Garima Singh