【CV入門】画像から3次元形状を推定しよう！【照度差ステレオ】

こんにちは。ももやまです。

突然ですが、下の2枚の画像をご覧ください。

どちらも丸が描かれていますが、右側の丸は外側が暗くなっていますね。

おそらく、右側の画像は見ただけで、「なんとなく球っぽい形をしているな」というのがわかると思います。

言い換えると、我々人間の目は画像の明暗からどんな形状をしているかを認識できるということですね。

実は、人間の目ではなくとも、コンピューターを用いて物体の形状を認識することができるのです。

そこで今回は、コンピューター上で「画像から立体形状（3次元形状）を求める方法」の1つである照度差ステレオについて見ていきましょう。

スポンサードリンク

１．明るさの関係と反射

(1) 明るさと反射

最初に、人間の目は「明るさの差」から形状を想像すると説明しましたね。

なので、3次元形状を求める前にまずは「明るさ（輝度）の関係」について見ていきましょう。

ある物体（の表面）に光を当てると、反射が起こりますね。

反射が起こることによって我々は物体を見ることができます。

反射には大きく分けて

• 拡散反射
• 鏡面反射

の2つがあります。

今回説明する照度差ステレオでは、2つの反射のうち、拡散反射しか起こらない（と仮定した）ときに適用できる方法です。

照度差ステレオでは、拡散反射の中でも最も単純であるランバートモデルを仮定して物体の3次元形状を推定します。

それでは、ランバートモデルと仮定したときに物体の明るさがどう決まるかを見ていきましょう。

※ 拡散反射、鏡面反射については下の記事で詳しくまとめているので興味がある人はこちらもご覧ください。

www.momoyama-usagi.com

(2) 拡散反射と明るさ

拡散反射による輝度 \( L_{d} \) は、光源の明るさ \( I \), 拡散反射率*1 \( K_{d} \), 物体表面の法線と光源方向のなす角度 \( \theta \) を用いて、以下の式(1)で表されます。\[
L_{d} = I \ K_{d} \cos \theta \tag{1}
\]

※ ただし、\( \theta > 90^{\circ} \) のときは、物体に光が全く当たらないため、\( L_{d} = 0 \) として考えます。

さらに物体表面の法線を \( \vec{n} \), 光源方向を \( \vec{s} \) とし、どちらのベクトルの大きさも1とすると、\( \cos \theta = \vec{s} \cdot \vec{n} \) となるため、式(1)を以下のように書き換えることができます。\[\begin{align}
L_{d} & = I \ K_{d} \cos \theta
\\ & = I \ K_{d} ( \vec{s} \cdot \vec{n} ) \tag{2}
\end{align}\]

スポンサードリンク

２．照度差ステレオの仕組み

(1) 3次元形状と法線

物体の表面は基本的になめらかですよね。

ここで、物体の表面を細かく細かく刻んでいくことを考えましょう。

すると、それぞれの表面は平面と仮定することができます。

この物体の表面を細かく細かく刻んだものを1画素*2だと考えてください。

ある平面の法線 \( \vec{n} \) というのは、平面から垂直方向に出るベクトルを指します。

この垂直方向に出るベクトルのことを、平面が向いている方向だと考えてください。

ここで、先程細かく刻んだ平面（つまり画素）それぞれに対して、法線を求めると、それぞれの画素ごとに平面がどこを向いているのかがわかりますね。

面の傾き具合がわかってしまえば、あとは積分をすることで、物体の3次元形状を求めることができるのです！

（積分する部分は、数2Bや数3で習う積分よりは難しい積分なので、今回は理論部分の説明を省略します。）

(2)では、3次元形状を求めるために必要な法線 \( \vec{n} \) の求め方について説明していきましょう。

(2) 法線の求め方

3次元の世界での法線ベクトル \( \vec{n} \) は、\( x \) 成分、\( y \) 成分、\( z \) 成分の3つからなる3次元ベクトル\[
\vec{n} = \left( \begin{array}{ccc} x \\ y \\ z  \end{array} \right)
\]ですね。

法線を求めるということは、3つの変数 \( x \), \( y \), \( z \) の値をすべて求めればOKということですね…！

例えば、物体表面上のある点を、強さが1の光源方向\[
\vec{s} = \frac{1}{3} \left( \begin{array}{ccc} 2 \\ 2 \\ 1  \end{array} \right)
\]から照らしたときの輝度値（明るさ）が60だったとします。

すると、\( \[\begin{align*}
K_{d} ( \vec{s} \cdot \vec{n} ) & =  \frac{1}{3} K_{d} (2x+2y+z)
\\ & = 60 
\end{align*}\]という式が成り立ちますね。

この式には \( K_{d} \), \( x \), \( y \), \( z \) の4つの未知数がありますね。

ここで、拡散反射率 \( K_{d} \) は、同じ点であればどの光源方向から照らそうが変わりませんよね。

また、法線ベクトルは向きさえあっていれば大きさはどうでもいいので、\[
K_{d} ( \vec{s} \cdot \vec{n} ) = \vec{s} \cdot ( K_{d} \ \vec{n} )
\]とし、法線ベクトル \( \vec{ \tilde{n} } \) を\[
 \vec{ \tilde{n} } = \frac{ K_{d} \ \vec{n} }{ | K_{d} \ \vec{n} | }
\]のように大きさ1に正規化したものとします。

すると、 \( K_{d} \) の値はどうでもよくなるので、未知数は \( x \), \( y \), \( z \) の3つになりますね。

なので、この連立方程式を解くためには3つの式が必要になってきますね。

なので、最低でも3枚の画像、言い換えると3つの異なる方向から撮影した画像があればOKですね。

それぞれの光源方向で強度が違う場合は…？

上の例では \( I = 1 \) と仮定しましたが、常にどの光源も同じ強度になるとは限りませんよね。

なので、各光源方向で強度が違う場合を考えて見ましょう。

光源方向を、光源の強度×光源方向（つまり \( I \vec{s} \) として考えます。

さらに、光源強度×光源方向を表すベクトルを \[
\vec{l} = \left( \begin{array}{ccc} l_x \\ l_y \\ l_z  \end{array} \right)
\] とすると、拡散反射の式(1)を\[\begin{align*}
L_{d} & = I K_{d} ( \vec{s} \cdot \vec{n} ) \\ & = (I \vec{s}) \cdot ( K_{d} \ \vec{n})
\\ & = \vec{l} \cdot \vec{ \tilde{n} } \tag{3}
\\ & = l_x \cdot x + l_y \cdot y + l_z \cdot z  
\end{align*}\]と変形することができます*3。

最後に、\[
 \vec{ n } = \frac{ \vec{ \tilde{n} } }{ | \vec{ \tilde{n} } | }
\]と法線ベクトル \( \vec{ \tilde{n} } \) の大きさを1に正規化することで、強度が異なる場合でも法線を推定することができますね。

(3) 法線を手計算してみよう

照度差ステレオの仕組みを理解したかを確認するために、1問例題を解いてみましょう。

例題

物体表面上のある点Xを3つの光源\[
\vec{l}_1 = \frac{1}{3} \left( \begin{array}{ccc} 2 \\ 2 \\ 1  \end{array} \right) , \ \ \ 
\vec{l}_2 = \left( \begin{array}{ccc} 1 \\ 0 \\ 1  \end{array} \right) , \ \ \ 
\vec{l}_3 = \left( \begin{array}{ccc} 0 \\ 0 \\ 1  \end{array} \right) 
\]で照らしたとき、輝度値は60, 90, 40となった。

点Xの大きさ1の法線ベクトルを求めなさい。

解説

\( \vec{l}_1 \) より、\[
\frac{1}{3} (2x+2y+z) = 60 \]\[
2x+2y+z = 180 \tag{1-1}
\]の式が得られる。また、\( \vec{l}_2 \) より、\[
x+z=90
\]の式が得られる。さらに \( \vec{l}_3 \) より、\[
z=40
\]の式が得られる。式(1-1)～(1-3)より、\[
\left\{\begin{array}{l} 2x+2y+z= 180 \\ x+z=90 \\ z=40  \end{array}\right.
\qquad \therefore \quad \left\{\begin{array}{l} \displaystyle x = 50  \\ \displaystyle y = 20 \\ \displaystyle z = 40   \end{array}\right.
\]

よって、\( \vec{ \tilde{n} } \) は\[
 \vec{ \tilde{n} } = \left( \begin{array}{ccc} 50 \\ 20 \\ 40  \end{array} \right)
\]となるので、大きさを1に正規化した法線ベクトルは\[
 \vec{ n } = \frac{1}{ 3 \sqrt{5} } \left( \begin{array}{ccc} 5 \\ 2 \\ 4  \end{array} \right)
\]となる。

スポンサードリンク

３．連立方程式と行列

理系大学だと、連立方程式を行列を用いて解く方法を習うと思います。

忘れてしまった人は下の記事にて復習をしましょう。

www.momoyama-usagi.com

先程の連立方程式で表された式(3)を3つ集めた式\[ 
\vec{l}_1 \cdot \vec{ \tilde{n} } = i_{1} \]\[
\vec{l}_2 \cdot \vec{ \tilde{n} } = i_{2} \]\[
\vec{l}_3 \cdot \vec{ \tilde{n} } = i_{3} 
\]を行列を用いて表してみましょう。

（※ \( i_1 \), \( i_2 \), \(i_3 \) はそれぞれある点Xを1, 2, 3番目の光源で照らしたときの輝度値を表す。）

光源行列 \( L \)、および輝度値を並べたベクトル \( \vec{i} \) を\[
L = \left( \begin{array}{ccc} {}^t\! \ \vec{l}_1 \\ {}^t\! \ \vec{l}_2 \\ {}^t\! \ \vec{l}_3  \end{array} \right) , \ \ \ \vec{i} = \left( \begin{array}{ccc} i_1 \\ i_2 \\ i_3  \end{array} \right)
\]で定義します。

（※ \( {}^t\! \ \vec{l} \) は、ベクトル \( \vec{l} \) の転置を表します。）

すると、連立方程式を\[
L \vec{ \tilde{n} }  = \vec{i}
\]と1つの式で表すことができます。 

線形代数の力を借りると、下のように逆行列を使うことで\[
\vec{ \tilde{n} } = L^{-1} \vec{i}
\]とすることで、法線ベクトル \( \vec{ \tilde{n} } \) を求めることができます。

（※ 最後に正規化 \( \vec{n} = \vec{ \tilde{n} } / |\vec{ \tilde{n} }| \) を忘れずに）

４．画像が4枚以上の場合は…？

実際に照度差ステレオのプログラムを組んでみるとわかるのですが、画像が3枚の場合だとあまりいい精度が出ません。

そこで、もっと精度を上げるために4枚（4光源で撮影）以上の照度差ステレオについて考えてみましょう。

4光源以上あるということは、4つ以上の式があるということですね。

未知数が3つしかないのに4つ以上の式があると、基本的にはすべての式を満たすような解を求めることができません。例えば\[
\left\{ \begin{array}{l} 2x+2y+z= 180 \\ x+z=90 \\ z=40 \\ 4x + 3y + 5z = 550 \end{array} \right.
\]を満たすような \( x \), \( y \), \( z \) は存在しないことがわかりますね。

すべての式を満たすような解を求めることができないならダメじゃんとおもう人もいるかもしれません。しかし、完全に満たすような解は求められなくても、答えに近い解を求めることはできますね。

この答えに最も近い解を求める方法のことを最小2乗法といいます。

(i) 最小2乗法とは？

最小2乗法の中で一番有名なのは、2つ以上の \( n \) 個の点 \( (x_1,y_1) \), \( (x_2,y_2) \), …, \( (x_n,y_n) \) が与えられたときにもっと近い直線\[
y = ax + b
\]を求める手法ですね。

この手法の場合、\( k\) 番目による誤差 \( r_k \) を\[
r_k = y_k - a x_k - b
\]を2乗して全部足したもの、つまり\[
\sum^{n}_{k = 1} r_k^2
\]が最小となるような \( a \), \( b \) を求めることで最も近い直線を求めることができます。

最小2乗法については、こちらの記事にて説明をしているので、もし忘れてしまった人は確認をおすすめします。

www.momoyama-usagi.com

(ii) 照度差ステレオで使う最小2乗法

(i)の最小2乗法は、誤差の2乗を最も小さくするように \( a \), \( b \) という2つの未知数を求めるという方法でしたね。

照度差ステレオで使う最小2乗法では、法線の \( x \), \( y \), \( z \) 成分、つまり3つの未知数を求める手法を使います。

3つの未知数を求める場合でも、(i)と同じ考え方で誤差の2乗を小さくします。

例えば、下のように \( n \) 個の式が与えられていると考えましょう。\[\begin{align*}
\left\{ \begin{array}{l} 
\alpha_1 x + \beta_1 y + \gamma_1 z & = i_1 \\
\alpha_2 x + \beta_2 y + \gamma_2 z & = i_2 \\
\alpha_3 x + \beta_3 y + \gamma_3 z & = i_3 \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots \\ 
\alpha_n x + \beta_n y + \gamma_n z & = i_n \\
\end{array} \right. \end{align*} \]

まず、\( k \) 番目の式における誤差 \( r_k \) を下のように定義しましょう。\[
r_k = i_k - \alpha_k x - \beta_k y - \gamma_k z 
\]

つぎに \( k\) 番目による2乗誤差 \( r_k^2 \) を全部足したもの、つまり\[
\sum^{n}_{k = 1} r_k^2
\]が最小となるような \( x \), \( y \), \( z \) を求めることで、答えに最も近くなるような解を求めることができます。

(iii) 最小2乗法を行列で表してみよう

未知数の数に限らず、最小2乗法をコンピュータ上で解く際には、下のように行列（擬似逆行列・一般化逆行列）を用いた方法がよく使われます。

うさぎでもわかる線形代数　応用編第9羽　行列を使って最小2乗法を解いてみよう！

2022.5.25

そもそも解をもたないような連立方程式 Ax = b に対して、一番それっぽい解(最小2乗解) x を計算する行列の最小2乗法について仕組みや計算例などを詳しくまとめています。...

照度差ステレオでは、\( n \) 個の式（つまり \( n \) 個の光源・画像）の最適解（つまり法線）を行列を用いた最小2乗法で計算していきます。

まず、光源行列 \( L \)、および輝度値を並べたベクトル \( \vec{i} \) を\[
L = \left( \begin{array}{ccc} {}^t\! \ \vec{l}_1 \\ {}^t\! \ \vec{l}_2 \\ \vdots \\ {}^t\! \ \vec{l}_n  \end{array} \right) = \left( \begin{array}{ccc} \alpha_1 & \beta_1 & \gamma_1 \\ \alpha_2 & \beta_2 & \gamma_2 \\ \vdots & \vdots & \vdots  \\ \alpha_n & \beta_n & \gamma_n  \end{array} \right)
\]\[ \vec{i} = \left( \begin{array}{ccc} i_1 \\ i_2 \\ \vdots \\ i_n  \end{array} \right)
\]で定義します。

（※1 \( \vec{l}_k \), \( i_k \) はそれぞれ \( k \) 番目の光源方向、画素値を表す）
（※2 \( \alpha_k \), \( \beta_k \), \( \gamma_k \) は \( k \) 番目の光源方向の \( x \), \( y \), \( z \) 成分を表す。(ii)の変数と同じ。）

すると、3つの式のときと同じように法線 \( \vec{ \tilde{n} } \) を求めるための連立方程式を\[
L \vec{ \tilde{n} }  = \vec{i}
\]と1つの式で表すことができます。

ここで、2乗誤差 \( r_k^2 \) を全部足したものは行列を用いて\[
| L \vec{ \tilde{n} } - \vec{i} |
\]と書き換えることができます。

上の式を最小にするような \( \vec{ \tilde{n} } \) は、疑似逆行列\[
L^{+} = ( {}^t\! L L )^{-1} \  {}^t\! L
\]を用いて、\[
\vec{ \tilde{n} } = L^{+} \vec{i}
\]と解くことができます。

※1 最後に \( \vec{n} = \vec{ \tilde{n} } / | \vec{ \tilde{n} }| \) を忘れずに。
※2 擬似逆行列 \( L^{+} \) を求める過程を知りたい方は、こちらをご覧ください。

５．法線推定結果の評価方法

法線推定を行っても、結果が正しいかどうかを確認できなければ意味がありません。

そこで、結果が正しいかを確認する2つの方法について説明しましょう。

(1) 画像出力による評価（定性的評価）

それぞれの画素の法線が具体的に何の値であったかを数値で確認するのはキリがありません。

そこで、画素ごとの法線がどの向き（値）だったかを色で表現します。この色をつけた画像のことを法線マップと呼びます。

具体的には、法線の \( x \), \( y \), \( z \) 成分を下のようにカラー成分のR, G, Bに対応させます。

• \( x \) 成分の -1～1 → カラー画像の赤(R)成分の 0～255
• \( y \) 成分の -1～1 → カラー画像の緑(G)成分の 0～255
• \( z \) 成分の -1～1 → カラー画像の青(B)成分の 0～255

数式で表すと、\[
\left( \begin{array}{ccc} R \\ G \\ B  \end{array} \right) = \frac{255}{2} \left( \begin{array}{ccc} x+1 \\ y+1 \\ z+1  \end{array} \right)
\]の変換となります。

(2) 数値による評価（定量的評価）

(1)の画像による評価は、結果が正しいかどうかを直感的に確かめることができますが、推定の精度として使うには大雑把です。

そこで使われるのが誤差の度合いを数値により評価する方法です。定量的評価とも呼ばれます。

法線の推定誤差を評価する際には、平均角度誤差というのがよく使われます。

平均角度誤差というのは、法線の正しい値 \( \vec{n}_t \) と推定結果 \( \vec{n}_e \) の2つのベクトルが何度離れているかを画素ごとに求め、それらを平均したものを表します。

ここで、ある画素 \( p \) の正しい値 \( \vec{n}_{t_p} \) と推定結果 \( \vec{n}_{e_p} \) の角度誤差を \( \theta_p \) とします。

法線の大きさは必ず1なので、式(4)が成立します。\[
\cos \theta_p = \vec{n}_{t_p} \cdot  \vec{n}_{e_p} \tag{4}
\]

なので、角度誤差 \( \theta_p \) を\[
\theta_p = \cos^{-1} (  \vec{n}_{t_p} \cdot  \vec{n}_{e_p} ) \tag{5}
\]と求めることができます。

各画素の角度誤差さえ求まってしまえば、あとは画素ごとの平均を取るだけでOKです。

式で書くと、下のようになります。\[\begin{align*}
\theta & = \frac{1}{n} \sum^{n}_{p = 1} \theta_p
\\ & = \frac{1}{n} \sum^{n}_{p = 1} \cos^{-1}  (  \vec{n}_{t_p} \cdot  \vec{n}_{e_p} ) 
\end{align*}\]

（ \( \theta \) が平均角度誤差、\( n \) は画素数を表す）

６．かげに要注意！

実際に照度差ステレオのアルゴリズムを組んでみるとわかるのですが、何も工夫しないと何枚画像を使っても推定結果があまりよくなりません。

角度誤差はもちろんのこと、法線マップも真値とはちょっと違う色合いになっていますね。

（10枚画像を使った場合でも平均角度誤差が6°～12°程度になってしまいます。）

輝度を求める式を思い出してみましょう。

拡散反射による輝度値は、物体に光が当たらない（物体表面の法線と光源方向のなす角度 \( \theta \) が90°を超えるとき）ときは 0 になるのでしたね。

しかし、実際は \( \theta > 90^{\circ} \) のとき、\( \cos \theta < 0 \) となるため、\( I \ K_{d} \cos \theta  \) は必ず負になります。

つまり、輝度値が0となる画素からは、正しい式を得ることができません。

どんなに画像の数（式の数）を増やしても、その中に正しくない式が混ざっていると推定結果はよくなりませんよね。

そこで、輝度値が0となる画素の式を除き、残りの式を使って最小2乗法を使って法線推定を行うようにアルゴリズムを変えてみましょう。

すると、どうでしょう。

誤差は歴然と減りましたね！

７．照度差ステレオを体験してみよう！

github内に照度差ステレオを体験するためのプログラムなどを用意しました。

もし体験したい人や、照度差ステレオのプログラムコードを見てみたい人は下のURLをご覧ください。

github.com

※ Python, MATLABのどちらかがつかえる環境であれば体験することができます。

８．さいごに

今回は、コンピューター上で「画像から立体形状（3次元形状）を求める方法」の1つである照度差ステレオについて説明しました。

実は画像から3次元形状を求めることを業界用語（専門用語）で「形状復元」と呼びます。

この形状復元はコンピュータビジョン(CV)の中でも重要な研究テーマとなっており、2021年になろうとしている現在でも、研究が行われています。

では、2020年良いお年を。

９．参考文献

R. Woodham, “Photometric method for determining surface orientation from multiple images”, Optical Engineering, Vol.19, No.1, pp.139–144, 1980.