空間図形の可視化（その1）球と平面の交線：3D空間の円

はじめに

プログラミングは、アイディアそれ自体を深め、考えるためのツールとしても役に立ちます。その一つが空間図形の可視化による理解（＝数式の理解）です。MATLABは可視化機能が充実していますし、Symbolic Math Toolboxと組み合わせて使うことで、数式との対応関係が付きやすいので、プログラムに慣れていない中高生や、大学の新入生にも数式が持つ意味を理解するのに役立つと思います。

実行条件

• MATLABとSymbolic Math Toolboxのみ

前回ご紹介したように、MathWorksアカウントを作成してMATLAB Mobileアプリを使う事で無償で何方でも試せます。受験生の皆さんも遊んでみてくださいね。

問題設定

空間に配置された球と平面の交線、円を求める
平面上の円の方程式は良く知られていますが、3次元空間の円の表現はあまり見ません。そのためか、”空間の円の方程式”等のキーワード(英語でも）で良く検索されているようです。

• ベクトルを使った定式化
  • $d$を通り、法線ベクトル$n$を持つ平面上の点$x$は$n\cdot \left(x-d\right)=0$を満たす
  • 中心$c$、半径$r$の球上の点$x$は、$\left|x-c\right|=r$を満たす

この両式を満たす点$x$の集まりが交線となります。

まず球と平面を透過表示してみましょう。
可視化には無名関数と、陰関数（$f(x,y,z)=0$の形式）を表示する関数fimplicit3を使います。$f$と$g$は各々平面と球を表す無名関数です。fimplicit3のオプションはとりあえずおまじないとして実行してみます。

close all; clear

n = [1 3 4]';    % Direction vector of the plane
d = [2 3 2.5]';  % Fix point on the plane
c = [1 1 0]';    % Center of the Sphere
r = 4;           % Radius of the Sphere

f = @(x,y,z) n(1).*(x-d(1))+n(2).*(y-d(2))+n(3).*(z-d(3));
g = @(x,y,z) (x-c(1)).^2+(y-c(2)).^2+(z-c(3)).^2-r^2;

fimplicit3(f,'FaceColor','c','EdgeColor','none','FaceAlpha',0.2);
hold on
fimplicit3(g,'FaceColor','b','EdgeColor','none','FaceAlpha',0.2);
axis equal;grid on
xlabel('X'),ylabel('Y'),zlabel('Z')
hold off
view(136,30)

実行結果

円らしい交線が見えます。ここで、交線自体の計算は全く行なっていないことに注意しましょう。fimplicit3を使うと（関数ハンドルを使うので）、ほとんど数式そのままの記述で、空間図形を簡単に表示できます。さらに透過オプションを使って複数の図形を重ねて表示すると、交線が自然に見えてくるので便利です。

implicit3のオプション設定（透過表示の設定）を簡単にするための関数を定義しましょう。

function h = my_implicit3(f,color,alpha)
h = fimplicit3(f);
% Set ImplicitFunctionSurface properties
h.FaceColor = color;
h.EdgeColor ='none';
h.FaceAlpha = alpha;
end

また後ほど使いますので、準備をかねて、球を定義する関数を関数ハンドルで定義します。無名関数では難しい、ベクトル表記も簡単です。

function g = sphere_implicit(r,c)
syms x y z real
X=[x y z]';
f(x,y,z) = norm(X-c)^2-r^2;
g = matlabFunction(f);
end

ここでは、シンボリック関数を定義した上で、matlabFunctionを使って、関数ハンドルに変換しています。Symbolic Math ToolboxはMATLABの数値計算と文法の記述方法が同じで、MATLABの関数にそのまま変換できるので便利です。

close all; clear

n = [1 3 4]';    % Direction vector of the plane
d = [2 3 2.5]';  % Fix point on the plane
c = [1 1 0]';    % Center of the Sphere
r = 4;           % Radius of the Sphere

f = @(x,y,z) n(1).*(x-d(1))+n(2).*(y-d(2))+n(3).*(z-d(3));
%g = @(x,y,z) (x-c(1)).^2+(y-c(2)).^2+(z-c(3)).^2-r^2;
g = sphere_implicit(r,c);

my_implicit3(f,'c',0.2);
hold on
my_implicit3(g,'b',0.2);
axis equal;grid on
xlabel('X'),ylabel('Y'),zlabel('Z')
hold off
view(136,30)

function h = my_implicit3(f,color,alpha)
h = fimplicit3(f);
% Set ImplicitFunctionSurface properties
h.FaceColor = color;
h.EdgeColor ='none';
h.FaceAlpha = alpha;
end

function g = sphere_implicit(r,c)
syms x y z real
X=[x y z]';
f(x,y,z) = norm(X-c)^2-r^2;
g = matlabFunction(f);
end

交線を求める

さて、先ほどの表示結果から、交線が円だとするとどのように表現できるでしょうか。
円を定めて作図するために必要な情報は、中心と半径です。この円は平面上にあり、中心を$p_c$とすると、$p_c$と球の中心$c$を結ぶベクトルは、平面の法線ベクトルと平行です。したがって、
　$n\cdot \left(p_c-d\right)=0$
　$p_c -c=t\cdot n$
を満たします。この2式から$t$を求めると、$p_c$が得られます。
次に、この円の半径を$r_c$とすると、ピタゴラスの定理
　$r^2 ={r_c }^2 +{\left|p_c -c\right|}^2$
から$r_c$が求められます。

function [pc,rc] = CenterCircle(d,n,c,r)
t = n'*(d-c)/(n'*n);
pc = c+t*n;
rc = sqrt(r^2-norm(pc-c)^2); % Radius of the circle
end

では、円自体はどのようになるでしょう。3次元空間の表示と言っても、この平面上では2次元の円にすぎません。この円（すなわち交線）上の点$p$は、$\theta$をパラメータとして、以下のように表されます¹。

ここで、$u$と$v$は互いに直行するこの平面上の単位ベクトルで、法線ベクトル$n$とも直交しています。どのように定めても良いのですが、

　$u={\left(n_y ,{-n}_x ,0\right)}^T$
　$v=n\times u$

として定めましょう。ここで$\times$は外積を表します。

function p = circle3d(n,c,r,th)
u = [n(2);-n(1);0];
u = u/norm(u);
v = cross(n,u); v = v/norm(v);
p = c + r*cos(th).*u + r*sin(th).*v;
end

最後にパラメータ表現された空間の円を、plot3を使ってプロットします。全体のプログラムと実行結果は以下のようになります。

全体のプログラム

close all; clear

n = [1 3 4]';    % Direction vector of the plane
d = [2 3 2.5]';  % Fix point on the plane
c = [1 1 0]';    % Center of the Sphere
r = 4;           % Radius of the Sphere
th = 0:pi/100:2*pi;

f = @(x,y,z) n(1).*(x-d(1))+n(2).*(y-d(2))+n(3).*(z-d(3));
g = sphere_implicit(r,c);

[pc,rc] = CenterCircle(d,n,c,r);
p = circle3d(n,pc,rc,th);

my_implicit3(f,'c',0.2);
hold on
my_implicit3(g,'b',0.2);
plot3(p(1,:),p(2,:),p(3,:),'r');
axis equal;grid on
xlabel('X'),ylabel('Y'),zlabel('Z')
hold off
view(136,30)

function g = sphere_implicit(r,c)
syms x y z real
X=[x y z]';
f(x,y,z) = norm(X-c)^2-r^2; % Express x,y,z explicitly
g = matlabFunction(f);
end

function [pc,rc] = CenterCircle(d,n,c,r)
t = n'*(d-c)/(n'*n);
pc = c+t*n;
rc = sqrt(r^2-norm(pc-c)^2); % Radius of the circle
end

function [p,u,v] = circle3d(n,c,r,th)
u = [n(2);-n(1);0];
u = u/norm(u);
v = cross(n,u); v = v/norm(v);
p = c + r*cos(th).*u + r*sin(th).*v;
end

function h = my_implicit3(f,color,alpha)
h = fimplicit3(f);
% Set ImplicitFunctionSurface properties
h.FaceColor = color;
h.EdgeColor ='none';
h.FaceAlpha = alpha;
end

実行結果

ちゃんと重なりました！
#まとめ

• 球と平面の交線を可視化した上で仮説を立て、数式使ってプロットすることで仮説を検証した。
  • fimplicit3を使うと、ほとんど数式の通りに陰関数で表された空間図形を可視化でき、重なった図形の交線は数式表現しなくても自然に見える。
  • シンボリック関数はMATLABの文法がそのまま使え、matlabFunctionを使うことで、関数への変換ができる。
  • plot3を使い、パラメータ表示の関数をプロットすることで、正しく交線が求められたか確認できる

1. https://en.wikipedia.org/wiki/Sphere
   　$p\ ＝;p_c+r_c ;\mathrm{cos}\theta ;u+r_c ;\mathrm{sin}\theta ;v$ ↩