An Area-Preserving Parametrization for Spherical Rectangles を読む
An Area-Preserving Parametrization for Spherical Rectangles
Carlos Ureña, Marcos Fajardo, Alan King , “An Area-Preserving Parametrization for Spherical Rectangles” は、 レイトレ合宿6 に提出しました MoodyRender で組み込んだものの、あまりしっかりとしたまとめを書いていませんでしたので、こちらにまとめておこうとおもいます。 こちらの論文は、SolidAngle社が EGSR 2013で発表した論文となっております。この論文では、単位面積 から、任意の長方形の領域を単位球面に射影してできた領域への写像を、面積を保存する形で提示します。文章で書くとややこしい印象になってしまいますが、イメージは以下のような感じです。
これは一度定義してしまえば、射影立体角内で、立体角一様なサンプルを生成したり、準乱数によるサンプルを生成したりも簡単にできます。このテクニックは特に直接照明の計算、いわゆるNext Event Estimationに便利です。
立体角でサンプリングしたくなる理由
Next Event Estimationを行う際、最も簡単なのは面積に一様なサンプルを生成することでしょう。しかしなぜわざわざややこしいことをして、立体角でサンプリングを生成するのでしょうか? 例えば以下のようなシチュエーションを考えてみましょう
これは結構横長な光源があり、そして手前側の半球の中心が現在光線追跡中の点で、ここの直接照明を考えるとしましょう。 面積に一様なサンプリングをした場合には赤い点のように生成されるでしょう。なんだか良さげですが、これを立体角に投射してみましょう。
...かなり偏っています。これが運良く、重点サンプリングのように機能するなら良いでしょうが、多くのシチュエーションではそんなことはないでしょう。これでは分散が増加してしまうのが目に見えています。なのでせめて立体角には一様なサンプリングをしたくなります。例えば以下のように。
そんなわけで、これをどのように達成しているのかを見ていきましょう。
写像
まず最初に面積を保存する写像とは何かを定義します。論文 3. Properties of map M では、以下のような式があります。
この式の意味を考えます。 ここではまず、単位面積の中の任意の領域 T を考えて、その面積を とします。
そして、写像する立体角の全領域をQ として、その立体角を とします。
そして T の写像を とすると、
全立体角あたりの、写像してきた立体角は、やはりTと一致する、これが面積を保存する写像というわけです。
座標系
次に座標系についてです。この論文では、もろもろ後の問題を単純にするために、長方形領域 P を基準にローカル座標系を定義します。 これは論文のFigure 1にあたります。
※ 元論文 Figure 1
原点oを基準に、長方形 P 上にサンプルを生成するのをイメージします。 この座標系の重要なところは 、
- の向きが x軸
- の向きが y軸
- 長方形 P の原点o側の法線がz軸
となるような座標系であるというところです。 それぞれの軸の単位ベクトルは、長方形Pと、原点oが決まれば自動的に決まってくることになります。 ここでローカル座標系での そして もやはり一意に決まりますので、ここで整理します。
sのローカル座標は、
であるので、座標軸に単に投射すれば良いため、 は、原点o から 長方形Pへのベクトルを d 、 を基底ベクトルとすると、
そしてそこから長方形のサイズ分ずらせば良いので は、
となります。 長方形の4点については重要な頂点になりますので、名前を と名前をつけておくことにします。
球面四角形の面積
ところで長方形の領域が射影された領域の立体角すなわち面積はどのように求めることができるでしょうか? 一般に単位球面三角形について、3つの角の角度から求めることができることが知られています。(証明は今回すっ飛ばします)
平面上の通常の三角形の内角の和はもちろん180° () でありますが、球面三角形の内角の和は膨らみがある分 以上になることは想像に難くないでしょう。 一方でこの式は、球面三角形の内角の和が のときどうなるかというと、
面積が0になってしまいました。どういうことかといえば、球面三角形の内角の和が に近づいていくということは、それだけ平面上の三角形に近づいていくことであり、結局それは以下の図のようにどんどん小さくなっていってしまうわけです。
では逆に球面三角形の最大値はどんな値になるでしょうか?これを考えるために、上の図と"反対側"の球面三角形を考えてみます。
この場合、それぞれの外側の角度は、 となるわけなので、面積は、
となり、さらに の場合は、
のように小さくなっていき、式としても、
なるほど、単位球の表面積と一致しています。どうやらこれが最大値のようです。だいぶイメージも掴めてきました。 ところで、球面四角形の場合はどうかといえば、二つの球面三角形に分割して考えることができます。
したがって、
この時2つの角については、一つにまとめることができるので、
と4つの角の角度で考えることができます。 この角度で を再定義すると、
なんとシンプルでしょうか。これで立体角について簡単に考えていくことができそうです。 余談なのですが、ここまでの話は、自然に任意の球面N角形に拡張することができ、それはもちろん、
のようになります。 この拡張、以下のようにどこか適当な1点は実は球面N角形の1頂点だったと考えてみても、ちゃんと成り立ちます。
式の左の項を見ると、 N が1増えると面積は 増えそうですが、それは右側の項でちゃんと相殺されますね。 なかなか面白いものです。
を求める
と立体角には密接な関係性があることがわかりました。そこで をそれぞれ求めたいわけです。座標系のセクションですでに定義した、原点o と、 の4つの頂点で作られる平面の法線から考えます。
ベクトルの外積からそれぞれ、
のように求めることができると、図からわかります。 法線がわかってしまえば、 の値は、以下の図をみながら、
となります。
ここでついでにひとつ重要な話題があって、今回はローカル座標系としてうまいこと軸が長方形に整列した座標系を定義しました。
これを踏まえて面について後ろから観察すると、
と、このようになっています。このことから、
- は Y軸にいつも直交しており、ローカル座標系でのy 成分は常に0
- は X軸にいつも直交しており、ローカル座標系でのx 成分は常に0
であることがわかります。これは実装の最適化にも寄与します。
このことから、それぞれの法線はそれぞれ1つの角度の値で表現できることもわかります。
※元論文では、 だけなぜかマイナスで定義されています。理由がちょっと定かではないのですが、この部分はさほど重要な点ではないように思いましたので、ここではシンプルにすべて+での定義を使いました。
u 方向の累積分布関数
いよいよコア部分に入っていきます。単位面積におけるx座標を u 、y座標をvとパラメータ化することを考えます。
ここでは u について考えます。 uは当然 0 ~ 1 の間を行き来することになります。そのようなシチュエーションを考えた時、
のとき、
そして、
のとき、
もっと進んで、任意のuで、
が成り立つような を考えます。ここではこれは ローカル座標系でのxに対応させて考えると、
とこのような領域を定義したいです。 uについてアニメーションさせると、以下のようになるでしょう。
ただ多くのシチュエーションでは上の図のように綺麗な形をしておらず、例えば以下のように歪んでいるということは、一応確認しておきたいところです。例えば以下のように。
なのでuの値が線形に変化するとしても、当然面積が線形にはならないということです。
では の面積について考えてみます。元論文 Figure 2を見てみましょう。 今までの議論でイメージができていれば、とても立体的な図として見えるはずです。
※元論文 Figure 2
この中で、多くのものは定義済みですが、以下青で丸をつけた部分については未定義でした。そして uに依存して変化する部分でもあります。
を計算したいということから考えると、 というのはやはり必要な値であり、それは から求めます。 とは何か?といえば、uの値に応じて、 から へ回転する法線であり、境界線の法線でもあります。 上から見るとわかりやすいです。
以上より、 の面積は以下のように書き下せることがわかります。
の角度パラメータ1つによる表現はすでに上の方で触れました。
これと同様に、 についてもやはり同様の表記が可能で、
と、このようになります。ここで、内積 と、について考えてみます。 は今回のローカル座標系においては、
と、このような値になるでしょう。 同様に、 もローカル座標系においては、
となるわけです。つまり結局のところ、内積 と、についてはローカル座標系でz成分のみの掛け算により、内積が得られます。 したがって、 は、
と記述することができ、元論文では、
とおいており、
と短く記述でき、さらに
から、
これが元論文の式 (8) になります。これはまるで逆関数法における、累積分布関数のようにも感じられるでしょう。 ちなみに、
- は、 のローカル座標系のz成分
- は、 のローカル座標系のz成分
であるとも整理できます。
次に は、まだこの段階ではどういった表現で得られるかはわかりませんが、仮に都合良く が求まりさえすれば、射影する前の長方形の座標(ローカル座標系)の x成分 を得ることができます。
は、上から見ると、次の図のようになります。
したがって、これを90度回転させると、原点から境界線へと向かうベクトルが手に入ります。
ここで、長方形までの距離 と、 求めたい は、次の図のような配置になっており、
したがって、
これで を得ることができました。 ただ、
については、 や、 では符号が逆になってしまうので、注意したいところではあります。
とはいえ、 は、明らかに、
なので、心配する必要はありませんでした。
というわけで、最終的には都合の良い を見つけることに焦点があたります。 ※ と書いていますが、実際には とuの関数として記述したほうが直観的かもしれません。
を求める
こちらは元論文の Appendices, "A. Expressing cu as a function of u" の部分になります。 は定数であり、 だけが未知のため、これを求めることが目標となります。 論文の説明が十分わかりやすいので、ここはかなりストレートに追いかけます。
まず定数をまとめて表記を簡単にするため、関数 を導入します。
これにより、
はともにコサインであり、範囲は -1 ~ 1をはみ出ることはないので、
の中の引き算を外に引きずり出すため、ここで加法定理、
で、
すこし複雑化してしまいましたが、ここで、 は、
コサインからサインを求める形になるため、単位円の上半分を表しており、
から、
であり、とすると、
これを少し変形して、
こうすると、かなりの部分が から分離でき、
とおいてしまうと、
とかなり簡単になります。 あとは について解きます。ここも逆関数法みたいな香りがしますね。
おっと、符号が行方不明になってしまいました。 ただ、
をよく見ると、
なのは明らかです。したがって、 と は同じ符号でなければなりません。
加えて、
であるため、以上をふまえると、 は、
であるべきでしょう。 長くなってしまったため、以上をまとめると、
となります。
の除算について
上記の式変形の中で、 で除算する部分があります。元論文中では、 は常に 0より小さいことが書かれています。しかしこれは事実ではありません。というのも、これの根拠にしているのが、 がそれぞれ よりも大きいこと、と記述されておりますが、これが真実ではありません。反例として、例えば以下のようなシチュエーションがあります。
これは明らかに です。ということは自然と も となり得るということになります。
ただ、論文中で を示したかった理由は単にゼロ除算を避けるためにあります。 では実際にゼロ除算は実装上問題になるのでしょうか?
であるため、まず
になります。次に、
なので、
になります。最後に、
より、
そう、これは結局単に の地点であり、図示すると、
となり、異常な特異点というわけではないことがわかります。 ただ、ゼロ除算に近づいていく過程で0付近は数値的に誤差が生じやすいかもしれません。
この点については、知人の紹介で、Carlos Ureña氏とAlan King氏にメールで確認させて頂いたところ、確かにerrorであると確認していただけました。
なにが解決したか?
の へのマッピングが以上で明らかになりました。これで何が実現できたのかを再確認します。 は連続的であり、イメージがややしにくいため、一旦均等に極めて細かく面積を分割することを考えます。 例えば極めて細かく、均等に4つに分割しますと、
この時立体角上に写像した a, b, c, dの面積(=立体角) も均等となります。
このようなマッピングが成立したことになります。
ということは、残る 方向については、以下のように経度方向に細かく分割した中で残りの縦方向の写像を考えれば良いことになります。
のマッピングどうするか?
ではここで、θに対して をマッピングするのは正しいでしょうか?
これは正しくありません。 なぜならθに均一に分割しても、それぞれの面積(=立体角)は一様ではないためです。 これは、
のような立体角積分でも として出現することからも明らかです。 これをなんとかするために、球とそれをすっぽりと覆う円柱の側面積の面積の関係性を使います。
ここでいう面積の関係性は、 横方向にスライスした部分の面積は、水平に切る限りどこでも一致する、というものです。(例えば赤い部分) ※これについて、単位球面に対してランダムに点を生成するでは、より詳細な説明をしています。
Y軸に対して面積は線形であり、そしてこの関係性は、経度方向に薄くスライスしても、同様の性質が保持されることは明らかです。さらに言えば、今回のケースは上下に切られている状況ではありますが、これも関係がありません。 ということは結局のところY軸に対して を線形にマッピングするだけで良いということになります。
の詳細
マッピングのアイディアがわかってしまえば、あとは幾何学的な部分を整理するだけになります。 元論文のFigure 3 を見てみましょう。
※元論文のFigure 3
d は、xz平面上のp までの距離であり、
で求められます。
と はどんな値かというと、oから および に向かうベクトルを正規化したあとのy成分で求めることができます。
ここの間を線形にマッピングすれば良いため、
と が 確定します。 ここから、 のときとほぼ同様に、
と を求めることができました。
グローバル座標へ
はローカル座標系での値ですので、最後にグローバル座標へ戻して終了です。
実装
以上を踏まえての実装です。
まずはナイーヴなほうです。
次に最適化したバージョンです。まだできることが残っていたらすみません。。
PathTracing + Next Event Estimation 比較
ナイーヴなPathTracing + Next Event Estimation の実装で、 面積に一様な場合 ( 64 spp ) と、立体角に一様な場合 ( 64 spp ) の比較です。
PathTracing + Next Event Estimation + MIS(パワーヒューリスティック) 比較
BRDFのサンプリングと光源のサンプリングの2戦略でのMISをした場合、 面積に一様な場合 ( 32 spp ) と、立体角に一様な場合 ( 32 spp ) の比較です。
MISなしのに比べると、差は少ないですが、確かにこの場合でも分散が低くなっています。特に光源付近では、MISウエイトによって分散の高いサンプルを十分に潰しきれていないことがわかります。
サンプリング生成コストが少し増加するので、同じsppで比較するのは完全にはフェアではありません。しかし単純なシーンならともかく、複雑なシーンではサンプリング生成コストの増加は相対的に影響度は少なく、この分散の低下のメリットが遥かに大きいと思います。 論文では、極めて最適化された実装で、プロダクションのシーンで、数パーセントのパフォーマンスの低下が見られたが、本手法のノイズ低減効果を考えると十分回収できるものだったとありました。
サンプリング手法と組み合わせる
また、今回の手法では、他のサンプリング手法と組み合わせるのが容易です。 例えばBlueNoiseのようなサンプリングからマッピングすると以下のような結果が得られます。
こちらのサンプルデモは以下にあります。
まとめ
この手法は、法線や頂点の数がそこそこ多く出てくる上、3Dなのでイメージするがやや難しいですが、実際にプロットしつつぐるぐる3Dビューを回して確認さえできれば、あとはそれほど高度な数学も登場してこないので、それほど難易度も高くなく、実装もシンプルで最適化もすれば十分効率的です。座標系のとりかたがうまいことがとても良く効いているように感じます。また分散低減は効果もわかりやすく、まだやったことがない方は試してみる価値は十分にあるのではないでしょうか。ただいったん三角形でいいや、という方は先行研究である、James Arvo, "Stratified sampling of spherical triangles" のほうを先に試したほうが良いかもしれません。ただ長方形の場合にそれを三角形二枚にしてサンプリングするのは、本手法と等価ではなく、論文によるとわずかに本手法のほうが分散が少ないようですね。
サンプリング対象が長方形限定であるのは、純粋にリミテーションであり、今後拡張される余地もあるでしょう。(もちろん三角形ポリゴンのバージョンはすでにありますが)また、cos項まで考慮したり、放射輝度が均一でない光源や、BSDFなども同時に考慮に入れれたら、さらにすごいかもしれません。
参考文献
Carlos Ureña, Marcos Fajardo, Alan King , “An Area-Preserving Parametrization for Spherical Rectangles”