Real-Time Skin Rendering

Wrap Lighting

最简单的Trick就是这个方法。假设uniform光照,把原本的 N ⋅ L N\cdot L NL变成 N ⋅ L + w r a p 1 + w r a p \frac{N\cdot L + wrap}{1+wrap} 1+wrapNL+wrap,这样原本变成0的cos就会被推迟,于是就在着色阴影边界就会出现一个类似于SSS的东西。

基于BSSRDF的卷积方法

主要是基于【6】【7】的理论来做,而且Skin渲染这个工作有一系列的研究。此文更加侧重于实时方法。

BRDF的可测性

首先定义S:
d L o ( x o , x i , ω o ) [ W ⋅ m − 2 ⋅ s r − 1 ] = S ( x i , ω i , x o , ω o ) d Φ i ( x i , ω i ) dL_o(x_o,x_i,\omega_o)[W\cdot m^{-2}\cdot sr^{-1}]=S(x_i,\omega_i,x_o,\omega_o)d\Phi_i(x_i,\omega_i) dLo(xo,xi,ωo)[Wm2sr1]=S(xi,ωi,xo,ωo)dΦi(xi,ωi)
其中 d Φ i ( x i , ω i ) [ W ] d\Phi_i(x_i,\omega_i)[W] dΦi(xi,ωi)[W]代表单位面积 d A ( x i ) [ m 2 ] dA(x_i)[m^2] dA(xi)[m2]上的入射微分辐射流 d Φ i ( x i , ω i ) = L i ( x i , ω i ) cos ⁡ < ω i , n > d ω i d A ( x i ) [ = d E i ⋅ d A i ] d\Phi_i(x_i,\omega_i)=L_i(x_i,\omega_i)\cos<\omega_i,n>d\omega_idA(x_i)[=dE_i\cdot dA_i] dΦi(xi,ωi)=Li(xi,ωi)cos<ωi,n>dωidA(xi)[=dEidAi].
注意:此处的 d L o dL_o dLo显然是和入射点的位置还有关系的。
于是将入射位置和方向积分,显然就有:
L o ( x o , ω o ) = ∫ A ∫ Ω S ( x i , ω i , x o , ω o ) L i ( x i , ω i ) ( ω i ⋅ n ) d ω i d A ( x i ) L_o(x_o,\omega_o)=\int_A\int_{\Omega}S(x_i,\omega_i,x_o,\omega_o)L_i(x_i,\omega_i)(\omega_i\cdot n)d\omega_idA(x_i) Lo(xo,ωo)=AΩS(xi,ωi,xo,ωo)Li(xi,ωi)(ωin)dωidA(xi)
通过S的定义
d L r = S ⋅ d Φ i = S d E i ⋅ d A i dL_r=S\cdot d\Phi_i = SdE_i\cdot dA_i dLr=SdΦi=SdEidAi
如果假设散射均匀且相对于参考表面具有各向同性,可以两端对入射面积A积分,求得入射位置的均值( d L o ( x o , x i , ω o ) → d L o ( x o , ω o ) dL_o(x_o,x_i,\omega_o) \rightarrow dL_o(x_o,\omega_o) dLo(xo,xi,ωo)dLo(xo,ωo)),自然地得到:
∫ A i d L r ( x o , x i , ω o ) = d L r ( x o , ω o ) = d E i ⋅ ∫ A i S ⋅ d A i \int_{A_i}dL_r(x_o,x_i,\omega_o)=dL_r(x_o,\omega_o)=dE_i\cdot \int_{A_i}S\cdot dA_i AidLr(xo,xi,ωo)=dLr(xo,ωo)=dEiAiSdAi
于是就可以得到 f r = ∫ A i S ⋅ d A i f_r=\int_{A_i}S\cdot dA_i fr=AiSdAi.
于是BRDF的定义 d L r d E i [ s r − 1 ] \frac{dL_r}{dE_i}[sr^{-1}] dEidLr[sr1].

介于BSSRDF和BRDF的定义都来自于无穷小量的推导,辐射度引入无穷小量分析详见【2】【3】,所以他们终究也不能被直接测量到(因为实际测量终究无法得到无穷小量,只能得到极小间隔的平均值,无穷小的立体角不会包含任何辐射流)。更多内容参考下一节后半部分

Reflectance

反射流(flux)和入射流的比值。
复习一下各个radiometric quantities:
Energy:当个光子携带能量 Q = h c λ [ J ] Q=\frac{hc}{\lambda}[J] Q=λhc[J]
面积相关
Flux:单位时间内通过一个任意表面或者区域的能量 Φ = d Q d t [ J / s ( W ) ] \Phi=\frac{dQ}{dt}[J/s(W)] Φ=dtdQ[J/s(W)]
Irradiance&Radiant Exitance:单位面积的Flux, E ( p ) = d Φ ( p ) d A [ W / m 2 ] E(p)=\frac{d\Phi(p)}{dA}[W/m^2] E(p)=dAdΦ(p)[W/m2],来的叫做Irradiance,离开的叫做Radiant Exitance.
立体角相关
Intensity I = d Φ d ω [ W / s r ] I=\frac{d\Phi}{d\omega}[W/sr] I=dωdΦ[W/sr]
面积角全部相关
Radiance:单位垂直面积单位立体角单位时间的能量 L ( p , ω ) = d E ω d ω = d Φ d ω d A p e r p L(p,\omega)=\frac{dE_{\omega}}{d\omega}=\frac{d\Phi}{d\omega dA_{perp}} L(p,ω)=dωdEω=dωdAperpdΦ,注意$ dA_{perp} 表 示 垂 直 于 表示垂直于 \omega$的投影面积

因为BRDF的定义来自于无穷小量定义,真实测量中,无穷小的立体角不会包含任何辐射流,所以没有办法直接测定,所以需要更加一般的反射表达以便测定反射属性,Reflectance就是这样一个量,测试他可以使用任何配置的beam,同时还能被BRDF表达。
注意下面的推导都是在讨论如何测定反射属性!
通过一个非无穷小立体角抵达/离开一个面积元(dA)的辐射流(Flux)可以表达为:
d Φ i = d A ⋅ ∫ ω i L i ( ω i ) ⋅ c o s θ d ω i [ W ] d\Phi_i=dA\cdot \int_{\omega_i}L_i(\omega_i)\cdot cos\theta d\omega_i[W] dΦi=dAωiLi(ωi)cosθdωi[W]
d Φ r = d A ⋅ ∫ ω r L r ( ω r ) ⋅ c o s θ d ω r [ W ] d\Phi_r=dA\cdot \int_{\omega_r}L_r(\omega_r)\cdot cos\theta d\omega_r[W] dΦr=dAωrLr(ωr)cosθdωr[W]
其中带入反射方程可以展开 d Φ r d\Phi_r dΦr:
d Φ r = d A ⋅ ∫ ω r ∫ ω i f r ( ω i , ω r ) ⋅ L i ( ω i ) cos ⁡ θ d ω i cos ⁡ γ d ω r [ W ] d\Phi_r=dA\cdot \int_{\omega_r}\int_{\omega_i}f_r(\omega_i,\omega_r)\cdot L_i(\omega_i)\cos\theta d\omega_i \cos\gamma d\omega_r[W] dΦr=dAωrωifr(ωi,ωr)Li(ωi)cosθdωicosγdωr[W]
于是我们可以得到一般的Reflectance, ρ = d Φ r d Φ i \rho = \frac{d\Phi_r}{d\Phi_i} ρ=dΦidΦr
ρ ( ω i , ω r , L i ) = ∫ ω r ∫ ω i f r ( ω i , ω r ) ⋅ L i ( ω i ) cos ⁡ θ d ω i cos ⁡ γ d ω r ∫ ω i L i ( ω i ) ⋅ c o s θ d ω i [ d i m e n s i o n l e s s ] \rho(\omega_i,\omega_r,L_i)=\frac{\int_{\omega_r}\int_{\omega_i}f_r(\omega_i,\omega_r)\cdot L_i(\omega_i)\cos\theta d\omega_i \cos\gamma d\omega_r}{\int_{\omega_i}L_i(\omega_i)\cdot cos\theta d\omega_i}[dimensionless] ρ(ωi,ωr,Li)=ωiLi(ωi)cosθdωiωrωifr(ωi,ωr)Li(ωi)cosθdωicosγdωr[dimensionless]
此时我们假设 L i L_i Li位置无关,仅仅和方向有关。如果我们进一步假设其与方向也无关,那么他就是一个常量,于是上式化为:
R = ρ ( ω i , ω r ) = ∫ ω r ∫ ω i f r ( ω i , ω r ) ⋅ cos ⁡ θ d ω i cos ⁡ γ d ω r ∫ ω i c o s θ d ω i [ d i m e n s i o n l e s s ] R = \rho(\omega_i,\omega_r)=\frac{\int_{\omega_r}\int_{\omega_i}f_r(\omega_i,\omega_r)\cdot \cos\theta d\omega_i \cos\gamma d\omega_r}{\int_{\omega_i} cos\theta d\omega_i}[dimensionless] R=ρ(ωi,ωr)=ωicosθdωiωrωifr(ωi,ωr)cosθdωicosγdωr[dimensionless]
使用此表达式测定反射属性时,需要保证 L i L_i Li在各个空间和方向上都是相同的。
由此我们可以知道,理想漫反射表面的R如何测定,测定到R之后有: R = f r ∗ π ∗ π π → f r = R π R=\frac{f_r*\pi*\pi}{\pi}\rightarrow f_r =\frac{R}{\pi} R=πfrππfr=πR

关于BSSRDF:
根据【7】, R d R_d Rd是指单位反射面积上的diffuse reflectance,所以在单位反射面积上依然有: R d = π 2 ∗ S / π → S = R d / π R_d=\pi^2*S/\pi\rightarrow S = R_d/\pi Rd=π2S/πS=Rd/π(*)

完整的BSSRDF再考虑一下菲涅尔:
S d ( x i , ω i . x o , ω o ) = 1 π F t ( x i , ω i ) R ( ∣ x i − x o ∣ 2 ) F t ( x o , ω o ) S_d(x_i,\omega_i.x_o,\omega_o)=\frac{1}{\pi}F_t(x_i,\omega_i)R(|x_i-x_o|_2)F_t(x_o,\omega_o) Sd(xi,ωi.xo,ωo)=π1Ft(xi,ωi)R(xixo2)Ft(xo,ωo)
其中 F t F_t Ft代表菲涅尔系数。

使用 R d R_d Rd渲染

根据【6】(【6】中方程8参考【7】中方程4以获得更加详细的推导),有:
d M o ( x o ) = R d ⋅ d Φ i ( x i ) = R d L i ( x i , ω ) cos ⁡ θ ⋅ d ω d A → dM_o(x_o) = R_d\cdot d\Phi_i(x_i)= R_d L_i(x_i,\omega) \cos\theta \cdot d\omega dA\rightarrow dMo(xo)=RddΦi(xi)=RdLi(xi,ω)cosθdωdA

M o ( x o ) = ∫ Ω ∫ A R d L i ( x i , ω ) cos ⁡ θ ⋅ d ω d A M_o(x_o) = \int_{\Omega}\int_A R_d L_i(x_i,\omega) \cos\theta \cdot d\omega dA Mo(xo)=ΩARdLi(xi,ω)cosθdωdA

= ∫ A R d E i d A =\int_A R_dE_i dA =ARdEidA

如果 R d R_d Rd是对称的,那么对任意表面点计算radiant exitance有卷积:
M o ( x , y ) = ∫ R 2 E ( x ′ , y ′ ) R d ( x − x ′ , y − y ′ ) d x ′ d y ′ M_o(x,y) = \int_{R^2}^{} E(x',y')R_d(x-x',y-y')dx'dy' Mo(x,y)=R2E(x,y)Rd(xx,yy)dxdy

如果假设纯漫反射,可以丢掉所有的菲涅尔项,进而可以把radiant exitance转换为radiance:

L o = M o π L_o=\frac{M_o}{\pi} Lo=πMo

上式可以考虑一下 R d R_d Rd与BSSRDF的关系 S = R d / π S = R_d/\pi S=Rd/π

纹理空间

这种方法用到纹理表面,先计算纹理表面的irradiance map,然后直接计算 M o M_o Mo,细节和注意事项(stretch map等)可以参考GPU Gem 3. 这种方法如果用在离线渲染上还是比较精确的,但是用在实时上效率就不好,因为每个觉得都将会要求来计算一次纹理空间的高斯卷积,对于GPU不是特别适应(顶点着色器并不擅长把顶点展开到clip空间,而是赋予每个顶点一个uv值,这些缺陷有后续论文做改进),对于角色比较多的场景也相当慢。

屏幕空间

【8】指出这种使用在屏幕空间的算法只是一种看起来比较好方法,也就是说不是完全的基于物理,他们对比了

  • (a) no simulation of subsurface scattering,
  • (b) a state-of-the-art, texture-space algorithm from nVIDIA,
  • © naive mapping of an existing algorithm into screen space, and
  • (d) screen-space algorithm in 【8】

并且指出这种【8】的SS的方法独立于surface距离和灯光位置,对比于广为接受的纹理空间的方法在感知上更加贴近,同时更快,也更好地能用在多角色的场景下。
因为是在屏幕空间,所以需要通过考虑深度来决定blur的宽度,

  • A pixel representing a further away object should use a narrower kernel.
  • Greater gradients in the depth map should also use narrower kernels. This is similar in spirit to using stretch maps, without the need to actually calculate them.

具体的处理方法很简单就是在kernel的宽度上乘以一个factor:

s x = α d ( x , y ) + β ⋅ a b s ( d d x ( d ( x , y ) ) s_x =\frac{\alpha}{d(x,y)+\beta\cdot abs(ddx(d(x,y))} sx=d(x,y)+βabs(ddx(d(x,y))α
s y = α d ( x , y ) + β ⋅ a b s ( d d y ( d ( x , y ) ) s_y =\frac{\alpha}{d(x,y)+\beta\cdot abs(ddy(d(x,y))} sy=d(x,y)+βabs(ddy(d(x,y))α

这里 α \alpha α控制全局次表面散射等级, β \beta β控制次表面程度与深度值梯度的关系,d是线性深度,ddx/ddy很方便地对应于屏幕空间导数。参考值 α = 11 , β = 800 \alpha = 11,\beta=800 α=11β=800.

缺陷:

  • Adjacent points in 2D screen space may not be adjacent in 3D world space. This produces artifacts in the form of small haloes around some areas of the model such as the ears or nose.
  • It fails to simulate the light transmitted through high-curvature features, since we lack lighting information from the back of the objects.
  • It requires the usage of additional textures in order to store the specular channel and object’s matte.

可分离核方法(此方法可以用在任何SS和TS)

【5】在前面SS的基础上提到了一个kernel可分离的idea.

大致思想就是把原本的高斯拟合的diffusion profile再利用高斯核的可分离性把N个二维高斯卷积简化为2N个单方向的一维卷积:

( E ∗ R d ) ( x , y ) ≈ ( E ∗ A ) ( x , y ) = ∑ i = 1 N ( ( E ∗ a i ) ∗ a i ) ( x , y ) (E*R_d)(x,y)\approx(E*A)(x,y)=\sum_{i=1}^N((E*a_i)*a_i)(x,y) (ERd)(x,y)(EA)(x,y)=i=1N((Eai)ai)(x,y)

w i t h . A ( x , y ) = ∑ i = 1 N a i ( x ) a i ( y ) with. A(x,y)=\sum_{i=1}^Na_i(x)a_i(y) with.A(x,y)=i=1Nai(x)ai(y)

【10】指出这种近似可以使用均值为0的高斯函数来表达:

R d = ( x , y ) ≈ A g ( x , y ) = ∑ i = 1 N w i G ( x , y ; σ i ) R_d=(x,y)\approx A_g(x,y)=\sum_{i=1}^N w_iG(x,y;\sigma_i) Rd=(x,y)Ag(x,y)=i=1NwiG(x,y;σi)

介于高斯核的可分离性,和自然地就可以转换为2N个一维高斯卷积。

高斯核的 separability 属性

G 2 D ( x , y ; σ 1 2 + σ 2 2 ) = G 1 D ( x ; σ 1 2 ) G 1 D ( y ; σ 2 2 ) G_{2D}(x,y;\sigma_1^2+\sigma_2^2)=G_{1D}(x;\sigma_1^2)G_{1D}(y;\sigma_2^2) G2D(x,y;σ12+σ22)=G1D(x;σ12)G1D(y;σ22)
介于高斯核的此属性,我们显然有卷积:
∫ Y ∫ X G 2 D ( x , y ; σ 1 2 + σ 2 2 ) f ( x 0 − x , y 0 − y ) d x d y \int_Y\int_XG_{2D}(x,y;\sigma_1^2+\sigma_2^2)f(x_0-x,y_0-y)dxdy YXG2D(x,y;σ12+σ22)f(x0x,y0y)dxdy
= ∫ Y ∫ X G 1 D ( x ; σ 1 2 ) G 1 D ( y ; σ 2 2 ) f ( x 0 − x , y 0 − y ) d x d y =\int_Y\int_XG_{1D}(x;\sigma_1^2)G_{1D}(y;\sigma_2^2)f(x_0-x,y_0-y)dxdy =YXG1D(x;σ12)G1D(y;σ22)f(x0x,y0y)dxdy
= ∫ Y G 1 D ( y ; σ 2 2 ) ∫ X G 1 D ( x ; σ 1 2 ) f ( x 0 − x , y 0 − y ) d x d y =\int_YG_{1D}(y;\sigma_2^2)\int_XG_{1D}(x;\sigma_1^2)f(x_0-x,y_0-y)dxdy =YG1D(y;σ22)XG1D(x;σ12)f(x0x,y0y)dxdy
如此一个二维卷积可以化为两个分离的一维卷积,实现上,先卷x方向,再对卷好x方向的卷y方向。
对于任何一个n维高斯核,也都可以类似地实现为2N个一维卷积。

但是对于实时来说2N个高斯卷积依然是消耗很大的,但是在离散的像素平面上高斯核就变成了2D矩阵,于是整个问题就变成了如何来分解整个矩阵,使得它能够尽可能地简单又保留更多的diffusion profile的信息。基本的思想是对矩阵做奇异值分解,得到原来kernel对应离散表达矩阵的低秩近似。

Matrix factorization

由于大的奇异值对应着矩阵中的主要信息,因此可以运用奇异值分解进行数据分析,提取矩阵的主要信息。
现在的目标是要最 R d R_d Rd做一个separable kernel近似:
R d ( x , y ) ≈ A ( x , y ) = ∑ i = 1 N a i ( x ) a i ( y ) . . . . . ( 1 ) R_d(x,y)\approx A(x,y)=\sum_{i=1}^Na_i(x)a_i(y).....(1) Rd(x,y)A(x,y)=i=1Nai(x)ai(y).....(1)
考虑求解一个最优化问题:
给定一个值N, min ⁡ A ∣ R d − A ∣ F \min_A{|R_d-A|_F} minARdAF,s.t. r a n k ( A ) ≤ N rank(A)\le N rank(A)N
Eckart–Young–Mirsky theorem【11】指出一个矩阵的SVD分解能够得出F范数最佳的低秩近似。于是:
A s = U Σ N V T A_s=U\Sigma_NV^T As=UΣNVT
注意A的F范数是上述对角矩阵元素向量的2范数。关于分解具体可以参考高等工程数学。
于是 A s A_s As可以写为 A s = ∑ i = 1 N u ( i ) σ i v ( i ) T A_s=\sum_{i=1}^N u^{(i)} \sigma_i v^{(i)T} As=i=1Nu(i)σiv(i)T,其中 σ i \sigma_i σi Σ N \Sigma_N ΣN的对应对角元素, u / v u/v u/v是对应列。注意这个式子的意义,两个向量相乘得到的矩阵为秩一矩阵,如果把 σ \sigma σ从大到小排列,取得的前若干项就能很好的近似原来的矩阵。
由于 R d R_d Rd具有对称性,所以 u ( i ) = v ( i ) u^{(i)}=v^{(i)} u(i)=v(i)我们可以得到最初式(1)要求的 a i = σ i u ( i ) a_i = \sqrt{\sigma_i}u^{(i)} ai=σi u(i).
注意它不球对称,不能量守恒;但是当N大于3的情况下,这些限制将几乎不可感知。参考【4】的图S6.
另外我们也可以强制对其进行归一化 A s ′ = A s ∣ R d ∣ 1 ∣ A s ∣ 1 A_s^{'}=A_s\frac{|R_d|_1}{|A_s|_1} As=AsAs1Rd1,虽然归一化后不再满足2范数最优,但是结果看起来更好。

文章指出仅使用SVD的第1项还是不够好,一般使用前2~6项就比较接近原始kernel了,但是计算时间依旧很耗。最好是使用一个秩1的近似仅仅做两次一维卷积就能得到较好的结果,当然前提就是需要做一些无上大雅的看起来比较正确的假设。

但是但是,这些方法都不是很好,**首先真实的diffusion profile是不可分离的,所以完全重建是不可能的。而且它分离近似也不径向对称,可以参考【4】的图5. **

【5】又提出了一种更好的方法,做一个假设 E ( x , y ) = E 1 ( x ) + E 2 ( y ) E(x,y)=E_1(x)+E_2(y) E(x,y)=E1(x)+E2(y)注意这个假设比较适用于轴对齐的阴影边界,但是对于一般的函数,实验表明依然具有较好的结果,大部分真实世界的E都能被这个假设局部近似【5】),然后带入到计算 M e M_e Me的公式:

M e ( x , y ) = ∬ E ( x ′ , y ′ ) R d ( x − x ′ , y − y ′ ) d x ′ d y ′ M_{e}(x, y)=\iint E\left(x^{\prime}, y^{\prime}\right) R_{d}\left(x-x^{\prime}, y-y^{\prime}\right) d x^{\prime} d y^{\prime} Me(x,y)=E(x,y)Rd(xx,yy)dxdy

= ∬ ( E 1 ( x ′ ) + E 2 ( y ′ ) ) R d ( x − x ′ , y − y ′ ) d x ′ d y ′ =\iint\left(E_{1}\left(x^{\prime}\right)+E_{2}\left(y^{\prime}\right)\right) R_{d}\left(x-x^{\prime}, y-y^{\prime}\right) d x^{\prime} d y^{\prime} =(E1(x)+E2(y))Rd(xx,yy)dxdy

= ∫ E 1 ( x ′ ) ∫ R d ( x − x ′ , y − y ′ ) d y ′ d x ′ + ∫ E 2 ( y ′ ) ∫ R d ( x − x ′ , y − y ′ ) d x ′ d y ′ =\int E_{1}\left(x^{\prime}\right) \int R_{d}\left(x-x^{\prime}, y-y^{\prime}\right) d y^{\prime} d x^{\prime} +\int E_{2}\left(y^{\prime}\right) \int R_{d}\left(x-x^{\prime}, y-y^{\prime}\right) d x^{\prime} d y^{\prime} =E1(x)Rd(xx,yy)dydx+E2(y)Rd(xx,yy)dxdy

令:

a p ( x − x ′ ) = ∫ R d ( x − x ′ , y − y ′ ) d y ′ ; a p ( y − y ′ ) = ∫ R d ( x − x ′ , y − y ′ ) d x ′ a_p(x-x')=\int R_d(x-x',y-y')dy';a_p(y-y')=\int R_d(x-x',y-y')dx' ap(xx)=Rd(xx,yy)dy;ap(yy)=Rd(xx,yy)dx

则:

= ∫ E 1 ( x ′ ) a P ( x − x ′ ) [ 1 ∥ a p ∥ 1 ∫ a p ( y − y ′ ) d y ′ ] d x ′ + ∫ E 2 ( y ′ ) a P ( y − y ′ ) [ 1 ∥ a p ∥ 1 ∫ a p ( x − x ′ ) d x ′ ] d y ′ =\int E_{1}\left(x^{\prime}\right) a_{P}\left(x-x^{\prime}\right) [\frac{1}{\left\|a_{p}\right\|_{1}} \int a_{p}\left(y-y^{\prime}\right) d y^{\prime}] d x^{\prime}+\int E_{2}\left(y^{\prime}\right) a_{P}\left(y-y^{\prime}\right) [\frac{1}{\left\|a_{p}\right\|_{1}} \int a_{p}\left(x-x^{\prime}\right) d x^{\prime}] d y^{\prime} =E1(x)aP(xx)[ap11ap(yy)dy]dx+E2(y)aP(yy)[ap11ap(xx)dx]dy

注意中括号中的部分为1.上式最终等于:

∬ E ( x ′ , y ′ ) 1 ∥ a p ∥ 1 a p ( x − x ′ ) a p ( y − y ′ ) d x ′ d y ′ \iint E\left(x^{\prime}, y^{\prime}\right) \frac{1}{\left\|a_{p}\right\|_{1}} a_{p}\left(x-x^{\prime}\right) a_{p}\left(y-y^{\prime}\right) d x^{\prime} d y^{\prime} E(x,y)ap11ap(xx)ap(yy)dxdy

此时, a p a_p ap就代表沿着对应的轴预积分的 R d R_d Rd.因为 R d R_d Rd的径向对称属性,所以 a p ( x ) = a p ( y ) a_p(x) = a_p(y) ap(x)=ap(y).另外 ∣ ∣ a p ∣ ∣ 1 = ∫ a p ( y − y ′ ) d y ′ = ∫ ∫ R d ( x − x ′ , y − y ′ ) d x ′ d y ′ = ∣ ∣ R d ∣ ∣ 1 ||a_p||_1 = \int a_p(y-y')dy'=\int \int R_d(x-x',y-y')dx' dy' = ||R_d||_1 ap1=ap(yy)dy=Rd(xx,yy)dxdy=Rd1

为了可以艺术调参,对于式子:

A ( x , y ) = ∑ i = 1 N a i ( x ) a i ( y ) A(x,y)=\sum_{i=1}^Na_i(x)a_i(y) A(x,y)=i=1Nai(x)ai(y)

a i a_i ai可以使用更加一般的函数,这样每个波长都能有2N中调节方法。

也可以使用作者提出的一个分离散射模型,分别在1维度的情况下控制近远场的1D高斯函数 G ( x , σ ) G(x,\sigma) G(x,σ),于是有了N=1时候的新可分离kernel:

a m ( x ) = w G ( x , σ n ) + ( 1 − w ) G ( x , σ f ) A m ( x , y ) = a m ( x ) a m ( y ) \begin{aligned} a_{m}(x)=& w G\left(x, \sigma_{n}\right)+(1-w) G\left(x, \sigma_{f}\right) \\ & A_{m}(x, y)=a_{m}(x) a_{m}(y) \end{aligned} am(x)=wG(x,σn)+(1w)G(x,σf)Am(x,y)=am(x)am(y)

渲染与实现

高频照明走样,因为近似的kernel本身不是径向对称的,这个通过随机旋转kernel卷积的local坐标轴来克服,好处在于使用原来的预积分kernel就好,不需要再单独地处理。同时只对那些比较接近摄像机的像素进行旋转。距离小于30%的kernel宽度就可以触发随机旋转。

另外,再计算 R d ( d ) R_d(d) Rd(d)的d,这个d是在世界坐标下的,所以要把像素距离转换为对应的世界距离,可以根据深度和SS坐标位置来重建世界空间距离。对于SS上的卷积方法,也可以使用【8】中提出的临时factor来乘以kernel的宽度,对于上面指出的双一维高斯kernel可以使用作者在补充材料【4】中提出的方法,这种方法把d展开为 d = d x y 2 + d z 2 d=\sqrt{d_{xy}^2+d_z^2} d=dxy2+dz2 ,然后展开高斯,把关于 d z d_z dz的部分单独提取出来,如此就能对SS距离和深度单独处理,深度就可以做预计算。

皮肤的实时Translucency渲染【12】

离线方法暂时就不说了,可以直接参考【10】。

1、对于比较薄的物体,假设光线出射位置的法线和入射位置的法线相反(入射位置的法线有了)
2、观察者不知道背面的精确E分布(我们可以根据这一边的E推测假设出来)
3、dense的材质光线入射后再出射已经diffuse化,没有明显的高频现象(高频现象被移除,当使用法线贴图上的高频法线的时候也没有高频现象,甚至可以使用低频法线(顶点法线)来计算,同时假设那个地方周围都是常数E)
4、人类皮肤的albedo在表面上没有剧烈变化(假设入射和出射点的颜色都一样)

先计算入射面的E,再利用homogeneous材质中的transmittance计算公式进行透射。代码如下(论文拷贝):

float distance(float3 pos, float3 N, int i)
{
    //使用最普通的shadow map
    //把世界坐标转换到灯光view下
    //-0.005*N是防止shadow map走样
    float4 shrinkedpos = float4(pos - 0.005 ∗ N, 1.0);
    float4 shwpos = mul(shrinkedpos, lights[i].viewproj);
    float d1 = shwmaps[i].Sample(sampler, shwpos.xy/shwpos.w);
    float d2 = shwpos.z;
    return abs(d1 - d2);
}
// This function can be precomputed for efficiency
float3 T(float s) 
{
    return float3(0.233, 0.455, 0.649)exp(-s*s/0.0064) +
    float3(0.1, 0.336, 0.344)exp(-s*s/0.0484) +
    float3(0.118, 0.198, 0.0)exp(-s*s/0.187) +
    float3(0.113, 0.007, 0.007)exp(-s*s/0.567) +
    float3(0.358, 0.004, 0.0)exp(-s*s/1.99) +
    float3(0.078, 0.0, 0.0)exp(-s*s/7.41);
}
float s = scale ∗ distance(pos, Nvertex, i);
//+0.3防止 Nvertex*L=0 造成不连续
float E = max(0.3 + dot(-Nvertex, L), 0.0);
float3 transmittance = T(s) ∗ lights[i].color ∗
attenuation ∗ spot ∗ albedo.rgb ∗ E;
// We add the contribution of this light
M += transmittance + reflectance;

另外一种预积分方法

这种方法做起来非常简单,先预计算,然后渲染的时候查找LUT,比较麻烦的是需要烘焙一个曲率纹理,不过这个都是可以trick的,另外如果进一步降低精度,此模型可以进一步简化到不需要曲率纹理。

首先一块平整的光滑的皮肤在uniform的灯光(四面八方过来的光都是均匀的)照射下,各个地方都是一样的亮,只有三个情况可以引起亮度的变化,表面曲率/法线变化/阴影
假设uniform灯光,那么渲染方程的结果主要就是被 N ⋅ L N\cdot L NL决定的(Wrap Lighting就是如此思想),但是这种调节不是基于物理测量的。

从完整的次表面渲染方程开始:

L o ( x o , ω o ) = ∫ A ∫ Ω S ( x i , ω i , x o , ω o ) L i ( x i , ω i ) ( ω i ⋅ n ) d ω i d A ( x i ) L_o(x_o,\omega_o)=\int_A\int_{\Omega}S(x_i,\omega_i,x_o,\omega_o)L_i(x_i,\omega_i)(\omega_i \cdot n)d\omega_idA(x_i) Lo(xo,ωo)=AΩS(xi,ωi,xo,ωo)Li(xi,ωi)(ωin)dωidA(xi)

我们假设光照强度 L i L_i Li均匀分布在上半球:

L o ( x o , ω o ) = L i ∫ A ∫ Ω S ( x i , x o ) ( ω i ⋅ n ) d ω i d A ( x i ) L_o(x_o,\omega_o)=L_i\int_A\int_{\Omega}S(x_i,x_o)(\omega_i \cdot n)d\omega_idA(x_i) Lo(xo,ωo)=LiAΩS(xi,xo)(ωin)dωidA(xi)

L o ( x o , ω o ) = L i ∫ A S ( x i , x o ) ∫ Ω ( ω i ⋅ n ) d ω i ⋅ d A ( x i ) L_o(x_o,\omega_o)=L_i\int_AS(x_i,x_o)\int_{\Omega}(\omega_i \cdot n)d\omega_i \cdot dA(x_i) Lo(xo,ωo)=LiAS(xi,xo)Ω(ωin)dωidA(xi)

由于点光源加上狄拉克函数直接变为:

L o ( x o , ω o ) = L i ∫ A S ( x i , x o ) ( ω i ⋅ n ) ⋅ d A ( x i ) L_o(x_o,\omega_o)=L_i\int_AS(x_i,x_o)(\omega_i \cdot n)\cdot dA(x_i) Lo(xo,ωo)=LiAS(xi,xo)(ωin)dA(xi)

对于曲率因素以及比较光滑的部分,我们考虑是在不同半径的球上做积分,也就是上述方程对A的积分是在球面上,另外,而且球形积分为了简单起见被简化为在圆环上积分。固定一个 L i L_i Li,对于每一根法线,都要在环上做积分,由于各种简化,这个模型本身具有对称性,所以把法线作为参数不如直接把 cos ⁡ θ = L i ⋅ N \cos\theta = L_i\cdot N cosθ=LiN作为参数,此时在圆环上做积分变成了(从 L L L方向开始旋转):

L o ( x o , ω o ) → L i ∫ − π π C ⋅ R d ( x i , x o ) cos ⁡ ( θ + ϕ ) ⋅ d ϕ L_o(x_o,\omega_o)\rightarrow L_i\int_{-\pi}^{\pi}C\cdot R_d(x_i,x_o)\cos(\theta + \phi) \cdot d\phi Lo(xo,ωo)LiππCRd(xi,xo)cos(θ+ϕ)dϕ

注意上式的不完整性,因为我们把原本定义在面积上的 R d R_d Rd(它本身具有能量守恒的属性)简化到了环上来积分,找个时候就不一定能量守恒了,所以我们需要归一化这个东西:

L o ( x o , ω o ) ≈ L i ∫ − π π R d ( ∣ x i − x o ∣ ) ∫ − π π R d ( ∣ x i − x o ∣ ) cos ⁡ ( θ + ϕ ) ⋅ d ϕ L_o(x_o,\omega_o)\approx L_i\int_{-\pi}^{\pi}\frac{ R_d(|x_i-x_o|)}{\int_{-\pi}^{\pi} R_d(|x_i-x_o|)}\cos(\theta + \phi) \cdot d\phi Lo(xo,ωo)LiππππRd(xixo)Rd(xixo)cos(θ+ϕ)dϕ

∣ x i − x o ∣ |x_i-x_o| xixo按照 R d R_d Rd的意义就是直线距离 ∣ x i − x o ∣ = 2 r ⋅ sin ⁡ ( ϕ 2 ) |x_i-x_o|=2r\cdot\sin(\frac{\phi}{2}) xixo=2rsin(2ϕ)

另外对于每一个曲率计算一个圆,也就是半径的倒数 1 r \frac{1}{r} r1.

最后做一张2dLUT,两个参数 1 r \frac{1}{r} r1 cos ⁡ θ \cos\theta cosθ.

缺点就是,每个点都当成了球来算,但实际上一个着色点周围的那些点并不在同一个球上,所以曲率肯定是不一样的,有一定的差异。另外这个模型在平滑的地方很给力,但是如上所述因为是基于球的,遇到曲率变化过快的地方就不行了。

所以那些细节变化大的地方就需单独考虑,通常这些变化是由法线贴图带来的。

主要是改法线,有了法线贴图我们可以随意访问一个法线周围的所有法线,比如:采样很多法线/直接执行散射光照计算(也就是每个点算漫反射然后用skin profile加权求和)。但这不是很效率,不能预计算。所以我们要预过滤法线贴图,这样我们在渲染的时候就只需要采样一个法线采样计算了。但光照过程对于法线本身不是线性过程,所以直接过滤法线是有问题的。

但其实对于漫反射来说(SSS在前面已经被预积分为了一个近似的漫反射BRDF), 1 n ∑ i = 1 n ( f r L ⋅ N i ) = f r L i ⋅ ( 1 n ∑ i = 1 n N i ) \frac{1}{n}\sum_{i=1}^n(f_r L \cdot N_i)=f_r L_i \cdot (\frac{1}{n}\sum_{i=1}^n N_i) n1i=1n(frLNi)=frLi(n1i=1nNi),也就是说我们在漫反射BRDF下可以直接过滤法线。但是对于 L ⋅ N L\cdot N LN是负数的情况下,可能会出现一些问题,但是都是非常接近的。于是我们使用skin profile来预过滤法线贴图,这时候我们法线包含高光,RGB(每个通道profile都不同)一共需要四个法线贴图,也就是要采样四次,太cost了,于是作者发现这样一情况(通过一个捕捉真实物体法线的reseach作者发现的)

figure

从上到下越来越糊,于是我们就可以只过滤而且采样第一个和最后一个(通过一个可以调节miplevel的sampler),然后通过blend得到中间两个法线。

关于阴影边界的散射效果,这个效果还是蛮重要的,这里就是要利用shadow map,对于过滤过的shadow map,将边界看成一个过度函数,根据阴影的值反算出位置,然后再来和skin profile做过滤,过程和平滑表面的类似,不过个人感觉这么搞有点复杂,也没有太仔细的看,还要看看现在游戏中流行的做法是怎么搞的。

来龙去脉

关于BRDF和测量,【1】【2】都是重要参考。
关于DA,可以参考:https://dannyteng.bitcron.com/post/ti-chuan-shu
想要透彻理解皮肤渲染这一块,【9】,【6】【7】都是相当重要的工作
Advanced techniques for realistic real-time skin rendering(section 7). In GPU Gems 3, 以及Efficient Rendering of Human Skin使用高斯卷积计算
论文 Screen-space perceptual rendering of human skin初步把高斯卷积用到SS.
这几篇文章有助于解释为甚可以使用卷积方程来进行渲染,以及为何在SS方程不变
GPU Pro 2 预积分皮肤渲染是另外一种廉价的SSS渲染方法

实现

To be done…

Ref

【1】Geometrical Considerations and Nomenclature for Reflectance
【2】Part I. Concepts, Chapters 1-3, of [28], Nat. Bur. Stand. (U.S.), NBS Tech. Note 910-1, 93 pages (March 1976).
【3】 Jones, R. Clark, Appendix [on Radiometry] in: I. F. Spiro, R. Clark Jones, and D. Q. Wark, Atmospheric Transmission: Concepts, Symbols, Units and Nomenclature, Infrared Physics, 5, No. 1, March 1965, pp. 11- 36. (With the help of Dr. Jones, a number of unfortunate typographical errors were corrected for its reproduction in C37l.l
【4】Separable Subsurface Scattering Supplementary Materials
【5】Separable Subsurface Scattering
【6】A Rapid Hierarchical Rendering Technique for Translucent Materials
【7】A Practical Model for Subsurface Light Transport
【8】Screen-Space Perceptual Rendering of Human Skin
【9】Donner’s PhD thesis 涵盖了分层SSS材质的一系列方法,包括比较精确的三层皮肤模型等等
【10】Efficient rendering of human skin
【11】https://en.wikipedia.org/wiki/Low-rank_approximation#Basic_low-rank_approximation_problem
【12】Real-Time Realistic Skin Translucency

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值