Spherical-Harmonic

最新推荐文章于 2025-01-18 22:31:22 发布

「已注销」

最新推荐文章于 2025-01-18 22:31:22 发布

阅读量623

点赞数

分类专栏： Unity 文章标签： 3d渲染图形学 shader unity

本文链接：https://blog.youkuaiyun.com/qq_33967521/article/details/103656859

版权

Unity 专栏收录该内容

33 篇文章

订阅专栏

本文深入探讨了基于物理的渲染方程，特别是球谐函数在全局照明中的作用。通过球谐函数，可以近似还原复杂的球面光照方程。在Unity中，3阶球谐被用于实现光照探针，简化了光照计算并提供了基本的全局光照效果。文章还介绍了球谐函数的数学原理，并展示了其在实时渲染中的应用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

图片来自 Snaps Art HD | School
读本文前请确认掌握了足够的基础
至少达到了入门水平，或者像我一样踏进来半只脚，我是菜逼顺带一提 =_=
上篇(本文)主要讲解基于与计算辐射度传输技术的全局照明效果，也就是SH
这篇文章很蹩脚，但比上篇好一点，给美术写教程让我感觉我是傻逼
中、下篇有时间会写，写什么再说，但早晚会出
感谢熊大lao

基于物理的渲染方程

祖传的渲染方程：
$L_0 = L_e + \int_\Omega f_r \cdot L_i \cdot (w_i \cdot n) \cdot dw_i$
$L_0$ 、 $L_e$ 、 $f_r$ 、 $L_i$ 分别是p点的出射光亮度、发出的光亮度、入射方向到出射方向的反射比例、入射光亮度；
$(w_i \cdot n)$ 是入射角带了的入射光衰减；
$\int_\Omega ...dw_i$ 是入射方向半球的积分；
现在，以辐射度的角度重新看一下这个方程
$L(x,w_o) = L_e(x,w_o) + \int_s f_r(x,w_i \rarr w_o) \cdot L(x \rq,w_i) \cdot G(x,x \rq) \cdot V(x,x \rq) \cdot dw_i$
$L(x,w_o)$ 表示在位置 $x$ 上，方向为 $w_o$ 的光线的辐射度；
$L_e(x,w_o)$ 表示位置 $x$ 上，自发射防线为 $w_o$ 的光线的辐射度；
$f_r(x,w_r \rarr w_o)$ 表示位置 $x$ 上， $w_i$ （入射方向）到 $w_o$ （出射方向）的反射比例；
$\rq,w_i)$ 表示其他物体位置 $\rq$ 处，沿 $w_i$ 方向照射过来的光线辐射度；
$\rq)$ 表示位置 $x$ 与 $\rq$ 之间的几何关系；
$\rq)$ 表示位置 $x$ 与 $x\rq$ 之间是否存在遮挡。
那么，只要能解出入射辐射度值的立体角微分在半球球面上的积分的结果，就可以得到出射辐射度值。
但是，实时算球面每个点的出射辐射度现在的硬件水平几乎不可能。

球谐函数

因为光具有波粒二象性的特点，可以将光视为一种电磁波。
对于各向同性的在均匀介质中传递的电磁波，可以麦克斯韦方程组中推导出波速、电场强度、磁场强度三者建立联系的波动方程。
$\begin{cases} v = \cfrac{1}{\sqrt {u \varepsilon}} \\ \nabla ^2 E - \cfrac{1}{v^2} \cfrac{\partial ^2 E}{\partial t^2} = 0 \\ \nabla ^2 H - \cfrac{1}{v^2} \cfrac{\partial ^2 H}{\partial t^2} = 0 \\ \end{cases}$
因为是近似，再加上磁场的作用比电场弱的多，可以进一步砍成：
$\begin{cases} v = \cfrac{1}{\sqrt {u \varepsilon}} \\ \nabla ^2 E - \cfrac{1}{v^2} \cfrac{\partial ^2 E}{\partial t^2} = 0 \\ \end{cases}$
现在，电场强度矢量 $E$ 用符号 $f$ 表示；因为咱们是在球面坐标系算，所以先将拉普拉斯方程转换到球面坐标。
$\cfrac{1}{r^2} \cfrac{\partial }{\partial r}(r^2 \cfrac{\partial f}{\partial r}) + \cfrac{1}{r^2 sin \theta} \cfrac{\partial }{\partial \theta}(sin \theta \cfrac{\partial f}{\partial \theta}) + \cfrac{1}{r^2 sin ^2 \theta} \cfrac{\partial ^ 2 f}{\partial \varphi ^ 2} = 0$
利用分离变量法可以分解成多个常微分方程，即使 $f(r,\theta,\varphi) = R(r)\Phi(\theta)\Phi(\varphi)$
两边同时乘 $\cfrac{r^2}{YR}$ ，可以得到
$\cfrac{1}{R} \cfrac{d }{d r}(r^2 \cfrac{d R}{d r}) = - \cfrac{1}{Y sin \theta} \cfrac{\partial }{\partial \theta}(sin \theta \cfrac{\partial Y}{\partial \theta}) - \cfrac{1}{Y sin ^2 \theta} \cfrac{\partial ^ 2 Y}{\partial \varphi ^ 2} = l(l+1)$
我们利用常数项将左右分离为两个式子，再对同时包含 $\theta$ 与 $\varphi$ 的一边进一步分离
$\cfrac{sin \theta}{\Theta} \cfrac{d}{d \theta}(sin \theta \cfrac{d \Theta}{d \theta}) + l(l+1)sin ^2 \theta = - \cfrac{1}{\Theta} \cfrac{d^2 \Theta}{d \varphi ^2} = m^2$
最后我们可以拆成
$\begin{cases} \cfrac{d}{dr} (r^2 \cfrac{dR}{dr}) - l(l + 1)R = 0 \\ \cfrac{d ^2 \Phi}{d \varphi ^2} + m^2 \Phi = 0\\ sin \theta \cfrac{d}{d \theta}(sin \theta \cfrac{d \Theta}{d \theta}) + [l(l + 1)sin^2 \theta - m^2] \Theta = 0 \\ \end{cases}$
前两个可以说是白给，可以很容易算出解为：
$Cr^l + \cfrac{D}{r ^ {l + 1}}$ , $l = 0, 1, 2, . . ., C 、 D 为常数$
$\Phi (\varphi)= e^{im\varphi}$
第三个比较麻烦，需要利用到勒让德，因为我觉得大家也应该没学过吧，反正我以前都没听说过。。关于勒让德可以在下面的文章中了解

https://ccjou.wordpress.com/2012/02/20/legendre-%E5%A4%9A%E9%A0%85%E5%BC%8F/
https://en.wikipedia.org/wiki/Associated_Legendre_polynomials
http://silviojemma.com/public/papers/lighting/spherical-harmonic-lighting.pdf

首先令 $\theta$ , $\Theta(\theta)$ ，则有
$\cfrac{d \Theta(\theta)}{d \theta} = \cfrac{d y}{dx}\cfrac{dx}{d \theta} = -sin \theta \cfrac{dy}{dx} = (1-x^2) \cfrac{d}{dx}[(1-x^2) \cfrac{dy}{dx}]$
带入到第三个式子中，可以得到 $m$ 阶连带勒让德方程
$\cfrac {d}{dx}[(1-x^2) \cfrac{dy}{dx}] + [1(l+1) - \cfrac{m^2}{(1-x^2)}]y = 0$
当 $m = 0$ 时，此方程退化为勒让德方程
$\cfrac {d}{dx}[(1-x^2) \cfrac{dy}{dx}] + [1(l+1)]y = 0$
将 $y$ 写为 $P^m_l(x)$ ，上式可解为
$P^0_l(x) = \cfrac{1}{2^ll!} \cfrac{d^l}{dx^l} (x^2 - 1)^l,l \in [0,\infin]$
那么上文中 $m$ 阶连带勒让德方程可以解为
$\begin{cases} \tag{1} P^m_l(x) = \cfrac{(-1)^m}{2^ll!} (1 - x^2)^\cfrac{m}{2} \cfrac{d^{l + m}}{dx^{l + m}} (x^2 - 1)^l, l \in N,l \in [0,\infin] \\ P^{-m}_l(x) = (-1)^m \cfrac{(l - m)!}{(l + m)!} P^m_l(x), m<0,l \in[0,\infin] \end{cases}$
用伴随勒让德多项式定义函数 $\Theta (\theta)$ 为
$\tag{2} \Theta(\theta) = P^m _l(cos \theta), l \in N，l\geq |m|$
所以球谐表达式克表示为
$\tag{3} Y^m_l(\theta, \varphi) = K^m_l \Phi(\varphi) \Theta(\theta) = K^m_l e^{im \varphi} P^m_l (cos \theta), l \in N,m=0, \pm 1,\pm 2, ... , \ \pm l$
$K^m_l$ 为归一化因子，归一化后，关于l和m的球谐函数 $Y^m_l$ 可以表示为：
$Y^m_l(\theta, \varphi) = (-1)^m \sqrt{ \cfrac{(2l + 1)}{4 \pi} \cfrac{(l-|m|)!}{(l+|m|)!}}P^m_l (cos \theta)e^{im \varphi}$

球谐光照

众所周知，一个原始信号波可以由一系列带有缩放因子的简谐波叠加而成，而球谐光照的核心就是球面亮度信号编码和重建。
$\tag{4} O(x) \approx \tilde{O}(x) = \sum ^n_{i=1} c_iB_i(x)$
为此，我们需要基函数和缩放因子。
球谐光照中使用的球谐函数序列是一组正交基函数，而且在应用于渲染时，通常采用实数值的球谐，结合（1）(2)，基于球面坐标系下的实数球谐函数可定义如下：
$\tag{基函数} Y^m_l(\theta,\varphi) = \begin{cases} \sqrt{2} K^m_l P^{-m}_l (cos \theta) sin(-m \varphi) , &当 m < 0时\\ K^m_l P^m_l (cos \theta) , &当 m = 0时\\ \sqrt{2} K^m_l P^m_l (cos \theta) sin(m \varphi) , &当 m > 0时\\ \end{cases}$
有了正交基函数，便可以将原始的函数改用球谐基函数进行近似替代。
假设原始函数为 $g (s)$ , $s$ 是关于球面坐标系的天顶角 $\theta$ 和方位角 $\varphi$ 的函数，作为基函数 $Y^m_l(\theta,\varphi)$ 的可以改成 $Y^m_l(s)$ ，投影后的近似函数为 $\tilde{g}(s)$ 。
因为正交基函数的对偶函数便是它自身，所以加权系数 $c^m_l$ 为
$\tag{加权系数} c^m_l = \int_S g(s) Y^m_l (s)ds$
结合 $(4)$ ，可得到投影后的近似函数 $\tilde{g}(s)$ 为
$\tilde{g}(s) = \sum ^{lmax}_{l = 0} \sum ^{l}_{m=-l} c^m_l Y^m_l(s)$
现在通过 $i = l^2 + l + m$ ，映射 $l, m$
$c_i = \int_S g(s) Y_i (s)ds$
通过概率密度函数和随机样本进一步估算，也就是蒙特卡洛估算法
上式可近似为
$c_i = \cfrac{4 \pi}{N} \sum ^{N}_{j = 0} g(s_j) Y_i (s_j)$

还原球面上的光照方程

$L(x,w_o) = L_e(x,w_o) + \int_s f_r(x,w_i \rarr w_o) \cdot L(x \rq,w_i) \cdot G(x,x \rq) \cdot V(x,x \rq) \cdot dw_i$
1. 对于漫反射表面，各个方向上反射光线相同，所以可以吧 $f_r(x,w_i \rarr w_o)$ ，替换为常数 $k$ .
2. 如果现在不关心AO的话， $\rq)$ ，可替换为1。
3. $\rq,w_i) \cdot G(x,x \rq)$ 则用球谐函数替代。
4. 假设当前物体不存在自发光。
那么现在，我们可以把光照方程写为：
$L(x,w_o) \approx\sum ^{lmax}_{l = 0} \sum ^{l}_{m=-l} \stackrel{L}{C^m_l} \stackrel{G}{C^m_l} Y^m_l$

Unity中的球谐

在Lit中
因为 $i = l^2 + l + m$ ，所以 $l$ 值越大，咱们需要的系数和基函数数量越大。
unity中使用了3阶球谐，即 $l$ 取值为{0,1,2}。
所以球谐函数的系数共有9个，Unity用7个float4存储9*3(RGB)个数字。

下面的代码是简化过的，实际buffer里绝不可能就这几行
```
CBUFFER_START(UnityPerDraw)
float4 unity_SHAr, unity_SHAg, unity_SHAb;
float4 unity_SHBr, unity_SHBg, unity_SHBb;
float4 unity_SHC;
CBUFFER_END
```
假设当前片段受到光照探针的影响
```
CBUFFER_START(UnityProbeVolume)
float4 unity_ProbeVolumeParams;
CBUFFER_END
```
如果unity_ProbeVolumeParams.x == 0，那么当前采用的是常规探针而不是LPPV
现在，我们的光照探针采样函数可以写为
```
float3 SampleLightProbes (LitSurface s) {
	if (unity_ProbeVolumeParams.x) {
		//这种情况，涉及到一些其他因素，在这里先不做解释
		return LPPVSampleValue;
	}
	else {
		float4 coefficients[7];
		coefficients[0] = unity_SHAr;
		coefficients[1] = unity_SHAg;
		coefficients[2] = unity_SHAb;
		coefficients[3] = unity_SHBr;
		coefficients[4] = unity_SHBg;
		coefficients[5] = unity_SHBb;
		coefficients[6] = unity_SHC;
		return max(0.0, SampleSH9(coefficients, s.normal));
	}
}
```
将采样结果*diffuse叠加到color上
```
 color.rgb += SampleLightProbes (LitSurface s) * surface.diffuse;
```

SampleSH9
实际上，就是在做分量和系数相乘之后的叠加

float3 SampleSH9(float4 SHCoefficients[7], float3 N)
{
    float4 shAr = SHCoefficients[0];
    float4 shAg = SHCoefficients[1];
    float4 shAb = SHCoefficients[2];
    float4 shBr = SHCoefficients[3];
    float4 shBg = SHCoefficients[4];
    float4 shBb = SHCoefficients[5];
    float4 shCr = SHCoefficients[6];

    float3 res = SHEvalLinearL0L1(N, shAr, shAg, shAb);
    res += SHEvalLinearL2(N, shBr, shBg, shBb, shCr);

    return res;
}
real3 SHEvalLinearL0L1(real3 N, real4 shAr, real4 shAg, real4 shAb)
{
    real4 vA = real4(N, 1.0);
    real3 x1;
    
    x1.r = dot(shAr, vA);
    x1.g = dot(shAg, vA);
    x1.b = dot(shAb, vA);

    return x1;
}
real3 SHEvalLinearL2(real3 N, real4 shBr, real4 shBg, real4 shBb, real4 shC)
{
    real3 x2;
    real4 vB = N.xyzz * N.yzzx;
    
    x2.r = dot(shBr, vB);
    x2.g = dot(shBg, vB);
    x2.b = dot(shBb, vB);

    real vC = N.x * N.x - N.y * N.y;
    real3 x3 = shC.rgb * vC;

    return x2 + x3;
}

SHEvalLinearL0L1比较简单， $l_0,L_1$ 一共四项，换换系数顺序，rgb一一对应，妙啊。
SHEvalLinearL2中
real4 vB = N.xyzz * N.yzzx;构造出了 $l = 2$ 时，前四项多项式中的xy，yz，z^2，xz；
real vC = N.x * N.x - N.y * N.y; 则构造了 $l = 2, m = 2$ 时球谐系数的 x^2 - y^2，最后一项单独算。

还记得上文中启用LPPV的情况吗
在这种情况下，我们需要预计算计算产生的Texture3D中取得各个分量的球谐系数

float3 SampleProbeVolumeSH4(
TEXTURE3D_PARAM(SHVolumeTexture, SHVolumeSampler), 
float3 positionWS, float3 normalWS,
float4x4 WorldToTexture,float transformToLocal, 
float texelSizeX, float3 probeVolumeMin, 
float3 probeVolumeSizeInv)
{
    float3 position = (transformToLocal == 1.0) ? mul(WorldToTexture, float4(positionWS, 1.0)).xyz : positionWS;
    
    float3 texCoord = (position - probeVolumeMin) * probeVolumeSizeInv.xyz;
    texCoord.x = clamp(texCoord.x * 0.25, 0.5 * texelSizeX, 0.25 - 0.5 * texelSizeX);

    float4 shAr = SAMPLE_TEXTURE3D_LOD(SHVolumeTexture, SHVolumeSampler, texCoord, 0);
    texCoord.x += 0.25;
    float4 shAg = SAMPLE_TEXTURE3D_LOD(SHVolumeTexture, SHVolumeSampler, texCoord, 0);
    texCoord.x += 0.25;
    float4 shAb = SAMPLE_TEXTURE3D_LOD(SHVolumeTexture, SHVolumeSampler, texCoord, 0);
    
    return SHEvalLinearL0L1(normalWS, shAr, shAg, shAb);
}