获取属于凸包的点

Posted 2023-02-14

【中文标题】获取属于凸包的点

【英文标题】：Getting the points belonging to the convex hull

【发布时间】：2016-01-02 16:33:31

【问题描述】：

我在 Matlab 中有一个颗粒的二进制图像。我可以通过以下函数找到颗粒的凸包：

[K, V] = convhull(granule);

如何找到原始图像中属于凸包但不属于颗粒的所有像素？我的意思是我想做这样的事情：

granule2 = zeros(size(granule));
granule2(K == 1 & granule == 0) = 2;

它不起作用，因为 K 的大小是 x 乘 3，其中 x 是凸包中三角形的数量。

编辑：根据文档，凸包应该是一个数组，其中点的索引构成每行中的凸包。那么如何找到由这些点确定的体积内的所有点。

Edit2：让我换句话说：我有一个图像，它是一个 3D 点数组。它不是一个球体，它有一些凹痕（所以凸包不会位于我的图像表面上）。

我想找到凸包，然后找到所有在凸包内但在颗粒之外的点。这是它在 2D 中的样子（我想找到红色像素）：

Edit3：NicolaSysnet，您的算法应该返回我图片中所有红色的像素（它们的索引）（图片是二维的，因为它更容易绘制）。

【问题讨论】：

granule 的大小是多少？

您使用的 [K,V] 与 documentation 所说的完全相反。我认为这在您的代码中也是错误的，因为 K==1 在这里是一个相当晦涩的分配

@user2738748 你的问题很混乱。您似乎想要凸包和原始形状之间的像素？答案完全取决于图形的屏幕分辨率、仰角和方位角等因素。由于像素始终是 2D，因此您对 3D 的要求毫无意义。如果您实际上想要船体和您的形状之间的 3D 点，答案是inf。如果您指定点之间的最小间隙，即 3D 网格大小，则答案是有界的。

@user2738748 你没有澄清任何事情。在凸包包围的表面和实际表面之间有无数个点。你要求一个无限的答案。以您的 edit3 为例。有哪些指标？这些点不在颗粒上，所以它们没有索引。它们是船体和空地中的颗粒之间的点。您是否要询问颗粒上不在船体上的所有点？您需要进一步澄清。

@Matt，我的输入是一个 3D 数组，比方说，大小为 1024 x 1024 x 1024。每个点都有三个索引：x、y、z。索引的范围从 1 到 1024。每个点也有一个值 - 如果它是 0，它不属于颗粒，如果它是 1，它。所有点都有三个索引，不仅是颗粒上（或内部）的那些。当我说“点”时，我指的不是平面上的点（数学点），而是原始数组中的一个点（我有超过十亿个这样的点——这是一个有限的数字）。我认为这对任何了解 Matlab 的人来说都是清楚的。

【参考方案1】：

[K, V] = convhull(granule);
granule2 = zeros(size(granule));
tmp = granule(K,:)==0; %// all points on the convex hull but not in the granule
tmp2 = tmp(:,1)==tmp(:,2); %// both indices of your granule are zero
granule2(tmp(tmp2)) = 2;

K 是与凸包上的点相对应的点的行号，V 只是该凸包跨越的体积。因此，您可以使用此行索引来查找粒度中的零索引。

使用下面的例子：

granule = rand(1e3,2);
[K, V] = convhull(granule);
granule2 = zeros(size(granule));
tmp = granule(K,:)<0.1; %// all points on the convex hull but not in the granule
tmp2 = tmp(:,1)==tmp(:,2); %// both indices of your granule are below 0.1
granule2(tmp(tmp2)) = 2;

导致sum(tmp2)=11，因此在这种情况下有 11 个点位于凸包上，并且两个索引都低于 0.1（我无法使用 ==0，因为我的数组中没有零）。

您可能希望根据您的实际需要为tmp2 切换条件。

不出所料，越来越多的人为此苦苦挣扎，并为此编写了代码，请参阅 MathWorks Central，John D'Errico 的代码。

【讨论】：

感谢您的回答。但是你能更详细地解释为什么我需要 tmp2 变量吗？ tmp 变量似乎完全符合我的要求。

tmp 是一个二维逻辑数组，tmp2 是一个一维数组，其中tmp 中每一行的both 项都为真。由于您没有发布任何数据供我使用，我不确定您需要哪个，因此我将其用作更多功能，并添加了您可能希望将 tmp2 切换为您需要的任何内容的注释。

好的，但是为什么我需要执行这一行：tmp2 = tmp(:,1)==tmp(:,2);？这不是多余的吗？我不明白这一点。

该行确保tmp2 填充满足两列条件的值，从而使其成为一维数组。 tmp 是一个二维数组，满足每个值而不是每个行的条件。所以它可能是多余的，这取决于你的数据和条件是什么。由于您仍未发布这些内容，因此我不知道您对数据的期望以及您在这种情况下实际想要实现的目标。我添加这个纯粹是因为答案的多功能性。

我的数据大小为 (1024)^3（它是 3D 格式）。我怎样才能发布它？我的颗粒或多或少是椭圆形的，有一些孔和凹痕。它大约占变量“颗粒”的 5%。

【参考方案2】：

这应该可行：

ix=find(granule);
[x,y,z]=ind2sub(size(granule),ix);
K=convhulln([x,y,z]);
% take the index of the points on the frontier
front=unique([K(:,1);K(:,2);K(:,3)]);

% calculate the minimum box containing the whole granule
minx=min(x(front));
maxx=max(x(front));
miny=min(y(front));
maxy=max(y(front));
minz=min(z(front));
maxz=max(z(front));
% Make a mesh for the points of this box
[X, Y, Z]=meshgrid(minx:maxx,miny:maxy,minz:maxz);
X=X(:);Y=Y(:);Z=Z(:);
% Remove the points of the granule from the box
ind=find(granule(sub2ind(size(granule),X,Y,Z))==0);
X=X(ind);Y=Y(ind);Z=Z(ind);

% For every face
for face=1:length(K),
    % take the coordinates of the vertices of the face
    P1=[x(K(face,1)), y(K(face,1)), z(K(face,1))];
    P2=[x(K(face,2)), y(K(face,2)), z(K(face,2))];
    P3=[x(K(face,3)), y(K(face,3)), z(K(face,3))];

    % calculate the normal to the face (croos product of two vectors on the face, 
    % we take the two sides of the face using the first point as pivot)
    normal=cross(P2-P1,P3-P1);
    % take another point on the convex hull
    inner=1;
    direction=0;
    while direction==0,
        P4=[x(inner), y(inner), z(inner)];

        % calculate the projection of the other point on the normal (always using P1 as pivot) 
        % to verify which is the sign for the projection of the points inside the convex hull
        direction=dot(P4-P1,normal);
        inner=inner+1;
    end
    % do the same for all the remaining points of the mesh
    project=sum(([X,Y,Z]-repmat(P1,size(X,1),1)).*repmat(normal,size(X,1),1),2);
    % if the sign of the projection is different to the one of the test point, 
    % the point is outside the convex hull so remove it
    ind=find(sign(project)*sign(direction)>=1);
    X=X(ind);Y=Y(ind);Z=Z(ind);
    plot3(X,Y,Z,'.')
end

此代码的作用是为每个面验证矩形网格中的每个点是否位于凸包的该面的同一侧。

为了加快计算速度，矩形网格仅针对三个坐标平面上接触颗粒投影的坐标，这样，如果颗粒占据了你所需要的面积的 5%可以裁剪90-95%的点进行检查，无需计算。

【讨论】：

谢谢，这是个好主意。但是颗粒内部而不是表面上的点呢？

在[X, Y, Z]=meshgrid(minx:maxx,miny:maxy,minz:maxz); 中生成的网格包括表面内的点和表面外的所有点，这些点包含在包含整个颗粒的盒子中。因此，颗粒内的点被考虑在内。

@user2738748 我在一个由一个顶点切割的立方体上尝试了这个算法，结果它给了我切割sfere与立方体的交点（没有四面体开始从切出的顶点）。我以为那是你要求的。

交叉点是凹的吗？你能不能发一些照片。在您的情况下，我希望得到不属于切割立方体但位于切割立方体和创建凸包的平面之间的点。你不应该得到属于原始立方体的所有点。

我发布了一张我对你的立方体的期望的图片。

【参考方案3】：

您需要做的是使用 convhull 生成船体的防水表面，然后我看到两个选项：

对于 1024^3 网格中的每个点，检查它是在船体内部还是外部 - 可能使用“多面体测试”（例如 this）

或者，您可以对凸包进行体素化（即光栅化），而不是对每个体素（3d 像素）进行测试 - 即将这个防水表面转换为一组体素。 here 提供了一个选项。

一旦您知道凸包内有哪些点，您就可以轻松地从颗粒中移除这些点。

【讨论】：

检查每个体素是否都在我的凸包内不是一个好主意——我有超过 500 亿个这样的体素。而且您的第二个链接不起作用。

抱歉，应该修复断开的链接。

你也不需要测试属于应该消除相当数量的颗粒的点......很多点......洪水填充应该可以工作。有兴趣了解您的解决方案！

洪水填充是什么意思？我不能使用 imfill 函数，因为凸包是一组未连接的点。

@user2738748 在使用洪水填充之前，您需要先将凸包转换（光栅化）为一组连接点 - 其输出将是逻辑值的 3d 网格（onhull/outside船体）。第二个问题是 imfill 只能在 2d 中工作。看起来我在回答中指出的“voxeliser”可以处理这两个问题。

【参考方案4】：

目前正在研究一个类似的问题，我在 Matlab Exchange 平台上发现了一个名为 inhull 的函数，它可以确定某些 测试点 是在给定凸包内部还是外部。

效果很好，但是如果要测试的点数很大，处理时间可能会很长......