2016-11-10 23 views
3

この質問はいくつかの場所で聞かれたようです(including on SO)。私は最近、trilateration問題の結果を視覚化するときに、これが必要であることを知りました。どのように球の交点をMATLABで可視化できますか?

ほぼすべての場合、答えはWolfram for the mathを見るように照会を指示しますが、コードはすべて除外します。数学は本当に参考になりますが、プログラミングに関する質問をする場合は、いくつかのコードも役立ちます。 (コードの質問への回答が、「コードを書くのは簡単です」といった臆病なコメントは避けてください)

どのようにして球の交点をMATLABで可視化できますか?私は以下の簡単な解決策を持っています。

答えて

4

私はちょうどこれを行う小さなスクリプトを書いた。お気軽にご意見をお寄せください。これは、各球の表面が他のすべての球の体積内に入るかどうかをチェックすることによって機能します。

球の交差については、sphere()関数呼び出しでより多くの数の面を使用する方が良い(ただし遅くなります)。これにより、視覚化の結果がより高密度になるはずです。球のみの視覚化では、より小さい数(〜50)で十分です。それぞれを視覚化する方法のコメントを参照してください。

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機能は、これらの例で使用される6球については

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 

サンプル視覚化

あり、それが実行に1.934秒の平均を取りました私のラップトップで組み合わせた視覚化。

6球の交差点: Intersection of 6 spheres

実際の6球:以下 not intersection

、私は2つを組み合わせてきたので、あなたは、球のビューで交差点を見ることができます。これらの例については both

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 

私がこの効果を視覚化することを望むかもしれない他の誰に役立ちます願っています。

+0

非常にクールです。共有してくれてありがとう。 – rayryeng

+0

私はあなたの視覚化がポイントの "削除"のために壊れていると言っても過言であり、あなたは小さな数値エラーを考慮しません。残りは素晴らしいです。 –

+0

@AnderBiguri:ポイントの削除が厳しすぎるということはどういう意味ですか?さらに説明できますか?良いコードを得るためには、このコードが密集した球体を必要とするのは確かです。さもなければ、あなたは正しい、ポイントは抜け出すことができます。 – marcman

関連する問題