如何将仿射变换（4x4矩阵）应用于ndgrid/meshgrid结果？

Question

M 是一种仿射变换矩阵，可以将坐标从一个坐标系转换为另一个坐标系，具体如下：

v_old = [ x_old y_old z_old 1];
v_new = M * v_old;
% v_new contains [ x_new y_new z_new 1 ]

现在，我有ndgrid/meshgrid格式的坐标：

[ X_old Y_old Z_old ] = ndgrid( 0:15 ); % for instance

我该如何将它们转换为X_new等？

我可以使用三个 for 循环（X_old(i,j,k) 对应上面的 x_old），但肯定有更好的解决方案。

- Anna

Answer 1

您只需重新排列值，使每个点的四个坐标在一行中：

result = [X_old(:) Y_old(:) Z_old(:) ones(numel(X_old),1)] * M;

每行 result 给出每个点的新坐标。

请注意，这是因为将矩阵 A 乘以矩阵 B 相当于将 A 的每一行分别乘以 B。

如果由于内存限制无法构建上述矩阵，则可使用：

result = X_old(:)*M(1,:) + Y_old(:)*M(2,:) + Z_old(:)*M(3,:) + ones(numel(X_old),1)*M(4,:);

无论哪种情况，如果您想重新排列结果，使其大小与X_old、Y_old、Z_old相同，请将此用作最后一步：

s = size(X_old);
X_new = reshape(result(:,1), s);
Y_new = reshape(result(:,2), s);
Z_new = reshape(result(:,3), s);

- Luis Mendo

---

不错！不过我有一个非常大的矩阵（X_old 是 512x512x300），所以在我还拥有 ndgrid 的结果时，我无法承担创建该矩阵的成本。 - Anna

"afford"（指“在RAM中保留”） - Anna

@Anna，在这种情况下，您可能需要一个循环来构建矩阵的每一行：result(ii,:) = [X_old(ii) Y_old(ii) Z_old(ii) 1] * M;，其中ii是循环变量（从1到numel(X)）。或者使用尽可能多的行来分块构建矩阵。 - Luis Mendo

@Anna 我想我找到了一种方法。请查看更新的答案。 - Luis Mendo

谢谢。我尝试了你的更新，但结果与其他方法不同。第四列甚至不全是1。 - Anna

啊，M 必须被转置。谢谢！ - Anna