球面調和関数変換の理論と実装

こんにちは。球面調和関数とはいったい何なのか。ふんわりとした理論なら分かるという人は多いと思います。今回はもう少し理論的な部分や実装を踏まえた理解をしたいという方向けに記事を書いてみました。

#目次
1.はじめに
2.球面調和関数
3.球面調和関数変換
4.補足と用いたパラメータ
5.おわりに

#1. はじめに
球面調和関数は主に量子力学の世界で使われている関数です。しかし、数値計算の高速化やデータ量の圧縮を目的としてスペクトル法を用いた大気シミュレーションやコンピュータグラフィックスのPRT法を用いたレンダリングにも使用されています。今回は私の専門に合わせて、CGを意識したモンテカルロ法による球面調和関数変換についてご紹介します。

球面調和関数はよく上図右側のようにウニのような形の関数だと説明されます。しかし、この説明はCGの世界で生きる人たちにとっては少しわかりにくいと感じると思うので、別の表現を用いて説明します。球面調和関数は上図の左側のように塗り分けられた球面用テクスチャだと解釈できます。

レイトレーシングを用いたレンダリングでは、ある一点の色を計算する時、膨大な計算量になってしまいます。映画などの映像作品では問題ないですが、ゲームなど１フレームの計算時間がミリ秒単位になってしまう場合には、実用的な手法とは言えなくなってしまいます。そこで、ある一点についてのレイトレーシング結果を事前に下図の青線で描かれたような半球テクスチャに保存しておく手法が開発されました。ただし、あらゆる位置xについてテクスチャを保存しておくと膨大なデータ量になってしまいます。球面調和関数の展開係数として保存しておくことでデータ量を数百分の一に圧縮できるため、有用な手法だと言えます。

#2. 球面調和関数
この節では、実装を見据えて数式を用いながら球面調和関数について解説します。

緯度$\theta$と経度$\phi$は以上のような定義を与えています。

用いる数式群

P_n(x) = \frac{1}{2^n n!} \bigg(\frac{d}{dx}\bigg)^n (x^2 - 1)^n \tag{1}

P^m_n (x) = (1-x^2)^{\frac{m}{2}}\bigg(\frac{d}{dx}\bigg)^m P_n(x) \tag{2}

X_{nm} (\theta) =  \biggl[ \frac{2n+1}{4\pi}\frac{(n-|m|)!}{(n+|m|)!} \biggr]^{\frac{1}{2}}P^{|m|}_n(\cos\theta)\tag{3}

Y_{nm}(\theta,\phi)\equiv
\begin{cases}
     \sqrt{2}X_{n|m|}(\theta) ~\cos (m\phi)~~~~(-n<m<0)\\
     X_{n|m|}(\theta)~~~~~~~~~~~~~~~~~~~~~~~~~(m=0)\\
     \sqrt{2}X_{nm}(\theta) ~\sin(m\phi)~~~~(-n<m<0)  \tag{4}
\end{cases}

P^{m+1}_{n+1}(x) = x P^{m+1}_{n}(x)~+~(n+m+1)(1-x^2)^{\frac{1}{2}} P^{m}_{n}(x) \tag{5}

P^{n}_{n}(x) = \frac{(2n)!}{2^n n!}(1-x^2)^{\frac{n}{2}} \tag{6}

P^{m}_{n+1}(x) = \frac{1}{n-m+1} \Big\{(2n+1)xP^m_n(x) -(n+m)P^{m}_{n-1}(x) \Big\} \tag{7}

P^0_0(x) = 1,~~P^0_1(x)=x,~~P^1_1(x)=\sqrt{1-x^2},~~P^1_2(x)=3x\sqrt{1-x^2} \tag{8}

以上の(1)～(8)の式、漸化式を利用すれば、球面調和関数を実装できる。各式について簡単に説明する。
・(1)：ルジャンドル多項式
・(2)：ルジャンドル培関数
・(3)(4)：球面調和関数
・(5)(6)(7)(8)：低波数成分から高波数成分を導出するための漸化式と条件

//Factorial
float Factorial(int n) {
	int nm = n;
	float Out = n;
	if (nm == 0) {
		return 1;
	}
	else if (nm == 1) {
		return 1;
	}
	while (nm > 1) {
		nm = nm - 1;
		Out *= nm;
	}
	return Out;
}

//(8),(7)
float Pnpm(int n, int m, float ct) {
	if (n == 0 && m==0) {
		return 1.0;
	}
	else if (n == 1 && m == 0) {
		return ct;
	}
	else if (n == 1 && m == 1) {
		return sqrt(1 - ct * ct);
	}
	else if (n == 2 && m == 1) {
		return 3 * ct * sqrt(1 - ct * ct);
	}
	else {
		int nm = n - 1;
		return ((2 * nm + 1) *ct* Pnpm(nm, m, ct) - (nm + m) * Pnpm(nm - 1,m, ct)) /float(nm - m + 1);
	}
}

//(6)
float Pnn(int n, float ct) {
	return Factorial(2 * n) * pow(1 - ct * ct, 0.5 * n) / (pow(2, n) * Factorial(n));
}

//(5)
float Pnmp(int n, int m, float ct) {
	if (m == 0) {
		return Pnpm(n, 0, ct);
	}
	else if (m == 1) {
		return Pnpm(n, 1, ct);
	}
	else if (n == 0) {
		return Pnpm(0, 0, ct);
	}
	else if (n == m) {
		return Pnn(n, ct);
	}
	else {
		int mm = m - 1;
		int nm = n - 1;
		return ct * Pnmp(nm, m, ct) + (nm + mm + 1) * sqrt(1 - ct * ct) * Pnmp(nm, mm, ct);
	}
}

//(3)
float Xnm(int n, int m, float ct) {
	float S;
	S = (2 * n + 1) / (4 * acos(-1))*Factorial(n-abs(m))/float(Factorial(n+abs(m)));
	S = sqrt(S);

	float P;
	P = Pnmp(n, m, ct);
	
	if (m >= 0) {
		return S * P * pow(-1, m);
	}
	else {
		return S * P;
	}
}

//(4)
float Ynm(int n, int m, float theta, float phi) {
	float ct = cos(theta);
	if (m < 0) {
		return sqrt(2) * Xnm(n, -m, ct) * cos(m * phi);
	}
	else if(m==0){
		return Xnm(n, m, ct);
	}
	else {
		return sqrt(2) * Xnm(n, m, ct) * sin(m * phi);
	}
}

データの配列順序

球面調和関数の展開係数は以下のような配列順で格納する。(変数xの関数であることを省略する。)

P^0_0,~P^0_1,\cdots,P^0_N,

P^1_1,~P^{-1}_1,~P^1_2,~P^{-1}_2,~\cdots,P^1_N,~P^{-1}_N,

\vdots

P^M_M,~P^{-M}_M,\cdots,~P^M_N,~P^{-M}_N

配列番号から($n,m$)を求める関数は以下のように書ける。

void knm(int k, int& n, int& m,int MMax) {
	int l = MMax + 1 - 1;
	if (k <= l) {
		m = 0;
		n = k;
		return;
	}

	for (int i = 1; i <= MMax; i++) {
		int lp = l;
		l += 2 * (MMax+1 - i);
		if (k <= l) {
			if (k % 2 == 0) {
				m = i;
			}
			else m = -i;

			n = (k - lp + 1) / 2 + abs(m)-1;
			return;
		}
	}

}

#3. 球面調和関数変換
以下の天球画像を正変換することで球面調和展開係数に変換した後で、逆変換して画像に戻す過程を実装した。

モンテカルロ法を用いた正変換

手順
➀天球上にランダムにで一様な($\theta, \phi$)をサンプリングする。以下のリンクが参考になります。
https://qiita.com/aa_debdeb/items/e416ae8a018692fc07eb

u=\frac{\phi}{2\pi},~~v=\frac{\theta}{\pi}

➁画像のu,vに対応したピクセルからカラー(C_i)を取得する。

\hat{C}_{nm} = \frac{2\pi}{MN} \sum^{MN}_{i=1} C_i Y_{nm}(\theta, \phi)

➂すべての展開係数に対して、上記のようにサンプリング計算する。

float u = phi / (2 * PI);
float v = theta/PI; 

for (int k = 0; k < NumMM; k++) {
     int n=0;
     int m=0;
     knm(k, n, m,MMax);
     vec3 sampleCol = color_data->value(u, v);
     Chat[k] += 2 * PI * sampleCol * Ynm(n, m, theta, phi) /(MX*NY);
}

逆変換による画像生成

手順
➀各ピクセル位置($u,v$)から($\theta, \phi$)を計算
➁そのピクセル位置における全スペクトルからの寄与の総和を計算。

C_i = \sum^{N}_{n=0} \sum^{n}_{m=-|n|} \hat{C}_{nm}~ Y_{nm}(\theta, \phi)

Ci=0
for (int k = 0; k < NumMM; k++) {
    knm(k, n, m,MMax);
    Ci += Chat[k] * Ynm(n, m, theta, phi);
}

#4. 補足と用いたパラメータ
・Mは切断波数と呼ばれており、CG界の大域照明ではM=9が一般的であると言われている。
・一般的な手法ではM=Nであり、展開係数の数は全部で$(M+1)^2$個

#5. おわりに
今回の記事で球面調和関数の正変換と逆変換を扱った。CGではこれを発展させた様々な応用があるので、今後紹介していきたい。