はじめに

最近、3D業界で大きな衝撃を与えた「3D Gaussian Splatting」¹について、ご存知でしょうか？数少ない写真から、目を奪われるほど美しい三次元シーンを再構成できるデモを見て私も大感動しました。なぜこんなに美しいのか、どんな技術で実現したのか、興味が湧いています！

"普通の3D物体ではなく、カメラの移動に合わせて、水面に映る景色も正確に表現しています。これはなかなか凄い..."

私も時間をかけて論文や公開されたコード²を勉強しました。本家の実装はCUDA化されており、難解な部分が多く、論文に書かれていないこともあります。そのため、「3D Gaussian Splatting」を勉強したい人にむけ、わかりやすい解説記事を書こうと思いました。単に概念や考え方だけでなく、ゼロから再実装できるように、すべてのロジックを数式として整理し、徹底的に解説しようと思います。

「3D Gaussian Splatting」の概説

「3D Gaussian Splatting」を簡単に説明すると：

3D Gaussian Splattingは、その名の通り、複数の元写真から3D Gaussianを生成し、これらの3D Gaussianをカメラの姿勢に合わせて2D画面にsplatすることで、任意の視点でも2D画像を生成できる技術です。

わざと、変な単語が使われていますので、普通に読んでも全然わからないと思いますが、おそらく以下の３つの質問が出てくると思います。

１. 3D Gaussianとは？（3D物体の表現方法）
２. splatとは？（2D画像の生成方法）
３. 元の写真からどうやって、3D Gaussianを作るの？（学習方法）

これらの３つの質問に答えたら、「3D Gaussian Splatting」の本質が完全に理解できるかと思いますので、本記事では、順番に質問を答えようと思います。

3D Gaussianとは？（3D物体の表現方法）

これまで、3D物体を表現する方法はたくさんありますが、ゲームやCGなどでは、主に三角形（メッシュ）を使っています。三角形の数が多いほど、より詳細な3D物体が表現できます。一方、ロボット業界では、点群や3Dボクセルを使って3D物体を表現する方法が一般的です。例えば、SLAMで使われるNDTという技術は、3D空間を3Dボクセルで分割し、ボクセル内の物体を正規分布（Normal Distribution）で表現します。

3D Gaussian Splattingは、三角形や点群、ボクセルではなく、3D Gaussianを使用します。この方法は、NDTと似ていますが、NDTのように固定なボクセルで空間を分割せず、かなり高密度なガウス分布に従う単位を使うことで、細かい物体の形状を表現できます。

以下のデモでは、4つの3D Gaussianが可視化されています。それぞれの3D Gaussianの共分散が異なるため、サイズと形が異なります。ガウス分布に従って、中心から離れるほど透明になります。私たちの世界では、多くのものがガウス分布に従っているため、3D Gaussianを基本単元として世界を表現すると、三角メッシュよりも表現力が高くなるでしょう。

3D Gaussianは、物体の形だけでなく、視点方向に依存する色情報も持っています。上のデモでは、カメラ位置を移動させると、3D Gaussianの色も少しずつ変わります。この能力を利用することで、複雑な光の反射特性を近似できるため、実物に近い色合いを表現することができます。

なお、この文章で3D Gaussianを使用する際には、単なる3次元のガウス分布を指すのではなく、中心点の位置、共分散、透明度、視点方向依存の色情報を持つ3次元表現の基本単位を指します。

また、「ガウス分布」という言葉を使う時には、中心点と共分散から定義される確率分布のみを指します。

Splattingとは？（2D画像の生成方法）

次は、これらの3D Gaussianがどのように2D画像を生成するのかについて説明します。"Splat"という単語が理解するヒントになると思います。

chatGPTさんより：

"Splat"は何かが強く衝突して、平らに広がる様子を表現する言葉です。

3D GaussianをSplatすると、3D Gaussianが2D画面にぶつかって（投影されて）広がる様子を想像できます。3D Gaussianは形状を持っているため、2D画面に投影されると、2Dガウス分布として平らに広がります。以下の図では、異なる2つのカメラ視点から見た3D Gaussianの投影の様子が示されています。

ちなみに、任天堂の「Splatoon」もsplat + platoonの造語であり、インクを壁にぶつかるというイメージです。もしかしたら、芸術的な才能を持つイカさんたちは、3D Gaussianの特性を持つインクを使って、驚くほど美しい絵を描けるかもしれませんね。

2D画像の生成方法の定式化

あやふやなイメージだけでなく、数学的な原理を正しく理解するために、本記事ではすべてのロジックを数式として整理します。
まず、2D画像を生成するために必要な入力データを、パラメータと設定値の2種類に分けて定義します。

パラメータ:画像から学習できるデータ
設定値:既知のデータで、学習する必要がないデータ

パラメータ

3D Gaussian Splattingのパラメータは、物体もしくは空間の性質を記述する3D Gaussianに関する情報が含まれます。少ない元画像から学習することができます。

${p_w}_i$: 3D Gaussian分布の中心点（世界座標における）
$q_i$: 3D Gaussianの回転（クォータニオン）
$s_i$: 3D Gaussianのスケール
$h_i$: 3D Gaussianの球面調和パラメータ（視点方向依存の色情報）
$\alpha_i$: 3D Gaussianの不透明度

iは3D Gaussianのインデックスです。以下では、簡略化のためiを省略します。

設定値

設定値は、仮想カメラの設定に関する値です。これらの値を変化させると、投影された画像が変わりますが、3D世界の本質には関係せず、学習の対象ではありません。以下に、カメラの設定値を定義します。

$R_{cw}$: カメラの回転（ローカル座標系）
$t_{cw}$: カメラの並進（ローカル座標系）
$K$: カメラの内部パラメータ (i.e. $f_x, f_y, c_x, c_y$)

全体のパイプライン

3D Gaussianから2D画像を生成する全体的なパイプラインは以下です。

1.　3D Gaussianの中心点を2D画面に投影するときの座標値を計算する。
2.　3D Gaussianの3D共分散を計算する。
3.　3D Gaussianの3D共分散を2D画面に投影するときの2D共分散を計算する。
4.　カメラの姿勢に応じて、見えた3D Gaussianの色を計算する。
5.　1~4に計算された情報を用いて、2Dの画像を生成する。

以下、これらの処理を順番に定式化します。

※注意: より理解しやすくするため、本記事の数式は、以下のルールを用いて表示します。そのため、元の論文とは異なる記号を使用することがあります。

大文字：行列
子文字：ベクトル or スカラー量
ローマン体：関数名

(1).2D画面に投影するときの座標値

3D Gaussianの中心点($p_w$)をカメラの回転($R_{cw}$)と並進($t_{cw}$)を用いて、カメラ座標系の点($p_c$)に変換します。この処理は式(1)で示します。

$$
\begin{aligned}
p_c & = \mathrm{T_{cw}}(p_w) \\
& = R_{cw} p_w + t_{cw}
\end{aligned}
\tag{1.1}
$$

そして、ピンホールモデルを用いて、カメラ座標系の点($p_c$)を画面に落とすときのピクセル座標値を計算します。

$$
\begin{aligned}
\mathrm{u}(p_c) &=
\begin{bmatrix} x f_x /z + c_x \\ y f_y /z + c_y \end{bmatrix}
\end{aligned}
\tag{1.2}
$$

ここで、$p_c$の各成分をそれぞれx,y,zと表します。焦点距離($f_x$、$f_y$)と画像センサーの中心位置($c_x$、$c_y$)は既知のカメラの内部パラメータです。この部分の詳細については、以前の記事で説明したことがありますので、詳しくは「バンドル調整」を参考にしてください。

(2).3D Gaussianの3D共分散を計算

3D Gaussianの3D共分散は3x3の半正定値行列ですが、これを直接に行列として表現すると、勾配降下法などの最適化アルゴリズムで更新する際に、生成された行列が半正定値行列であることを保証できません。そのため、3Dガウス分布を式(2)のように、結果の正定値性を保証できる回転$q$と各軸のスケール$s$の合成で表現します。

$$
\begin{aligned}
\mathrm{\Sigma}(q, s) &= RSS^TR^T \\
\mathrm{\sigma}(q, s) & = \mathrm{upper\_triangular}(\Sigma) \\
\end{aligned}
\tag{2}
$$

ここで、$R$はクォータニオン$q$の回転行列表現です。$S$はベクトル$s$から作られた対角行列です。

$$
\begin{aligned}
R &=
\begin{bmatrix}
1 - 2(q_y^2 + q_z^2) & 2(q_xq_y - q_zq_w) & 2(q_xq_z + q_yq_w) \\
2(q_xq_y + q_zq_w) & 1 - 2(q_x^2 + q_z^2) & 2(q_yq_z - q_xq_w) \\
2(q_xq_z - q_yq_w) & 2(q_yq_z + q_xq_w) & 1 - 2(q_x^2 + q_y^2)
\end{bmatrix}
\end{aligned}
\tag{2.1}
$$

$$
\begin{aligned}
S &=
\begin{bmatrix}
s_0 & 0 & 0 \\
0 & s_1 &0 \\
0 & 0 & s_2
\end{bmatrix}
\end{aligned}
\tag{2.2}
$$

式(2)の計算は、共分散行列の固有値分解としても解釈できます。分解から得られた固有ベクトルは、データの最大分散の方向(主方向)を示します。固有値はその方向における分散の大きさを表します。言い換えれば、固有ベクトル行列はデータの回転を表し、固有値はスケーリングを示します。

3D共分散行列 $\Sigma$ は対称行列なので、そのまま3x3の行列を保存するのではなく、上三角行列のデータを6次元のベクトル $\sigma$ として保存します。本文では、$\mathrm{upper\_triangular}$ は対称行列から上三角部分を取り出してベクトルに変換する操作を示し、$\mathrm{symmetric\_matrix}$ はその逆操作であり、ベクトルから対称行列を再構築する操作を示します。

$$
\Sigma = \mathrm{symmetric\_matrix}(\sigma)
$$

(3). 3D Gaussianを2次元画像に投影したときの2次元共分散行列

2次元の画像に投影されると、3次元のガウス分布は2次元ガウス分布になります。この2次元ガウス分布の共分散行列は、式(3)を用いて計算できます。

$$
\begin{aligned}
\Sigma^{\prime}(\sigma, p_c) &= J R_{cw}\Sigma R_{cw}^T J^T \\
\sigma^{\prime}(\sigma, p_c) &= \mathrm{upper\_triangular}(\Sigma^{\prime})
\end{aligned}
\tag{3}
$$

元の論文にはこの式の導出が説明されていませんが、少し変わった式ですので、ここで式の導出を試みます。

式(3)の証明

$\mathbb{p}_w$を世界座標3Dの点の集合、$m_w$を$\mathbb{p}_w$の平均とすると、共分散行列の定義により、$\mathbb{p}_w$の共分散行列は以下のように計算されます。

$$
\Sigma = \mathrm{E}[(\mathbb{p}_w-m_w)(\mathbb{p}_w-m_w)^T]
\tag{3.1}
$$

※$\mathrm{E}$は平均値を取る計算の記号

式(1.1)を用いて、これらの点をカメラ座標に変換すると、カメラ座標の共分散行列$\Sigma_c$を計算できます。

$$
\begin{aligned}
\Sigma_c
&= \mathrm{E}[(\mathbb{p} _c-m_c)(\mathbb{p} _c-m_c)^T] \\
&= \mathrm{E}[(\mathrm{T}\mathrm{_cw}(\mathbb{p}_w)- \mathrm{T}\mathrm{_cw} (m_c))(\mathrm{T}\mathrm{_cw}(\mathbb{p}_w)-\mathrm{T}\mathrm{_cw}(m_w))^T] \\
&= {R} _{cw}\mathrm{E}[(\mathbb{p}_w - m_w)(\mathbb{p}_w-m_w)^T] R _{cw}^T \\
&= R _{cw} \Sigma R _{cw}^T \\
\end{aligned}
\tag{3.2}
$$

式(1.2)を用いて、カメラ座標系の点を画面に投影し、投影点の共分散は以下のように計算できます。

$$
\begin{aligned}
\Sigma^{\prime}
&= \mathrm{E}[(\mathrm{u}(\mathbb{p}_c)-\mathrm{u}(m_c))(\mathrm{u}(\mathbb{p}_c)-\mathrm{u}(m_c))^T]
\end{aligned}
\tag{3.3}
$$

ここで、$m_c$は$\mathbb{p}_c$の平均であり、$\mathbb{p}_c - m_c$は比較的に小さい値ですので、$\delta$として定義します。式(1.2)の$\mathrm{u}$は非線形関数ですが、$\mathrm{u}$のヤコビ行列$J$を用いた以下のような近似計算が存在します。

$$
\begin{aligned}
&\mathrm{u}(m_c + \delta) = \mathrm{u}(m_c) + J \delta　\\
&\mathrm{u}(\mathbb{p}_c) - \mathrm{u}(m_c) = J(\mathbb{p}_c - m_c)
\end{aligned}
\tag{3.4}
$$

式(3.4)と(3.2)を(3.3)に代入すると、2Dの共分散は以下のように計算できます。

$$
\begin{aligned}
\Sigma^{\prime}
&= \mathrm{E}[J(\mathbb{p}_c-m_c)(\mathbb{p}_c-m_c)^TJ^T] \\
&= J\Sigma_c J^T \\
&= J R _{cw}\Sigma R _{cw}^T J^T \\
\end{aligned}
\tag{3.5}
$$

以上より、式(3)を導出しました。

ここのヤコビ行列$J$の詳細について、以前の記事で議論したことがあります。(バンドル調整の式(8))。

$$
J = \begin{bmatrix}
f_x/x & 0 & -{f_x x}/{z^2} \\
0 & f_y/z & -{f_y y}/{z^2}
\end{bmatrix}
\tag{3.6}
$$

(4). カメラの姿勢に応じて、見えた3D Gaussianの色を計算

自然界では、異なる光源から放射された光は、物体の表面に当たると、特定の波長の光が吸収されたり、または屈折や反射によって方向を変えて、別の物体にもう一度当たったりします。これらの光が最終的に融合し、目に入り、私たちが多彩な現実世界を見ることができます。

レイトレーシングという技術は、その名の通り、光の経路を一本一本追跡し、現実の光をシミュレーションすることで、よりリアルな世界を描写することができますが、計算量が非常に多くなります。

3D Gaussian Splattingでは、レイトレーシングを直接に利用しませんが、球面調和関数（spherical harmonics）を使用して複雑な光の反射特性を近似できる能力を持っています。少ない計算量で、高いリアリティな光表現ができます。

この球面調和関数を簡単に説明すると、球面版のフーリエ級数です。フーリエ級数は、無限個の三角関数の基底とそれぞれの係数の線形結合で、任意の関数を再現することができます。一方、球面調和関数は、無限個の球面基底とそれぞれの係数の線形結合で、視点方向に依存し、色が変化できる球を表現することができます。

以下は、球面調和関数を用いて、色の計算式です。

$$
\mathrm{c}(r, h) = \sum_{l=0}^{l_{max}}{\sum_{m=-l}^{l}{h_{lm}\mathrm{Y}_{l}^{m}(r)}}
\tag{4}
$$

ここで、$r$はカメラと3D Gaussianの間の単位方向ベクトルです。$Y_l^m$はl次元m番の球面基底です。$Y_l^m$の詳細についてはreal spherical harmonics tableを参照してください。

$$
r = \frac{p_w - t_{wc}}{\lVert p_w - t_{wc} \rVert }
\tag{4.1}
$$

球面調和関数を視覚的に理解するために、地球の表面を近似するデモを作成しました。球面調和関数の低い次元は低周波数成分を表し、$l_{max}$が0のとき、単色な球体しか表現できませんが、$l_{max}$が高くなるほど、高周波成分に対する表現能力が高くなり、地球らしさが増していきます。ただし、$l_{max}$が大きくなると必要なパラメータも増えるため、品質とパラメータ数のバランスを考慮して、3D Gaussian Splattingの実装では$l_{max}$を3に制限しています。

"The ground truth Earth image is modified from URL. By Solar System Scope. Licensed under CC-BY-4.0"

注意：このデモでは、球面調和関数を視覚的に理解しやすくするために、各部分の色が異なる球体で表現されていますが、実際の3D Gaussianでは、全体の色は一つだけであり、視点の変化によってその色が変化します。

(5). 2Dの画像を生成

上記の1から4の計算から、以下の情報が得られました。このステップでは、これらの情報を使用して、2D画像のすべてのピクセルの最終色（RGBの値）を計算し、画像を生成します。

$u$: 3D Gaussianの中心を2D画面に投影するときの座標値(by. 式(1.1)(1.2))
$\Sigma^{\prime}$: 3D Gaussianをカメラ視点に投影したときに見えた2D共分散(by. 式(3))
$c$: カメラの姿勢に応じて、見えた3D Gaussianの色 (by. 式(4))
$\alpha$: この3D Gaussian固有の不透明度 (by. 入力パラメータ)

2D画面上の各ピクセルのRGB値は、そのピクセルに投影されたすべての3D Gaussianの重なり合いによって決まります。上図のように、後ろの物体から発射された光が手前の物体を透過すると、その光は弱くなり、ピクセルのRGB値への寄与が減少します。

j番目ピクセルのRGB値$\gamma_{j}$ は式(5)で計算できます。3D Gaussianの前後関係を正しく処理するため、式(5)の前処理として、カメラの方向から見た深度(Depth)に基づいて3D Gaussianをソートします。つまり、3D Gaussianのインデックスiが小さいほど、カメラとの距離が近くなります。

$$
\begin{aligned}
\gamma_{j}
&= \sum_{i \in N} \alpha_{ij}^{\prime}(x_{j}, \alpha_i, \sigma_i^{\prime}) c_i \tau_{ij}
\end{aligned}
\tag{5}
$$

$\alpha_{ij}^{\prime}$は、i番目の3D Gaussianがピクセルjにおける不透明度を表し、以下の式で計算します。

$$
\alpha_{ij}^{\prime}(x_{j}, \alpha_i, \sigma_i^{\prime})= \begin{aligned}
\exp(- 0.5 (u_{i}-x_{j}) \Sigma_i^{\prime-1} (u_{i}-x_{j})^T)\alpha_i
\end{aligned}
\tag{5.1}
$$

この計算は、i番目の3D Gaussianの固有の$\alpha$に基づいて、Gaussianの中心のピクセル値($u_{i}$)と描画したいピクセル($x_{j}$)とのマハラノビス距離が遠いほど、このピクセルでの不透明度が小さくなることを示しています。

$\tau_{ij}$は、i番目の3D Gaussianが前方の3D Gaussianを透過して、ピクセルjに到達した際の減衰を示す係数です。この値が1の場合、減衰がなくそのまま表示されます。一方、0の場合は完全に減衰し、画面から完全に見えなくなります。

$$
\begin{aligned}
&\tau_{ij} = \prod_{k=1}^{i-1} (1 - \alpha_{kj}^{\prime}) \\
&\tau_{1j}= 1
\end{aligned}
\tag{5.2}
$$

これにより、前後関係による減衰を考慮した上で、ピクセルに投影したすべての3D Gaussianの色を線形結合して、ピクセルの色値を計算できるようになりました。すべてのピクセルを処理すると、最終的に2D画像を生成することができました。

デモ

検証のために、式(1)から(5)の計算をPythonで実装し、学習済みの3D Gaussianモデルを2D画面に投影するデモを作成しました。

このデモでは、深度順に3D Gaussianを2D画面に投影し、少しずつ2D画像を生成しました。

まとめ

本文では、「3D Gaussian Splatting」の基本的な考え方を説明し、学習済みの3D Gaussianを使用して2D画像を生成する手法について議論しました。残る課題は2つあります。

画面のピクセル数が多く、3D Gaussianの数も多いため、そのままでは処理が遅くなります。この処理をGPUが得意とするラスタライズ法（rasterize）によって、大幅に加速することができます。
「元の写真からどうやって、3D Gaussianを作るの？(学習方法)」という質問まだ回答していません。

本文が長くなりすぎたので、次回は上記の課題について解説します。

「3D Gaussian Splatting」を徹底的に理解するため、数式ドキュメント、それぞれアルゴリズムのデモ、専用ビューア、およびゼロからの再実装を含むプロジェクトをGithubで公開しています。ぜひご参考にしてください。

驚くほどキレイな三次元シーン復元、「3D Gaussian Splatting」を徹底的に解説する