球面均勻采樣

寫games101作業七的時候，計算間接光照的時候需要對着色點所在的半球面進行一次均勻采樣來獲得其他物體到着色點的入射光的方向。但是框架中的采樣函數我實在看不懂

Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){
	switch(m_type){
		case DIFFUSE:
		{
			// uniform sample on the hemisphere
			float x_1 = get_random_float(), x_2 = get_random_float();
			float z = std::fabs(1.0f - 2.0f * x_1);
			float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;
			Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);
			return toWorld(localRay, N);
			break;
		}
	}
}

為了搞清楚這里采樣的原理，搜索了一下相關知識。
本文參考了球諧光照與PRT學習筆記（二）：蒙特卡洛積分與球面上的均勻采樣和球面采樣 Spherical Sampling以及逆變換采樣 (inverse transform sampling) 的原理這三篇博客

對球面的均勻采樣

由於我們采樣的結果是着色點的入射光方向，因此我們可以通過球坐標系的兩個分量：$ (\theta,\phi) $來表示采樣結果。

直接對\((\theta,\phi)\)進行均勻采樣

直接在\((\theta,\phi)\)進行均勻采樣即：

\[\begin{cases} \theta=\xi_x,&\xi_x \sim U(0，1)\\ \phi=\xi_y,&\xi_y \sim U (0,1) \end{cases} \]

但是這樣得到的采樣結果不是均勻分布

為什么采樣結果不均勻呢？原作者提供了一種直觀的理解方式：球面可以看作垂直於\(z\)軸的圓的累加，對\(\phi\)采樣即是對這些圓的采樣，由於這些圓的半徑不同，根據圓弧長度與半徑和角度之間的關系(如下圖所示)圓上采樣點之間的距離也不同，因此球面上的采樣並不均勻

面積微元 \(dA=r^2sin\theta d\theta d\phi\)，當\(sin\theta\)小的時候立體角也小，面積微元也更小，因此采樣點更加密集，隨着\(sin\theta\)的增大，立體角也不斷增大，面積微元也更大，因此采樣點更加稀疏。所以采樣結果在球面的兩極更加密集。因此直接對\((\theta,\phi)\)均勻采樣無法得到對球面均勻采樣的結果。

將均勻隨機變量映射到\((\theta,\phi)\)進行采樣

我們已知由於立體角在不同的\(\theta\)下大小並不相同，這就造成了不同角度下面積微元內的樣本密度不同，我們需要構造新的映射方式：

\[\begin{cases} \theta=f(\xi_x),&\xi_x \sim U(0,1)\\ \phi=g(\xi_y),&\xi_y \sim U(0,1) \end{cases} \]

使得在不同角度下的面積微元的樣本密度保持一致。
考慮對單位球面采樣的情況下，設\(p(v)\)為采樣的概率密度函數，\(v\)為采樣方向，理想狀態下即均勻采樣下\(p(v)\)是一個常數，由於單位球面立體角為\(4\pi\)，即\(p(v)=1/4\pi\)。單位球面下面積微元內的采樣概率與單位立體角內的采樣概率相等，即:

\[P(A)=\int p(v)dw=1/4\pi dw=1/4\pi sin\theta d\theta d\phi \]

用\((\theta,\phi)\)表達的概率密度函數為\(p(\theta,\phi)=dP(A)/d\theta d\phi=1/4\pi sin\theta.\)
現在我們有了理想狀態下的概率分布函數\(p(\theta,\phi)\)，還需要求出\((\theta,\phi)\)與\((\xi_x,\xi_y)\)之間的映射關系，我們可以通過Inverse Transform Sampling來得到二者之間的映射關系，即通過均勻分布變量\((\xi_x,\xi_y)\),構造出\((\theta,\phi)\leftarrow(F^{-1}_\theta(\xi_x),F^{-1}_\phi(\xi_y))\),使得\(\theta \sim F_\theta(\theta),\phi \sim F_\phi(\phi).\)因此只要求出\(F^{-1}_\theta(\theta)\)與\(F^{-1}_\phi(\phi)\)即可得到映射關系。我們求出\(\theta\)與\(\phi\)的邊緣概率密度：

\[\begin{cases} f(\theta)=\int_{0}^{2\pi} p(\theta,\phi)d\phi=sin\theta/2\\ f(\phi)=\int _{0}^\pi p(\theta,\pi)d\theta=1/2\pi \end{cases} \]

根據邊緣概率分布函數我們可以求出邊緣分布函數：

\[\begin{cases} F_\theta(\theta)=\int_{0}^\theta f(\hat \theta)d\hat \theta={1-cos\theta\over 2}\\ F_\phi(\phi)=\int _{0}^\phi f(\hat \phi)d\hat \phi=\phi/2 \end{cases} \]

這樣就得到了\((\theta,\phi)\)與\((\xi_x,\xi_y)\)之間的映射關系：

\[\begin{cases} \theta=arccos(1-2\xi_x)\\ \phi=2\pi\xi_y \end{cases} \]

將球面換成半球面求得的結果即為框架中的采樣函數的映射方式。

Inverse Transform Sampling的證明

給定隨機變量\(X \sim P_X(x),Y=f(X)\),當f(X)為遞增函數時,隨機變量\(Y\)的分布函數\(F_Y(y)\)為：

\[F_Y(y)=Pr(Y\le y)=Pr(X\le f^{-1}(y))=F_X(f^{-1}(y)). \]

隨機變量\(Y\)的概率密度函數為：

\[P_Y(y)={dF_Y(y)\over dy}={dF_X(f^{-1}(y))\over dy}=P_X(f^{-1}(y)){df^{-1}(y)\over dy}=P_X(x){dx\over dy}. \]

當\(f(x)\)為遞減函數時,隨機變量\(Y\)的分布函數為：

\[F_Y(y)=Pr(Y\le y)=Pr(X\ge f^{-1}(y))=1-Pr(X\le f^{-1}(y))=1-F_X(f^{-1}(y)). \]

概率密度函數為：

\[P_Y(y)={dF_Y(y)\over dy}={-dF_X(f^{-1}(y))\over dy}=-P_X(f^{-1}(y)){df^{-1}(y)\over dy}=-P_X(x){dx\over dy}. \]

綜上：

\[P_Y(y)=P_X(x)|{dx\over dy}| \]

現在給定隨機變量\(X \sim U(0,1),Y \sim F_Y(y),Y=f(X)\),令:

\[f(x)=F^{-1}(x) \]

\(F\)為任意分布函數，則\(Y\)的概率密度函數為：

\[P_Y(y)=P_X(x)|{dx\over dy}|=1*|{dF(y)\over dy}|=P(y) \]

因此

\[F_Y(y)=F(y) \]