duadua blog
Toggle navigation
duadua blog
Home
Study
About Me
Archives
Tags
LTC(Linearly Transformed Cosines)
2024-05-04 15:52:17
80
0
0
duadua
--- # abstract > - Linear Transformed Spherical Distribution(LTSD, 线性变换球面分布)。 > - Linear Transformed Cosines(LTC,线性余弦变换)。 > - M矩阵(逆矩阵/归一化/bake lut)。**描述bsdf形状**。 > - roughness / anisotropy / skewness。 > - Nelder–Mead method(downhill simplex method),下降单纯形法。 > - Shading -- 多边形边积分(非面积积分)。 > - Stoke’s theorem。 > - 水平面裁剪。 > - if-else check。 > - proxy sphere to lut by vector form factor(vector irradiance)。 > - Fresnel修正(bake lut2)。**描述bsdf大小**。 > - like split sum dfg lut。 --- # LTSD 对于一个球面分布,可以使用一个3x3的矩阵对它的方向向量(即该分布的自变量参数)进行线性变换,从而将其转为另一个球面分布。 ![title](/api/file/getImage?fileId=66358e5ec6875d1c89a1d829) 设 > - $\omega_{0}$。原分布的方向向量(输入)。 > - $M$。线性变换矩阵。 > - $\omega$。变换后的方向向量。 则有 $$ \omega = \frac{M \omega_{0}}{|| M \omega_{0} || } \\ \omega_0 = \frac{M^{-1} \omega}{|| M^{-1} \omega || } $$ 又设 > - $D_0$。原分布。 > - $D$。新分布。 则又有 $$ \begin{aligned} D(\omega) &= D_0(\omega_0) \frac{\partial \omega_0}{\partial \omega} \\ &= D_0(\frac{M^{-1} \omega}{|| M^{-1} \omega ||}) \frac{M^{-1}}{|| M^{-1} \omega || ^{3}} \end{aligned} \tag{1.1} $$ 其中$\frac{\partial \omega_0}{\partial \omega} = \frac{M^{-1}}{|| M^{-1} \omega || ^{3}}$ 为M变换的雅可比行列式(Jacobian)。 ## 性质 > 1. 旋转or缩放不改变形状。M为旋转or缩放变换时,雅可比项为1,分布的形状不改变。 > 2. 积分形式不变性。 $$ \int_{\Omega} D(\omega) \mathrm{d}\omega = \int_{\Omega} D_0(\omega_0) \frac{\partial \omega_0}{\partial \omega} \mathrm{d}\omega = \int_{\Omega} D(\omega_0) \mathrm{d}\omega_0 $$ > 3. 多边形积分结果不变。 在分布$D$上对一个多边形P积分,和在分布$D_0$上对多边形$P_0$积分的结果是一样的,如果$P_0$上的每一个顶点都是由M矩阵的逆矩阵变换而来的话,即$P_0 = M^{-1}P$。 > 4. 重要性采样。 如果$D_0$可以使用IS,则$D$也可以。 而且因为我们知道了Jacobian,所以$D$的pdf也可以求出。 --- # LTC ## 分布 如果$D_0$的分布为`clamped cosine`(切线),即 $$ D_0(\omega_0 = (x, y, z)) = \frac{1}{\pi} max(0, \vec{n} \cdot \vec{\omega_0}) $$ 又切线空间,$n = (0, 0, 1)$,所以有 $$ D_0(\omega_0 = (x, y, z)) = \frac{1}{\pi} max(0, z) $$ 而$D$则为基于各种bxdf所组成的分布,用$\rho$表示,有 $$ D \approx \rho(\omega_v, \omega_l) \cos\theta_l $$ 如果我们能找到一个变换矩阵$M$,它可以将分布$D_0$变换到分布$D$上,那么**在bxdf分布上对多边形积分**就可以转换为**在$D_0$分布上对多边形积分**的问题。 --- ## 变换矩阵的形式 ![title](/api/file/getImage?fileId=66359ddec6875d1c89a1d82a) $$ M = \begin{bmatrix} m00 & 0 & m02 \\ 0 & m11 & 0 \\ m20 & 0 & m22 \end{bmatrix} $$ 由于其缩放不改变形状,所以我们可以将主对角线上的任意一个数值归一化为1而不影响形状,只有两个值便可表示。这里将右下角的值归一化为一。 最终版本的M矩阵会选择将中间值归一化为一。这样会比右下角归一化出来的5个值更加平缓,数值精度更高。 > - m00和m11模拟材质的粗糙度,其非一致缩放模拟其各向异性。 > - m20和m02为错切变换,模拟偏斜度(skewness)。 --- 对于右下角为一的情况。 $$ M = \begin{bmatrix} m00 & 0 & m02 \\ 0 & m11 & 0 \\ m20 & 0 & m22 \end{bmatrix} \overset{normalize}{=\!=\!=\!=\!=} \begin{bmatrix} a & 0 & b \\ 0 & c & 0 \\ d & 0 & 1 \end{bmatrix} $$ $$ \det(M) = ac - bcd, \\ M^{-1} = \frac{ \begin{bmatrix} c & 0 & -bc \\ 0 & a-bd & 0 \\ -cd & 0 & ac \end{bmatrix} }{ ad - bc } $$ ![title](/api/file/getImage?fileId=663ee394c6875d1c89a1d8c9) ``` /* * a 0 b * 0 c 0 * d 0 1 * det = ac - bcd * adj = c 0 -bc * 0 a-bd 0 * -cd 0 ac */ float I_ltc_line_width_factor(float3 ortho, float3x3 m_inv) { // transpose(inverse(M)) = (1.0 / determinant(M)) * cofactor(M). // Take into account that m12 = m21 = m23 = m32 = 0 and m33 = 1. float a = m_inv._11; float b = m_inv._31; float c = m_inv._22; float d = m_inv._13; float det = a*c - b*c*d; float3x3 cof = { c, 0.0, -b*d, 0.0, a - b*d, 0.0, -c*d, 0.0, a*c }; // 1.0 / length(mul(V, (1.0 / s * M))) = abs(s) / length(mul(V, M)). return abs(det) / length(mul(ortho, cof)); } ``` --- 对于中间为一的情况。 $$ M = \begin{bmatrix} m00 & 0 & m02 \\ 0 & m11 & 0 \\ m20 & 0 & m22 \end{bmatrix} \overset{normalize}{=\!=\!=\!=\!=} \begin{bmatrix} a & 0 & b \\ 0 & 1 & 0 \\ c & 0 & d \end{bmatrix} $$ $$ \det(M) = ad - bc, \\ M^{-1} = \frac{ \begin{bmatrix} d & 0 & -b \\ 0 & ad-bc & 0 \\ -c & 0 & a \end{bmatrix} }{ ad - bc } $$ ``` /* * a 0 b * 0 1 0 * c 0 d * det = ad - bc * adj = d 0 -b * 0 ad-bc 0 * -c 0 a */ float I_ltc_line_width_factor_mid(float3 ortho, float3x3 m_inv) { // transpose(inverse(M)) = (1.0 / determinant(M)) * cofactor(M). // Take into account that m12 = m21 = m23 = m32 = 0 and m33 = 1. float a = m_inv._11; float b = m_inv._31; float c = m_inv._13; float d = m_inv._33; float det = a*d - b*c; float3x3 cof = { d, 0.0, -b, 0.0, a*d - b*c, 0.0, -c, 0.0, a }; // 1.0 / length(mul(V, (1.0 / s * M))) = abs(s) / length(mul(V, M)). return abs(det) / length(mul(ortho, cof)); } ``` ![title](/api/file/getImage?fileId=663ee36bc6875d1c89a1d8c8) > - 注意,此处的M为从LUT里读取的最终的矩阵,一般为$M_Inv$,而其逆矩阵(即原始的M矩阵)通常在计算线光源时会用到。 --- ## 变换矩阵的用法 我们需要根据不同粗糙度$\alpha$和观察方向$\theta_v$分别计算出一个变换矩阵$M^{-1}$,并将其存储在一张2D纹理中使用,并以$(\alpha, \theta_v)$作为uv采样。 > - 事实上,在真正实现的时候,并不需要真地计算角度之后作为采样uv,而是直接将其cos值作为uv更高效。 ![title](/api/file/getImage?fileId=6635c766c6875d1c89a1d830) 我们需要将当前多边形的顶点$P$通过$M^{-1}$逆变换为$P_0$,便可在简单的$D_0$分布下对多边形$P_0$进行积分,等价于在bxdf分布(即$D$分布)下对$P$积分。 ## 变换矩阵的求法 该变换矩阵使用下降单纯形法(downhill simplex method, Nelder–Mead method)拟合求得。 该方法使用MIS计算拟合的结果和原始结果的误差,并不断缩小这个误差。 https://blog.csdn.net/myf_666/article/details/129072785 ``` // compute the error between the BRDF and the LTC // using Multiple Importance Sampling float computeError(const LTC& ltc, const Brdf& brdf, const vec3& V, const float alpha) { double error = 0.0; for (int j = 0; j < Nsample; ++j) for (int i = 0; i < Nsample; ++i) { const float U1 = (i + 0.5f)/Nsample; const float U2 = (j + 0.5f)/Nsample; // importance sample LTC { // sample const vec3 L = ltc.sample(U1, U2); float pdf_brdf; float eval_brdf = brdf.eval(V, L, alpha, pdf_brdf); float eval_ltc = ltc.eval(L); float pdf_ltc = eval_ltc/ltc.magnitude; // error with MIS weight double error_ = fabsf(eval_brdf - eval_ltc); error_ = error_*error_*error_; error += error_/(pdf_ltc + pdf_brdf); } // importance sample BRDF { // sample const vec3 L = brdf.sample(V, alpha, U1, U2); float pdf_brdf; float eval_brdf = brdf.eval(V, L, alpha, pdf_brdf); float eval_ltc = ltc.eval(L); float pdf_ltc = eval_ltc/ltc.magnitude; // error with MIS weight double error_ = fabsf(eval_brdf - eval_ltc); error_ = error_*error_*error_; error += error_/(pdf_ltc + pdf_brdf); } } return (float)error / (float)(Nsample*Nsample); } ``` ![title](/api/file/getImage?fileId=6635e746c6875d1c89a1d836) --- # Integrate 在使用$M^{-1}$将多边形顶点转到cosine分布上之后,便可以在$D_0$上对该多边形进行积分。 但是这里的积分并不是对面积积分,而是对**边**进行积分。 > 这个积分的等价转换有点类似于**斯托克斯定理(Stoke's theorem)**。斯托克斯定理讲的是:在向量场中,场中一个闭合曲线上的线积分等于其所围面积中散度(Curl)的通量(Flux)之和。 $$ \int_{\Omega} \approx \frac{1}{\pi} \cdot \frac{1}{2} \sum_{i=1}^{n}\arccos(\vec{v_i} \cdot \vec{v_j}) \frac{\vec{v_i} \times \vec{v_j}}{|| \vec{v_i} \times \vec{v_j} || } \cdot \vec{n} $$ 其中i为当前边的起点,j为当前边的终点。如果多边形的边集是以极角排序的,那么将有$j = i \mod n + 1$。 其中$\arccos(\vec{v_i} \cdot \vec{v_j})$为当前边所对应的角度(单位半径下即为弧长)。 ![title](/api/file/getImage?fileId=6635bcd5c6875d1c89a1d82b) 其中$\frac{\vec{v_i} \times \vec{v_j}}{|| \vec{v_i} \times \vec{v_j} || }$为当前边与圆心连线组成的扇形面的法线(与该扇形面垂直)。 ![title](/api/file/getImage?fileId=6635bd1ac6875d1c89a1d82c) 其中$\vec{n}$则为当前着色点的法线,与扇形面的法线的点乘则表示其在当前着色点所在平面的投影,便是最终的radiance。 ![title](/api/file/getImage?fileId=6635c0f6c6875d1c89a1d82d) 除此之外,在实现的时候,计算$arccos$会有精度丢失从而导致比较明显的瑕疵。论文中的解决方法是通过一个三次函数来近似,当然也有许多其它近似方法。 ![origin](/api/file/getImage?fileId=6635c6d1c6875d1c89a1d82f) ![approximation](/api/file/getImage?fileId=6635c618c6875d1c89a1d82e) > - 此处将$\theta$和$\sin(\theta)$一起近似。 > - 这里的$\sin(\theta)$其实是$cross$的长度,即$|| \vec{v_i} \times \vec{v_j} ||$。 --- # clip to horizon(upper hemisphere) ## if-else check 一种方案是每条边与水平面求交,使用交点将水平面以下的边裁掉。这将会产生大量的if-else语句。 ![title](/api/file/getImage?fileId=6635c7f9c6875d1c89a1d831) ## proxy sphere to lut by vector form factor(vector irradiance)。 另一种方案是使用代理体球(proxy sphere)。 该方案重新考察边积分函数,删除与表面法线的点积,从而不投影到平面上,即 $$ \vec{F} = \arccos(\vec{v_i} \cdot \vec{v_j}) \frac{\vec{v_i} \times \vec{v_j}}{|| \vec{v_i} \times \vec{v_j} || } $$ 由此可得一个向量,称为向量辐射度(vector irradiance)$\vec{F}$。 $\vec{F}$的模长但是光源在$\vec{F}$方向上的irradiance大小。 ![title](/api/file/getImage?fileId=6635d874c6875d1c89a1d832) 由于我们可以假设irradiance大小为$||\vec{F}||$的光源来自一个球体。且可以通过$\vec{F}$获得一个该球体对于当前cosine分布的张角(angular extent)和倾斜度(elevation angle)。 $$ \left\{ \begin{array} \ angular\_extent = \arcsin\sqrt{|| \vec{F} || } \\ elvation\_angle = \vec{n} \cdot \frac{\vec{F}}{|| \vec{F} ||} \end{array} \right. $$ ![title](/api/file/getImage?fileId=6635de24c6875d1c89a1d833) ![title](/api/file/getImage?fileId=6635de29c6875d1c89a1d834) 对于不同的张角$(\sigma \in (0, \frac{\pi}{2}))$,和不同的倾斜度$(\omega)$。我们可以离线求出裁剪后的radiance结果大小,从而获得裁剪后与未裁剪的比例,并将其存储到一张2D Lut贴图上,并在运行时根据张角和倾斜度计算uv进行采样即可。 其中计算裁剪后的radiance结果可参考下图。 ![title](/api/file/getImage?fileId=6635e1c7c6875d1c89a1d835) 该方法根据倾斜度和张角确定的范围与半球的关系分为4类。 > - 最下方边缘处的倾斜角也小于90度,说明没有在水平面以下的情况。 > - 即$\omega + \sigma < \frac{\pi}{2}$。即$\omega \in [0, \frac{\pi}{2} - \sigma]$。 > - 最下方边缘处的倾斜角大于90度,但中间的倾斜角(即$\omega$)小于90度。 > - 此时有一半以上在水平面以上,但有一部分在水平面以下。 > - 即$\omega + \sigma > \frac{\pi}{2}$且$\omega < \frac{\pi}{2}$。 > - 即$\omega \in [\frac{\pi}{2} - \sigma, \frac{\pi}{2}]$。 > - 中间倾斜角大于90度,但最上方边缘处的倾斜角小于90度。 > - 此时有小于一半在水平面以上,有一半以上在水平面以下。 > - 即$\omega > \frac{\pi}{2}$且$\omega - \sigma < \frac{\pi}{2}$。 > - 即$\omega \in [\frac{\pi}{2},\frac{\pi}{2} + \sigma]$。 > - 最上方边缘处的倾斜角也大于90度,中间倾斜角小于180度。 > - 此时所有都在水平面以下。 > - 即$\omega - \sigma > \frac{\pi}{2}$。 > - 即$\omega \in [\frac{\pi}{2} + \sigma, \pi]$。 ``` float ihemi(float w, float s) { float g = asinf(cosf(s)/sinf(w)); float sinsSq = sqr(sinf(s)); if (w >= 0.0f && w <= (pi/2.0f - s)) return pi*cosf(w)*sinsSq; if (w >= (pi/2.0f - s) && w < pi/2.0f) return pi*cosf(w)*sinsSq + G(w, s, g) - H(w, s, g); if (w >= pi/2.0f && w < (pi/2.0f + s)) return G(w, s, g) + H(w, s, g); return 0.0f; } ``` --- # fresnel approximation 由于没有考虑能量的遮蔽和吸收,上面的LTC积分的结果始终为1。 也就是说,LTC积分只考虑了bxdf的形状,但是并没有考虑其大小,我们需要考虑Fresnel项求其大小。 这个过程跟IBL里求DFG那张LUT有些类似。最终的$n_D$和$f_D$项也序列化到一张2D Lut上,根据roughness和NoV采样即可。 ![title](/api/file/getImage?fileId=6635e859c6875d1c89a1d837) > - 这里也可以存储$n_D-f_D$和$f_D$,从而减少一次乘法。 > - 这张LUT和水平裁剪的那张LUT可以合并在一起。 --- # shading with textured polygonal lights todo... --- # ref > - Real-Time Area Lighting https://blog.selfshadow.com/publications/s2016-advances/s2016_ltc_rnd.pdf > - https://zhuanlan.zhihu.com/p/345364404 > - https://zhuanlan.zhihu.com/p/426217725 > - https://zhuanlan.zhihu.com/p/350694154 > - Projective Texture Mapping https://developer.download.nvidia.cn/assets/gamedev/docs/projective_texture_mapping.pdf > - learn opengl area light https://learnopengl-cn.github.io/08%20Guest%20Articles/2022/03%20Area%20Lights/#_9 > - self shadow ltc github code https://github.com/selfshadow/ltc_code > - LTC Fresnel Approximation https://blog.selfshadow.com/publications/s2016-advances/s2016_ltc_fresnel.pdf > - sphere / disk light by ltc https://blog.selfshadow.com/publications/s2017-shading-course/heitz/s2017_pbs_ltc_lines_disks.pdf [sig17][arealight][disk light]Real-Time Line- and Disk-Light Shading with Linearly Transformed Cosines.pdf
Pre: No Post
Next:
Test
0
likes
80
Weibo
Wechat
Tencent Weibo
QQ Zone
RenRen
Submit
Sign in
to leave a comment.
No Leanote account?
Sign up now.
0
comments
More...
Table of content
No Leanote account? Sign up now.