はじめに
以前の記事で3Dモデリングにおける「法線」というものの意味と大切さについて語りました。
その法線はMATLABの中で自動的に定義される場合が多くて便利ですが、それだけでは足りなくて自分で計算する必要がある場合もあります。
この記事は法線の計算に関する話をしていきます。
規定のサーフェスの法線
MATLABでサーフェスオブジェクトを作ったら法線が勝手に定義されて、それがVertexNormalsというプロパティに保存されて、私達はそれを利用することができます。
例えばこのようにランダムで凸凹にされた表面を作って法線を矢印として表示してみます。
[mx,my] = meshgrid(-1:0.4:1,-1:0.5:1);
mz = rand(size(mx))*0.5;
hold on
sf = surf(mx,my,mz,FaceColor=[0.9,0.5,0.6],FaceLighting="gouraud");
mn = sf.VertexNormals;
quiver3(mx,my,mz,mn(:,:,1),mn(:,:,2),mn(:,:,3),0,LineWidth=1.5)
light(Position=[-2 1 10])
view(30,45)
axis equal; grid
矢印を見て、凸凹になっている表面でも法線がちゃんと適切に決められているとわかりますね。光の反射もこの法線によって計算されるからリアルっぽく見えます。
でもそのVertexNormalsは自分で設定できるプロパティだから、自動計算を使わずに自分で値を与えることもできます。例えばそもそも滑らかな表面でもランダムで法線を与えたら凸凹な表面に見えます。
[mx,my] = meshgrid(-1:0.2:1,-1:0.2:1);
mz = -mx/2;
mn = rand([size(mx) 3]); % ランダムの法線
sf = surf(mx,my,mz,VertexNormals=mn,FaceColor=[0.95,0.6,0.7],FaceLighting="gouraud");
light(Position=[-3 2 8])
view(30,45)
axis equal
だから法線は光の反射を決める大事な要因だとわかります。
自動的に計算されるのはいいことですが、それがどうやって計算されているのか、気になりますね。だから私はこれについて色々調べてテストして結論ができたのでこの記事を書くこになりました。
三角形の場合から考える
まず一番簡単な場合から理解していきましょう。
一番簡単なモデルは一つの三角形しかない形です。その場合法線はその三角形の表面と垂直するベクトルとなります。
描写するためにこのように画像を描きます。
ar_v = [0 0 0.1
0.7 1.6 0.1
1.6 0 0.1];
ar_f = [1 2 3];
hold on
patch(Faces=ar_f,Vertices=ar_v,FaceColor=[0.8 0.9 0.6]);
quiver3(ar_v(:,1),ar_v(:,2),ar_v(:,3),[0;0;0],[0;0;0],[1;1;1],0,LineWidth=1.5)
view(15,30)
axis equal; grid
この画像で三角形の表面はxy面にあるから法線は[0,0,1]、つまりz軸に沿うものだと簡単にわかるでしょう。
しかし一般的にどうやってその垂直のベクトルを求めればいいでしょう？答えは「クロス積」です。
三角形の頂点が$v_1$、$v_2$、$v_3$とします。ここで$v_1$を基準として、そうしたら$v_{1→2}$と$v_{1→3}$のクロス積は法線と同じ方向となります。
では表面と点とクロス積から計算された法線を描いてみます。
v1 = [0 0 0];
v2 = [1 0.5 0.5];
v3 = [0 0.8 0.5];
vvv = [v1;v2;v3];
cv = cross(v2-v1,v3-v1);
hold on
patch(Faces=[1 2 3],Vertices=vvv,FaceAlpha=0.1);
scatter3(vvv(:,1),vvv(:,2),vvv(:,3),"filled")
quiver3(v1(1),v1(2),v1(3),v2(1)-v1(1),v2(2)-v1(2),v2(3)-v1(3),0,LineWidth=1.5)
quiver3(v1(1),v1(2),v1(3),v3(1)-v1(1),v3(2)-v1(2),v3(3)-v1(3),0,LineWidth=1.5)
quiver3(v1(1),v1(2),v1(3),cv(1),cv(2),cv(3),0,LineWidth=1.5,LineStyle="--")
text(v1(1),v1(2),v1(3),"v_1",VerticalAlignment="top",FontSize=15)
text(v2(1),v2(2),v2(3),"v_2",VerticalAlignment="top",FontSize=15)
text(v3(1),v3(2),v3(3),"v_3",VerticalAlignment="top",FontSize=15)
text(cv(1),cv(2),cv(3)," 法線",HorizontalAlignment="left",FontSize=20)
view(30,20)
axis equal; grid
ただ、ここでできたクロス積はあくまで法線と同じ方向であり、法線は普段単位ベクトルなので、実際に使う時はvecnorm関数とかを使って正則化しておいた方がいいでしょう。
なお$v_2$又は$v_3$から考えても同じ方向のベクトルができるので、$v_1$の場合だけ計算してもいいです。
ここではどの点を使っても同じ法線ができると描写する画像を描いておきましょう。
ar_x = [0 1 0];
ar_y = [0 0.6 1];
ar_z = [0 0.8 0.5];
ar_v = [ar_x;ar_y;ar_z]';
ar_f = [1 2 3];
cv1 = cross(ar_v(2,:)-ar_v(1,:),ar_v(3,:)-ar_v(1,:));
cv1 = cv1/vecnorm(cv1);
cv2 = cross(ar_v(3,:)-ar_v(2,:),ar_v(1,:)-ar_v(2,:));
cv2 = cv2/vecnorm(cv2);
cv3 = cross(ar_v(1,:)-ar_v(3,:),ar_v(2,:)-ar_v(3,:));
cv3 = cv3/vecnorm(cv3);
hold on
patch(Faces=ar_f,Vertices=ar_v,FaceColor=[0.8 0.7 0.9]);
quiver3(ar_x,ar_y,ar_z,[cv1(1) cv2(1) cv3(1)],[cv1(2) cv2(2) cv3(2)],[cv1(3) cv2(3) cv3(3)],0,LineWidth=1.5)
view(60,40)
axis equal; grid
一つの四角形
3Dモデリングでは複数の三角形からなされる表面が基本ですが、四角形も一般に使われています。普段三角形よりも四角形の方がわかりやすいから。MATLABのサーフェスでも四角形からなされています。
ただし厳密にいうと四角形の表面は実際に2つの三角形からなされるものです。だからその四角形の法線はその2つの三角形の足し算が使われます。
まず一つの四角形ないサーフェスを作成して、自動的に計算された法線を見てみましょう。
mx = [0 1
0 1];
my = [0 0
1 1];
mz = [0 0
0 1];
hold on
sf = surf(mx,my,mz,FaceColor=[0.9 0.5 0.8]);
mn = sf.VertexNormals;
quiver3(mx,my,mz,mn(:,:,1),mn(:,:,2),mn(:,:,3),0,LineWidth=1.5)
view(45,15)
axis equal
grid
これを見て、四角形の表面が2つの三角形からなされるとはっきりとわかりますね。それでも4つの頂点の法線が同じ方向に向いているものです。
では次に自分で計算して結果を比べてみます。
mx = [0 1
0 1];
my = [0 0
1 1];
mz = [0 0
0 1];
mxyz = cat(3,mx,my,mz);
% 1つ目の三角形のクロス積を計算する
mc1 = cross(mxyz(1,2,:)-mxyz(1,1,:),mxyz(2,2,:)-mxyz(1,1,:));
% 2つ目の三角形
mc2 = cross(mxyz(2,2,:)-mxyz(1,1,:),mxyz(2,1,:)-mxyz(1,1,:));
% 2つのクロス積の結果を合わせる
mn = mc1+mc2;
% 単位ベクトルへ正則化
mn = mn/vecnorm(mn);
hold on
sf = surf(mx,my,mz,FaceColor=[0.9 0.5 0.8]);
quiver3(mx(1,1),my(1,1),mz(1,1),mn(1),mn(2),mn(3),0,LineWidth=1.5,Color=[0.1 0.6 0.3])
for i = 1:2
for j = 1:2
text(mx(i,j),my(i,j),mz(i,j),sprintf("v_{%d,%d}",i,j),VerticalAlignment="bottom",HorizontalAlignment="right",FontSize=15)
end
end
view(45,15)
axis equal; grid
矢印一つしか描いていないが、上の画像と比べたら同じだとわかります。他の頂点も同じです。
これで一つの四角形の法線の求め方はわかりましたね。
連結の四角形の面の場合
サーフェスオブジェクトは基本的に複数の四角形の面からなされます。その場合各頂点の法線は、その頂点が使われている四角形の面の平均となります。
例えば2つの四角形からなされるサーフェスの場合を見てみます。
mx = [0 1 2
0 1 2];
my = [0 0 0
1 1 1];
mz = [0 0 0
0 1 0];
hold on
sf = surf(mx,my,mz,FaceColor=[0.9 0.8 0.6],FaceLighting="gouraud");
mn = sf.VertexNormals;
quiver3(mx,my,mz,mn(:,:,1),mn(:,:,2),mn(:,:,3),0)
for i = 1:2
for j = 1:3
text(mx(i,j),my(i,j),mz(i,j),sprintf("v_{%d,%d}\n",i,j),VerticalAlignment="bottom",FontSize=15)
end
end
view(15,20)
axis equal; grid
ここで$v_{1,2}$と$v_{2,2}$は2つの面によって使われるので、法線はその2つの面の平均から計算されます。
では実際に$v_{1,2}$だけ自分で計算してみます。（勿論$v_{2,2}$も同じ）
mx = [0 1 2
0 1 2];
my = [0 0 0
1 1 1];
mz = [0 0 0
0 1 0];
mxyz = cat(3,mx,my,mz);
% 左の面の法線
mc11 = cross(mxyz(1,2,:)-mxyz(1,1,:),mxyz(2,2,:)-mxyz(1,1,:));
mc12 = cross(mxyz(2,2,:)-mxyz(1,1,:),mxyz(2,1,:)-mxyz(1,1,:));
mn1 = mc11+mc12;
mn1 = mn1/vecnorm(mn1);
% 右の面の法線
mc21 = cross(mxyz(1,3,:)-mxyz(1,2,:),mxyz(2,3,:)-mxyz(1,2,:));
mc22 = cross(mxyz(2,3,:)-mxyz(1,2,:),mxyz(2,2,:)-mxyz(1,2,:));
mn2 = mc21+mc22;
mn2 = mn2/vecnorm(mn2);
% 両方合わせる
mn = mn1+mn2;
mn = mn/vecnorm(mn);
hold on
sf = surf(mx,my,mz,FaceColor=[0.9 0.8 0.6],FaceLighting="gouraud");
quiver3(mx(1,2),my(1,2),mz(1,2),mn(1),mn(2),mn(3),0,LineWidth=1.5,Color=[0.9 0.3 0.5])
for i = 1:2
for j = 1:3
text(mx(i,j),my(i,j),mz(i,j),sprintf("v_{%d,%d}\n",i,j),VerticalAlignment="bottom",Color=[0.4 0.1 0.6],FontSize=15)
end
end
view(15,20)
axis equal; grid
これで複数の四角形から構成されたサーフェスの各頂点の法線の計算はわかりました。次は全部の頂点にこんな計算を適用するだけです。
自分で計算した法線を適用する
計算の方法がわかったからここで最初の画像と同じようなものを自分で計算した法線で再現してみます。
[mx,my] = meshgrid(-1:0.4:1,-1:0.5:1);
mz = rand(size(mx))*0.5;
mxyz = cat(3,mx,my,mz);
k1 = mxyz(1:end-1,2:end ,:) - mxyz(1:end-1,1:end-1,:);
k2 = mxyz(2:end ,2:end ,:) - mxyz(1:end-1,1:end-1,:);
k3 = mxyz(2:end ,1:end-1,:) - mxyz(1:end-1,1:end-1,:);
mc = cross(k1,k2,3) + cross(k2,k3,3);
mc = mc./vecnorm(mc,2,3);
mn = zeros([size(mx) 3]);
mn(1:end-1,1:end-1,:) = mc;
mn(1:end-1,2:end ,:) = mn(1:end-1,2:end ,:) + mc;
mn(2:end ,1:end-1,:) = mn(2:end ,1:end-1,:) + mc;
mn(2:end ,2:end ,:) = mn(2:end ,2:end ,:) + mc;
mn = mn./vecnorm(mn,2,3);
hold on
sf = surf(mx,my,mz,VertexNormals=mn,FaceColor=[0.9,0.5,0.6],FaceLighting="gouraud");
quiver3(mx,my,mz,mn(:,:,1),mn(:,:,2),mn(:,:,3),0,LineWidth=1.5)
light(Position=[-2 1 10])
view(30,45)
axis equal; grid
表面がランダムなので結果は同じになるわけがないが、同じようにできているとわかるでしょう。
円柱の場合
ただの長方形のサーフェスだったらそのまま自動的に計算してくれた法線を使っても基本的に問題ないのですが、サーフェスから円柱みたいな形を作りたい場合はちょっと問題があります。
それを見せるために円柱らしいサーフェスを作成して普通に法線を表示してみます。
[mtheta,mz] = meshgrid(0:30:360,-1:0.4:1);
mr = (4-mz.^2)/4;
mx = mr.*cosd(mtheta);
my = mr.*sind(mtheta);
hold on
sf = surf(mx,my,mz,FaceColor=[0.9,0.6,0.8],FaceLighting="gouraud");
mn = sf.VertexNormals;
quiver3(mx,my,mz,mn(:,:,1),mn(:,:,2),mn(:,:,3),0,LineWidth=1.5)
light(Position=[2 -1 3])
view(85,15)
axis equal; grid
これを見て正面にあるところの法線がおかしいとわかりますね。そして法線がおかしいせいで反射の様子も違和感がる気がしますね。
それはここが境界線だからです。本来長方形であるはずの紙を巻いて縁を合わせるだけのものです。そのためその縁にある頂点の法線は別々に計算されて、違う方向に向いてしまうのです。
それを解決するためちょっと計算し直す必要があります。
[mtheta,mz] = meshgrid(0:30:360,-1:0.4:1);
mr = (4-mz.^2)/4;
mx = mr.*cosd(mtheta);
my = mr.*sind(mtheta);
hold on
sf = surf(mx,my,mz,FaceColor=[0.9,0.6,0.8],FaceLighting="gouraud");
% この部分だげ追加
mn = sf.VertexNormals;
mn(:,1,:) = mn(:,1,:) + mn(:,end,:);
mn(:,end,:) = mn(:,1,:);
mn = mn./vecnorm(mn,2,3);
sf.VertexNormals = mn;
quiver3(mx,my,mz,mn(:,:,1),mn(:,:,2),mn(:,:,3),0,LineWidth=1.5)
light(Position=[2 -1 3])
view(85,15)
axis equal; grid
これで法線が合わせられて綺麗に反射される表面になります。
又は、自分で全部の法線を計算してパッチオブジェクトに作ることもできます。
[mtheta,mz] = meshgrid(0:30:330,-1:0.4:1);
mr = (4-mz.^2)/4;
mx = mr.*cosd(mtheta);
my = mr.*sind(mtheta);
mxyz = cat(3,mx,my,mz);
mxyzc = circshift(mxyz,-1,2);
k1 = mxyzc(1:end-1,:,:) - mxyz(1:end-1,:,:);
k2 = mxyzc(2:end ,:,:) - mxyz(1:end-1,:,:);
k3 = mxyz(2:end ,:,:) - mxyz(1:end-1,:,:);
mc = cross(k1,k2,3) + cross(k2,k3,3);
mc = mc./vecnorm(mc,2,3);
ar_vn = zeros([size(mx) 3]);
ar_vn(1:end-1,:,:) = mc;
ar_vn(2:end ,:,:) = ar_vn(2:end,:,:) + mc;
ar_vn = ar_vn + circshift(ar_vn,1,2);
ar_vn = ar_vn./vecnorm(ar_vn,2,3);
ar_vn = reshape(permute(ar_vn,[2 1 3]),[],3);
ar_v = reshape(permute(cat(3,mx,my,mz),[2 1 3]),[],3);
sx = size(mz,2);
sy = size(mz,1);
fv1 = (1:sx) + (0:sx:sx*(sy-1)-1)';
fv2 = (1:sx) + (sx:sx:sx*sy-1)';
fv3 = [2:sx 1] + (sx:sx:sx*sy-1)';
fv4 = [2:sx 1] + (0:sx:sx*(sy-1)-1)';
ar_f = cat(3,fv1,fv2,fv3,fv4);
ar_f = reshape(permute(ar_f,[2 1 3]),[],4);
hold on
patch( ...
Faces=ar_f, ...
Vertices=ar_v, ...
VertexNormals=ar_vn, ...
FaceColor=[0.9,0.6,0.8], ...
FaceLighting="gouraud");
quiver3(ar_v(:,1),ar_v(:,2),ar_v(:,3),ar_vn(:,1),ar_vn(:,2),ar_vn(:,3),0,LineWidth=1.5)
light(Position=[2 -1 3])
view(85,15)
axis equal; grid
（結果は上の画像のサーフェスと同じなのでここに載せません）
ドーナツ
長方形を巻いて循環にしたら円柱になりますが、更に円柱を巻いて循環にたらドーナツの形になります。（一般に数学では「トーラス」と呼ばれますが、「ドーナツ」という呼び方が美味しそうだからここでそう呼びます）
円柱では一つの軸が循環になりますが、ドーナツになると両方循環になるので、もしサーフェスから作ったら縁のところに法線のズレが見えて、円柱よりも違和感が大きいでしょう。
だから自分で法線を計算してパッチを作った方がいいと思います。では法線を計算して適用して表示してみます。
r = 3;
sr = 1;
sx = 18;
sy = 9;
[ar_theta,ar_phi] = meshgrid( ...
linspace(0,360-(360/sx),sx), ...
linspace(0,360-(360/sy),sy));
ar_a = sr.*cosd(ar_phi)+r;
ar_x = ar_a.*cosd(ar_theta);
ar_y = ar_a.*sind(ar_theta);
ar_z = sr.*sind(ar_phi);
mxyz = cat(3,ar_x,ar_y,ar_z);
mxyzc1 = circshift(mxyz,-1,2);
mxyzc2 = circshift(mxyzc1,-1,1);
mxyzc3 = circshift(mxyz,-1,1);
k1 = mxyzc1 - mxyz;
k2 = mxyzc2 - mxyz;
k3 = mxyzc3 - mxyz;
mc = cross(k1,k2,3) + cross(k2,k3,3);
mc = mc./vecnorm(mc,2,3);
mcc1 = circshift(mc,1,2);
mcc2 = circshift(mcc1,1,1);
mcc3 = circshift(mc,1,1);
ar_vn = mc + mcc1 + mcc2 + mcc3;
ar_vn = ar_vn./vecnorm(ar_vn,2,3);
ar_vn = reshape(ar_vn,[],3);
ar_v = reshape(mxyz,[],3);
fv1 = (1:sy) + (0:sy:sy*(sx-1))';
fv2 = [2:sy,1] + (0:sy:sy*(sx-1))';
fv3 = [2:sy,1] + [sy:sy:sy*sx-1,0]';
fv4 = (1:sy) + [sy:sy:sy*sx-1,0]';
ar_f = cat(3,fv1.',fv2.',fv3.',fv4.');
ar_f = reshape(ar_f,[],4);
hold on
patch( ...
Faces=ar_f, ...
Vertices=ar_v, ...
VertexNormals=ar_vn, ...
FaceColor=[1 0.8 0.4], ...
FaceLighting="gouraud");
quiver3(ar_v(:,1),ar_v(:,2),ar_v(:,3),ar_vn(:,1),ar_vn(:,2),ar_vn(:,3),0,LineWidth=1.5)
light(Position=[2 1 2])
view(85,15)
axis equal; grid
綺麗な違和感がないドーナツができています。
映り込み
最後に、前回の記事と同じような映り込みのドーナツを作ってみたいです。ただし今回は自分で計算した法線を使います。
パッチオブジェクトには頂点の法線が自動的に計算それず、面の法線しかないので、前回は面ごとに色を塗るという形にしましたが、本来はサーフェスと同じように頂点に色を設定したいです。
やっと自分で頂点の法線を計算することができるようになったので、今回はこの頂点の法線を使って各頂点の映り込みの色を計算します。
今回使う画像はこれにします。（前回と同じシーンで、ただカメラの位置が違うだけ）
|
|
ではコードと結果を載せてこの記事はこれで終了します。
r = 6;
sr = 2.5;
sx = 2880;
sy = 1440;
capo = [40 0 20];
[ar_theta,ar_phi] = meshgrid( ...
linspace(0,360-(360/sx),sx), ...
linspace(0,360-(360/sy),sy));
ar_s = sr.*(1-0.25.*cosd(ar_phi/2*12-ar_theta/2*12).^16);
ar_a = ar_s.*cosd(ar_phi) + r + randn(size(ar_s))*0.0002;
ar_x = ar_a.*cosd(ar_theta);
ar_y = ar_a.*sind(ar_theta);
ar_z = ar_s.*sind(ar_phi);
mxyz = cat(3,ar_x,ar_y,ar_z);
% 法線の計算
mxyzc1 = circshift(mxyz,-1,2);
mxyzc2 = circshift(mxyzc1,-1,1);
mxyzc3 = circshift(mxyz,-1,1);
k1 = mxyzc1 - mxyz;
k2 = mxyzc2 - mxyz;
k3 = mxyzc3 - mxyz;
mc = cross(k1,k2,3) + cross(k2,k3,3);
mc = mc./vecnorm(mc,2,3);
mcc1 = circshift(mc,1,2);
mcc2 = circshift(mcc1,1,1);
mcc3 = circshift(mc,1,1);
ar_vn = mc + mcc1 + mcc2 + mcc3;
ar_vn = ar_vn./vecnorm(ar_vn,2,3);
ar_vn = reshape(ar_vn,[],3);
ar_v = reshape(mxyz,[],3);
% 映り込みの計算
ma = ar_v - capo;
ma = ma./vecnorm(ma,2,2);
mb = ma - 2.*dot(ma,ar_vn,2).*ar_vn;
mb_phi = atan2d(mb(:,2),mb(:,1));
mb_theta = acosd(mb(:,3));
sphimg = double(imread("eliiro.jpg"))/255;
sphimg = [sphimg sphimg(:,1,:)];
[sphphi,sphtheta] = meshgrid( ...
linspace(180,-180,size(sphimg,2)), ...
linspace(0,180,size(sphimg,1)));
ar_c = [
interp2(sphphi,sphtheta,sphimg(:,:,1),mb_phi,mb_theta), ...
interp2(sphphi,sphtheta,sphimg(:,:,2),mb_phi,mb_theta), ...
interp2(sphphi,sphtheta,sphimg(:,:,3),mb_phi,mb_theta)];
fv1 = (1:sy) + (0:sy:sy*(sx-1))';
fv2 = [2:sy,1] + (0:sy:sy*(sx-1))';
fv3 = [2:sy,1] + [sy:sy:sy*sx-1,0]';
fv4 = (1:sy) + [sy:sy:sy*sx-1,0]';
ar_f = cat(3,fv1.',fv2.',fv3.',fv4.');
ar_f = reshape(ar_f,[],4);
patch( ...
Faces=ar_f, ...
Vertices=ar_v, ...
VertexNormals=ar_vn, ...
FaceVertexCData=ar_c, ...
FaceColor="interp", ...
LineStyle="none");
camproj("perspective")
set(gca,"CameraPosition",capo)
set(gca,"CameraTarget",[0 0 0])
set(gca,"CameraViewAngle",15)
lightangle(105,70)
axis equal off
