Bruneton 预计算大气渲染 Unity版代码阅读与原理理解#1
[{“source”:{“position”:0,“lines”:[“
开始之前
本人对各种渲染技术很感兴趣，但由于现实原因不能深入探索。对Unity中预计算大气渲染技术研究是出于爱好，在课余业余时间自学的，仅仅是阅读到刚刚能看懂个大概的水平。我觉得从理论/物理，到数学公式，再到代码的过程十分神奇，我的学习过程就是通过代码的实现，反推其中的几何/物理原理。尽管我是自学，并没有深厚的学术背景，但通过反向工程和参考各种资源，我逐渐拼凑出了大气散射的复杂机制及其实现方式。
参考的Unity项目仓库：https://github.com/Scrawk/Brunetons-Improved-Atmospheric-Scattering
参考的散射原理（后期看到多次散射才找到ebruneton的原文）：
各常量解析
CIE_2_DEG_COLOR_MATCHING_FUNCTIONS
// Values from \"CIE (1931) 2-deg color matching functions\""]},"target":{"position":0,"lines":["开始之前本人对各种渲染技术很感兴趣，但由于现实原因不能深入探索。对Unity中预计算大气渲染技术研究是出于爱好，在课余业余时间自学的，仅仅是阅读到刚刚能看懂个大概的水平。我觉得从理论/物理，到数学公式，再到代码的过程十分神奇，我的学习过程就是通过代码的实现，反推其中的几何/物理原理。尽管我是自学，并没有深厚的学术背景，但通过反向工程和参考各种资源，我逐渐拼凑出了大气散射的复杂机制及其实现方式。参考的Unity项目仓库：https://github.com/Scrawk/Brunetons-Improved-Atmospheric-Scattering参考的散射原理（后期看到多次散射才找到ebruneton的原文）：https://ebruneton.github.io/precomputed_atmospheric_scattering/atmosphere/functions.glsl.htmlhttps://youtu.be/DxfEbulyFcY?si=rCxfRidhp1PmpGh2本篇内容// 计算透射度材质"," compute.Dispatch(compute_transmittance, CONSTANTS.TRANSMITTANCE_WIDTH / NUM_THREADS, CONSTANTS.TRANSMITTANCE_HEIGHT / NUM_THREADS, 1);"," // 计算没有经过任何散射的地面直接辐照度"," compute.Dispatch(compute_direct_irradiance, CONSTANTS.IRRADIANCE_WIDTH / NUM_THREADS, CONSTANTS.IRRADIANCE_HEIGHT / NUM_THREADS, 1);"," // 计算3D单次散射纹理"," for (int layer = 0; layer < CONSTANTS.SCATTERING_DEPTH; ++layer) "," {"," // 计算3D texture ScatteringArray"," compute.SetInt(\"layer\", layer);"," compute.Dispatch(compute_single_scattering, CONSTANTS.SCATTERING_WIDTH / NUM_THREADS, CONSTANTS.SCATTERING_HEIGHT / NUM_THREADS, 1);"," }"," // 循环计算多次散射"," for (int scattering_order = 2; scattering_order <= num_scattering_orders; ++scattering_order) "," {"," // 计算大气辐射度"," for (int layer = 0; layer < CONSTANTS.SCATTERING_DEPTH; ++layer) "," {"," compute.SetInt(\"layer\", layer);"," compute.Dispatch(compute_scattering_density, CONSTANTS.SCATTERING_WIDTH / NUM_THREADS, CONSTANTS.SCATTERING_HEIGHT / NUM_THREADS, 1);"," }"," // 计算来自地面的辐照度"," compute.Dispatch(compute_indirect_irradiance, CONSTANTS.IRRADIANCE_WIDTH / NUM_THREADS, CONSTANTS.IRRADIANCE_HEIGHT / NUM_THREADS, 1);"," // 利用上面两项计算出本阶的多次散射"," for (int layer = 0; layer < CONSTANTS.SCATTERING_DEPTH; ++layer) "," {"," compute.SetInt(\"layer\", layer);"," compute.Dispatch(compute_multiple_scattering, CONSTANTS.SCATTERING_WIDTH / NUM_THREADS, CONSTANTS.SCATTERING_HEIGHT / NUM_THREADS, 1);"," }"," }本篇主要讲解各种常量上面代码中的第一阶段——计算透射度材质各常量解析CIE_2_DEG_COLOR_MATCHING_FUNCTIONS// Values from \"CIE (1931) 2-deg color matching functions\""]},"type":"CHANGE"},{"source":{"position":25,"lines":["};这个3×3矩阵用于将线性 XYZ 颜色空间的值转换为线性 sRGB 颜色空间，转换矩阵的数值取自 Wikipedia 上关于 sRGB 的标准描述。转换过程中可能还会加入伽玛校正以得到最终的显示结果。距离相关static readonly float kSunAngularRadius = 0.00935f / 2.0f;"]},"target":{"position":53,"lines":["};这个3×3矩阵用于将线性 XYZ 颜色空间的值转换为线性 sRGB 颜色空间，转换矩阵的数值取自 Wikipedia 上关于 sRGB 的标准描述。转换过程中可能还会加入伽玛校正以得到最终的显示结果。https://en.wikipedia.org/wiki/SRGB距离/角度相关static readonly float kSunAngularRadius = 0.00935f / 2.0f;"]},"type":"CHANGE"},{"source":{"position":27,"lines":["static readonly float kLengthUnitInMeters = 1000.0f;kSunAngularRadius : 太阳的角半径（除以2后），用于描述太阳在天空中的视角大小。角直径与角半径天文中描述一个天体在天空中的大小时，常用到“角直径”这个概念。角直径指的是从观察者视点看到天体两个边缘之间形成的角度。例如，太阳的角直径大约为 0.00935 弧度（约 0.53°）。角半径 则是角直径的一半，表示从天体中心到边缘所占用的角度。kBottomRadius : 地球的半径，也是大气层下界kLengthUnitInMeters : 1 单位 = 1000 米，便于缩放计算太阳辐照度数据double[] kSolarIrradiance = new double[]"]},"target":{"position":55,"lines":["static readonly float kLengthUnitInMeters = 1000.0f;","double max_sun_zenith_angle = (UseHalfPrecision ? 102.0 : 120.0) / 180.0 * Mathf.PI;kSunAngularRadius : 太阳的角半径（除以2后），用于描述太阳在天空中的视角大小。角直径与角半径天文中描述一个天体在天空中的大小时，常用到“角直径”这个概念。角直径指的是从观察者视点看到天体两个边缘之间形成的角度。例如，太阳的角直径大约为 0.00935 弧度（约 0.53°）。角半径 则是角直径的一半，表示从天体中心到边缘所占用的角度。kBottomRadius : 地球的半径，也是大气层下界kLengthUnitInMeters : 1 单位 = 1000 米，便于缩放计算max_sun_zenith_angle : 太阳光顶角的最大余弦值，在此角度下大气散射仍然有效。必须预先计算（为了获得最大的精度，需要使用最合适的最小的太阳高度角）。例如，对于地球情况，102° 是一个不错的选择，此时 mu_s_min = -0.2）。太阳辐照度数据double[] kSolarIrradiance = new double[]"]},"type":"CHANGE"},{"source":{"position":49,"lines":["double kMaxOzoneNumberDensity = 300.0 * kDobsonUnit / 15000.0;DU 等于 2.687×10^20 个臭氧分子每平方米，即如果将臭氧在一个垂直大气柱中全部压缩到标准状况下，厚度大约为 0.01 毫米。在代码中，kDobsonUnit 被定义为 2.687e20，表示每 1 DU 的臭氧分子的数量。MaxOzoneNumberDensity : 臭氧的最大分子数密度。表示假设 300 DU 均匀分布在 15 公里（15000 米）的厚度中，得到平均臭氧分子的数密度。预计算常量List<double> wavelengths = new List<double>();"]},"target":{"position":78,"lines":["double kMaxOzoneNumberDensity = 300.0 * kDobsonUnit / 15000.0;DU 等于 2.687×10^20 个臭氧分子每平方米，即如果将臭氧在一个垂直大气柱中全部压缩到标准状况下，厚度大约为 0.01 毫米。在代码中，kDobsonUnit 被定义为 2.687e20，表示每 1 DU 的臭氧分子的数量。MaxOzoneNumberDensity : 臭氧的最大分子数密度。表示假设 300 DU 均匀分布在 15 公里（15000 米）的厚度中，得到平均臭氧分子的数密度。大气密度剖面// An atmosphere layer of width 'width', and whose density is defined as","// 'exp_term' * exp('exp_scale' * h) + 'linear_term' * h + 'constant_term',","// clamped to [0,1], and where h is the altitude.","struct DensityProfileLayer {"," Length width;"," Number exp_term;"," InverseLength exp_scale;"," InverseLength linear_term;"," Number constant_term;","};一个宽度为'width'的大气层，其密度定义为'exp_term' * exp('exp_scale' * h) + 'linear_term' * h + 'constant_term'，并且被限制在[0,1]范围内，其中 h 是海拔。下面是各个密度层的定义：double kRayleighScaleHeight = 8000.0;","double kMieScaleHeight = 1200.0;DensityProfileLayer rayleigh_layer = new DensityProfileLayer(\"rayleigh\", 0.0, 1.0, -1.0 / kRayleighScaleHeight, 0.0, 0.0);","DensityProfileLayer mie_layer = new DensityProfileLayer(\"mie\", 0.0, 1.0, -1.0 / kMieScaleHeight, 0.0, 0.0);","","// Density profile increasing linearly from 0 to 1 between 10 and 25km, and","// decreasing linearly from 1 to 0 between 25 and 40km. This is an approximate","// profile from http://www.kln.ac.lk/science/Chemistry/Teaching_Resources/","// Documents/Introduction%20to%20atmospheric%20chemistry.pdf (page 10).","List<DensityProfileLayer> ozone_density = new List<DensityProfileLayer>();","ozone_density.Add(new DensityProfileLayer(\"absorption0\", 25000.0, 0.0, 0.0, 1.0 / 15000.0, -2.0 / 3.0));","ozone_density.Add(new DensityProfileLayer(\"absorption1\", 0.0, 0.0, 0.0, -1.0 / 15000.0, 8.0 / 3.0));预计算常量double kTopRadius = 6420000.0;","double kRayleigh = 1.24062e-6;","double kMieAngstromAlpha = 0.0;","double kMieAngstromBeta = 5.328e-3;","double kMieSingleScatteringAlbedo = 0.9;","double kMiePhaseFunctionG = 0.8;","double kGroundAlbedo = 0.1;List<double> wavelengths = new List<double>();"]},"type":"CHANGE"},{"source":{"position":64,"lines":[" if (UseConstantSolarSpectrum)"," solar_irradiance.Add(kConstantSolarIrradiance);"," else"," solar_irradiance.Add(kSolarIrradiance[(l - kLambdaMin) / 10]);"]},"target":{"position":118,"lines":[" solar_irradiance.Add(kSolarIrradiance[(l - kLambdaMin) / 10]);"]},"type":"CHANGE"},{"source":{"position":94,"lines":["m_model.LengthUnitInMeters = kLengthUnitInMeters;Note： kMiePhaseFunctionG 是在 Mie 散射相函数中用来描述散射方向性（或称为各向异性）的一个参数，通常称为“渐近因子”或“偏斜因子\"。Henyey-Greenstein (HG) 相函数（入射光线方向和散射光线方向有夹角θ，散射值与这个夹角有关，这部分被提取出来叫做相函数）中，θ 是散射角，而 g 就是 kMiePhaseFunctionG，表示平均散射角的余弦值。预计算Luminance？？private static void ComputeSpectralRadianceToLuminanceFactors(IList<double> wavelengths, IList<double> solar_irradiance, double lambda_power, out double k_r, out double k_g, out double k_b) "]},"target":{"position":145,"lines":["m_model.LengthUnitInMeters = kLengthUnitInMeters;Note： kMiePhaseFunctionG 是在 Mie 散射相函数中用来描述散射方向性（或称为各向异性）的一个参数，通常称为“渐近因子”或“偏斜因子\"。Henyey-Greenstein (HG) 相函数（入射光线方向和散射光线方向有夹角θ，散射值与这个夹角有关，这部分被提取出来叫做相函数）中，θ 是散射角，而 g 就是 kMiePhaseFunctionG，表示平均散射角的余弦值。预计算天空和太阳的IrradianceToLuminanceprivate static void ComputeSpectralRadianceToLuminanceFactors(IList<double> wavelengths, IList<double> solar_irradiance, double lambda_power, out double k_r, out double k_g, out double k_b) "]},"type":"CHANGE"},{"source":{"position":125,"lines":["}在函数开头，通过 Interpolate 函数分别获得在预定义的参考波长（kLambdaR、kLambdaG、kLambdaB）上的太阳辐照度。然后通过一个 for 循环，从 kLambdaMin 到 kLambdaMax（波长范围）以每1nm（dlambda = 1）的步长进行积分，累加每个波长对亮度转换因子贡献的部分。对于每个波长 lambda，调用 SampleCIEColorMatchingFunction(lambda, index) 得到 CIE 1931 标准观察者的颜色匹配函数：x_bar, y_bar, z_bar 分别对应不同颜色通道在 XYZ 三刺激值中的响应。再将 (x̅, y̅, z̅) （x_bar, y_bar, z_bar）转换成 (r_bar, g_bar, b_bar)，这三个值给出了在该波长下，对应 RGB 颜色分量的权重。用 Interpolate 获取当前波长 lambda 对应的太阳辐照度 irradiance。irradiance / solar_channel是为了归一化当前波长与参考波长下的太阳辐照度。 Math.Pow(lambda / kLambdaX, lambda_power); 是额外的权重项，用 lambda_power 调整每个波长的贡献，可能用来控制预计算时的数值范围或其他校正目的。？？在积分完成后，将每个颜色通道的累积分数乘以：CONSTANTS.MAX_LUMINOUS_EFFICACY×dlambdaprivate void SkySunRadianceToLuminance(out Vector3 skySpectralRadianceToLuminance, out Vector3 sunSpectralRadianceToLuminance)"]},"target":{"position":176,"lines":["}在函数开头，通过 Interpolate 函数分别获得在预定义的参考波长（kLambdaR、kLambdaG、kLambdaB）上的太阳辐照度。然后通过一个 for 循环，从 kLambdaMin 到 kLambdaMax（波长范围）以每1nm（dlambda = 1）的步长进行积分，累加每个波长对亮度转换因子贡献的部分。对于每个波长 lambda，调用 SampleCIEColorMatchingFunction(lambda, index) 得到 CIE 1931 标准观察者的颜色匹配函数：x_bar, y_bar, z_bar 分别对应不同颜色通道在 XYZ 三刺激值中的响应。再将 (x̅, y̅, z̅) （x_bar, y_bar, z_bar）转换成 (r_bar, g_bar, b_bar)，这三个值给出了在该波长下，对应 RGB 颜色分量的权重。用 Interpolate 获取当前波长 lambda 对应的太阳辐照度 irradiance。irradiance / solar_channel是为了归一化当前波长与参考波长下的太阳辐照度。 Math.Pow(lambda / kLambdaX, lambda_power); 是额外的权重项，用 lambda_power 调整每个波长的贡献，可能用来控制预计算时的数值范围或其他校正目的。在积分完成后，将每个颜色通道的累积分数乘以：CONSTANTS.MAX_LUMINOUS_EFFICACY×dlambdaprivate void SkySunRadianceToLuminance(out Vector3 skySpectralRadianceToLuminance, out Vector3 sunSpectralRadianceToLuminance)"]},"type":"CHANGE"},{"source":{"position":141,"lines":["}然后通过这个函数获取到 skySpectralRadianceToLuminance 和 sunSpectralRadianceToLuminance ？？第一阶段：预计算大气透射度 Transmittance下面是预计算大气渲染的第一阶段，计算 Transmittance Texture#pragma kernel ComputeTransmittance"]},"target":{"position":192,"lines":["}然后通过这个函数获取到 skySpectralRadianceToLuminance 和 sunSpectralRadianceToLuminanceWikipedia里有完整的介绍：https://en.wikipedia.org/wiki/Luminous_efficiency_function但代码中实际上的计算公式是：也许两者之间有什么相关性第一阶段：预计算大气透射度 Transmittance下面是预计算大气渲染的第一阶段，计算 Transmittance Texture#pragma kernel ComputeTransmittance"]},"type":"CHANGE"}]