この前の記事ではMATLABにおける綺麗に3Dレンダリングする機能について説明しました。
ただし複雑な表現ができないという欠点があるということも挙げました。その中の一つは「鏡みたいな映り込み」のことです。これができたらもっと綺麗な画像が作れるのに、と残念に思ってしまいました。
ないなら自分で作るしかないですね。幸い映り込みの原理はそこまで複雑ではないのでそれらしいものはできなくもないのです。
ということで、この記事では映り込みを作る考え方と簡単な実装について説明します。
鏡面反射
MATLABには鏡面反射を再現する機能が一応持っています。ただし鏡面反射というのは大まかに2つ：
- 鏡面ハイライト
- 周りの物の映り込み
鏡面ハイライトというのは光源の映り込みです。
MATLABで特に何をしなくても表面と照明があれば鏡面ハイライトが描写されます。
例えばこのように少し凸凹ある球体の表面のレンダリングです。
[phi,theta] = meshgrid(deg2rad(-90:0.5:90),deg2rad(180:-0.5:0));
r = 1+randn(size(phi))*0.002;
mx = r.*cos(phi).*sin(theta);
my = r.*sin(phi).*sin(theta);
mz = r.*cos(theta);
surf(mx,my,mz,EdgeColor="none",FaceColor=[0.8 0.8 0.9],FaceLighting="gouraud");
camproj("perspective")
view(90,0)
lightangle(135,45)
axis equal off
このように鏡面ハイライトを表現して綺麗な表面ができますが、もしこれが鏡のようにピカピカな表面だったら周りの景色がこの表面に映るはずですね。
残念ながらMATLABではそんな映り込みを簡単に表現できる関数を持っているわけではないのです。
それでもちょっと工夫したら方法があります。
反射に関する考え方
まずは反射について説明します。
解説するためにとりあえずこのような絵を描きます。
view(-90,90)
lightangle(135,135)
axis equal
xlim([-2 0])
ylim([-1.2 1.2])
grid
hold on
[phi,theta] = meshgrid(deg2rad(90:0.5:270),deg2rad(180:-0.5:0));
r = 1+randn(size(phi))*0.002;
surf(r.*cos(phi).*sin(theta),r.*sin(phi).*sin(theta),r.*cos(theta), ...
EdgeColor="none",FaceColor=[0.8 0.8 0.9],FaceLighting="gouraud");
% 進行ベクトル
qx = zeros(1,10)-2;
qy = zeros(1,10);
qtheta = deg2rad(-45:10:45);
dqx = 2-cos(qtheta);
dqy = sin(qtheta);
atheta = asin(dqy./sqrt(dqx.^2+dqy.^2));
quiver(qx,qy,dqx,dqy,0,LineWidth=1.5)
% 法線ベクトル
qx = -cos(qtheta);
qy = sin(qtheta);
dqx = qx*0.5;
dqy = qy*0.5;
quiver(qx,qy,dqx,dqy,0,LineWidth=2,LineStyle=":",ShowArrowHead="off")
% 反射ベクトル
dqx = -cos(qtheta*2+atheta)*0.5;
dqy = sin(qtheta*2+atheta)*0.5;
quiver(qx,qy,dqx,dqy,0,LineWidth=1.5)
簡単のために2次元で描いています。私達が鏡の表面に映り込んでいる何かを眺める時はこんな感じですね。
ここでは
| 呼び方 | 説明 | |
|---|---|---|
| 青の線 | 進行ベクトル | カメラから鏡面までのベクトル |
| 赤の点線 | 法線ベクトル | 鏡面から直交方向を指すベクトル |
| 橙色の線 | 反射ベクトル | 鏡面から映り込み対象へ向かうベクトル |
実際に光が映り込み対象から鏡面にぶつかって私達の目（カメラ）に入るので、この矢の方向は本来逆ですが、3Dではカメラから考えることが多いのでこのような矢の方向にします。
自然では「反射ベクトルと法線ベクトルとの角度」と「進行ベクトルと法線ベクトルとの角度」は必ず同等なので、「鏡面でどの方向の物が写っているのか」は、簡単な計算で求められます。
ベクトルの方程式に書くとこうなります。
\vec{b} = \vec{a}-(\vec{a}\cdot\vec{n})\vec{n}
ただし$\vec{a}$は進行ベクトルで、$\vec{n}$は法線ベクトルで、$\vec{b}$は反射ベクトルで、$\cdot$はベクトルの内積です。こうやって進行ベクトルと法線ベクトルから反射ベクトルが求められます。
考え方について詳しくはこの記事を参考に。
view(0,45)
lightangle(135,135)
axis equal off
hold on
[phi,theta] = meshgrid(deg2rad(90:0.5:270),deg2rad(180:-0.5:0));
r = 1+randn(size(phi))*0.001;
surf(r.*cos(phi).*sin(theta),r.*sin(phi).*sin(theta),r.*cos(theta), ...
EdgeColor="none",FaceColor=[0.8 0.8 0.9],FaceLighting="gouraud");
quiver(-1.5,1,0.5,-1,0,LineWidth=1.5)
text(-1.25,0.5,0,"$\vec{a}$",Interpreter="latex",FontSize=20,VerticalAlignment="bottom")
quiver(-1,0,-0.8,0,0,LineWidth=2,LineStyle=":",ShowArrowHead="off")
text(-1.4,0,0,"$\vec{n}$",Interpreter="latex",FontSize=20,VerticalAlignment="bottom",HorizontalAlignment="center")
quiver(-1,0,-0.5,-1,0,LineWidth=1.5)
text(-1.25,-0.5,0,"$\vec{b}$",Interpreter="latex",FontSize=20,VerticalAlignment="bottom",HorizontalAlignment="right")
法線ベクトル
反射を表現するのに法線ベクトルと進行ベクトルが必要だとわかりました。では次の問題はその「法線ベクトルがどうやってわかるの？」です。
普段サーフェスやパッチを作る時に必要なのは頂点のx,y,zの位置の他に法線の情報を入れることもできますが、そうしなくても法線は頂点から自動時に決められるから便利です。
サーフェスオブジェクトが作成された後、法線の情報はVertexNormalsというプロパティに入ってあるから、これを使うことができます。
では球体のサーフェスを作って法線を描いてみます。
[phi,theta] = meshgrid(deg2rad(-90:5:90),deg2rad(180:-5:0));
r = 1+randn(size(phi))*0.01;
mx = r.*cos(phi).*sin(theta);
my = r.*sin(phi).*sin(theta);
mz = r.*cos(theta);
hold on
sf = surf(mx,my,mz,FaceColor=[0.8 0.8 0.9],FaceLighting="gouraud");
mn = sf.VertexNormals; % 法線情報を取得
for i=1:size(mx,1) % 一線ずつ法線を描いていく
for j=1:size(my,2)
plot3([mx(i,j) mx(i,j)+mn(i,j,1)*0.1], ...
[my(i,j) my(i,j)+mn(i,j,2)*0.1], ...
[mz(i,j) mz(i,j)+mn(i,j,3)*0.1])
end
end
camproj("perspective")
light(Position=[2 2 2])
set(gca,"CameraTarget",[0 0 0])
set(gca,"CameraPosition",[10 0 0])
set(gca,"CameraViewAngle",11)
axis equal off
今回はわかりやすいようにメッシュの線も描いて、メッシュの間隔も大きくしています。
このように法線がちゃんと球面から出ています。凸凹になっているのもちゃんと法線がわかります。
進行ベクトル
次は進行ベクトルも描いておきます。これは簡単です。カメラを各頂点の位置はら引くだけ。
では進行ベクトルを描いてみます。ただし今回は別のカメラからの視点で描きます。
[phi,theta] = meshgrid(deg2rad(-90:5:90),deg2rad(180:-5:0));
r = 1+randn(size(phi))*0.01;
mx = r.*cos(phi).*sin(theta);
my = r.*sin(phi).*sin(theta);
mz = r.*cos(theta);
capo = [2 0 0];
hold on
sf = surf(mx,my,mz,FaceColor=[0.8 0.8 0.9],FaceLighting="gouraud");
ma = cat(3,mx,my,mz)-reshape(capo,1,1,3); % カメラの位置を引く
% カメラの位置を示す点を描く
scatter3(capo(1),capo(2),capo(3),100,[0.6 0 0.8],MarkerFaceColor="flat")
% 進行ベクトルを描く
for i=1:size(mx,1)
for j=1:size(my,2)
plot3([mx(i,j) mx(i,j)-ma(i,j,1)], ...
[my(i,j) my(i,j)-ma(i,j,2)], ...
[mz(i,j) mz(i,j)-ma(i,j,3)])
end
end
camproj("perspective")
light(Position=[2 2 2])
set(gca,"CameraTarget",[0.5 0 0])
set(gca,"CameraPosition",[10 7 -4])
set(gca,"CameraViewAngle",9)
axis equal off
反射ベクトル
次は法線ベクトルと進行ベクトルから反射ベクトルを計算して描いてみます。
[phi,theta] = meshgrid(deg2rad(-90:5:90),deg2rad(180:-5:0));
r = 1+randn(size(phi))*0.01;
mx = r.*cos(phi).*sin(theta);
my = r.*sin(phi).*sin(theta);
mz = r.*cos(theta);
capo = [10 0 0];
hold on
sf = surf(mx,my,mz,FaceColor=[0.8 0.8 0.9],FaceLighting="gouraud");
mn = sf.VertexNormals; % 法線ベクトル
ma = cat(3,mx,my,mz)-reshape(capo,1,1,3); % 進行ベクトル
ma = ma./vecnorm(ma,2,3); % 進行ベクトルを単位ベクトルにする
mb = ma-2.*dot(ma,mn,3).*mn; % 反射ベクトルの計算
% 反射ベクトルを描く
for i=1:size(mx,1)
for j=1:size(my,2)
plot3([mx(i,j) mx(i,j)+mb(i,j,1)*0.1], ...
[my(i,j) my(i,j)+mb(i,j,2)*0.1], ...
[mz(i,j) mz(i,j)+mb(i,j,3)*0.1])
end
end
camproj("perspective")
light(Position=[2 2 2])
set(gca,"CameraTarget",[0 0 0])
set(gca,"CameraPosition",capo)
set(gca,"CameraViewAngle",11)
axis equal off
これで原理として整いました。次は本番です。
球体パノラマ画像による映り込み
例として映り込みに使う球体パノラマの画像を準備しておいています。
|
|
これをinterp2関数で内挿を行うことで映り込みの色を計算できます。
[phi,theta] = meshgrid(-90:0.1:90,180:-0.1:0);
r = 1+randn(size(phi))*0.00005;
mx = r.*cosd(phi).*sind(theta);
my = r.*sind(phi).*sind(theta);
mz = r.*cosd(theta);
capo = [10 0 0];
figure(Position=[300 100 600 600])
sf = surf(mx,my,mz,EdgeColor="none",FaceColor="interp",FaceLighting="gouraud");
mn = sf.VertexNormals; % 法線ベクトル
ma = cat(3,mx,my,mz) - reshape(capo,1,1,3);
ma = ma./vecnorm(ma,2,3); % 進行ベクトル
mb = ma - 2.*dot(ma,mn,3).*mn; % 反射ベクトル
mb_phi = atan2d(mb(:,:,2),mb(:,:,1)); % 反射φ
mb_theta = acosd(mb(:,:,3)); % 反射θ
sphimg = double(imread("iroeli.jpg"))/255; % 球体パノラマ画像を読み込む
sphimg = [sphimg sphimg(:,1,:)];
[sphphi,sphtheta] = meshgrid( ...
linspace(180,-180,size(sphimg,2)), ...
linspace(0,180,size(sphimg,1)));
% 計算した映り込みで表面の色を設定する
sf.CData = cat(3, ...
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));
material shiny
camproj("perspective")
set(gca,"CameraTarget",[0 0 0])
set(gca,"CameraPosition",capo)
set(gca,"CameraViewAngle",10)
light(Position=[2 2 2])
axis equal off
これで出来上がり！
本当にザラザラの球体金属を見ているような間隔です。
因みにザラザラなく滑らかな表面にしたらこうなります。ちゃんと上手くいって嬉しいです。
アニメーションを作る場合視点が変わったり鏡面が動いたりすると当然毎回計算し直して表面の色を入れ替える必要があります。動いても映り込みがそのままだと不自然になるでしょう。
ではどんどん球体に近づいていくアニメーションを作ってみます。
sphimg = double(imread("iroeli.jpg"))/255;
sphimg = [sphimg sphimg(:,1,:)];
[sphphi,sphtheta] = meshgrid( ...
linspace(180,-180,size(sphimg,2)), ...
linspace(0,180,size(sphimg,1)));
[phi,theta] = meshgrid(0:0.1:180,180:-0.1:0);
r = imresize(1-abs(randn(ceil(size(phi)/50)))*0.01,size(phi),"nearest");
mx = r.*cosd(phi).*sind(theta);
my = r.*sind(phi).*sind(theta);
mz = r.*cosd(theta);
figure(Position=[300 100 600 600])
sf = surf(mx,my,mz,EdgeColor="none",FaceColor="interp",FaceLighting="gouraud");
mn = sf.VertexNormals; % 法線ベクトル
material shiny
camproj("perspective")
set(gca,"CameraTarget",[0 0 0])
light(Position=[2 2 2])
axis equal off
for i = 0:8
ky = -2-((8-i)^2); % 球体の中心からの距離
capo = [0 -ky 0]; % カメラの位置
set(gca,"CameraPosition",capo) % カメラを移す
set(gca,"CameraViewAngle",asind(-1/ky)*2) % 視野の角度もそれに従って変える
ma = cat(3,mx,my,mz) - reshape(capo,1,1,3);
ma = ma./vecnorm(ma,2,3); % 進行ベクトル
mb = ma - 2.*dot(ma,mn,3).*mn; % 反射ベクトル
mb_phi = atan2d(mb(:,:,2),mb(:,:,1));
mb_theta = acosd(mb(:,:,3));
sf.CData = cat(3, ...
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));
pause(0.2)
end
同じ角度から見ても距離が違うと球面の反射で見える景色は少し違いますね。
パッチに使う
最後にパッチにも適用してみます。ただしパッチはサーフェスと違って、作った時に頂点の法線ベクトルがVertexNormalsに入ったわけではないのです。その代わりに各面の法線ベクトルがFaceNormalsに入ります。
法線が頂点ではなく面のものなので、色を設定する時も面ごとの色を使うことになります。そのせいでもしポリゴンを細かく分割しないとギザギザって感じになりやすいです。
ではドーナツの形を作って映り込みさせてみます。
r = 6;
sr = 2;
sx = 2880;
sy = 1440;
[ar_theta,ar_phi] = meshgrid( ...
linspace(0,360-(360/sx),sx), ...
linspace(0,360-(360/sy),sy));
ar_s = sr+0.5.*cosd(ar_phi/2*3-ar_theta/2*12).^2;
ar_a = ar_s.*cosd(ar_phi)+r;
ar_x = ar_a.*cosd(ar_theta);
ar_y = ar_a.*sind(ar_theta);
ar_z = ar_s.*sind(ar_phi);
ar_v = [ar_x(:),ar_y(:),ar_z(:)];
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);
pt = patch(Faces=ar_f,Vertices=ar_v,FaceColor="flat",LineStyle="none",SpecularExponent=30);
camproj("perspective")
capo = [30 0 15];
set(gca,"CameraPosition",capo)
set(gca,"CameraTarget",[0 0 0])
set(gca,"CameraViewAngle",20)
lightangle(135,70)
axis equal off
drawnow
mn = pt.FaceNormals;
ma = [mean(pt.XData);mean(pt.YData);mean(pt.ZData)]' - capo;
ma = ma./vecnorm(ma,2,2);
mb = ma - 2.*dot(ma,mn,2).*mn;
mb_phi = atan2d(mb(:,2),mb(:,1));
mb_theta = acosd(mb(:,3));
sphimg = double(imread("iroeli.jpg"))/255;
sphimg = [sphimg sphimg(:,1,:)];
[sphphi,sphtheta] = meshgrid( ...
linspace(180,-180,size(sphimg,2)), ...
linspace(0,180,size(sphimg,1)));
pt.FaceVertexCData = [
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)];
こうやって綺麗に映り込んでいるドーナツができました。
ただ、勘がいい人はこれを見て何か違和感を感じるでしょう。それはドーナツ自身の姿はドーナツの表面に映っていないということです。
当然ですが、今の手法はあくまで予め準備しておいて球体パノラマ背景を映り込ませることだけで、自分を含め近くにあるものを映り込ませるのは更に難易度が高いのでまだそこまでできていません。
終わりに
こうやって映り込みを再現することができました。ちょっと手間がかかりますが、擬似的にこんなことができるだけで嬉しいことです。まだ物足りなくてもっと機能拡張が必要かもしれませんが、これだけでも可能性が広げています。
反射はできたけど、屈折や影などはもっと難しい課題です。更に現実っぽくレンダリングできるように、今後はレイトレーシングとかもやってみたいと思いますが、あまり簡単なことではないのでまだもっと勉強が必要です。
とりあえず、元々搭載していない機能であってももしアルゴリズムがわかっていたら、ちょっと面倒かもれしないけど自分で書いて実装していくということです。そうすることで理解が深まるというメリットもありますし。
追記
[2024/06/23]
法線の計算に関してはこの記事に纏めておきました。これでパッチの頂点の法線を計算して頂点に映り込みの色を与えることができます。