蒙特卡洛积分(Monte Carlo Integration)应用：利用蒙特卡洛积分生成 McBeth表

最新推荐文章于 2023-05-11 12:22:35 发布

0小龙虾0

最新推荐文章于 2023-05-11 12:22:35 发布

阅读量1.5k

点赞数 3

蒙特卡洛积分(Monte Carlo Integration)应用

蒙特卡洛积分

通常函数f(x)的积分：

可以解释为计算函数曲线下方的面积：

而我们的蒙特卡洛积分则是通过近似的方式来获取一个函数的积分值。现在想象一下，我们只选取一个随机值，例如[a，b]范围内的x，对x处的函数f(x)求值，然后将结果乘以（b-a）。下图显示了结果：

这是一个矩形（其中f(x)是该矩形的高度，而（ba）是其宽度），从某种意义上讲，这可以作为f(x)积分的近似值。当我们继续采样a和b之间不同随机点处的函数，将矩形的面积相加并求和平均时，所得的数字越来越接近积分的实际结果。

而当我们采样的点越来越多，获得的答案则越来越精确，这就是蒙特卡洛积分。我们上图的例子需要确保我们采样是完全平均随机的，在这里它的形式为：

但是，我们可以将蒙特卡洛积分扩展到具有任意PDF的随机变量。然后，更通用的公式是：

这意味着不只有均匀分布可以用蒙特卡洛积分，任何分布都可以使用蒙特卡洛积分，一些特殊的分布可以帮助我们更快的收敛。

颜色空间XYZ到RGB的转换与McBeth图

光谱功率分布（SPD）

通常，每个光源的特征在于它在可见光谱中的每个波长处产生的颜色的范围和强度。这可以看作是光的特征，称为光谱或光谱功率分布（SPD）。它表示为一条曲线，该曲线表示每个特定光源从可见光谱发出每种颜色的光的强度（该曲线是光波长的函数）。下图显示了典型白炽灯泡的光谱：

从该图看出它主要发出红光和红外光。

SPD的概念不仅限于光源，还可以用于定义对象的颜色，但是对于对象，我们所说的是光谱反射率曲线。光源和物体之间存在重要区别。对于光源，SPD定义了光源的颜色及其强度。对于物体，光谱反射率曲线仅代表物体表面反射的光在照亮物体的白光量中所占的比例。下图显示了三个物体（黄油，番茄和生菜）的光谱反射率曲线。黄油（A）主要反射绿色和红色光，这些光结合在一起会产生淡黄色。番茄（B）主要反射红光，而生菜（C）主要反射绿光：

XYZ颜色模型

XYZ的基色也是红，绿，蓝。我们获取XYZ需要用到CIE标准观察者颜色匹配函数的曲线，该曲线如下图：

值得注意的是绿色曲线正好对应于光度函数，它代表人眼对亮度的敏感性。从这些曲线，我们现在可以使用以下公式将光谱转换为三个值，分别称为X，Y，Z：

其中 $S_{e}(\lambda )$ 是我们上文提到的光谱功率分布的函数。我们XYZ颜色空间的模型如下：

其中Y值决定的是该颜色的亮度，因此可以直接将该颜色模型二维表示。

XYZ到RGB的转换

在计算机图形学中，RGB是我们主流使用的颜色模型。由于XYZ和RGB都可以用一个三维向量表示，因此我们直接可以使用一个三维矩阵来将一个向量转换为一个新的向量。要使用该矩阵只需要左乘XYZ向量即可。由于该变换的是固定的，我们可以直接写出该矩阵：

glm::mat3 XYZ_to_RGB = {
    { 2.3706743, -0.5138850,  0.0052982},
    {-0.9000405,  1.4253036,  -0.0146949},
    { -0.4706338, 0.0885814,  1.0093968}
};

这里要注意的是glm库中矩阵的构造是按列向量来的。所以上面的代码每一行是一个列向量。获得的RGB颜色空间的所有颜色实际上要少于XYZ空间的，下图中三角形部分则为RGB颜色空间中的颜色：

McBeth图

麦克白图表是一块硬纸板，上面有二十四个彩色贴片（六列四行）供您选择。这些颜色是固定的，并且对应于典型材料（例如人体皮肤，天空和树叶）的平均反射率：

无论在哪里生产Macbeth，这些颜色都应始终与图表相同，由于这二十四个色块的颜色始终是相同的，并且是在受控条件下产生的，因此我们可以使用分光光度计来测量这些正方形中每一个的光谱功率分布（SPD）（在恒定和受控的光照条件下）。您可以从下图的Macbeth图表中看到两种颜色（浅色皮肤和蓝天）的光谱曲线：

Macbeth图表光谱数据可以很容易地以表格形式在网上找到，该表格由24个条目组成，每个条目对应于图表上的一种颜色。根据测量的精度，可以在我们可见光（通常为380 nm至780 nm）范围内，每隔1、2、5或10纳米采样颜色的SPD。

用蒙特卡洛积分绘制McBeth图

为了尝试了解蒙特卡洛积分的应用，我们可以使用蒙特卡洛积分绘制McBeth图，绘制步骤如下：

1.获得给定的CIE标准观察者颜色匹配函数的曲线数据与Macbeth图中24种颜色的SPD函数数据。

2.利用1获得的数据通过蒙特卡洛积分计算XYZ值。

3.将XYZ转换为RGB。

4.利用获得的RGB绘制McBeth图。

我们的目的是了解蒙特卡洛积分如何应用，步骤1可以非常方便的在网上查找得出，我们数据代码为：

//CIE标准观察者颜色匹配数据
const float colorMatchingFunc[3][72] =
{
    {0.001368, 0.002236, 0.004243, 0.007650, 0.014310, 0.023190, 0.043510, 0.077630, 0.134380, 0.214770, 0.283900, 0.328500, 0.348280, 0.348060, 0.336200, 0.318700, 0.290800, 0.251100, 0.195360, 0.142100, 0.095640, 0.057950, 0.032010, 0.014700, 0.004900, 0.002400, 0.009300, 0.029100, 0.063270, 0.109600, 0.165500, 0.225750, 0.290400, 0.359700, 0.433450, 0.512050, 0.594500, 0.678400, 0.762100, 0.842500, 0.916300, 0.978600, 1.026300, 1.056700, 1.062200, 1.045600, 1.002600, 0.938400, 0.854450, 0.751400, 0.642400, 0.541900, 0.447900, 0.360800, 0.283500, 0.218700, 0.164900, 0.121200, 0.087400, 0.063600, 0.046770, 0.032900, 0.022700, 0.015840, 0.011359, 0.008111, 0.005790, 0.004106, 0.002899, 0.002049, 0.001440, 0.000000},
    {0.000039, 0.000064, 0.000120, 0.000217, 0.000396, 0.000640, 0.001210, 0.002180, 0.004000, 0.007300, 0.011600, 0.016840, 0.023000, 0.029800, 0.038000, 0.048000, 0.060000, 0.073900, 0.090980, 0.112600, 0.139020, 0.169300, 0.208020, 0.258600, 0.323000, 0.407300, 0.503000, 0.608200, 0.710000, 0.793200, 0.862000, 0.914850, 0.954000, 0.980300, 0.994950, 1.000000, 0.995000, 0.978600, 0.952000, 0.915400, 0.870000, 0.816300, 0.757000, 0.694900, 0.631000, 0.566800, 0.503000, 0.441200, 0.381000, 0.321000, 0.265000, 0.217000, 0.175000, 0.138200, 0.107000, 0.081600, 0.061000, 0.044580, 0.032000, 0.023200, 0.017000, 0.011920, 0.008210, 0.005723, 0.004102, 0.002929, 0.002091, 0.001484, 0.001047, 0.000740, 0.000520, 0.000000},
    {0.006450, 0.010550, 0.020050, 0.036210, 0.067850, 0.110200, 0.207400, 0.371300, 0.645600, 1.039050, 1.385600, 1.622960, 1.747060, 1.782600, 1.772110, 1.744100, 1.669200, 1.528100, 1.287640, 1.041900, 0.812950, 0.616200, 0.465180, 0.353300, 0.272000, 0.212300, 0.158200, 0.111700, 0.078250, 0.057250, 0.042160, 0.029840, 0.020300, 0.013400, 0.008750, 0.005750, 0.003900, 0.002750, 0.002100, 0.001800, 0.001650, 0.001400, 0.001100, 0.001000, 0.000800, 0.000600, 0.000340, 0.000240, 0.000190, 0.000100, 0.000050, 0.000030, 0.000020, 0.000010, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000},
};
 
// SPD函数数据
const float spectralData[24][36] = {
    {0.055, 0.058, 0.061, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.063, 0.065, 0.070, 0.076, 0.079, 0.081, 0.084, 0.091, 0.103, 0.119, 0.134, 0.143, 0.147, 0.151, 0.158, 0.168, 0.179, 0.188, 0.190, 0.186, 0.181, 0.182, 0.187, 0.196, 0.209},
    {0.117, 0.143, 0.175, 0.191, 0.196, 0.199, 0.204, 0.213, 0.228, 0.251, 0.280, 0.309, 0.329, 0.333, 0.315, 0.286, 0.273, 0.276, 0.277, 0.289, 0.339, 0.420, 0.488, 0.525, 0.546, 0.562, 0.578, 0.595, 0.612, 0.625, 0.638, 0.656, 0.678, 0.700, 0.717, 0.734},
    {0.130, 0.177, 0.251, 0.306, 0.324, 0.330, 0.333, 0.331, 0.323, 0.311, 0.298, 0.285, 0.269, 0.250, 0.231, 0.214, 0.199, 0.185, 0.169, 0.157, 0.149, 0.145, 0.142, 0.141, 0.141, 0.141, 0.143, 0.147, 0.152, 0.154, 0.150, 0.144, 0.136, 0.132, 0.135, 0.147},
    {0.051, 0.054, 0.056, 0.057, 0.058, 0.059, 0.060, 0.061, 0.062, 0.063, 0.065, 0.067, 0.075, 0.101, 0.145, 0.178, 0.184, 0.170, 0.149, 0.133, 0.122, 0.115, 0.109, 0.105, 0.104, 0.106, 0.109, 0.112, 0.114, 0.114, 0.112, 0.112, 0.115, 0.120, 0.125, 0.130},
    {0.144, 0.198, 0.294, 0.375, 0.408, 0.421, 0.426, 0.426, 0.419, 0.403, 0.379, 0.346, 0.311, 0.281, 0.254, 0.229, 0.214, 0.208, 0.202, 0.194, 0.193, 0.200, 0.214, 0.230, 0.241, 0.254, 0.279, 0.313, 0.348, 0.366, 0.366, 0.359, 0.358, 0.365, 0.377, 0.398},
    {0.136, 0.179, 0.247, 0.297, 0.320, 0.337, 0.355, 0.381, 0.419, 0.466, 0.510, 0.546, 0.567, 0.574, 0.569, 0.551, 0.524, 0.488, 0.445, 0.400, 0.350, 0.299, 0.252, 0.221, 0.204, 0.196, 0.191, 0.188, 0.191, 0.199, 0.212, 0.223, 0.232, 0.233, 0.229, 0.229},
    {0.054, 0.054, 0.053, 0.054, 0.054, 0.055, 0.055, 0.055, 0.056, 0.057, 0.058, 0.061, 0.068, 0.089, 0.125, 0.154, 0.174, 0.199, 0.248, 0.335, 0.444, 0.538, 0.587, 0.595, 0.591, 0.587, 0.584, 0.584, 0.590, 0.603, 0.620, 0.639, 0.655, 0.663, 0.663, 0.667},
    {0.122, 0.164, 0.229, 0.286, 0.327, 0.361, 0.388, 0.400, 0.392, 0.362, 0.316, 0.260, 0.209, 0.168, 0.138, 0.117, 0.104, 0.096, 0.090, 0.086, 0.084, 0.084, 0.084, 0.084, 0.084, 0.085, 0.090, 0.098, 0.109, 0.123, 0.143, 0.169, 0.205, 0.244, 0.287, 0.332},
    {0.096, 0.115, 0.131, 0.135, 0.133, 0.132, 0.130, 0.128, 0.125, 0.120, 0.115, 0.110, 0.105, 0.100, 0.095, 0.093, 0.092, 0.093, 0.096, 0.108, 0.156, 0.265, 0.399, 0.500, 0.556, 0.579, 0.588, 0.591, 0.593, 0.594, 0.598, 0.602, 0.607, 0.609, 0.609, 0.610},
    {0.092, 0.116, 0.146, 0.169, 0.178, 0.173, 0.158, 0.139, 0.119, 0.101, 0.087, 0.075, 0.066, 0.060, 0.056, 0.053, 0.051, 0.051, 0.052, 0.052, 0.051, 0.052, 0.058, 0.073, 0.096, 0.119, 0.141, 0.166, 0.194, 0.227, 0.265, 0.309, 0.355, 0.396, 0.436, 0.478},
    {0.061, 0.061, 0.062, 0.063, 0.064, 0.066, 0.069, 0.075, 0.085, 0.105, 0.139, 0.192, 0.271, 0.376, 0.476, 0.531, 0.549, 0.546, 0.528, 0.504, 0.471, 0.428, 0.381, 0.347, 0.327, 0.318, 0.312, 0.310, 0.314, 0.327, 0.345, 0.363, 0.376, 0.381, 0.378, 0.379},
    {0.063, 0.063, 0.063, 0.064, 0.064, 0.064, 0.065, 0.066, 0.067, 0.068, 0.071, 0.076, 0.087, 0.125, 0.206, 0.305, 0.383, 0.431, 0.469, 0.518, 0.568, 0.607, 0.628, 0.637, 0.640, 0.642, 0.645, 0.648, 0.651, 0.653, 0.657, 0.664, 0.673, 0.680, 0.684, 0.688},
    {0.066, 0.079, 0.102, 0.146, 0.200, 0.244, 0.282, 0.309, 0.308, 0.278, 0.231, 0.178, 0.130, 0.094, 0.070, 0.054, 0.046, 0.042, 0.039, 0.038, 0.038, 0.038, 0.038, 0.039, 0.039, 0.040, 0.041, 0.042, 0.044, 0.045, 0.046, 0.046, 0.048, 0.052, 0.057, 0.065},
    {0.052, 0.053, 0.054, 0.055, 0.057, 0.059, 0.061, 0.066, 0.075, 0.093, 0.125, 0.178, 0.246, 0.307, 0.337, 0.334, 0.317, 0.293, 0.262, 0.230, 0.198, 0.165, 0.135, 0.115, 0.104, 0.098, 0.094, 0.092, 0.093, 0.097, 0.102, 0.108, 0.113, 0.115, 0.114, 0.114},
    {0.050, 0.049, 0.048, 0.047, 0.047, 0.047, 0.047, 0.047, 0.046, 0.045, 0.044, 0.044, 0.045, 0.046, 0.047, 0.048, 0.049, 0.050, 0.054, 0.060, 0.072, 0.104, 0.178, 0.312, 0.467, 0.581, 0.644, 0.675, 0.690, 0.698, 0.706, 0.715, 0.724, 0.730, 0.734, 0.738},
    {0.058, 0.054, 0.052, 0.052, 0.053, 0.054, 0.056, 0.059, 0.067, 0.081, 0.107, 0.152, 0.225, 0.336, 0.462, 0.559, 0.616, 0.650, 0.672, 0.694, 0.710, 0.723, 0.731, 0.739, 0.746, 0.752, 0.758, 0.764, 0.769, 0.771, 0.776, 0.782, 0.790, 0.796, 0.799, 0.804},
    {0.145, 0.195, 0.283, 0.346, 0.362, 0.354, 0.334, 0.306, 0.276, 0.248, 0.218, 0.190, 0.168, 0.149, 0.127, 0.107, 0.100, 0.102, 0.104, 0.109, 0.137, 0.200, 0.290, 0.400, 0.516, 0.615, 0.687, 0.732, 0.760, 0.774, 0.783, 0.793, 0.803, 0.812, 0.817, 0.825},
    {0.108, 0.141, 0.192, 0.236, 0.261, 0.286, 0.317, 0.353, 0.390, 0.426, 0.446, 0.444, 0.423, 0.385, 0.337, 0.283, 0.231, 0.185, 0.146, 0.118, 0.101, 0.090, 0.082, 0.076, 0.074, 0.073, 0.073, 0.074, 0.076, 0.077, 0.076, 0.075, 0.073, 0.072, 0.074, 0.079},
    {0.189, 0.255, 0.423, 0.660, 0.811, 0.862, 0.877, 0.884, 0.891, 0.896, 0.899, 0.904, 0.907, 0.909, 0.911, 0.910, 0.911, 0.914, 0.913, 0.916, 0.915, 0.916, 0.914, 0.915, 0.918, 0.919, 0.921, 0.923, 0.924, 0.922, 0.922, 0.925, 0.927, 0.930, 0.930, 0.933},
    {0.171, 0.232, 0.365, 0.507, 0.567, 0.583, 0.588, 0.590, 0.591, 0.590, 0.588, 0.588, 0.589, 0.589, 0.591, 0.590, 0.590, 0.590, 0.589, 0.591, 0.590, 0.590, 0.587, 0.585, 0.583, 0.580, 0.578, 0.576, 0.574, 0.572, 0.571, 0.569, 0.568, 0.568, 0.566, 0.566},
    {0.144, 0.192, 0.272, 0.331, 0.350, 0.357, 0.361, 0.363, 0.363, 0.361, 0.359, 0.358, 0.358, 0.359, 0.360, 0.360, 0.361, 0.361, 0.360, 0.362, 0.362, 0.361, 0.359, 0.358, 0.355, 0.352, 0.350, 0.348, 0.345, 0.343, 0.340, 0.338, 0.335, 0.334, 0.332, 0.331},
    {0.105, 0.131, 0.163, 0.180, 0.186, 0.190, 0.193, 0.194, 0.194, 0.192, 0.191, 0.191, 0.191, 0.192, 0.192, 0.192, 0.192, 0.192, 0.192, 0.193, 0.192, 0.192, 0.191, 0.189, 0.188, 0.186, 0.184, 0.182, 0.181, 0.179, 0.178, 0.176, 0.174, 0.173, 0.172, 0.171},
    {0.068, 0.077, 0.084, 0.087, 0.089, 0.090, 0.092, 0.092, 0.091, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.089, 0.089, 0.088, 0.087, 0.086, 0.086, 0.085, 0.084, 0.084, 0.083, 0.083, 0.082, 0.081, 0.081, 0.081},
    {0.031, 0.032, 0.032, 0.033, 0.033, 0.033, 0.033, 0.033, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.032, 0.033}
};

为了实现步骤2，我们拿出蒙特卡洛积分估计器：

由于我们的X，Y，Z公式几乎一致的，因此我们拿X公式举例：

该式子中对 $S_{e}(\lambda )$ 以及 $\bar{X}(\lambda )$ 的积分可以利用蒙特卡洛积分来近似，我们取N个波长在380～730之间的样本，计算所有 $S_{e}(\lambda_{i} )$ 的和与 $\bar{X}(\lambda_{i})$ 的和，最后分别对其求平均值，所得的值就是这两个函数积分的近似解。我们分母对 $Y(\lambda )$ 积分是用来标准化的，依旧可以使用蒙特卡洛积分计算。同理Y和Z的值也是这样计算得到，代码如下：

//Monte Carlo integration(curveIndex 为 MacBeth表的颜色块索引)
void monteCarloIntegration(const int& curveIndex,float& X,float& Y,float& Z){
    float S=0;//用于归一化
    for(int i=0;i<N;i++){
        float lambda=drand48()*(lambdaMax-lambdaMin);//差值为350
        float b=lambda/10;//MacBeth光谱数据每10nm采样
        int j=b;
        float t=b-j;
        float fx=linerp(spectralData[curveIndex], j, t, nbins-1);//(spd(x))
        b=lambda/5;//XYZ颜色对照数据每5nm采样
        j=b;
        t=b-j;
        X+=linerp(colorMatchingFunc[0], j, t, 2*nbins-1)*fx;
        Y+=linerp(colorMatchingFunc[1], j, t, 2*nbins-1)*fx;
        Z+=linerp(colorMatchingFunc[2], j, t, 2*nbins-1)*fx;
        S+=linerp(colorMatchingFunc[1], j, t, 2*nbins-1);
    }
    // 积分 = (b-a) * 1/N * sum_{i=0}^{N-1} f(X_i) 然后标准化
    S*=(lambdaMax-lambdaMin)/N;
    X*=(lambdaMax-lambdaMin)/N/S;
    Y*=(lambdaMax-lambdaMin)/N/S;
    Z*=(lambdaMax-lambdaMin)/N/S;
}

由于我们给的数据是离散点，为了使结果接近函数连续的效果，我们还要对结果插值，插值代码：

//我们数据是离散的，我们为了拟合连续的函数，需要插值。
float linerp(const float *f, const int &i, const float &t, const int &max)
{
    return f[i] * (1 - t) + f[std::min(max, i + 1)] * t;
}

获取到McBeth图的XYZ值后，我们为了将其在OpenGL中展示，就使XYZ转换为RGB，其代码为：

glm::vec3 XYZtoRGB(const float &X, const float &Y, const float &Z)
{
    glm::vec3 rgb=glm::vec3(X,Y,Z);
    return XYZ_to_RGB*rgb;
}

最后绘制完成在N=32的效果为：

我们尝试提升N为64和256，其效果分别为：

可以看到N的提升会是效果也有显著的提升，但是收敛的速度确实非常的慢。

0小龙虾0

关注

3
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
蒙特卡洛积分(Monte Carlo Integration)应用：利用蒙特卡洛积分生成 McBeth表

蒙特卡洛积分(Monte Carlo Integration)应用蒙特卡洛积分通常函数f(x)的积分：可以解释为计算函数曲线下方的面积：而我们的蒙特卡洛积分则是通过近似的方式来获取一个函数的积分值。现在想象一下，我们只选取一个随机值，例如[a，b]范围内的x，对x处的函数f(x)求值，然后将结果乘以（b-a）。下图显示了结果：这是一个矩形（其中f(x)是...
复制链接

扫一扫