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で発表した論文となっております。この論文では、単位面積  \left[0\, 1 \right]^2 から、任意の長方形の領域を単位球面に射影してできた領域への写像を、面積を保存する形で提示します。文章で書くとややこしい印象になってしまいますが、イメージは以下のような感じです。

f:id:ushiostarfish:20180924165141p:plain

これは一度定義してしまえば、射影立体角内で、立体角一様なサンプルを生成したり、準乱数によるサンプルを生成したりも簡単にできます。このテクニックは特に直接照明の計算、いわゆるNext Event Estimationに便利です。

立体角でサンプリングしたくなる理由

Next Event Estimationを行う際、最も簡単なのは面積に一様なサンプルを生成することでしょう。しかしなぜわざわざややこしいことをして、立体角でサンプリングを生成するのでしょうか? 例えば以下のようなシチュエーションを考えてみましょう

f:id:ushiostarfish:20180919233606p:plain

これは結構横長な光源があり、そして手前側の半球の中心が現在光線追跡中の点で、ここの直接照明を考えるとしましょう。 面積に一様なサンプリングをした場合には赤い点のように生成されるでしょう。なんだか良さげですが、これを立体角に投射してみましょう。

f:id:ushiostarfish:20180919234011p:plain

...かなり偏っています。これが運良く、重点サンプリングのように機能するなら良いでしょうが、多くのシチュエーションではそんなことはないでしょう。これでは分散が増加してしまうのが目に見えています。なのでせめて立体角には一様なサンプリングをしたくなります。例えば以下のように。

f:id:ushiostarfish:20180919234504p:plain

そんなわけで、これをどのように達成しているのかを見ていきましょう。

写像

まず最初に面積を保存する写像とは何かを定義します。論文 3. Properties of map M では、以下のような式があります。

 
\displaystyle {
A\left( T \right) =\frac { 1 }{ \sigma \left( Q \right)  } \sigma \left( M\left( T \right)  \right) 
}

この式の意味を考えます。 ここではまず、単位面積の中の任意の領域 T を考えて、その面積を  A\left( T \right) とします。

f:id:ushiostarfish:20180919235500p:plain

そして、写像する立体角の全領域をQ として、その立体角を  \sigma \left( Q \right) とします。 f:id:ushiostarfish:20180919235750p:plain

そして T の写像 M \left( T \right) とすると、 f:id:ushiostarfish:20180919235950p:plain

全立体角あたりの、写像してきた立体角は、やはりTと一致する、これが面積を保存する写像というわけです。 f:id:ushiostarfish:20180920000128p:plain

座標系

次に座標系についてです。この論文では、もろもろ後の問題を単純にするために、長方形領域 P を基準にローカル座標系を定義します。 これは論文のFigure 1にあたります。

f:id:ushiostarfish:20180920232954p:plain ※ 元論文 Figure 1

原点oを基準に、長方形 P 上にサンプルを生成するのをイメージします。 この座標系の重要なところは 、

  •  e_x の向きが x軸
  •  e_y の向きが y軸
  • 長方形 P の原点o側の法線がz軸

となるような座標系であるというところです。 それぞれの軸の単位ベクトルは、長方形Pと、原点oが決まれば自動的に決まってくることになります。 ここでローカル座標系での  x_0, x_1, y_0, y_1 そして  z_0もやはり一意に決まりますので、ここで整理します。

sのローカル座標は、

 
\displaystyle {
s=\left( x_{ 0 },y_{ 0 },z_{ 0 } \right) 
}

であるので、座標軸に単に投射すれば良いため、  x_0, y_0, z_0 は、原点o から 長方形Pへのベクトルを d 、 x,y, z を基底ベクトルとすると、

 
\displaystyle {
x_{ 0 }=d\cdot x\\ y_{ 0 }=d\cdot y\\ z_{ 0 }=d\cdot z
}

そしてそこから長方形のサイズ分ずらせば良いので  x_1, y_1 は、

 
\displaystyle {
x_{ 1 }=x_{ 0 }+\left\| e_{ x } \right\| \\ y_{ 1 }=y_{ 0 }+\left\| e_{ y } \right\| 
}

となります。 長方形の4点については重要な頂点になりますので、名前を  v_{00}, v_{10}, v_{01}, v_{11} と名前をつけておくことにします。

f:id:ushiostarfish:20180921000625p:plain

球面四角形の面積

ところで長方形の領域が射影された領域の立体角すなわち面積はどのように求めることができるでしょうか? 一般に単位球面三角形について、3つの角の角度から求めることができることが知られています。(証明は今回すっ飛ばします)

f:id:ushiostarfish:20180922115249p:plain

 
\displaystyle {
A\left( R \right) =\gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 }-\pi 
}

平面上の通常の三角形の内角の和はもちろん180° ( =\pi) でありますが、球面三角形の内角の和は膨らみがある分 \pi 以上になることは想像に難くないでしょう。 一方でこの式は、球面三角形の内角の和が  \pi のときどうなるかというと、

 
\displaystyle {
A\left( R \right) =\pi -\pi =0
}

面積が0になってしまいました。どういうことかといえば、球面三角形の内角の和が  \pi に近づいていくということは、それだけ平面上の三角形に近づいていくことであり、結局それは以下の図のようにどんどん小さくなっていってしまうわけです。

f:id:ushiostarfish:20180922115959p:plain

では逆に球面三角形の最大値はどんな値になるでしょうか?これを考えるために、上の図と"反対側"の球面三角形を考えてみます。

f:id:ushiostarfish:20180922121025p:plain

この場合、それぞれの外側の角度は、 2\pi - \gamma となるわけなので、面積は、

 
\displaystyle {
A\left( R' \right) =\left( 2\pi -\gamma _{ 0 } \right) +\left( 2\pi -\gamma _{ 1 } \right) +\left( 2\pi -\gamma _{ 2 } \right) -\pi \\ =6\pi -\left( \gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 } \right) -\pi \\ =5\pi -\left( \gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 } \right) 
}

となり、さらに  \gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 } = \pi の場合は、

f:id:ushiostarfish:20180922121633p:plain

のように小さくなっていき、式としても、

 
\displaystyle {
A\left( R' \right) =5\pi -\pi =4\pi 
}

なるほど、単位球の表面積と一致しています。どうやらこれが最大値のようです。だいぶイメージも掴めてきました。 ところで、球面四角形の場合はどうかといえば、二つの球面三角形に分割して考えることができます。

f:id:ushiostarfish:20180922122550p:plain

したがって、

 
\displaystyle {
A\left( R_{ 1 } \right) +A\left( R_{ 2 } \right) \\ =\gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 }-\pi +\gamma '_{ 0 }+\gamma '_{ 1 }+\gamma '_{ 2 }-\pi \\ =\gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 }+\gamma '_{ 0 }+\gamma '_{ 1 }+\gamma '_{ 2 }-2\pi 
}

この時2つの角については、一つにまとめることができるので、

 
\displaystyle {
=\gamma _{ 0 }+\gamma '_{ 0 }+\left( \gamma _{ 1 }+\gamma '_{ 1 } \right) +\left( \gamma _{ 2 }+\gamma '_{ 2 } \right) -2\pi 
}

と4つの角の角度で考えることができます。 この角度で \gamma を再定義すると、

f:id:ushiostarfish:20180922112659p:plain

 
\displaystyle {
A\left( Q \right) =\gamma _{ 0 }+\gamma _{ 1 }+\gamma _{ 2 }+\gamma _{ 3 }-2\pi 
}

なんとシンプルでしょうか。これで立体角について簡単に考えていくことができそうです。 余談なのですが、ここまでの話は、自然に任意の球面N角形に拡張することができ、それはもちろん、

 
\displaystyle {
A_{ n }=\sum _{ i=0 }^{ N-1 }{ \gamma _{ i } } -\left( N-2 \right) \pi 
}

のようになります。 この拡張、以下のようにどこか適当な1点は実は球面N角形の1頂点だったと考えてみても、ちゃんと成り立ちます。

f:id:ushiostarfish:20180922125021p:plain

式の左の項を見ると、 N が1増えると面積は  \pi 増えそうですが、それは右側の項でちゃんと相殺されますね。 なかなか面白いものです。

 \gamma_i を求める

 \gamma_i と立体角には密接な関係性があることがわかりました。そこで  \gamma_i をそれぞれ求めたいわけです。座標系のセクションですでに定義した、原点o と、 v_{00}, v_{10}, v_{01}, v_{11} の4つの頂点で作られる平面の法線から考えます。

f:id:ushiostarfish:20180922132732p:plain

f:id:ushiostarfish:20180922133157p:plain

ベクトルの外積からそれぞれ、

 
\displaystyle {
n_{ 0 }=\frac { v_{ 00 }\times v_{ 10 } }{ \left\| v_{ 00 }\times v_{ 10 } \right\|  } \\ n_{ 1 }=\frac { v_{ 10 }\times v_{ 11 } }{ \left\| v_{ 10 }\times v_{ 11 } \right\|  } \\ n_{ 2 }=\frac { v_{ 11 }\times v_{ 01 } }{ \left\| v_{ 11 }\times v_{ 01 } \right\|  } \\ n_{ 3 }=\frac { v_{ 01 }\times v_{ 00 } }{ \left\| v_{ 01 }\times v_{ 00 } \right\|  } 
}

のように求めることができると、図からわかります。 法線がわかってしまえば、 \gamma_i の値は、以下の図をみながら、

f:id:ushiostarfish:20180922173231p:plain

 
\displaystyle {
\gamma _{ 0 }=\arccos { \left( -{ n }_{ { 0 } }\cdot { n }_{ { 1 } } \right)  } \\ \gamma _{ 1 }=\arccos { \left( -{ n }_{ { 1 } }\cdot { n }_{ { 2 } } \right)  } \\ \gamma _{ 2 }=\arccos { \left( -{ n }_{ { 2 } }\cdot { n }_{ { 3 } } \right)  } \\ \gamma _{ 3 }=\arccos { \left( -{ n }_{ { 3 } }\cdot { n }_{ { 0 } } \right)  } 
}

となります。

ここでついでにひとつ重要な話題があって、今回はローカル座標系としてうまいこと軸が長方形に整列した座標系を定義しました。

f:id:ushiostarfish:20180922134212p:plain

これを踏まえて面について後ろから観察すると、

f:id:ushiostarfish:20180922135806p:plain

と、このようになっています。このことから、

  •  n_1, n_3 は Y軸にいつも直交しており、ローカル座標系でのy 成分は常に0
  •  n_0, n_2 は X軸にいつも直交しており、ローカル座標系でのx 成分は常に0

であることがわかります。これは実装の最適化にも寄与します。

このことから、それぞれの法線はそれぞれ1つの角度の値で表現できることもわかります。

f:id:ushiostarfish:20180922143446p:plain

 
\displaystyle {
{ n }_{ { 3 } }=\cos  \varphi _{ { 0 } }{ z }+\sin  \varphi _{ { 0 } }{ x }\\ { n }_{ { 1 } }=\cos  \varphi _{ { 1 } }{ z }+\sin  \varphi _{ { 1 } }{ x }
}

f:id:ushiostarfish:20180922142752p:plain

 
\displaystyle {
{ n }_{ { 0 } }=\cos  \theta _{ { 0 } }{ z }+\sin  \theta _{ { 0 } }{ y }\\ { n }_{ { 2 } }{ =\cos  \theta _{ { 1 } }{ z }+\sin  \theta _{ { 1 } }{ y } }
}

※元論文では、 n_3 だけなぜかマイナスで定義されています。理由がちょっと定かではないのですが、この部分はさほど重要な点ではないように思いましたので、ここではシンプルにすべて+での定義を使いました。

u 方向の累積分布関数

いよいよコア部分に入っていきます。単位面積におけるx座標を u 、y座標をvとパラメータ化することを考えます。

f:id:ushiostarfish:20180922161555p:plain

ここでは u について考えます。 uは当然 0 ~ 1 の間を行き来することになります。そのようなシチュエーションを考えた時、

 u = 0 のとき、 A\left( Q_{ u } \right) =0

そして、

 u = 1 のとき、 A\left( Q_{ u } \right) =A\left( Q \right)

もっと進んで、任意のuで、

 
\displaystyle {
A\left( Q_{ u } \right) =A\left( Q \right) u
}

が成り立つような  Q_{ u } を考えます。ここではこれは ローカル座標系でのxに対応させて考えると、

f:id:ushiostarfish:20180922162545p:plain

とこのような領域を定義したいです。 uについてアニメーションさせると、以下のようになるでしょう。

f:id:ushiostarfish:20180922163018g:plain

ただ多くのシチュエーションでは上の図のように綺麗な形をしておらず、例えば以下のように歪んでいるということは、一応確認しておきたいところです。例えば以下のように。 f:id:ushiostarfish:20180922163501p:plain

なのでuの値が線形に変化するとしても、当然面積が線形にはならないということです。

では  Q_u の面積について考えてみます。元論文 Figure 2を見てみましょう。 今までの議論でイメージができていれば、とても立体的な図として見えるはずです。

f:id:ushiostarfish:20180922165047p:plain ※元論文 Figure 2

この中で、多くのものは定義済みですが、以下青で丸をつけた部分については未定義でした。そして uに依存して変化する部分でもあります。

f:id:ushiostarfish:20180922165345p:plain

 Q_u を計算したいということから考えると、 \gamma '_{ 0 },\gamma '_{ 1 } というのはやはり必要な値であり、それは  m_u から求めます。  m_u とは何か?といえば、uの値に応じて、 -n_3 から  n_1 へ回転する法線であり、境界線の法線でもあります。 上から見るとわかりやすいです。

f:id:ushiostarfish:20180922170435p:plain

以上より、 Q_u の面積は以下のように書き下せることがわかります。

 
\displaystyle {
A\left( Q_{ u } \right) =\gamma '_{ 0 }+\gamma '_{ 1 }+\gamma _{ 2 }+\gamma _{ 3 }-2\pi \\ A\left( Q_{ u } \right) =\arccos { \left( -{ n }_{ { 0 } }\cdot { m }_{ { u } } \right)  } +\arccos { \left( -{ n }_{ { 2 } }\cdot { m }_{ { u } } \right)  } +\gamma _{ 2 }+\gamma _{ 3 }-2\pi 
}

 n_1, n_3 の角度パラメータ1つによる表現はすでに上の方で触れました。

 
\displaystyle {
{ n }_{ { 0 } }=\cos  \theta _{ { 0 } }{ z }+\sin  \theta _{ { 0 } }{ y }\\ { n }_{ { 2 } }{ =\cos  \theta _{ { 1 } }{ z }+\sin  \theta _{ { 1 } }{ y } }
}

これと同様に、 m_u についてもやはり同様の表記が可能で、

 
\displaystyle {
m_{ { u } }=\cos  \varphi _{ { u } }{ z }+\sin  \varphi _{ { u } }{ x }
}

と、このようになります。ここで、内積  { n }_{ { 0 } }\cdot { m }_{ { u } } と、{ n }_{ { 2 } }\cdot { m }_{ { u } }について考えてみます。  n_0, n_1 は今回のローカル座標系においては、

 
\displaystyle {
{ n }_{ { 0 } }=\left( 0,a,b \right) \\ { n }_{ { 2 } }=\left( 0,c,d \right) \\ 
}

と、このような値になるでしょう。 同様に、 m_u もローカル座標系においては、

 
\displaystyle {
{ m }_{ { u } }=\left( g,0,h \right) 
}

となるわけです。つまり結局のところ、内積  { n }_{ { 0 } }\cdot { m }_{ { u } } と、{ n }_{ { 2 } }\cdot { m }_{ { u } }についてはローカル座標系でz成分のみの掛け算により、内積が得られます。 したがって、 A \left( Q_u \right) は、

 
\displaystyle {
A\left( Q_{ u } \right) =\arccos { \left( -\cos  \theta _{ { 0 } }\cos  \varphi _{ { u } } \right)  } +\arccos { \left( -\cos  \theta _{ { 1 } }\cos  \varphi _{ { u } } \right)  } +\gamma _{ 2 }+\gamma _{ 3 }-2\pi 
}

と記述することができ、元論文では、

 
\displaystyle {
b_{ 0 }=\cos  \theta _{ { 0 } }\\ b_{ 1 }=\cos  \theta _{ { 1 } }\\ c_{ u }=\cos  \varphi _{ { u } }
}

とおいており、

 
\displaystyle {
A\left( Q_{ u } \right) =\arccos { \left( -b_{ 0 }c_{ u } \right)  } +\arccos { \left( -b_{ 1 }c_{ u } \right)  } +\gamma _{ 2 }+\gamma _{ 3 }-2\pi 
}

と短く記述でき、さらに

 
\displaystyle {
A\left( Q_{ u } \right) =A\left( Q \right) u
}

から、

 
\displaystyle {

A\left( Q \right) u=\arccos { \left( -b_{ 0 }c_{ u } \right)  } +\arccos { \left( -b_{ 1 }c_{ u } \right)  } +\gamma _{ 2 }+\gamma _{ 3 }-2\pi 
}

これが元論文の式 (8) になります。これはまるで逆関数法における、累積分布関数のようにも感じられるでしょう。 ちなみに、

  •  b_0 は、 n_0 のローカル座標系のz成分
  •  b_1 は、 n_2 のローカル座標系のz成分

であるとも整理できます。

次に  c_u は、まだこの段階ではどういった表現で得られるかはわかりませんが、仮に都合良く  c_u が求まりさえすれば、射影する前の長方形の座標(ローカル座標系)の x成分  x_u を得ることができます。

 m_u は、上から見ると、次の図のようになります。 f:id:ushiostarfish:20180922184543p:plain

したがって、これを90度回転させると、原点から境界線へと向かうベクトルが手に入ります。

f:id:ushiostarfish:20180922185920p:plain

ここで、長方形までの距離  z_0 と、 求めたい x_u は、次の図のような配置になっており、

f:id:ushiostarfish:20180922190628p:plain

したがって、

 
\displaystyle {
\begin{align} x_{ u } & =z_{ 0 }\tan { \left( \varphi _{ { u } }+\frac { \pi  }{ 2 }  \right)  }  \\  & =-z_{ 0 }\frac { 1 }{ \tan { \left( \varphi _{ { u } } \right)  }  }  \\  & =-z_{ 0 }\frac { \cos { \left( \varphi _{ { u } } \right)  }  }{ \sin { \left( \varphi _{ { u } } \right)  }  }  \\  & =-z_{ 0 }\frac { \cos { \left( \varphi _{ { u } } \right)  }  }{ \sqrt { 1-\cos ^{ 2 }{ \left( \varphi _{ { u } } \right)  }  }  }  \\  & =-z_{ 0 }\frac { c_{ u } }{ \sqrt { 1-{ c_{ u } }^{ 2 } }  }  \end{align}
}

これで  x_u を得ることができました。 ただ、

 
\displaystyle {
\sin { \left( \varphi _{ { u } } \right)  } =\sqrt { 1-\cos ^{ 2 }{ \left( \varphi _{ { u } } \right)  }  } 
}

については、 -\pi \lt \varphi _{ { u } }\lt0 や、 \pi \lt\varphi _{ { u } }\lt2\pi  では符号が逆になってしまうので、注意したいところではあります。

f:id:ushiostarfish:20180922192857p:plain

とはいえ、 \varphi _{ { u } }+\frac { \pi  }{ 2 } は、明らかに、

 
\displaystyle {
\frac { \pi  }{ 2 } <\varphi _{ { u } }+\frac { \pi  }{ 2 } <\frac { \pi  }{ 2 } +\pi \\ 0<\varphi _{ { u } }<\pi 
}

なので、心配する必要はありませんでした。

というわけで、最終的には都合の良い  c_u を見つけることに焦点があたります。 ※  c_u と書いていますが、実際には  c_u \left(u \right) とuの関数として記述したほうが直観的かもしれません。

 
\displaystyle {
A\left( Q \right) u=\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  } +\arccos { \left( -b_{ 1 }c_{ u }\left( u \right)  \right)  } +\gamma _{ 2 }+\gamma _{ 3 }-2\pi 
}

 c_u \left(u \right) を求める

こちらは元論文の Appendices, "A. Expressing cu as a function of u" の部分になります。  A\left( Q \right), b_0, b_1, \gamma_2, \gamma_3 は定数であり、 c_u \left(u \right) だけが未知のため、これを求めることが目標となります。 論文の説明が十分わかりやすいので、ここはかなりストレートに追いかけます。

まず定数をまとめて表記を簡単にするため、関数  \phi \left( x \right) を導入します。

 
\displaystyle {
\phi \left( x \right) =xA\left( Q \right) -\gamma _{ 2 }-\gamma _{ 3 }+2\pi 
}

これにより、

 
\displaystyle {
A\left( Q \right) u=\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  } +\arccos { \left( -b_{ 1 }c_{ u }\left( u \right)  \right)  } +\gamma _{ 2 }+\gamma _{ 3 }-2\pi \\ A\left( Q \right) u-\gamma _{ 2 }-\gamma _{ 3 }+2\pi -\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  } =\arccos { \left( -b_{ 1 }c_{ u }\left( u \right)  \right)  } \\ \phi \left( u \right) -\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  } =\arccos { \left( -b_{ 1 }c_{ u }\left( u \right)  \right)  } 
}

 b_0, c_u はともにコサインであり、範囲は -1 ~ 1をはみ出ることはないので、

f:id:ushiostarfish:20180923114946p:plain

 
\displaystyle {
\cos { \left( \phi \left( u \right) -\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  } =\cos { \left( \arccos { \left( -b_{ 1 }c_{ u }\left( u \right)  \right)  }  \right)  } \\ \cos { \left( \phi \left( u \right) -\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  } =-b_{ 1 }c_{ u }\left( u \right) 
}

 \cos の中の引き算を外に引きずり出すため、ここで加法定理、

 
\displaystyle {
\cos { \left( \alpha -\beta  \right)  } =\cos { \alpha  } \cos { \beta  } +\sin { \alpha  } \sin { \beta  } 
}

で、

 
\displaystyle {
\begin{align} \cos { \left( \phi \left( u \right) -\arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  }  & =-b_{ 1 }c_{ u }\left( u \right)  \\ \cos { \left( \phi \left( u \right)  \right)  } \cos { \left( \arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  } +\sin { \left( \phi \left( u \right)  \right)  } \sin { \left( \arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  }  & =-b_{ 1 }c_{ u }\left( u \right)  \\ \cos { \left( \phi \left( u \right)  \right)  } \left( -b_{ 0 }c_{ u }\left( u \right)  \right) +\sin { \left( \phi \left( u \right)  \right)  } \sin { \left( \arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  }  & =-b_{ 1 }c_{ u }\left( u \right)  \\ -\cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) +\sin { \left( \phi \left( u \right)  \right)  } \sin { \left( \arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  }  & =-b_{ 1 }c_{ u }\left( u \right)  \end{align}
}

すこし複雑化してしまいましたが、ここで、 \sin\left(\arccos\left(x\right)\right) は、

f:id:ushiostarfish:20180923115955p:plain

コサインからサインを求める形になるため、単位円の上半分を表しており、

 
\displaystyle {
{ x }^{ 2 }+{ y }^{ 2 }=1
}

から、

 
\displaystyle {
\sin  \left( \arccos  \left( x \right)  \right) =\sqrt { 1-x^{ 2 } } 
}

f:id:ushiostarfish:20180923120709p:plain

であり、とすると、

 
\displaystyle {
\begin{align} -\cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) +\sin { \left( \phi \left( u \right)  \right)  } \sin { \left( \arccos { \left( -b_{ 0 }c_{ u }\left( u \right)  \right)  }  \right)  }  & =-b_{ 1 }c_{ u }\left( u \right)  \\ -\cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) +\sin { \left( \phi \left( u \right)  \right)  } \sqrt { 1-\left( -b_{ 0 }c_{ u }\left( u \right)  \right) ^{ 2 } }  & =-b_{ 1 }c_{ u }\left( u \right)  \\ -\cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) +\sin { \left( \phi \left( u \right)  \right)  } \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & =-b_{ 1 }c_{ u }\left( u \right)  \end{align}
}

これを少し変形して、

 
\displaystyle {
\begin{align} -\cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) +\sin { \left( \phi \left( u \right)  \right)  } \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & =-b_{ 1 }c_{ u }\left( u \right)  \\ \sin { \left( \phi \left( u \right)  \right)  } \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & =\cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) -b_{ 1 }c_{ u }\left( u \right)  \\ \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & =\frac { \cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }c_{ u }\left( u \right) -b_{ 1 }c_{ u }\left( u \right)  }{ \sin { \left( \phi \left( u \right)  \right)  }  }  \\ \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & =c_{ u }\left( u \right) \frac { \cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }-b_{ 1 } }{ \sin { \left( \phi \left( u \right)  \right)  }  }  \end{align}
}

こうすると、かなりの部分が  c_u \left(u \right) から分離でき、

 
\displaystyle {
f\left( u \right) =\frac { \cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }-b_{ 1 } }{ \sin { \left( \phi \left( u \right)  \right)  }  } 
}

とおいてしまうと、

 
\displaystyle {
\sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } } =c_{ u }\left( u \right) f\left( u \right) 
}

とかなり簡単になります。 あとは  c_u \left(u \right) について解きます。ここも逆関数法みたいな香りがしますね。

 
\displaystyle {
\begin{align} \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & =c_{ u }\left( u \right) f\left( u \right)  \\ 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } & =\left( c_{ u }\left( u \right)  \right) ^{ 2 }{ \left( f\left( u \right)  \right)  }^{ 2 } \\ \frac { 1 }{ \left( c_{ u }\left( u \right)  \right) ^{ 2 } } -b_{ 0 }^{ 2 } & ={ \left( f\left( u \right)  \right)  }^{ 2 } \\ \frac { 1 }{ \left( c_{ u }\left( u \right)  \right) ^{ 2 } }  & ={ \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } \\ \left( c_{ u }\left( u \right)  \right) ^{ 2 } & =\frac { 1 }{ { \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } }  \\ c_{ u }\left( u \right)  & =\pm \frac { 1 }{ \sqrt { { \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } }  }  \end{align}
}

おっと、符号が行方不明になってしまいました。 ただ、

 
\displaystyle {
\sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } } =c_{ u }\left( u \right) f\left( u \right) 
}

をよく見ると、

 
\displaystyle {
0\le \sqrt { 1-b_{ 0 }^{ 2 }\left( c_{ u }\left( u \right)  \right) ^{ 2 } } 
}

なのは明らかです。したがって、 c_{ u } \left( u \right)   f\left( u \right) は同じ符号でなければなりません。

加えて、

 
\displaystyle {
0\le \frac { 1 }{ \sqrt { { \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } }  } 
}

であるため、以上をふまえると、 c_{ u } \left( u \right)  は、

 
\displaystyle {
c_{ u }\left( u \right) =\frac { sign\left( f\left( u \right)  \right)  }{ \sqrt { { \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } }  } 
}

であるべきでしょう。 長くなってしまったため、以上をまとめると、

 
\displaystyle {
\begin{align} \phi \left( x \right)  & =xA\left( Q \right) -\gamma _{ 2 }-\gamma _{ 3 }+2\pi  \\ f\left( u \right)  & =\frac { \cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }-b_{ 1 } }{ \sin { \left( \phi \left( u \right)  \right)  }  }  \\ c_{ u }\left( u \right)  & =\frac { sign\left( f\left( u \right)  \right)  }{ \sqrt { { \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } }  }  \end{align}
}

となります。

 \sin { \left( \phi \left( u \right)  \right)  } の除算について

上記の式変形の中で、 \sin { \left( \phi \left( u \right)  \right)  } で除算する部分があります。元論文中では、 \sin { \left( \phi \left( u \right)  \right)  } は常に 0より小さいことが書かれています。しかしこれは事実ではありません。というのも、これの根拠にしているのが、 \gamma '_0, \gamma '_1 がそれぞれ  \frac { \pi  }{ 2 }  よりも大きいこと、と記述されておりますが、これが真実ではありません。反例として、例えば以下のようなシチュエーションがあります。

f:id:ushiostarfish:20180923131654p:plain

これは明らかに  \gamma _0 \lt \frac { \pi  }{ 2 }  , \gamma _1 \lt \frac { \pi  }{ 2 } です。ということは自然と \gamma '_0, \gamma '_1 \gamma' _0 \lt \frac { \pi  }{ 2 }  , \gamma' _1 \lt \frac { \pi  }{ 2 } となり得るということになります。

ただ、論文中で \sin { \left( \phi \left( u \right)  \right)  } \lt 0 を示したかった理由は単にゼロ除算を避けるためにあります。 では実際にゼロ除算は実装上問題になるのでしょうか?

 
\displaystyle {
f\left( u \right) =\frac { \cos { \left( \phi \left( u \right)  \right)  } b_{ 0 }-b_{ 1 } }{ \sin { \left( \phi \left( u \right)  \right)  }  } 
}

であるため、まず

 
\displaystyle {
f\left( u \right) =+inf
}

になります。次に、

 
\displaystyle {
c_{ u }\left( u \right) =\frac { sign\left( f\left( u \right)  \right)  }{ \sqrt { { \left( f\left( u \right)  \right)  }^{ 2 }+b_{ 0 }^{ 2 } }  } 
}

なので、

 
\displaystyle {
c_{ u }\left( u \right) = 0
}

になります。最後に、

 
\displaystyle {
x_{ u }=-z_{ 0 }\frac { c_{ u } }{ \sqrt { 1-{ c_{ u } }^{ 2 } }  } 
}

より、

 
\displaystyle {
x_{ u }=0
}

そう、これは結局単に x = 0 の地点であり、図示すると、

f:id:ushiostarfish:20180923133526p:plain

となり、異常な特異点というわけではないことがわかります。 ただ、ゼロ除算に近づいていく過程で0付近は数値的に誤差が生じやすいかもしれません。

この点については、知人の紹介で、Carlos Ureña氏とAlan King氏にメールで確認させて頂いたところ、確かにerrorであると確認していただけました。

なにが解決したか?

 u x_u へのマッピングが以上で明らかになりました。これで何が実現できたのかを再確認します。  u は連続的であり、イメージがややしにくいため、一旦均等に極めて細かく面積を分割することを考えます。 例えば極めて細かく、均等に4つに分割しますと、

f:id:ushiostarfish:20180924185614p:plain

この時立体角上に写像した a, b, c, dの面積(=立体角) も均等となります。

 
\displaystyle {
a'=b'=c'=d'\Leftrightarrow a=b=c=d
}

このようなマッピングが成立したことになります。

ということは、残る  v 方向については、以下のように経度方向に細かく分割した中で残りの縦方向の写像を考えれば良いことになります。 f:id:ushiostarfish:20180923152144p:plain

 vマッピングどうするか?

ではここで、θに対して  vマッピングするのは正しいでしょうか?

f:id:ushiostarfish:20180923153341p:plain

これは正しくありません。 なぜならθに均一に分割しても、それぞれの面積(=立体角)は一様ではないためです。 これは、

 
\displaystyle {
\int _{ \Omega  }^{  }{ d\omega  } =\int _{ 0 }^{ 2\pi  }{ \int _{ 0 }^{ \frac { \pi  }{ 2 }  }{ \sin { \theta  } d\theta d\phi  }  } 
}

のような立体角積分でも  \sin \theta として出現することからも明らかです。 これをなんとかするために、球とそれをすっぽりと覆う円柱の側面積の面積の関係性を使います。

f:id:ushiostarfish:20180923154604p:plain

ここでいう面積の関係性は、 横方向にスライスした部分の面積は、水平に切る限りどこでも一致する、というものです。(例えば赤い部分) ※これについて、単位球面に対してランダムに点を生成するでは、より詳細な説明をしています。

Y軸に対して面積は線形であり、そしてこの関係性は、経度方向に薄くスライスしても、同様の性質が保持されることは明らかです。さらに言えば、今回のケースは上下に切られている状況ではありますが、これも関係がありません。 ということは結局のところY軸に対して  v を線形にマッピングするだけで良いということになります。

 v の詳細

マッピングのアイディアがわかってしまえば、あとは幾何学的な部分を整理するだけになります。 元論文のFigure 3 を見てみましょう。

f:id:ushiostarfish:20180923160526p:plain ※元論文のFigure 3

d は、xz平面上のp までの距離であり、

 
\displaystyle {
d=\sqrt { { x_{ u } }^{ 2 }+{ z_{ 0 } }^{ 2 } } 
}

で求められます。

 h_0 h_1 はどんな値かというと、oから  a_0 および a_1 に向かうベクトルを正規化したあとのy成分で求めることができます。

 
\displaystyle {
h_{ 0 }=\frac { y_{ 0 } }{ \left\| a_{ 0 } \right\|  } =\frac { y_{ 0 } }{ \sqrt { { d }^{ 2 }+{ y_{ 0 } }^{ 2 } }  } \\ h_{ 1 }=\frac { y_{ 1 } }{ \left\| a_{ 1 } \right\|  } =\frac { y_{ 1 } }{ \sqrt { { d }^{ 2 }+{ y_{ 1 } }^{ 2 } }  } 
}

ここの間を線形にマッピングすれば良いため、

 
\displaystyle {
h_{ v }=h_{ 0 }+v\left( h_{ 1 }-h_{ 0 } \right) 
}

 h_{ v } が 確定します。 ここから、 c_u のときとほぼ同様に、

f:id:ushiostarfish:20180923163328p:plain

 
\displaystyle {
y_{ v }=d\tan { \theta  } =d\frac { \sin { \theta  }  }{ \cos { \theta  }  } =d\frac { \sin { \theta  }  }{ \sqrt { 1-\sin ^{ 2 }{ \theta  }  }  } \\ y_{ v }=d\frac { h_{ v } }{ \sqrt { 1-{ h_{ v } }^{ 2 } }  } 
}

 y_{ v } を求めることができました。

グローバル座標へ

 x_u, y_v はローカル座標系での値ですので、最後にグローバル座標へ戻して終了です。

 
\displaystyle {
o+x_{ u }x+y_{ v }y+z_{ 0 }z
}

実装

以上を踏まえての実装です。

まずはナイーヴなほうです。

次に最適化したバージョンです。まだできることが残っていたらすみません。。

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のようなサンプリングからマッピングすると以下のような結果が得られます。

f:id:ushiostarfish:20180929161549p:plain

f:id:ushiostarfish:20180929161557p:plain

こちらのサンプルデモは以下にあります。

github.com

まとめ

この手法は、法線や頂点の数がそこそこ多く出てくる上、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”