除了PO,还有一种方法介于Wave Solvers和Geometrical Optics之间,称为Hybrid GO-PO,我个人觉得应该叫做几何光学-物理光学混合方法。统一衍射理论(Uniform Theory of Diffraction, UTD)将衍射效应纳入几何光学来计算高频条件下的电磁波传输。个人理解,UTD通过计算绕射系数来补偿几何光学射线边界条件的不足,也就是说几何光学的射线也可以转弯了。这种操作在雷达探测天线设计领域非常实用。除了UTD,Hybrid GO-PO还涉及一种叫射线发射和反弹方法(Shooting and Bouncing Rays, SBR)的技术。这种技术模拟射线在物体表面多次反射,同样基于几何光学。
在进一步讨论利用边界元方法(Boundary Element Method, BEM)和自适应积分法(Adaptive Integral Method, AIM)加速的全波参考模拟器之前,首先简单回顾下前文的内容。前文介绍了PO方法和SDTE/RCSD理论。这些方法用于不同的散射计算需求,但它们的基本理论和适用范围不同。本文将讨论一种通过BEM和AIM结合提供高精度的面散射模拟的方法。
等号左边描述电场和磁场在空间中的“旋转”程度。 $M$ 和 $J$ 是表面电流密度(假想的电流),分别表示磁流和电流的密度(electric and magnetic current densities)。这个公式可以理解为,电场在边界附近“旋转”时,产生磁流和磁场的变化;磁场的旋转也会产生电场和电流的变化。
where $\mathbf{f}{m x}, \mathbf{f}{m y}, \mathbf{f}_{m z}$ are the $x, y, z$ components of the vector basis function $\mathbf{f}_m$ . Similarly, we have:
$$ \begin{aligned} & A_{\mathrm{EM}, i}^{m n}=A_{\mathrm{HJ}, i}^{m n}=\int_{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m z}(\mathbf{r}) g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n y}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m y}(\mathbf{r}) g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n z}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m x}(\mathbf{r}) g_{3, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n z}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m z}(\mathbf{r}) g_{3, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n x}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m y}(\mathbf{r}) g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n x}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m x}(\mathbf{r}) g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}_{n y}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \end{aligned} \tag{S.19} $$
Lastly, we also have:
$$ \begin{aligned} A_{\mathrm{HM}, i}^{m n} & =-j \omega \varepsilon_i \int_{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m x}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n x}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_i \int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m y}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n y}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_i \int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m z}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n z}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{\omega \mu_i} \int{\mathbf{f}m} \int{\mathbf{f}n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) g{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \end{aligned} \tag{S.20} $$
Kai Yang and Ali E Yilmaz. 2011. Comparison of precorrected FFT/adaptive integral method matching schemes. Microwave and Optical Technology Letters 53, 6 (2011), 1368–1372.
The required spatial frequencies are not on the FFT grid but can be interpolated; we add zero padding prior to the FFT step, to ensure enough resolution in the frequency domain for the trilinear interpolation to be sufficiently accurate.
$$ a_{st} = e^{j k \mathbf{p}{st} \cdot \omega_i} \tag{39} $$
当目标入射方向 $-\omega_i$ 与基本方向 $\mathbf{u}$ 之间的夹角接近各高斯光束的发散角(divergence of the small beams)时,aliasing artifacts begin to appear. An example is shown in Fig. 8 (d).
In our framework, we decide on a primary waist w and choose a collection of basic incident directions. In general, 较小的腰宽意味着更大的发散角,这样可以从每个基本方向派生出更多的入射方向,从而减少所需的基本方向数量。较大的腰宽会降低每个高斯光束的发散角,使得组合后的总场发散更小,从而产生更精确的入射方向。
where $\mathbf{n}$ is the surface normal at the macro scale ( $+\mathbf{z}$ ). Note that Eq. 42 and Eq. 43 can also be applied in single simulations, where $\Phi_i$ is computed from a single Gaussian beam.
毛发渲染一直以来主要基于光线追踪技术,无法处理波动光学效应(Wave Optics),例如毛发表面强烈的前向散射和微妙的颜色变化。之前的研究【Xia et al. 2020】证明了衍射效应在纤维的颜色和散射方向上起到了关键作用。然而,这项研究并没有考虑如表面粗糙度和纤维表皮层的微观结构(例如倾斜的角质鳞片)。
法线分布函数预计算(precomputed Normal Distribution Functions, NDFs),使用四维高斯分布(4D Gaussians)预计算并存储小型样本的法线分布函数(NDFs)。存储在多尺度 NDF 图中(multi-scale NDF maps),并在渲染时通过简单的纹理采样调用。
现有方法例如递归计算法和转移矩阵法(Transfer Matrix Method, TMM)的局限性(limitations of existing methods)在于计算复杂度随着层数的增加而显著提高。简化方法则会忽略薄膜中的多重反射。
2. 动机
提高多层结构的渲染效率(improving efficiency for multilayer structures)。
应用于复杂材料的物理渲染(physical rendering of complex materials)。
3. 方法
提出一种基于 Huxley 方法的多层干涉模型(multilayer interference model based on Huxley’s approach)。高效计算周期性多层结构中的反射和透射系数。且支持多种材料和吸收效应(support for multiple materials and absorption effects)。
现有的光传输模型通常基于几何光学和辐射度(radiometry)忽略了光的波动特性,无法模拟干涉和衍射等现象。不能准确再现如彩虹色效应、光栅、薄膜干涉、光偏振等波动光学效应,也就是经典辐射度光传输模型的局限性(limitations of classical radiometric light transport models)。
目前的模型只能处理局部处理波动效应(local treatment of wave effects),对于全局场景中光的传输和相干性搞不了。
2. 动机
实现全局光传输中波动效应的物理一致性(achieving global wave-optics consistency in light transport),也就是结合麦克斯韦电磁理论。
结合波动光学与传统几何光学(integrating wave optics with classical geometric optics),处理光的波动效应,又能在短波长极限下与经典几何光学一致。
3. 方法
光的部分相干性建模(modelling partially-coherent light)分为两个部分。两点相干性描述(two-point coherence description)和光源模型(light source model)。与传统的辐射度不同,本文基于光的部分相干性引入了“交叉谱密度函数”(cross-spectral density function),能够捕捉光的干涉特性。自然光源的物理模型,基于量子力学中的自发辐射原理。
光传输方程的推广(generalizing the light transport equation)。利用谱密度传输方程(spectral-density transport equation)计算光在传播过程中的干涉和衍射效应。并且本文证明了该框架在短波长极限下简化为经典的几何光学,因而可以与现有的光传输方法无缝衔接。
Mengqi Xia, Bruce Walter, Christophe Hery, Olivier Maury, Eric Michielssen, and Steve Marschner, “A Practical Wave Optics Reflection Model for Hair and Fur,” ACM Transactions on Graphics (TOG), vol. 42, no. 4, article 39, pp. 1-15, Jul. 2023.
这里作者上来就提到完全发展的散斑 (Fully Developed Speckle)这个概念。以下是个人理解。在光线照射到一个粗糙的表面(比如纤维/毛发表面)时,表面上的每个小区域都会散射光线。由于表面的微小不规则性,散射的光线之间会相互干涉,产生一个复杂的光强分布。这种分布表现为一系列亮点和暗点,我们称之为散斑。表面的微小特征(比如粗糙度)在整个照射区域(相对于光线的波长来说)变得足够不规则,导致散射光线在各个点上的相位和强度是随机的,这就是所谓的完全发展的散斑。
Improved Accuracy and Practicality(提高的精度和实用性)高精度模型(如Marschner模型、Yan[2015]等)需要复杂的数值计算和大量的预计算数据,因此难以在实时应用中实现。而低精度的经验模型则缺乏足够的物理真实性。于是在Yan[2017]中,做了很多改进。尽管模型简化了光散射路径,但通过结合物理现象(如毛小皮的多层反射结构、各向异性散射等),保持了毛发和毛皮的物理精度。通过对光的纵向和方位角散射进行合理简化,该模型在精确性上优于传统的经验模型。此外就是近场与远场的过渡处理。传统模型往往无法平滑地处理近场和远场之间的光学过渡问题。Yan等人引入了一个近场-远场的解析解决方案,在光线靠近毛发纤维时精确模拟反射,同时在远场时快速近似光线的整体反射行为。使得渲染效率可用于实时渲染。
Analytic Near/Far Field Solution(解析的近场/远场求解)近场(光线与单根毛发纤维之间的短距离交互,即散射行为)和远场(光线与大量毛发纤维之间的远距离集体效应)的处理存在巨大差异。为了实现近场和远场的无缝过渡,作者使用了解析积分方法,而非繁琐的数值积分。解析积分能够直接计算反射函数,不需要通过复杂的数值求解或预计算,极大地减少了计算时间。
简单总结一下,Yan[2017]提出的反射模型效果和性能都不错,通过统一皮层和髓质的IOR,使得模型只需要5个lobes就可以表示皮毛的复杂散射,利用张量近似最小化存储开销。在这个模型的基础上,提出远近场的解析积分,将模型拓展至多尺度渲染。想要在实时渲染中实现BCSDF模型是非常简单的,当前已经有不少实现方式了,并且已经应用于影视行业。[The Lion King (HD). 2017 movie] (2019 Oscar Nominee for Best Visual Effects)
早在XIA[2020]就已经提出了利用波动光学准确描述光与纤维的相互作用,用边界元法(Boundary Element Method, BEM)模拟光线在任意截面的纤维散射。并且XIA[2020]指出由于衍射效应,纤维表现出极强的前向散射效应。因此波动光学效应应该让光线专注集中前向散射的方向。还指出了,小的纤维散射效应显著依赖光的波长,导致强烈的波长散射。此外波动场带来的奇异性软化现象也是决定真实焦散效果的关键。为了BEM模拟的计算量可控,纤维的形状理想为具有规则的横截面形状。但是Marschner[2003]指出毛发表面的不规则对毛发外观有重要的影响,在波动光学中这样的效应是否显著仍然是一个需要探讨和解决的问题。
之前的研究已经探索了如何通过蒙特卡洛方法来模拟体积散射中的散斑效应Bar[2019, 2020]。然而,这些模型主要适用于均匀介质(homogeneous media),不适用于毛发纤维等异质结构。Steinberg, Yan [2022]研究了平面粗糙表面的散斑渲染。然而作者指出,纤维表面的散斑效应与平面表面不同,表现出不同的统计特性。
[1] J. T. Kajiya and T. L. Kay, “RENDERING FUR WITIt THREE DIMENSIONAL TEXTURES,” 1989.
[2] S. R. Marschner, H. W. Jensen, and M. Cammarano, “Light Scattering from Human Hair Fibers,” 2003.
[3] A. Zinke and A. Weber, “Light Scattering from Filaments,” IEEE Trans. Visual. Comput. Graphics, vol. 13, no. 2, pp. 342–356, Mar. 2007.
[4] L.-Q. Yan, C.-W. Tseng, H. W. Jensen, and R. Ramamoorthi, “Physically-accurate fur reflectance: modeling, measurement and rendering,” ACM Trans. Graph., vol. 34, no. 6, pp. 1–13, Nov. 2015.
[5] L.-Q. Yan, H. W. Jensen, and R. Ramamoorthi, “An efficient and practical near and far field fur reflectance model,” ACM Trans. Graph., vol. 36, no. 4, pp. 1–13, Aug. 2017.
[6] M. Xia, B. Walter, C. Hery, O. Maury, E. Michielssen, and S. Marschner, “A Practical Wave Optics Reflection Model for Hair and Fur,” ACM Trans. Graph., vol. 42, no. 4, pp. 1–15, Aug. 2023.
[1] L.-Q. Yan, M. Hašan, W. Jakob, J. Lawrence, S. Marschner, and R. Ramamoorthi, “Rendering glints on high-resolution normal-mapped specular surfaces,” ACM Trans. Graph., vol. 33, no. 4, pp. 1–9, Jul. 2014.
[2] L.-Q. Yan, M. Hašan, S. Marschner, and R. Ramamoorthi, “Position-normal distributions for efficient rendering of specular microstructure,” ACM Trans. Graph., vol. 35, no. 4, pp. 1–9, Jul. 2016.
对于散射场,用边界积分公式(Boundary Integral Formulation)将电磁波的散射问题转化为仅在表面边界上求解的积分方程,关键实现方法是边界元方法(BEM)。再用基于三维快速傅里叶变换(3D FFT)的自适应积分方法(Adaptive Integral Method, AIM)加速边界积分的计算过程。
并且利用GPU加速并行处理大规模的表面散射问题。
并且paper中采用了一种组合小规模模拟结果来表征表面双向散射行为。
相关工作
基于波动光学的反射模型
老生常谈的几何光学 vs 波动光学。本文主要对比表面散射模型。几何光学的经典模型包括:Cook-Torrance模型 [Cook and Torrance 1982]、Oren-Nayar模型[Michael 1994]。波动光学这边,主要是用物理光学近似(Physical Optics Approximations)来简化全波方程。也就是黑毛狗中的一阶近似(单次散射)来估计表面反射。经典模型包括Beckmann-Kirchoff理论和Harvey-Shack模型,它们使用标量形式的近似方程来模拟波动光学效应。它们被广泛应用于各种表面类型的反射估计,比如高斯随机表面、周期性表面等。但是这些方法的计算结果常常是空间上的平均结果,没办法进行高分辨率的细节反射。
A Full-Wave Reference Simulator for Computing Surface Reflectance
Petr Beckmann and Andre Spizzichino. 1987. The scattering of electromagnetic waves from rough surfaces. Artech House.
Andrey Krywonos. 2006. Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction. Ph. D. Dissertation. University of Central Florida.
Ling-Qi Yan, Miloš Hašan, Bruce Walter, Steve Marschner, and Ravi Ramamoorthi. 2018. Rendering Specular Microgeometry with Wave Optics. ACM Trans. Graph. 37, 4 (2018).
Ling-Qi Yan, Miloš Hašan, Steve Marschner, and Ravi Ramamoorthi. 2016. Positionnormal distributions for efficient rendering of specular microstructure. ACM Transactions on Graphics (TOG) 35, 4 (2016), 1–9.
R. L. Cook and K. E. Torrance. 1982. A Reflectance Model for Computer Graphics. ACM Trans. Graph. 1, 1 (jan 1982). https://doi.org/10.1145/357290.357293
Phong细分不需要知道相邻的拓扑信息,仅仅用插值计算,比PN triangles等算法效率更高。GAMES101上提到的Loop and Schaefer利用低度数四边形曲面近似Catmull-Clark曲面,这些方法输入的多边形都被一个多项式曲面替代。而本文的Phong细分不需要任何修正额外的几何区域的操作。
structTessellationFactors{floatedge[3] : SV_TessFactor;float inside :SV_InsideTessFactor;};// The patch constant function runs once per triangle, or "patch"// It runs in parallel to the hull functionTessellationFactorsPatchConstantFunction(InputPatch<TessellationControlPoint,3>patch) {UNITY_SETUP_INSTANCE_ID(patch[0]);// Set up instancing// Calculate tessellation factorsTessellationFactorsf;f.edge[0] =_FactorEdge1.x;f.edge[1] =_FactorEdge1.y;f.edge[2] =_FactorEdge1.z;f.inside=_FactorInside;returnf;}
structInterpolators{float3 normalWS :TEXCOORD0;float3 positionWS :TEXCOORD1;float4 positionCS :SV_POSITION;};// Call this macro to interpolate between a triangle patch, passing the field name#defineBARYCENTRIC_INTERPOLATE(fieldName) \patch[0].fieldName*barycentricCoordinates.x+ \patch[1].fieldName*barycentricCoordinates.y+ \patch[2].fieldName*barycentricCoordinates.z// The domain function runs once per vertex in the final, tessellated mesh// Use it to reposition vertices and prepare for the fragment stage[domain("tri")] // Signal we're inputting trianglesInterpolatorsDomain(TessellationFactorsfactors,// The output of the patch constant functionOutputPatch<TessellationControlPoint,3>patch,// The Input trianglefloat3barycentricCoordinates : SV_DomainLocation) {// The barycentric coordinates of the vertex on the triangleInterpolatorsoutput;// Setup instancing and stereo support (for VR)UNITY_SETUP_INSTANCE_ID(patch[0]);UNITY_TRANSFER_INSTANCE_ID(patch[0],output);UNITY_INITIALIZE_VERTEX_OUTPUT_STEREO(output);float3positionWS=BARYCENTRIC_INTERPOLATE(positionWS);float3normalWS=BARYCENTRIC_INTERPOLATE(normalWS);output.positionCS=TransformWorldToHClip(positionWS);output.normalWS=normalWS;output.positionWS=positionWS;returnoutput;}
// The patch constant function runs once per triangle, or "patch"// It runs in parallel to the hull functionTessellationFactorsPatchConstantFunction(InputPatch<TessellationControlPoint,3>patch) {UNITY_SETUP_INSTANCE_ID(patch[0]);// Set up instancing// Calculate tessellation factorsTessellationFactorsf;f.edge[0] =_FactorEdge1.x;f.edge[1] =_FactorEdge2;// -- Edited --f.edge[2] =_FactorEdge3;// -- Edited --f.inside=_FactorInside;returnf;}_FactorEdge2("Edge 2 factor",Float) =1// -- Edited --_FactorEdge3("Edge 3 factor",Float) =1// -- Edited --
// Returns true if it should be clipped due to frustum or winding cullingboolShouldClipPatch(float4p0PositionCS,float4p1PositionCS,float4p2PositionCS) {returnfalse;}
// Returns true if the point is outside the bounds set by lower and higherboolIsOutOfBounds(float3p,float3lower,float3higher) {returnp.x<lower.x||p.x>higher.x||p.y<lower.y||p.y>higher.y||p.z<lower.z||p.z>higher.z;}// Returns true if the given vertex is outside the camera fustum and should be culledboolIsPointOutOfFrustum(float4positionCS) {float3culling=positionCS.xyz;floatw=positionCS.w;// UNITY_RAW_FAR_CLIP_VALUE is either 0 or 1, depending on graphics API// Most use 0, however OpenGL uses 1float3lowerBounds=float3(-w,-w,-w*UNITY_RAW_FAR_CLIP_VALUE);float3higherBounds=float3(w,w,w);returnIsOutOfBounds(culling,lowerBounds,higherBounds);}
// Returns true if it should be clipped due to frustum or winding cullingboolShouldClipPatch(float4p0PositionCS,float4p1PositionCS,float4p2PositionCS) {boolallOutside=IsPointOutOfFrustum(p0PositionCS) &&IsPointOutOfFrustum(p1PositionCS) &&IsPointOutOfFrustum(p2PositionCS);// -- Edited -- returnallOutside;// -- Edited -- }
// Returns true if the points in this triangle are wound counter-clockwiseboolShouldBackFaceCull(float4p0PositionCS,float4p1PositionCS,float4p2PositionCS) {float3point0=p0PositionCS.xyz/p0PositionCS.w;float3point1=p1PositionCS.xyz/p1PositionCS.w;float3point2=p2PositionCS.xyz/p2PositionCS.w;float3normal=cross(point1-point0,point2-point0);returndot(normal,float3(0,0,1)) <0;}
上面的代码还存在一个跨平台问题。观察方向在不同API的朝向是不同的,因此修改一下代码。
// In clip space, the view direction is float3(0, 0, 1), so we can just test the z coord#ifUNITY_REVERSED_Zreturncross(point1-point0,point2-point0).z<0;#else// In OpenGL, the test is reversedreturncross(point1-point0,point2-point0).z>0;#endif
最后的最后,在 ShouldClipPatch 中添加刚写好的函数用于判断背面剔除。
// Returns true if it should be clipped due to frustum or winding cullingboolShouldClipPatch(float4p0PositionCS,float4p1PositionCS,float4p2PositionCS) {boolallOutside=IsPointOutOfFrustum(p0PositionCS) &&IsPointOutOfFrustum(p1PositionCS) &&IsPointOutOfFrustum(p2PositionCS);returnallOutside||ShouldBackFaceCull(p0PositionCS,p1PositionCS,p2PositionCS);// -- Edited -- }
然后在 PatchConstantFunction 中将需要剔除的Patch的顶点因子设置为0 。
...if (ShouldClipPatch(patch[0].positionCS,patch[1].positionCS,patch[2].positionCS)) {f.edge[0] =f.edge[1] =f.edge[2] =f.inside=0;// Cull the patch}...
// Returns true if the given vertex is outside the camera fustum and should be culledboolIsPointOutOfFrustum(float4positionCS,floattolerance) {float3culling=positionCS.xyz;floatw=positionCS.w;// UNITY_RAW_FAR_CLIP_VALUE is either 0 or 1, depending on graphics API// Most use 0, however OpenGL uses 1float3lowerBounds=float3(-w-tolerance,-w-tolerance,-w*UNITY_RAW_FAR_CLIP_VALUE-tolerance);float3higherBounds=float3(w+tolerance,w+tolerance,w+tolerance);returnIsOutOfBounds(culling,lowerBounds,higherBounds);}
// Returns true if the points in this triangle are wound counter-clockwiseboolShouldBackFaceCull(float4p0PositionCS,float4p1PositionCS,float4p2PositionCS,floattolerance) {float3point0=p0PositionCS.xyz/p0PositionCS.w;float3point1=p1PositionCS.xyz/p1PositionCS.w;float3point2=p2PositionCS.xyz/p2PositionCS.w;// In clip space, the view direction is float3(0, 0, 1), so we can just test the z coord#ifUNITY_REVERSED_Zreturncross(point1-point0,point2-point0).z<-tolerance;#else// In OpenGL, the test is reversedreturncross(point1-point0,point2-point0).z>tolerance;#endif}
可以在材质面板中暴露一个Range。
// .shaderProperties{_tolerance("_tolerance",Range(-0.002,0.001)) =0...}// .hlslfloat_tolerance;...// Returns true if it should be clipped due to frustum or winding cullingboolShouldClipPatch(float4p0PositionCS,float4p1PositionCS,float4p2PositionCS) {boolallOutside=IsPointOutOfFrustum(p0PositionCS,_tolerance) &&IsPointOutOfFrustum(p1PositionCS,_tolerance) &&IsPointOutOfFrustum(p2PositionCS,_tolerance);// -- Edited -- returnallOutside||ShouldBackFaceCull(p0PositionCS,p1PositionCS,p2PositionCS,_tolerance);// -- Edited -- }
// Calculate the tessellation factor for an edgefloatEdgeTessellationFactor(floatscale,floatbias,float3p0PositionWS,float3p1PositionWS) {floatfactor=distance(p0PositionWS,p1PositionWS) /scale;returnmax(1,factor+bias);}
f.inside= ( // If the compiler doesn't play nice...EdgeTessellationFactor(_TessellationFactor,_TessellationBias,patch[1].positionWS,patch[2].positionWS) +EdgeTessellationFactor(_TessellationFactor,_TessellationBias,patch[2].positionWS,patch[0].positionWS) +EdgeTessellationFactor(_TessellationFactor,_TessellationBias,patch[0].positionWS,patch[1].positionWS) ) /3.0;
// Calculate Phong projection offsetfloat3PhongProjectedPosition(float3flatPositionWS,float3cornerPositionWS,float3normalWS) {returnflatPositionWS-dot(flatPositionWS-cornerPositionWS,normalWS) *normalWS;}// Apply Phong smoothingfloat3CalculatePhongPosition(float3bary,float3p0PositionWS,float3p0NormalWS,float3p1PositionWS,float3p1NormalWS,float3p2PositionWS,float3p2NormalWS) {float3smoothedPositionWS=bary.x*PhongProjectedPosition(flatPositionWS,p0PositionWS,p0NormalWS) +bary.y*PhongProjectedPosition(flatPositionWS,p1PositionWS,p1NormalWS) +bary.z*PhongProjectedPosition(flatPositionWS,p2PositionWS,p2NormalWS);returnsmoothedPositionWS;}// The domain function runs once per vertex in the final, tessellated mesh// Use it to reposition vertices and prepare for the fragment stage[domain("tri")] // Signal we're inputting trianglesInterpolatorsDomain(TessellationFactorsfactors,// The output of the patch constant functionOutputPatch<TessellationControlPoint,3>patch,// The Input trianglefloat3barycentricCoordinates : SV_DomainLocation) {// The barycentric coordinates of the vertex on the triangleInterpolatorsoutput;...float3positionWS=CalculatePhongPosition(barycentricCoordinates,patch[0].positionWS,patch[0].normalWS,patch[1].positionWS,patch[1].normalWS,patch[2].positionWS,patch[2].normalWS);float3normalWS=BARYCENTRIC_INTERPOLATE(normalWS);float3tangentWS=BARYCENTRIC_INTERPOLATE(tangentWS.xyz);...output.positionCS=TransformWorldToHClip(positionWS);output.normalWS=normalWS;output.positionWS=positionWS;output.tangentWS=float4(tangentWS,patch[0].tangentWS.w);...}
structTessellationFactors{floatedge[3] : SV_TessFactor;float inside :SV_InsideTessFactor;float3bezierPoints[7] : BEZIERPOS;};//Bezier control point calculationsfloat3CalculateBezierControlPoint(float3p0PositionWS,float3aNormalWS,float3p1PositionWS,float3bNormalWS) {floatw=dot(p1PositionWS-p0PositionWS,aNormalWS);return (p0PositionWS*2+p1PositionWS-w*aNormalWS) /3.0;}voidCalculateBezierControlPoints(inoutfloat3bezierPoints[7],float3p0PositionWS,float3p0NormalWS,float3p1PositionWS,float3p1NormalWS,float3p2PositionWS,float3p2NormalWS) {bezierPoints[0] =CalculateBezierControlPoint(p0PositionWS,p0NormalWS,p1PositionWS,p1NormalWS);bezierPoints[1] =CalculateBezierControlPoint(p1PositionWS,p1NormalWS,p0PositionWS,p0NormalWS);bezierPoints[2] =CalculateBezierControlPoint(p1PositionWS,p1NormalWS,p2PositionWS,p2NormalWS);bezierPoints[3] =CalculateBezierControlPoint(p2PositionWS,p2NormalWS,p1PositionWS,p1NormalWS);bezierPoints[4] =CalculateBezierControlPoint(p2PositionWS,p2NormalWS,p0PositionWS,p0NormalWS);bezierPoints[5] =CalculateBezierControlPoint(p0PositionWS,p0NormalWS,p2PositionWS,p2NormalWS);float3avgBezier=0; [unroll] for (inti=0;i<6;i++) {avgBezier+=bezierPoints[i];}avgBezier/=6.0;float3avgControl= (p0PositionWS+p1PositionWS+p2PositionWS) /3.0;bezierPoints[6] =avgBezier+ (avgBezier-avgControl) /2.0;}// The patch constant function runs once per triangle, or "patch"// It runs in parallel to the hull functionTessellationFactorsPatchConstantFunction(InputPatch<TessellationControlPoint,3>patch) {...TessellationFactorsf= (TessellationFactors)0;// Check if this patch should be culled (it is out of view)if (ShouldClipPatch(...)) {...}else{...CalculateBezierControlPoints(f.bezierPoints,patch[0].positionWS,patch[0].normalWS,patch[1].positionWS,patch[1].normalWS,patch[2].positionWS,patch[2].normalWS);}returnf;}
$$ \begin{aligned} & b: \quad R^2 \mapsto R^3, \quad \text { for } w=1-u-v, \quad u, v, w \geq 0 \ & b(u, v)= \sum_{i+j+k=3} b_{i j k} \frac{3!}{i!j!k!} u^i v^j w^k \ &= b_{300} w^3+b_{030} u^3+b_{003} v^3 \ &+b_{210} 3 w^2 u+b_{120} 3 w u^2+b_{201} 3 w^2 v \ &+b_{021} 3 u^2 v+b_{102} 3 w v^2+b_{012} 3 u v^2 \ &+b_{111} 6 w u v . \end{aligned} $$
// Barycentric interpolation as a functionfloat3BarycentricInterpolate(float3bary,float3a,float3b,float3c) {returnbary.x*a+bary.y*b+bary.z*c;}float3CalculateBezierPosition(float3bary,floatsmoothing,float3bezierPoints[7],float3p0PositionWS,float3p1PositionWS,float3p2PositionWS) {float3flatPositionWS=BarycentricInterpolate(bary,p0PositionWS,p1PositionWS,p2PositionWS);float3smoothedPositionWS=p0PositionWS* (bary.x*bary.x*bary.x) +p1PositionWS* (bary.y*bary.y*bary.y) +p2PositionWS* (bary.z*bary.z*bary.z) +bezierPoints[0] * (3*bary.x*bary.x*bary.y) +bezierPoints[1] * (3*bary.y*bary.y*bary.x) +bezierPoints[2] * (3*bary.y*bary.y*bary.z) +bezierPoints[3] * (3*bary.z*bary.z*bary.y) +bezierPoints[4] * (3*bary.z*bary.z*bary.x) +bezierPoints[5] * (3*bary.x*bary.x*bary.z) +bezierPoints[6] * (6*bary.x*bary.y*bary.z);returnlerp(flatPositionWS,smoothedPositionWS,smoothing);}// The domain function runs once per vertex in the final, tessellated mesh// Use it to reposition vertices and prepare for the fragment stage[domain("tri")] // Signal we're inputting trianglesInterpolatorsDomain(TessellationFactorsfactors,// The output of the patch constant functionOutputPatch<TessellationControlPoint,3>patch,// The Input trianglefloat3barycentricCoordinates : SV_DomainLocation) {// The barycentric coordinates of the vertex on the triangleInterpolatorsoutput;...// Calculate tessellation smoothing multiplerfloatsmoothing=_TessellationSmoothing;#ifdef_TESSELLATION_SMOOTHING_VCOLORSsmoothing*=BARYCENTRIC_INTERPOLATE(color.r);// Multiply by the vertex's red channel#endiffloat3positionWS=CalculateBezierPosition(barycentricCoordinates,smoothing,factors.bezierPoints,patch[0].positionWS,patch[1].positionWS,patch[2].positionWS);float3normalWS=BARYCENTRIC_INTERPOLATE(normalWS);float3tangentWS=BARYCENTRIC_INTERPOLATE(tangentWS.xyz);...}
structTessellationFactors{floatedge[3] : SV_TessFactor;float inside :SV_InsideTessFactor;float3bezierPoints[10] : BEZIERPOS;};float3CalculateBezierControlNormal(float3p0PositionWS,float3aNormalWS,float3p1PositionWS,float3bNormalWS) {float3d=p1PositionWS-p0PositionWS;floatv=2*dot(d,aNormalWS+bNormalWS) /dot(d,d);returnnormalize(aNormalWS+bNormalWS-v*d);}voidCalculateBezierNormalPoints(inoutfloat3bezierPoints[10],float3p0PositionWS,float3p0NormalWS,float3p1PositionWS,float3p1NormalWS,float3p2PositionWS,float3p2NormalWS) {bezierPoints[7] =CalculateBezierControlNormal(p0PositionWS,p0NormalWS,p1PositionWS,p1NormalWS);bezierPoints[8] =CalculateBezierControlNormal(p1PositionWS,p1NormalWS,p2PositionWS,p2NormalWS);bezierPoints[9] =CalculateBezierControlNormal(p2PositionWS,p2NormalWS,p0PositionWS,p0NormalWS);}// The patch constant function runs once per triangle, or "patch"// It runs in parallel to the hull functionTessellationFactorsPatchConstantFunction(InputPatch<TessellationControlPoint,3>patch) {...TessellationFactorsf= (TessellationFactors)0;// Check if this patch should be culled (it is out of view)if (ShouldClipPatch(...)) {..}else{...CalculateBezierControlPoints(f.bezierPoints,patch[0].positionWS,patch[0].normalWS,patch[1].positionWS,patch[1].normalWS,patch[2].positionWS,patch[2].normalWS);CalculateBezierNormalPoints(f.bezierPoints,patch[0].positionWS,patch[0].normalWS,patch[1].positionWS,patch[1].normalWS,patch[2].positionWS,patch[2].normalWS);}returnf;}
并且需要注意,所有插值得到的法线向量都需要标准化。
float3CalculateBezierNormal(float3bary,float3bezierPoints[10],float3p0NormalWS,float3p1NormalWS,float3p2NormalWS) {returnp0NormalWS* (bary.x*bary.x) +p1NormalWS* (bary.y*bary.y) +p2NormalWS* (bary.z*bary.z) +bezierPoints[7] * (2*bary.x*bary.y) +bezierPoints[8] * (2*bary.y*bary.z) +bezierPoints[9] * (2*bary.z*bary.x);}float3CalculateBezierNormalWithSmoothFactor(float3bary,floatsmoothing,float3bezierPoints[10],float3p0NormalWS,float3p1NormalWS,float3p2NormalWS) {float3flatNormalWS=BarycentricInterpolate(bary,p0NormalWS,p1NormalWS,p2NormalWS);float3smoothedNormalWS=CalculateBezierNormal(bary,bezierPoints,p0NormalWS,p1NormalWS,p2NormalWS);returnnormalize(lerp(flatNormalWS,smoothedNormalWS,smoothing));}// The domain function runs once per vertex in the final, tessellated mesh// Use it to reposition vertices and prepare for the fragment stage[domain("tri")] // Signal we're inputting trianglesInterpolatorsDomain(TessellationFactorsfactors,// The output of the patch constant functionOutputPatch<TessellationControlPoint,3>patch,// The Input trianglefloat3barycentricCoordinates : SV_DomainLocation) {// The barycentric coordinates of the vertex on the triangleInterpolatorsoutput;...// Calculate tessellation smoothing multiplerfloatsmoothing=_TessellationSmoothing;float3positionWS=CalculateBezierPosition(barycentricCoordinates,smoothing,factors.bezierPoints,patch[0].positionWS,patch[1].positionWS,patch[2].positionWS);float3normalWS=CalculateBezierNormalWithSmoothFactor(barycentricCoordinates,smoothing,factors.bezierPoints,patch[0].normalWS,patch[1].normalWS,patch[2].normalWS);float3tangentWS=BARYCENTRIC_INTERPOLATE(tangentWS.xyz);...}
voidCalculateBezierNormalAndTangent(float3bary,floatsmoothing,float3bezierPoints[10],float3p0NormalWS,float3p0TangentWS,float3p1NormalWS,float3p1TangentWS,float3p2NormalWS,float3p2TangentWS,outfloat3normalWS,outfloat3tangentWS) {float3flatNormalWS=BarycentricInterpolate(bary,p0NormalWS,p1NormalWS,p2NormalWS);float3smoothedNormalWS=CalculateBezierNormal(bary,bezierPoints,p0NormalWS,p1NormalWS,p2NormalWS);normalWS=normalize(lerp(flatNormalWS,smoothedNormalWS,smoothing));float3flatTangentWS=BarycentricInterpolate(bary,p0TangentWS,p1TangentWS,p2TangentWS);float3flatBitangentWS=cross(flatNormalWS,flatTangentWS);tangentWS=normalize(cross(flatBitangentWS,normalWS));}[domain("tri")] // Signal we're inputting trianglesInterpolatorsDomain(TessellationFactorsfactors,// The output of the patch constant functionOutputPatch<TessellationControlPoint,3>patch,// The Input trianglefloat3barycentricCoordinates : SV_DomainLocation) {// The barycentric coordinates of the vertex on the triangle...float3normalWS,tangentWS;CalculateBezierNormalAndTangent(barycentricCoordinates,smoothing,factors.bezierPoints,patch[0].normalWS,patch[0].tangentWS.xyz,patch[1].normalWS,patch[1].tangentWS.xyz,patch[2].normalWS,patch[2].tangentWS.xyz,normalWS,tangentWS);...}