MATLABで鏡みたいな映り込みのレンダリングをする方法

この前の記事ではMATLABにおける綺麗に3Dレンダリングする機能について説明しました。

https://qiita.com/phyblas/items/150a34293df36b8da313

ただし複雑な表現ができないという欠点があるということも挙げました。その中の一つは「鏡みたいな映り込み」のことです。これができたらもっと綺麗な画像が作れるのに、と残念に思ってしまいました。

ないなら自分で作るしかないですね。幸い映り込みの原理はそこまで複雑ではないのでそれらしいものはできなくもないのです。

ということで、この記事では映り込みを作る考え方と簡単な実装について説明します。


# 鏡面反射

MATLABには**鏡面反射**を再現する機能が一応持っています。ただし**鏡面反射**というのは大まかに2つ：
- 鏡面ハイライト
- 周りの物の映り込み

**鏡面ハイライト**というのは光源の映り込みです。

MATLABで特に何をしなくても表面と照明があれば鏡面ハイライトが描写されます。

例えばこのように少し凸凹ある球体の表面のレンダリングです。
```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
```

![q01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/30b9fe8c-0f7e-6e64-1773-250343ea596c.jpeg)

このように鏡面ハイライトを表現して綺麗な表面ができますが、もしこれが鏡のようにピカピカな表面だったら周りの景色がこの表面に映るはずですね。

残念ながらMATLABではそんな映り込みを簡単に表現できる関数を持っているわけではないのです。

それでもちょっと工夫したら方法があります。


# 反射に関する考え方

まずは反射について説明します。

解説するためにとりあえずこのような絵を描きます。
```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)
```

![q02.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/ac137cd4-a3a8-308c-9f21-6308b43ddc3b.jpeg)

簡単のために2次元で描いています。私達が鏡の表面に映り込んでいる何かを眺める時はこんな感じですね。

ここでは
<table>
<tr><th style="text-align: center; border: 1.5px solid #8866AA;"></th><th style="text-align: center; border: 1.5px solid #8866AA;">呼び方</th><th style="text-align: center; border: 1.5px solid #8866AA;">説明</th></tr>
<tr></tr>
<tr><td style="border: 1.5px solid #8866AA;"><font color="#0a6ab6">青の線</font></td><td style="border: 1.5px solid #8866AA;">進行ベクトル</td><td style="border: 1.5px solid #8866AA;">カメラから鏡面までのベクトル</td></tr>
<tr><td style="border: 1.5px solid #8866AA;"><font color="#ff704b">赤の点線</font></td><td style="border: 1.5px solid #8866AA;">法線ベクトル</td><td style="border: 1.5px solid #8866AA;">鏡面から直交方向を指すベクトル</td></tr>
<tr><td style="border: 1.5px solid #8866AA;"><font color="#dab212">橙色の線</font></td><td style="border: 1.5px solid #8866AA;">反射ベクトル</td><td style="border: 1.5px solid #8866AA;">鏡面から映り込み対象へ向かうベクトル</td></tr>
</table>

実際に光が映り込み対象から鏡面にぶつかって私達の目（カメラ）に入るので、この矢の方向は本来逆ですが、3Dではカメラから考えることが多いのでこのような矢の方向にします。

自然では「**反射ベクトル**と**法線ベクトル**との角度」と「**進行ベクトル**と**法線ベクトル**との角度」は必ず同等なので、「鏡面でどの方向の物が写っているのか」は、簡単な計算で求められます。

ベクトルの方程式に書くとこうなります。

```math
\vec{b} = \vec{a}-(\vec{a}\cdot\vec{n})\vec{n}
```

ただし$\vec{a}$は進行ベクトルで、$\vec{n}$は法線ベクトルで、$\vec{b}$は反射ベクトルで、$\cdot$はベクトルの内積です。こうやって進行ベクトルと法線ベクトルから反射ベクトルが求められます。

考え方について詳しくはこの記事を参考に。

https://qiita.com/edo_m18/items/b145f2f5d2d05f0f29c9

```matlab
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")
```

![q02-.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/207853fd-660a-0c3a-0346-3d44fe8b8751.jpeg)


# 法線ベクトル

反射を表現するのに法線ベクトルと進行ベクトルが必要だとわかりました。では次の問題はその「法線ベクトルがどうやってわかるの？」です。

普段サーフェスやパッチを作る時に必要なのは頂点のx,y,zの位置の他に法線の情報を入れることもできますが、そうしなくても法線は頂点から自動時に決められるから便利です。

サーフェスオブジェクトが作成された後、法線の情報は`VertexNormals`というプロパティに入ってあるから、これを使うことができます。

では球体のサーフェスを作って法線を描いてみます。
```matlab
[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
```

![q03.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/15415272-7843-fc9e-ab1f-95d44c9c4500.jpeg)

今回はわかりやすいようにメッシュの線も描いて、メッシュの間隔も大きくしています。

このように法線がちゃんと球面から出ています。凸凹になっているのもちゃんと法線がわかります。


# 進行ベクトル

次は進行ベクトルも描いておきます。これは簡単です。カメラを各頂点の位置はら引くだけ。

では進行ベクトルを描いてみます。ただし今回は別のカメラからの視点で描きます。
```matlab
[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
```

![q04.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/3f74d68b-ea38-9e57-69cd-cf89df83b37f.jpeg)

# 反射ベクトル

次は法線ベクトルと進行ベクトルから反射ベクトルを計算して描いてみます。
```matlab
[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
```

![q05.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/eefbef5f-57a2-370e-d898-ed4b277148c2.jpeg)

これで原理として整いました。次は本番です。


# 球体パノラマ画像による映り込み

例として映り込みに使う球体パノラマの画像を準備しておいています。
<table><tr><td style="text-align: center">

![iroeli.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/26fe163c-e961-2ac9-db25-d3111db4e81a.jpeg)
<font color="red"><code>iroeli.jpg</code></font>
（[MMDモデルからMaya](https://qiita.com/phyblas/items/e5fe203c955e273b26a4)においてレンダリングしたもの）
</td></tr></table>

これを[`interp2`](https://jp.mathworks.com/help/matlab/ref/interp2.html)関数で内挿を行うことで映り込みの色を計算できます。

```matlab
[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
```

![q06.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/e6352917-a735-f70a-ba6a-2403fb16ed91.jpeg)

これで出来上がり！

本当にザラザラの球体金属を見ているような間隔です。

因みにザラザラなく滑らかな表面にしたらこうなります。ちゃんと上手くいって嬉しいです。

![q06-.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/6a39a996-47a6-4d84-be70-de5d0f31e56d.jpeg)

アニメーションを作る場合視点が変わったり鏡面が動いたりすると当然毎回計算し直して表面の色を入れ替える必要があります。動いても映り込みがそのままだと不自然になるでしょう。

ではどんどん球体に近づいていくアニメーションを作ってみます。
```matlab
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
```

![q07.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/d629e846-0722-5440-6bbe-577920334ee2.gif)

同じ角度から見ても距離が違うと球面の反射で見える景色は少し違いますね。


# パッチに使う

最後にパッチにも適用してみます。ただしパッチはサーフェスと違って、作った時に頂点の法線ベクトルが`VertexNormals`に入ったわけではないのです。その代わりに各面の法線ベクトルが`FaceNormals`に入ります。

法線が**頂点**ではなく**面**のものなので、色を設定する時も面ごとの色を使うことになります。そのせいでもしポリゴンを細かく分割しないとギザギザって感じになりやすいです。

ではドーナツの形を作って映り込みさせてみます。
```matlab
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)];
```

![q08.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/52c9fe5d-74f4-0d50-a88d-3c74c7f2143f.jpeg)

こうやって綺麗に映り込んでいるドーナツができました。

ただ、勘がいい人はこれを見て何か違和感を感じるでしょう。それはドーナツ自身の姿はドーナツの表面に映っていないということです。

当然ですが、今の手法はあくまで予め準備しておいて球体パノラマ背景を映り込ませることだけで、自分を含め近くにあるものを映り込ませるのは更に難易度が高いのでまだそこまでできていません。


# 終わりに

こうやって映り込みを再現することができました。ちょっと手間がかかりますが、擬似的にこんなことができるだけで嬉しいことです。まだ物足りなくてもっと機能拡張が必要かもしれませんが、これだけでも可能性が広げています。

反射はできたけど、屈折や影などはもっと難しい課題です。更に現実っぽくレンダリングできるように、今後はレイトレーシングとかもやってみたいと思いますが、あまり簡単なことではないのでまだもっと勉強が必要です。

とりあえず、元々搭載していない機能であってももしアルゴリズムがわかっていたら、ちょっと面倒かもれしないけど自分で書いて実装していくということです。そうすることで理解が深まるというメリットもありますし。


# 追記

*[2024/06/23]*
法線の計算に関してはこの記事に纏めておきました。これでパッチの頂点の法線を計算して頂点に映り込みの色を与えることができます。

https://qiita.com/phyblas/items/fba7523b95b96e3aeb43