HIYOKOのページ

球面調和関数展開

球面調和関数変換したキューブマップ画像を正距円筒図法に展開します
変換後データの確認にどうぞ
表示は700x300で行いますが、計算終了後に設定した大きさでダウンロードできます

アプリページ

高速化バージョン

球面調和関数計算部分を高速化しました
が、それによって誤差も生じているため高次になるとノイズが発生します
よって旧バージョンも残しておきます
高速バージョンのアプリページ

サンプルデータ

画像変換が面倒くさい！そもそもキューブマップ画像とかねーし！
でも遊んでみたい、試してみたい
そんなあなたへl=300まで変換したデータを用意しました
良ければどうぞ
非営利目的に限り展開以外の用途(研究など)に使ってみてもいいです
(使えるのか…?)

高速化について

今回行った球面調和関数、及び展開の高速化についてやや詳しく
正距円筒図法で表せば変換のほうにも応用可能なはず

　球面調和関数が遅い原因はルジャンドル陪関数の計算が大きなウェイトを占めています
漸化式によってルジャンドル陪関数を計算するため、次数が上がるほどループが大きくなります

ここで、球面調和関数の中でのルジャンドル陪関数の意味合いと言うものを考えて見ます
数式書くのが面倒なので概念的にざっくりと
球面調和関数では波の形を球座標のθとφ方向に分けて考えることが出来ます
地図みたいに考えると緯度方向にあたるのがθで経度にあたるのがφですね
そして球面調和関数では、θ方向の波をルジャンドル陪関数でφ方向の波を三角関数で表すことが出来ます
ルジャンドル陪関数部分はθ依存であり、θが同じ角度である(緯度が同じ)ところでは球面調和関数を求める時常にルジャンドル陪関数部分の結果が一定です

とすれば球面調和関数の各次数lにおいて、各オーダーmに対して予めθ方向の波であるルジャンドル陪関数を求めておけば、各々のオーダーにおける球面調和関数の計算時に逐一高負荷なルジャンドル陪関数を計算せずに済みます
つまるところ、ルジャンドル陪関数のルックアップテーブル化を図ります
これを用いる事で球面調和関数の計算は定数と三角関数とテーブル値の掛け算だけとなり非常に高速に計算できるようになります

ルックアップテーブルについては必要充分な分解能をもつ大きさを用意しないといけません
今回の実装では正距円筒図法のように展開するため、θは求める画像の高さに相当します
そこでルックアップテーブルの大きさは画像の高さの2倍の大きさで作ることにしました
(とはいえ、高次になるとノイズが目立ってきたためこの大きさでは充分では無いのかもしれません。というより、ルックアップテーブルの大きさはθ依存では無く次数lに依って決めるべきなのかも?）
そして、ただ参照しただけでは最大誤差が大きくなってしまうのでルックアップテーブルを補間して使います
今回の実装では高速性重視のため、また計算が簡便であるために線形補間を使いましたが、3次スプライン補間などを使うとより誤差を減らせるかもしれません
(3次スプライン補間ならば係数が4つとなるため、それぞれをRGBAに格納して使うと無駄なく高性能なはず)

　ここまではルジャンドル陪関数に対しての高速化でしたがここからは展開そのものの高速化に移ります
球面調和関数展開する際には変換によって得られた各オーダーに対するウェイトを球面調和関数に掛けて和を求めていきます
このとき、各次数lにおいてオーダーmの数は-l<m<lである事からl*2+1個となります
つまり各次数においてl*2+1回の球面調和関数の計算が必要です

ここで、今度はθでは無くφ方向についての球面調和関数の特性を利用します
球面調和関数はそのオーダーがm>0の時はφに対してcos(|m|*φ)という波を持ちます
対してm<0の時はsin(|m|*φ)という波を持ちます
どちらも三角関数で波長も同じです
つまりmと-mのφに対する波形は同じであり、位相がズレているだけなのです
そしてθに対してはルジャンドル陪関数に対して|m|として用いられるためmも-mも同じ特性を持ちます
と言うことはどちらか片方の球面調和関数を求めれば、φに対して位相をズラせば計算しなおさなくても答えが分かるのです

今回の実装ではGPGPUであるため球面調和関数の計算結果はテクスチャと言う形で保持されます
なのでテクスチャをリピートさせてm>0で求めたテクスチャを参照をズラして再利用すればm<0の計算にそのまま使えるのです

よって展開においてmと-mを交互に計算すれば、球面調和関数の計算は各次数につきl+1回の計算で済むことになり計算回数を約半分にすることが出来ます


　高速化については大体こんな感じです
ただ、展開自体の計算量は変わるものではなく高次元になると結局多大な時間がかかります
というわけで高速球面調和関数展開というよりは球面調和関数展開の最適化という表現の方が正しいかな?
とはいえ計算速度そのものは愚直な実装に比べて飛躍的に上がっているハズ
一応今のところ展開にだけこの高速化を適用しているが、変換の球面調和関数計算部分には共通点が多いので正距円筒図法にキューブマップを変換してしまえばほぼ同じ手法で高速化できるはず
いずれ暇があれば変換も高速化してみます
ただ変換については、キューブマップの方が空間的にサンプリング点が比較的均等に配置されていて、正距円筒図法だと極に集中して赤道付近が疎であるためテクスチャ解像度が無駄になるかも
あと、立体角もキューブマップはそれなりに均一であるのに対して、正距円筒図法だと極付近が極端に小さく赤道付近との差が大きくて、総和とるときに情報落ちしないか心配
スピードは上がっても精度が大きく下がりそうではある