如何在 MATLAB 中可视化球体的交集？

Question

似乎这个问题已经在几个地方被问过（including on SO）。我最近在可视化trilateration 问题的结果时遇到了这种需求。

在几乎所有情况下，答案都会引导查询查看Wolfram for the math，但不包括任何代码。数学确实是一个很好的参考，但如果我问一个关于编程的问题，一些代码也可能会有所帮助。 （如果对代码问题的回答避免像“编写代码是微不足道的”之类的精辟问题，这当然也值得赞赏。

那么如何在 MATLAB 中可视化球体的交集？我在下面有一个简单的解决方案。

Answer 1

我为此编写了一个小脚本。随时提出建议和编辑。它通过检查每个球体的表面是否落在所有其他球体的体积内来工作。

对于球体相交，在sphere() 函数调用中使用更多面会更好（但更慢）。这应该在可视化中给出更密集的结果。对于单独的球体可视化，较小的数字（~50）就足够了。请参阅 cmets 了解如何可视化每个。

close all
clear
clc

% centers   : 3 x N matrix of [X;Y;Z] coordinates
% dist      : 1 x N vector of sphere radii

%% Plot spheres (fewer faces)
figure, hold on % One figure to rule them all
[x,y,z] = sphere(50); % 50x50-face sphere
for i = 1 : size(centers,2)
    h = surfl(dist(i) * x + centers(1,i), dist(i) * y + centers(2,i), dist(i) * z + centers(3,i));
    set(h, 'FaceAlpha', 0.15)
    shading interp
end

%% Plot intersection (more faces)
% Create a 1000x1000-face sphere (bigger number = better visualization)
[x,y,z] = sphere(1000);

% Allocate space
xt = zeros([size(x), size(centers,2)]);
yt = zeros([size(y), size(centers,2)]);
zt = zeros([size(z), size(centers,2)]);
xm = zeros([size(x), size(centers,2), size(centers,2)]);
ym = zeros([size(y), size(centers,2), size(centers,2)]);
zm = zeros([size(z), size(centers,2), size(centers,2)]);

% Calculate each sphere
for i = 1 : size(centers, 2)
    xt(:,:,i) = dist(i) * x + centers(1,i);
    yt(:,:,i) = dist(i) * y + centers(2,i);
    zt(:,:,i) = dist(i) * z + centers(3,i);
end

% Determine whether the points of each sphere fall within another sphere
% Returns booleans
for i = 1 : size(centers, 2)
    [xm(:,:,:,i), ym(:,:,:,i), zm(:,:,:,i)] = insphere(xt, yt, zt, centers(1,i), centers(2,i), centers(3,i), dist(i)+0.001);
end

% Exclude values of x,y,z that don't fall in every sphere
xmsum = sum(xm,4);
ymsum = sum(ym,4);
zmsum = sum(zm,4);
xt(xmsum < size(centers,2)) = 0;
yt(ymsum < size(centers,2)) = 0;
zt(zmsum < size(centers,2)) = 0;

% Plot intersection
for i = 1 : size(centers,2)
    xp = xt(:,:,i);
    yp = yt(:,:,i);
    zp = zt(:,:,i);
    zp(~(xp & yp & zp)) = NaN;
    surf(xt(:,:,i), yt(:,:,i), zp, 'EdgeColor', 'none');
end

这里是insphere 函数

function [x_new,y_new,z_new] = insphere(x,y,z, x0, y0, z0, r)
    x_new = (x - x0).^2 + (y - y0).^2 + (z - z0).^2 <= r^2;
    y_new = (x - x0).^2 + (y - y0).^2 + (z - z0).^2 <= r^2;
    z_new = (x - x0).^2 + (y - y0).^2 + (z - z0).^2 <= r^2;
end

---

示例可视化

对于这些示例中使用的 6 个球体，在我的笔记本电脑上运行组合可视化平均需要 1.934 秒。

6 个球体的交点：

实际 6 个球体：

下面，我将两者结合起来，这样您就可以在球体视图中看到交叉点。

对于这些示例：

centers =

   -0.0065   -0.3383   -0.1738   -0.2513   -0.2268   -0.3115
    1.6521   -5.7721   -1.7783   -3.5578   -2.9894   -5.1412
    1.2947   -0.2749    0.6781    0.2438    0.4235   -0.1483

dist =

    5.8871    2.5280    2.7109    1.6833    1.9164    2.1231

我希望这可以帮助其他可能希望看到这种效果的人。