こんにちは！Iwaken Lab.のなおです。
現在は材料系の研究室でB4として卒論を頑張ってます。春からはRobotics Learning系の研究室に移る予定です。
本記事はIwaken Lab.アドベントカレンダー2023の7日目です。最近アルバイトなどで触っている 3D Gaussian Splatting について書いていきます。特にGitHubで公開されているコードを見ながら実装についても少し見ていきます。
正直CG系の論文を読むのは初めてで、NeRFはもちろんクォータニオンやSH関数、レイトレーシングなどの基礎知識もない状態でした。勉強とコードを読む作業を２ヶ月間繰り返してなんとか少しは分かったように思います。 非情報系だけど、理論的背景を少し見てみたい！ という方に向けた記事です。できるだけ元の論文を読むようにしていますが、間違いがあればご指摘いただけると助かります。
3D Gaussian Splattingって？
3D Gaussian Splattingは一文で表すと 「NeRFに代わる、空間のガウシアン群表現をベースにしたView Synthesis手法」 です。
概要を掴むだけならこちらの記事がとても分かりやすくておすすめです！
↓こちらが3D Gaussian Splattingの公式WebPageになります。
View Synthesis(視点合成)って？
View Synthesisとは、 対象物を(単一または複数の視点から)撮影した画像や動画を学習データとして、学習データにはない自由な視点からみた際の様子を推定するタスク です。 フォトグラメトリと近いものですが、強いて言えばView Synthesis画像データ、フォトグラメトリは3Dデータを出力することが主目的という違いがあります。基本的には静的(動きのない)シーンを扱いますが、最近は動的シーンを学習できるものも出てきています。
まずは具体例を見ましょう。私の机の上の画像100枚程度から3D Gaussian Splattingで学習したものです。
このように限られた視点からの画像データから、自由な視点からの様子を出力できます。
NeRF(Neural Radience Field)って？
NeRFは2020年のECCVで発表されたView Synthesis手法です。3次元シーンを各座標に対する放射輝度と密度を返す連続な5次元の関数で表現しています。この関数を近似するニューラルネットワークを学習することでView Synthesisを実現しています。
↓NeRFの公式WebPage
↓詳しくはAI-SCHOLARさんの記事が分かりやすいです。
2020年の発表を皮切りに、学習速度を高速化したInstant-NeRF(Instant-NGP)や、１枚の画像から学習できるpixelNeRF、NeRFとDreamFieldを組み合わせてtext2shapeを実現した論文など様々な改良、拡張版が発表されています。
3Dガウシアンってなに？
(1次元の)ガウシアン(ガウス分布、正規分布)とは以下の数式で表される分布です。分布と聞いて一番初めにイメージするものかもしれません。
\mathcal{N}(\mu, \sigma^2) = A \exp\{-(x-\mu)^2/2\sigma^2\}
*確率密度関数としては$A=\frac{1}{\sqrt{2\pi\sigma^2}}$ですが、3D Gaussian Splattingでのガウシアンはスケールを含むので$A$で表現しています。
これを3次元にしてみると
\mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}) = A \exp\{-{}^t(\boldsymbol{x}-\boldsymbol{\mu}) \boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})/2\}
ただし、
平均
\boldsymbol{\mu} = (\mu_x, \mu_y, \mu_z)
共分散行列
\boldsymbol{\Sigma} =
\begin{bmatrix}
\sigma_x^2 & \sigma_x\sigma_y & \sigma_x\sigma_z \\ \sigma_x\sigma_y & \sigma_y^2 & \sigma_y\sigma_z \\ \sigma_x\sigma_z & \sigma_y\sigma_z & \sigma_z^2
\end{bmatrix}
可視化は難しいですが、中央から離れるに従って値が小さくなる楕円体みたいな感じですね。
空間を3Dガウシアン群で表現するってどういうこと？
3D Gaussian Splattingでは
- 平均(座標, 3次元)
- 共分散行列(実装上は 回転のクォータニオン表現4次元 + スケール3次元)
に加えて
- 透明度(1次元)
- 色/放射照度(球面調和関数を用いる, 3+45次元?)
を各ガウシアンのパラメータとして定義しています。
「3次元シーンは上記のガウシアンの集合で表現できる」 という仮定のもと、以下の手順でガウシアンのパラメータを学習していきます。
- SfM手法(COLMAP)を使ってカメラ位置と特徴点をマッピングする。
- 各特徴点を座標は引き継ぎ、それ以外はランダムに初期化したガウシアンで置き換える。
- 学習データのカメラ姿勢についてレンダリングを行い、誤差を計算する。
- 誤差逆伝搬によりガウシアンのパラメータを更新する
- 学習の進み具合に応じて分散の大きなガウシアンを2つに分けたり、透明度の高いガウシアンを削除したりする。
- 指定のイテレーション数まで3~5を繰り返す。
※ SfM(Structure from Motion) : 複数の視点からのカメラ画像について特徴点検出・マッチングを用いて視差を求め、三角測量のような方法で各画像がどの座標からどの向きに撮影されたのかを推定する手法
3D Gaussian SplattingのContribution
私の認識では以下の3点です。
- 3次元空間をガウシアン群で表現するモデルを(微分可能な形で)提案した。
- ガウシアンの最適化の過程で、分散の大きなガウシアンを２つに分ける、透明度が高いガウシアンを削除するといった計算量を減らすための処理を導入した。
- ガウシアン群で表された空間を視点に応じて高速にレンダリングする手法を提案した。
実装を見てみる
最後に、出力されるpoint_cloud.plyからガウシアンモデルの実装を見ていきましょう。
学習された点群(iteration_30000/point_cloud.ply)をplyfileライブラリなどを用いて読み込むと、各ガウシアンが62次元のパラメータを持つnumpy.void形式で保存されていることが分かります。
import numpy as np
from plyfile import PlyData, PlyElement
# ファイルの指定
filepath = ${point_cloud.plyのパス}
# PLYファイルを読み込む
plydata = PlyData.read(filepath)
# カスタムプロパティを含むデータを取得
custom_properties = plydata['vertex'].data
print(plydata.__dict__) # 属性の表示
print(len(custom_properties)) # 点群の数を表示
print(len(custom_properties[0])) # １つのガウシアンの次元数を表示
{'_byte_order': '<', 'text': False, '_comments': [], '_obj_info': [], '_elements': (PlyElement('vertex', (PlyProperty('x', 'float'), PlyProperty('y', 'float'), PlyProperty('z', 'float'), PlyProperty('nx', 'float'), PlyProperty('ny', 'float'), PlyProperty('nz', 'float'), PlyProperty('f_dc_0', 'float'), PlyProperty('f_dc_1', 'float'), PlyProperty('f_dc_2', 'float'), PlyProperty('f_rest_0', 'float'), PlyProperty('f_rest_1', 'float'), PlyProperty('f_rest_2', 'float'), PlyProperty('f_rest_3', 'float'), PlyProperty('f_rest_4', 'float'), PlyProperty('f_rest_5', 'float'), PlyProperty('f_rest_6', 'float'), PlyProperty('f_rest_7', 'float'), PlyProperty('f_rest_8', 'float'), PlyProperty('f_rest_9', 'float'), PlyProperty('f_rest_10', 'float'), PlyProperty('f_rest_11', 'float'), PlyProperty('f_rest_12', 'float'), PlyProperty('f_rest_13', 'float'), PlyProperty('f_rest_14', 'float'), PlyProperty('f_rest_15', 'float'), PlyProperty('f_rest_16', 'float'), PlyProperty('f_rest_17', 'float'), PlyProperty('f_rest_18', 'float'), PlyProperty('f_rest_19', 'float'), PlyProperty('f_rest_20', 'float'), PlyProperty('f_rest_21', 'float'), PlyProperty('f_rest_22', 'float'), PlyProperty('f_rest_23', 'float'), PlyProperty('f_rest_24', 'float'), PlyProperty('f_rest_25', 'float'), PlyProperty('f_rest_26', 'float'), PlyProperty('f_rest_27', 'float'), PlyProperty('f_rest_28', 'float'), PlyProperty('f_rest_29', 'float'), PlyProperty('f_rest_30', 'float'), PlyProperty('f_rest_31', 'float'), PlyProperty('f_rest_32', 'float'), PlyProperty('f_rest_33', 'float'), PlyProperty('f_rest_34', 'float'), PlyProperty('f_rest_35', 'float'), PlyProperty('f_rest_36', 'float'), PlyProperty('f_rest_37', 'float'), PlyProperty('f_rest_38', 'float'), PlyProperty('f_rest_39', 'float'), PlyProperty('f_rest_40', 'float'), PlyProperty('f_rest_41', 'float'), PlyProperty('f_rest_42', 'float'), PlyProperty('f_rest_43', 'float'), PlyProperty('f_rest_44', 'float'), PlyProperty('opacity', 'float'), PlyProperty('scale_0', 'float'), PlyProperty('scale_1', 'float'), PlyProperty('scale_2', 'float'), PlyProperty('rot_0', 'float'), PlyProperty('rot_1', 'float'), PlyProperty('rot_2', 'float'), PlyProperty('rot_3', 'float')), count=883400, comments=[]),), '_element_lookup': {'vertex': PlyElement('vertex', (PlyProperty('x', 'float'), PlyProperty('y', 'float'), PlyProperty('z', 'float'), PlyProperty('nx', 'float'), PlyProperty('ny', 'float'), PlyProperty('nz', 'float'), PlyProperty('f_dc_0', 'float'), PlyProperty('f_dc_1', 'float'), PlyProperty('f_dc_2', 'float'), PlyProperty('f_rest_0', 'float'), PlyProperty('f_rest_1', 'float'), PlyProperty('f_rest_2', 'float'), PlyProperty('f_rest_3', 'float'), PlyProperty('f_rest_4', 'float'), PlyProperty('f_rest_5', 'float'), PlyProperty('f_rest_6', 'float'), PlyProperty('f_rest_7', 'float'), PlyProperty('f_rest_8', 'float'), PlyProperty('f_rest_9', 'float'), PlyProperty('f_rest_10', 'float'), PlyProperty('f_rest_11', 'float'), PlyProperty('f_rest_12', 'float'), PlyProperty('f_rest_13', 'float'), PlyProperty('f_rest_14', 'float'), PlyProperty('f_rest_15', 'float'), PlyProperty('f_rest_16', 'float'), PlyProperty('f_rest_17', 'float'), PlyProperty('f_rest_18', 'float'), PlyProperty('f_rest_19', 'float'), PlyProperty('f_rest_20', 'float'), PlyProperty('f_rest_21', 'float'), PlyProperty('f_rest_22', 'float'), PlyProperty('f_rest_23', 'float'), PlyProperty('f_rest_24', 'float'), PlyProperty('f_rest_25', 'float'), PlyProperty('f_rest_26', 'float'), PlyProperty('f_rest_27', 'float'), PlyProperty('f_rest_28', 'float'), PlyProperty('f_rest_29', 'float'), PlyProperty('f_rest_30', 'float'), PlyProperty('f_rest_31', 'float'), PlyProperty('f_rest_32', 'float'), PlyProperty('f_rest_33', 'float'), PlyProperty('f_rest_34', 'float'), PlyProperty('f_rest_35', 'float'), PlyProperty('f_rest_36', 'float'), PlyProperty('f_rest_37', 'float'), PlyProperty('f_rest_38', 'float'), PlyProperty('f_rest_39', 'float'), PlyProperty('f_rest_40', 'float'), PlyProperty('f_rest_41', 'float'), PlyProperty('f_rest_42', 'float'), PlyProperty('f_rest_43', 'float'), PlyProperty('f_rest_44', 'float'), PlyProperty('opacity', 'float'), PlyProperty('scale_0', 'float'), PlyProperty('scale_1', 'float'), PlyProperty('scale_2', 'float'), PlyProperty('rot_0', 'float'), PlyProperty('rot_1', 'float'), PlyProperty('rot_2', 'float'), PlyProperty('rot_3', 'float')), count=883400, comments=[])}}
883400
62
少し属性と次元数を整理すると、 以下のようになります。
- x,y,z: ガウシアンの中心座標(3次元)
- nx,ny,nz: 法線ベクトル(3次元) *3D Gaussian Splattingでは重要ではないので(0,0,0)が入っている
- f_dc_0 ~ 2: 色情報(?)(3次元)
- f_rest_0 ~ 44 : SH関数(?)(45次元)
- opacity: 透明度(1次元)
- scale_0 ~ 2: スケール(3次元)
- rot_0 ~ 3: 回転(4次元)
１つずつ見ていきましょう。
x,y,z
3Dガウシアンの中心座標(平均)を表しています。上で紹介した$(\mu_x,\mu_y,\mu_z)$に相当します。
f_dc_0~2, f_rest_0~44
これに関してはあまり自信がないのですが、3D Gaussian Splattingでは放射輝度を表現するために球面調和関数(Spherical Harmonics Function, SH関数)というモデルを用いています。
放射輝度とは「ある角度から視線が見られた時に、どれくらい明るく見えるか」です。視線の角度$(\theta,\phi)$に対して輝度を返す関数を作れればこれを求めることができます。そこで効率的かつ微分可能な形で放射輝度を表現するため、以下のような関数を用います。これを球面調和関数と呼びます。
$\displaystyle Y_{l, m}(\theta, \phi) = (-1)^{\frac{m+|m|}{2}} \sqrt{\frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!}} , P_l^{|m|}(\cos \theta), e^{im\phi} \tag{19}$
詳細は省きますが、SH関数は量子力学に由来します。電子の取りうる状態を表すシュレディンガー方程式を変数分離して解いて得られるのがSH関数です。なぜこれがレンダリングの世界で用いられているかというと、放射輝度を学習する際、初めは$l=1$から始め、学習が進むにつれて$l=2,3$と考慮する範囲を増やしていくことで、学習を安定化するからです。実際にSH関数の何が3+45次元の値に相当するのかは正直追い切れていないので割愛します。(間違っていたらすみません)
SH関数の導出については以下の記事がとても分かりやすいです。
opacity(透明度)
透明度は1次元で表現されています。ただし値をヒストグラムにしてみると0~1の間に収まっていません。どういう処理が行われているのかはまだよく分かってません。(すいません)
scale_0 ~ 2 と rot_0 ~ 3(スケールと回転)
上述のとおり、ガウシアンの形状的な特徴は共分散行列$\boldsymbol{\Sigma}$で表現されます。しかし$\boldsymbol{\Sigma}$は対称行列である必要があり、かつ１つのパラメータを変えると角度と大きさの両方が変わってしまうため学習が安定しにくいのではという懸念があります。
そこで3D Gaussian Splattingでは共分散行列をスケールと回転に分けて定義しています。
\boldsymbol{\Sigma} = \boldsymbol{R}\boldsymbol{S}\boldsymbol{S}^T\boldsymbol{R}^T
\boldsymbol{S} =
\begin{bmatrix} \sigma_{x'}^2 & 0 & 0\\ 0 & \sigma_{y'}^2 & 0 \\ 0 & 0 & \sigma_{z'}^2\end{bmatrix}
回転行列$\boldsymbol{R}$はクォータニオン(4元数)で表現されています。
\begin{align}
\boldsymbol{R_{xyz}} &=
\left(
\begin{array}{ccc}
\cos\theta_{y}\cos\theta_{z} & -\cos\theta_{y}\sin\theta_{z} & \sin\theta_{y} \\
\sin\theta_{x}\sin\theta_{y}\cos\theta_{z} + \cos\theta_{x}\sin\theta_{z} & -\sin\theta_{x}\sin\theta_{y}\sin\theta_{z} + \cos\theta_{x}\cos\theta_{z} & -\sin\theta_{x}\cos\theta_{y} \\
-\cos\theta_{x}\sin\theta_{y}\cos\theta_{z} + \sin\theta_{x}\sin\theta_{z} & \cos\theta_{x}\sin\theta_{y}\sin\theta_{z} + \sin\theta_{x}\cos\theta_{z} & \cos\theta_{x}\cos\theta_{y}
\end{array}
\right) \\
&= \left(
\begin{array}{ccc}
2q_w^2 + 2q_x^2 - 1 & 2q_xq_y - 2q_zq_w & 2q_xq_z + 2q_yq_w \\
2q_xq_y + 2q_zq_w & 2q_w^2 + 2q_y^2 - 1 & 2q_yq_z - 2q_xq_w \\
2q_xq_z - 2q_yq_w & 2q_yq_z + 2q_xq_w & 2q_w^2 + 2q_z^2 - 1
\end{array}
\right)
\end{align}
つまりscale_0~2が$(\sigma_{x'}^2, \sigma_{y'}^2, \sigma_{z'}^2, )$, rot_0~3が$(q_x,q_y, q_z, q_w)$に対応していると考えられます。
クォータニオンについては以下の記事が分かりやすいです。
最後に
最後まで読んでいただきありがとうございます。
正直まだ理解し切れていないところが多く、稚拙な文章となってしまいました。
ただ、分からない部分を調べていく中で、
- 共分散行列、院試で確率統計やったからわかるぞ！！
- SH関数の図、電子軌道っぽい...?よく調べたら量子論の授業で求めた1次元シュレディンガー方程式の解の3次元バージョンだ！本当に電子軌道では！？
- 回転行列、メカトロニクスの講義でやったぞ！けどクォータニオンって何や、、？
- SfMの特徴点検出、最近読んでる画像認識(MLPプロフェッショナルシリーズ)で出てきたやつだ！
などなど、今まで断片的に触れてきた様々な知識が繋がって理解を助けてくれました。来年からRobotics Learning系に進むので学部でやった内容の半分くらいは無駄になるかも、、と心配していたんですが、今はやって損はなかったかも！という気持ちです。
それでは少し早いですが良いクリスマスをお過ごしください。