Monte Carlo in Rendering (A Practical Example)

https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-rendering-practical-example

Rendering the McBeth Chart using Monte Carlo Integration
In lesson 5 (Colors and Digital Images) we expalined that colors could be represented as curve giving the amount of light at each wavelength (denoted with the greek letter lambda λ ) in the visible spectrum (these curves are called spectral power distribution or SPD).

Figure 1: SPD of three different objects.

examples of such curves for three different materials are shown in figure 1. as u can see, these curves can be interpreted as a function of wavelength: f(λ). generally, the visible spectrum is considered to be within the interval 380 to 730 nanometers. to convert this curve back to tristimulus 三色的 values X, Y and Z, we need to calculate three one-dimensional integrals:

Where the functions X¯,Y¯,Z¯X¯,Y¯,Z¯ are the CIE standard observer color matching functions. The resulting tristimulus values X, Y and Z can then be converted to RGB. We explained this process in details in the lesson on Color and Digital Images.

hoepfully, u have suggested that we were going to use Monte Carlo integration to solve these integrals. the SPDs we will be using are the spectral data for each of the 24 different colors used on a McBeth chart. Spectral data for these colors are readily available on the web. If you remember what we said about Monte Carlo integration in the previous chapter, this estimator is defined as:

in our example, the interval [a,b] is [380, 730] the range of wavelength defining the visible spectrum. to approximate the XYZ values for each color of this color checker, we need to choose a random wavelength within this interval, evaluate and multiply with each other the spd and the color matching functions at this random wavelength, repeat this process N times (where N is the number of samples), sum up all results, divide the final number by N, and multiply it by (b-a). in pseudo-code we get:

int N = 32; 
float lambdaMin = 380, lambdaMax = 730; 
float X = 0, Y = 0, Z = 0; 
for (int i = 0; i < N; ++i) { // for each sample 
    // select a random wavelength
    float lambda = lambdaMin + drand48() * (lambdaMax - lambdaMin); 
    // evaluate the sod and the color matching function at this wavelength,
    // multiply the results, and add the contribution of this samples to the others
    X += spd(lambda) * CIE_X(lambda); 
    Y += spd(lambda) * CIE_Y(lambda); 
    Z += spd(lambda) * CIE_Z(lambda); 
} 
// complete the integral, multiply by (b-a) and divide by N
X *= (lambdaMax - lambdaMin) / N; 
Y *= (lambdaMax - lambdaMin) / N; 
Z *= (lambdaMax - lambdaMin) / N; 

simple? indeede. to illustrate what variance is, we will actually render each color of the MacBeth chart 麦克白颜色表 as a 64x64 square using MC integration for each pixel of that square (each pixel of a color bucket, is obtained by running a new Monte Carlo simulation and therefore each pixel is likely to be “slightly” different than the others. how different depends mostly on the number of samples N used in the simulation). our final goal is to recreate something that looks like a real McBeth chart, but CG generated (and using Monte Carlo integration of course). in lesson 5, we provided the source code of a small program rendering a McBeth chart using a Riemann sum to compute the integrals. we will start from this program as it already contains the spectral data for the McBeth chart as well as the CIE color matching functions. the heart of the code is really the function in which the Monte Carlo integration is done.

unsigned int N = 32; 
 
inline float linerp(const float *f, const short &i, const float &t, const int &max) 
{ return f[i] * (1 - t) + f[std::min(max, i + 1)] * t; } 
 
void monteCarloIntegration(const short &curveIndex, float &X, float &Y, float &Z) 
{ 
    float S = 0; // sum, used to normalize XYZ values 
    for (int i = 0; i < N; ++i) { 
        float lambda = drand48() * (lambdaMax - lambdaMin); 
        float b = lambda / 10; 
        short j = (short)b; 
        float t = b - j; 
        // interpolate
        float fx = linerp(spd[curveIndex], j, t, nbins - 1); 
        b = lambda / 5; 
        j = (short)b; 
        t = b - j; 
        X += linerp(CIE_X, j, t, 2 * nbins - 1) * fx; 
        Y += linerp(CIE_Y, j, t, 2 * nbins - 1) * fx; 
        Z += linerp(CIE_Z, j, t, 2 * nbins - 1) * fx; 
        S += linerp(CIE_Y, j, t, 2 * nbins - 1); 
    } 
    // sum, normalization factor
    S *= (lambdaMax - lambdaMin) / N; 
    // integral = (b-a) * 1/N * sum_{i=0}^{N-1} f(X_i) and normalize
    X *= (lambdaMax - lambdaMin) / N / S; 
    Y *= (lambdaMax - lambdaMin) / N / S; 
    Z *= (lambdaMax - lambdaMin) / N / S; 
}