如何找到平面和3D矩阵的交点平面

Question

如何找到平面和3D矩阵的交点平面

4

如果我有一堆图片，它们的尺寸如下：

size(M) = [256 256 124];

我有三个点，它们的坐标分别是：

coor_a = [100,100,124];
coor_b = [256,156,0];
coor_c = [156,256,0];

如何从由3个点定义的平面与M相交的点中创建点图像？

- Q_1.0

你试过这个吗？http://www.mathworks.com/matlabcentral/fileexchange/32032-extract-slice-from-volume - Huá dé ní 華得尼

1

我把你的image-processing标签改成了geometry，这与图像处理没有任何关系。 - rayryeng

@rayryeng 如果矩阵M是医学图像（CT扫描切片）或高光谱图像（124个唯一波长的带），那么这很可能是一个图像处理问题。而且问题明确要求如何创建切片的图像，尽管在MATLAB中，图像只是一个矩阵。因此，这是一个图像处理和几何问题，尽管它可以被视为纯几何问题。只是随便提一下。 - user3834473

@TommyPKeane - OP没有明确说明这是医学图像切片，所以我删除了标签...虽然我在帖子中看到了“图像”关键字...这就是为什么对我来说有点不清楚。如果OP澄清并且这涉及图像体积，我很乐意把标签放回去...但现在，我会保持它的状态。感谢您的见解！顺便说一句，我+1了你的答案。非常详细和有趣。 - rayryeng

2个回答

1

我拟合了一个平面，然后找到了最近的点。

%This is the inputs you provided
M = floor(256*rand(256, 256, 124));

a = [100, 100, 124];
b = [256, 156, 0];
c = [156, 256, 0];

% change from a set of triplets to x, y and z vectors
points = [a', b', c'];
x = points(1,:)';
y = points(2,:)';
z = points(3,:)';

%Use Matlab's plane fit function
[ft, ~] = fit( [x, y], z, 'poly11' );

%Get the size of the input matrix
[XExt, YExt, ZExt] = size(M);

%These will be the inputs in the plane
X = 0:XExt-1;
Y = 0:YExt-1;
[X, Y] = meshgrid(X, Y);

%Find the points closest to the plane
Z = round(ft(X, Y));

%Get the x, y, z points of the points on the plane which are also in the extent of M
x = X((Z > 0) & (Z < ZExt)) + 1;
y = Y((Z > 0) & (Z < ZExt)) + 1;
z = Z((Z > 0) & (Z < ZExt));

%Preallocate for the points on the plane
MPrime = zeros(size(z));

%Get the points, need the loop because otherwise Matlab throws an error
for i = 1:length(z)
    MPrime(i) = M(x(i), y(i), z(i));
end

%Final output is a N x 3 vector with the row, column, value
NewImage = [x, y, MPrime];

我希望这可以帮助到您。

- user1543042

太好了。但是，你的NewImage是一个三维矩阵，对吧？你能解释一下如何显示它吗？谢谢。 - Huá dé ní 華得尼

2

请解释一下这段代码是如何工作的……如果您不想这样做，请描述预期的输入和输出应该是什么。我们在这里写答案，以便访问此问题的人们可以学习和提高自己。简单地倾泻代码与此目标相悖。 - rayryeng

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user3834473 · Accepted Answer

首先，你可以注意到如果你在一个（三维）欧几里得（代数）空间中，你要做的就是找到平面和格子之间的交点。

你拥有的是矩阵（M），其大小定义了格子。

然后，因为三个非共线点定义了一个唯一的平面（在三维空间中），你有一个唯一的平面，但你必须想出一个更好的描述平面的方法来找到它与格子的交点。

注意：在任何其他操作之前，你需要将你的点重新定义为适当的列向量。

coor_a = [100 ; 100 ; 124];
coor_b = [256 ; 156 ; 0];
coor_c = [156 ; 256 ; 0];

现在，首先，我们可以通过使用你提供的三个平面点中的两个点对来定义一个任意的二维平面向量基。

u = (corr_b - corr_a);
v = (corr_c - corr_b);

现在，u和v定义了平面流形空间的基向量。因此，从概念上讲，您正在尝试找到由u和v的线性组合“到达”的晶格点。这是因为平面实际上是一个基本的线性流形，也可以被视为一个函数：

P = f(u, v) = au + bv

请注意，此处的P是平面上的任意一点。

现在，您已经知道M是一个矩阵晶格，并具有以下定义域集合（集合符号）：

e1 = [0, 255]
e2 = [0, 255]
e3 = [0, 123]

您可以使用以下 MATLAB 代码获取晶格中的所有点：

X, Y, Z = meshgrid( 0:1:255, 0:1:255, 0:1:123 );

一旦您拥有了 `meshgrid` 的索引，您可以使用冒号索引将矩阵展平为向量，并且您可以创建一个新的矩阵 `A`，其大小为 `3 x (m • n • p)`，其中我们已经有了 `M` 大小为 `m x n x p`。

A = [ X(:); Y(:); Z(:) ];

因为您的矩阵M不是很大，所以这不应该造成内存问题。

现在我们可以回顾P，并将P定义为向量符号表示的线性方程（MATLAB符号）：

P = ( (a * u) + (b * v) );

我们想要的是所有满足以下条件的i的值：

A(i) = ( (a * u) + (b * v) );

但这是MATLAB，所以我们希望利用线性代数运算，并在可能的情况下避免迭代进行向量化。

我们正在构建的是平面代数方程的“Hesse Normal Form”。

http://en.wikipedia.org/wiki/Hesse_normal_form

我们可以利用基向量的事实来为平面找到一个法向量：

n = cross( (corr_b - corr_a), (corr_c - corr_b) );

现在缺少的只是到平面的距离d，可以通过原始点向量的平均值的范数来找到：

d = norm( mean( [ corr_a ; corr_b ; corr_c ] ), 2 );

从Hesse形式中，我们得到了n和d，现在缺少的只是r，即平面上格点的点，因此我们的r值是在平面上的A矩阵列。

因此，现在我们可以使用n作为行向量（列向量的转置）进行向量化计算，在MATLAB中非常快地获取每个格点到平面的距离，并利用矩阵乘法的优势。

L = ( ( (n.') * A ) - d );

因为我们使用矩阵乘法，L 变成了一个 1 x (m • n • p) 的数组，其中每个条目是 n 和 A 的列的点积。任何接近或等于 0 的值 - 取决于您想要使用的阈值 - 都是在平面上的格点。

您可以通过使用 find 方法并生成逻辑 NOT (~) 来从 l 中获取 A 的索引，将所有零变成 1，将所有非零值变成 0，其中 find 返回给定矩阵或向量的所有非零值的索引：

i = find( ~L );

有许多其他方法可以做到这一点，但我认为这种方法利用了代数和内置的MATLAB函数。如果你小心的话，实际上可以将所有先前的命令都放在一行MATLAB代码中完成。这将很难阅读，所以最好将它们保留为单独的行。

此外，这可以轻松设置为一个函数，使您输入3个点值和给定的晶格或晶格定义矩阵，然后您可以返回平面上的所有晶格点或返回向量L的距离，以便您可以选择“接近”0的程度是可接受的，具体取决于您手头的问题。

最后，为了明确起见，如果您想创建一个包含平面上所有点的1通道图像，则可以从以下位置获取点数组：

ImPoints = A(:,i);

一旦您这样做了，您可以使用基本的for循环从M中提取每个点，并将M的值放入一些新的图像数组I中，您需要预先用零构造它。

请注意，ImPoints的列是M的行-列-通道索引，这意味着您需要执行类似以下操作：

for j = [ 1 : 1 : NumPoints ]

    NewImIndex = SomeFunc( ImPoints, j );    

    I( NewImIndex ) = M( ImPoints(:,j) );

end;

您需要预定义I并想出您偏爱的方法，让SomeFunc将像素/数据从M放入您的新图像I中。