Because I am also a newbie, I can't ensure that everything is correct. I hope the experts can correct me.

Zhihu's formula is a bit ugly, you can go to:GitHub

Project source code:

https://github.com/Remyuu/GAMES202-Homeworkgithub.com/Remyuu/GAMES202-Homework

Precomputed spherical harmonic coefficients

The spherical harmonics coefficients are pre-computed using the framework nori.

Ambient lighting: Calculate the spherical harmonic coefficients for each pixel of the cubemap

ProjEnv::PrecomputeCubemapSH(images, width, height, channel); Use the Riemann integral method to calculate the coefficients of the ambient light spherical harmonics.

Complete code

// TODO: here you need to compute light sh of each face of cubemap of each pixel
// TODO: Here you need to calculate the spherical harmonic coefficients of a certain face of the cubemap for each pixel
Eigen::Vector3f dir = cubemapDirs[i * width * height + y * width + x];
int index = (y * width + x) * channel;
Eigen::Array3f Le(images[i][index + 0], images[i][index + 1],
                  images[i][index + 2]);
// Describe the current angle in spherical coordinates
double theta = acos(dir.z());
double phi = atan2(dir.y(), dir.x());
// Traverse each basis function of spherical harmonics
for (int l = 0; l <= SHOrder; l++){
    for (int m = -l; m <= l; m++){
        float sh = sh::EvalSH(l, m, phi, theta);
        float delta = CalcArea((float)x, (float)y, width, height);
        SHCoeffiecents[l*(l+1)+m] += Le * sh * delta;
    }
}

analyze

Spherical harmonic coefficientsIt is the projection of the spherical harmonic function on a sphere, which can be used to represent the distribution of the function on the sphere. Since we have three channels of RGB values, the spherical harmonic coefficients we will store as a three-dimensional vector. Parts that need to be improved:

/// prt.cpp - PrecomputeCubemapSH()
// TODO: here you need to compute light sh of each face of cubemap of each pixel
// TODO: Here you need to calculate the spherical harmonic coefficients of a certain face of the cubemap for each pixel
Eigen::Vector3f dir = cubemapDirs[i * width * height + y * width + x];
int index = (y * width + x) * channel;
Eigen::Array3f Le(images[i][index + 0], images[i][index + 1],
                  images[i][index + 2]);

First, we sample a direction (a 3D vector representing the direction from the center to the pixel) from each pixel of the six cubemaps (the images array) and convert the direction to spherical coordinates (theta and phi).

Then, each spherical coordinate is passed into sh::EvalSH() to calculate the real value sh of each spherical harmonic function (basis function) and the proportion delta of the spherical area occupied by each pixel in each cubemap is calculated.

Finally, we accumulate the spherical harmonic coefficients. In the code, we can accumulate all the pixels of the cubemap, which is similar to the original operation of calculating the integral of the spherical harmonic function.

$$
Ylm=∫ϕ=02π∫θ=0πf(θ,ϕ)Ylm(θ,ϕ)sin⁡(θ)dθdϕ
$$

in:

• θ is the zenith angle, ranging from 0 to π; ϕ is the azimuth angle, ranging from 0 to 2pi.
• f(θ,ϕ) is the value of the function at a point on the sphere.
• Ylm is a spherical harmonic function, which consists of the corresponding Legendre polynomials Plm and some trigonometric functions.
• l is the order of the spherical harmonics; m is the ordinal number of the spherical harmonics, ranging from −l to l.

In order to make the readers understand more specifically, here is the estimate of the discrete form of the spherical harmonics in the code, that is, the Riemann integral method for calculation.

$$
Ylm=∑i=1Nf(θi,ϕi)Ylm(θi,ϕi)Δωi
$$

in:

• f(θi,ϕi) is the value of the function at a point on the sphere.
• Ylm(θi,ϕi) is the value of the spherical harmonics at that point.
• Δωi is the tiny area or weight of the point on the sphere.
• N is the total number of discrete points.

Code Details

• Get RGB lighting information from cubemap

Eigen::Array3f Le(images[i][index + 0], images[i][index + 1],
                  images[i][index + 2]);

The value of channel is 3, corresponding to the three channels of RGB. Therefore, index points to the position of the red channel of a pixel, index + 1 points to the position of the green channel, and index + 2 points to the position of the blue channel.

• Convert direction vector to spherical coordinates

double theta = acos(dir.z());
double phi = atan2(dir.y(), dir.x());

theta is the angle from the positive z-axis to the direction of dir, and phi is the angle from the positive x-axis to the projection of dir on the xz plane.

• Traversing the basis functions of spherical harmonics

for (int l = 0; l <= SHOrder; l++){
    for (int m = -l; m <= l; m++){
        float sh = sh::EvalSH(l, m, phi, theta);
        float delta = CalcArea((float)x, (float)y, width, height);
        SHCoeffiecents[l*(l+1)+m] += Le * sh * delta;
    }
}

Unshadowed diffuse term

scene->getIntegrator()->preprocess(scene); calculation Diffuse Unshadowed Simplify the rendering equation and substitute the spherical harmonic function in the previous section to further calculate the coefficients of the spherical harmonic projection of the BRDF. The key function is ProjectFunction. We need to write a lambda expression for this function to calculate the transfer function term.

analyze

For the diffuse transmission term, we canThere are three situationsconsider:Shadowed,No shadowandMutually Reflective.

Let's first consider the simplest case without shadows. We have the rendering equation

in,

• is the incident radiance.
• It is a geometric function, and the microscopic properties of the surface are related to the direction of the incident light.
• is the incident light direction.

For a diffuse surface with equal reflection everywhere, we can simplify to Unshadowed Lighting equation

in:

• is the diffuse outgoing radiance of the point.
• is the surface normal.

The incident radiance and transfer function terms are independent of each other, as the former represents the contribution of the light sources in the scene, and the latter represents how the surface responds to the incident light. Therefore, these two components are treated independently.

Specifically, when using spherical harmonics approximation, we expand these two items separately. The input of the former is the incident direction of light, and the input of the latter is the reflection (or outgoing direction), and the expansion is two series of arrays, so we use a data structure called Look-Up Table (LUT).

auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);

Among them, the most important one is the function ProjectFunction above. We need to write a Lambda expression (shFunc) as a parameter for this function, which is used to calculate the transfer function term.

ProjectFunction function parameter passing:

• Spherical harmonic order
• Functions that need to be projected onto basis functions (that we need to write)
• Number of samples

This function will take the result returned by the Lambda function and project it onto the basis function to get the coefficient. Finally, it will add up the coefficients of each sample and multiply them by the weight to get the final coefficient of the vertex.

Complete code

Compute the geometric terms, i.e. the transfer function terms.

//prt.cpp
...
double H = wi.normalized().dot(n.normalized()) / M_PI;
if (m_Type == Type::Unshadowed){
    // TODO: here you need to calculate unshadowed transport term of a given direction
    // TODO: Here you need to calculate the unshadowed transmission term spherical harmonics value in a given direction
    return (H > 0.0) ? H : 0.0;
}

In short, remember to divide the final integral result by , and then pass it to m_TransportSHCoeffs.

Shadowed Diffuse Term

scene->getIntegrator()->preprocess(scene); calculation Diffuse Shadowed This item has an additional visible item.

analyze

The Visibility item () is a value that is either 1 or 0. The bool rayIntersect(const Ray3f &ray) function is used to reflect a ray from the vertex position to the sampling direction. If it hits the object, it is considered to be blocked and has a shadow, and 0 is returned; if the ray does not hit the object, it is still returned.

Complete code

//prt.cpp
...
double H = wi.normalized().dot(n.normalized()) / M_PI;
...
else{
    // TODO: here you need to calculate shadowed transport term of a given direction
    // TODO: Here you need to calculate the spherical harmonic value of the shadowed transmission term in a given direction
    if (H > 0.0 && !scene->rayIntersect(Ray3f(v, wi.normalized())))
        return H;
    return 0.0;
}

In short, remember to divide the final integral result by , and then pass it to m_TransportSHCoeffs.

Export calculation results

The nori framework will generate two pre-calculated result files.

Add run parameters:

./scenes/prt.xml

In prt.xml, you need to do the followingRevise, you can choose to render the ambient light cubemap. In addition, the model, camera parameters, etc. can also be modified by yourself.

//prt.xml

<!-- Render the visible surface normals -->
<integrator type="prt">
    <string name="type" value="unshadowed" />
    <integer name="bounce" value="1" />
    <integer name="PRTSampleCount" value="100" />
<!--        <string name="cubemap" value="cubemap/GraceCathedral" />-->
<!--        <string name="cubemap" value="cubemap/Indoor" />-->
<!--        <string name="cubemap" value="cubemap/Skybox" />-->
    <string name="cubemap" value="cubemap/CornellBox" />

</integrator>

Among them, the label optional value:

• type: unshadowed, shadowed, interreflection
• bounce: The number of light bounces under the interreflection type (not yet implemented)
• PRTSampleCount: The number of samples per vertex of the transmission item
• cubemap: cubemap/GraceCathedral, cubemap/Indoor, cubemap/Skybox, cubemap/CornellBox

The above pictures are the unshadowed rendering results of GraceCathedral, Indoor, Skybox and CornellBox, with a sampling number of 1.

Coloring using spherical harmonics

Manually drag the files generated by nori into the real-time rendering framework and make some changes to the real-time framework.

After the calculation in the previous chapter is completed, copy the light.txt and transport.txt in the corresponding cubemap path to the cubemap folder of the real-time rendering framework.

Precomputed data analysis

Cancel The comments on lines 88-114 in engine.js are used to parse the txt file just added.

// engine.js
// file parsing
... // Uncomment this code

Import model/create and use PRT material shader

In the materials folderEstablishFile PRTMaterial.js.

//PRTMaterial.js

class PRTMaterial extends Material {
    constructor(vertexShader, fragmentShader) {
        super({
            'uPrecomputeL[0]': { type: 'precomputeL', value: null},
            'uPrecomputeL[1]': { type: 'precomputeL', value: null},
            'uPrecomputeL[2]': { type: 'precomputeL', value: null},
        }, 
        ['aPrecomputeLT'], 
        vertexShader, fragmentShader, null);
    }
}

async Function buildPRTMaterial(vertexPath, fragmentPath) {
    let vertexShader = await getShaderString(vertexPath);
    let fragmentShader = await getShaderString(fragmentPath);

    return new PRTMaterial(vertexShader, fragmentShader);
}

Then import it in index.html.