Global Illumination and Path Tracing: a Practical Implementation——具体实现——重点参考

https://www.scratchapixel.com/lessons/3d-basic-rendering/global-illumination-path-tracing/global-illumination-path-tracing-practical-implementation

in the previous chapter, we introduced the principle of Monte Carlo integration to simulate global illumination. we also explained that uniformly sampling the hemisphere works well in practice for simulating indirect diffuse, but not well for simulating indirect specular or caustics. for this reason, in this lesson, we will limit ourselves to simulating indirect diffuse: that is, how diffuse objects illuminate each other as the “indirect” effect of these objects reflecting light themselves and acting somehow as a source of light.

we showed in the first chapter, how the monte carlo integration method wokred in 2D. extending the method to 3D is not “simple” and works on the same principle. here are the steps (which as u can see are very similar to the steps we followed in the 2D case).

concept: ideally we want to create samples on the hemisphere oriented about the normal N of the shaded point P. though as in the 2D case, it is easier to create these samples in a hemisphere oriented with the y-axis of some “canonical” Cartesian coordinate systems and then transform them into the shaded point local coordinate system (in which the hemisphere is oriented along the normal as shown in figure 1-bottom). the reason why we prefer to create samples in hemisphere whose up vector is oriented with the y-axis of the world Cartesian coordinate system, is because the coordinates of the sample in this coordinate system can be easily computed with the basic spherical-to-cartesian coordinates equations, which we know are:

and then transformed in the shaded point local coordinate system. We will explained in this lesson how samples are created, how to create a new coordinate system in which the shading normal NN defines the up axis of that coordinate system, and how the other axis of that local Cartesian coordinate system (which are technically tangent to the surface at PP and perpendicular to NN) are also computed. If you read the lesson on Geometry, you should know that once we have the axes of an orthogonal coordinate system CC, we can transform any point or vector (aka direction) from any Cartesian coordinate system to CC.

Step 1: create a Cartesian coordinate system in which the up vector is oriented along the shaded point normal NN (the shaded point normal NN and the up vector of the coordinate system are aligned).
Step2: create a sample using the spherical to Cartesian coordinates equations. We will show in this chapter how this can be done in practice.
Step 3: transform the sample direction from the original coordinate system to the shaded point coordinate system.
Step 4: trace a ray in the scene in the sampled direction.
Step 5: if the ray intersects an object, compute the color of that object at the intersection point and add this result to a temporary variable. Because the surface is diffuse, don’t forget to multiply the light intensity returned along each ray by the dot product between the ray direction (the light direction) and the shaded normal NN (see below).
Step 6: repeat step 2 to 5 N-times.
Step 7: divide the temporary variables that holds all the results of all the sampled rays by N, the total number of samples used. The final value is an approximation of the shaded point indirect diffuse illumination. The illumination of PP by other diffuse surfaces in the scene.

Figure 1: transforming sampling from the coordinate system in which they were created to world space (in the shaded point local coordinate system whose up vector is aligned with NN).

Figure 2: the principle of Monte Carlo integration applied to a 3D scene. Rays are randomly generated into a hemisphere of directions oriented around the shading normal NN. If these rays hit an object, the color at the intersection point is computed and returned to the shading point. All the colors are summed up and divided by the number N of cast samples. Monte Carlo theory also implies to divide each sample by the random variable PDF. We shouldn’t also forget to multiply the object’s color by the cosine of the angle between the shading normal and the ray direction (which in this case can be seen the light direction).

Technically, we need to divide each sample contribution by the random sample PDF. This will be explained later in this chapter.
In pseudo-code we get:

Vec3f castRay(Vec3f &orig, Vec3f &dir, const uint32_t &depth, ...) 
{ 
    if (depth > options.maxDepth) return 0; 
    Vec3f hitPointColor = 0; 
    // compute direct ligthing
    ... 
    // step1: compute shaded point coordinate system using normal N.
    ... 
    // number of samples N
    uint32_t N = 16; 
    Vec3f indirectDiffuse = 0; 
    for (uint32_t i = 0; i < N; ++i) { 
        // step 2: create sample in world space
        Vec3f sample = ...; 
        // step 3: transform sample from world space to shaded point local coordinate system
        sampleWorld = ...; 
        // step 4 & 5: cast a ray in this direction
        indirectDiffuse += N.dotProduct(sampleWorld) * castRay(P, sampleWorld, depth + 1, ...); 
    } 
    // step 7: divide the sum by the total number of samples N
    hitPointColor += (indirectDiffuse / N) * albedo; 
    ... 
    return hitPointColor; 
} 

What the castRay does, is essentially returning the amount of light that an object reflects towards the shaded point along the direction sampleWorld (if the ray intersected an object). The principle here to compute the contribution of that light ray to the color of PP is exactly the same than with standard light sources. In other words, if the surface is diffuse, we should also multiply the light intensity by the cosine of the angle between the light direction (sampleWorld) and the shading normal NN (line 18).
Let’s see how to implement these steps one by one.

Sep 1: Create a Local Coordinate System Oriented along the Shading Normal N
we want to build a Cartesian coordinate system in which the normal N is aligned with the up vector (or y-axis) of the coordinate system. what we know from the lesson on geometry, is that a vector can also be interpreted as the coefficients of a plane equation and that the plane in question is perpendicular to that vector:

in our case, we do not care about the variable D which set to 0 because, since we deal with vectors in this case, the plane passes through the world origin. since we are looking to build a Cartesian coordinate system, we can assume that any vector lying on that plane is perpendicular to N. Let’s call this vector NT. We could use NT as the second axis of our coordinate system. Now by taking the cross product of NT and N we would create a third vector NB perpendicular to both NT and NN. Et voila! N, NT and NBNB form a Cartesian coordinate system. So how do we find a vector lying in the plane anyway?

If you look at figure 3, you will realise the one of the vectors lying in the plane has actually its y-coordinate equals to 0. So we can assume that there is one such vector lying in this plane whose y-coordinate is actually 0 which reduces our plane equation to two unknowns:

Figure 3: two vectors lying in the plane perpendicular to the normal N exist.

Note that this equation is true if x=Nzx=Nz and z=−Nxz=−Nx or if x=−Nzx=−Nz and z=−Nxz=−Nx. So this is simple, if we take either NT=Nz,0,−NxNT=Nz,0,−Nx or NT=−Nz,0,NxNT=−Nz,0,Nx then we have a vector lying in the plane perpendicular to NN. You can choose either option as long as you take care of creating a right hand coordinate system which also depends on the order in which you order NN and NTNT in the final cross product. The only thing we need to take care of, is to normalise NTNT:

上面到底是啥意思呢？就是有了法线N，然后平面上的任意一个向量都垂直于N，假如N=(Nx,Ny,Nz).
又因为垂直于向量的平面上的向量y=0，所以有：
Nxx+Ny0+Nzz = 0
所以得到：
Nxx = -Nz*z
不妨设x=Nz，那么此时得到z=-Nx
那么此时这个平面上的向量NT=(Nz,0,-Nx)，当然还可以假设x=-Nz,那么z=Nx
都是可以的。

有了这两个向量，那么通过叉乘就得到了第三个向量，三个向量互相垂直，构成了一个坐标系。

void createCoordinateSystem(const Vec3f &N, Vec3f &Nt, Vec3f &Nb) 
{ 
    Nt = Vec3f(N.z, 0, -N.x) / sqrtf(N.x * N.x + N.z * N.z); 
    ... 
} 

All there is left now is to compute NBNB with a cross product:

void createCoordinateSystem(const Vec3f &N, Vec3f &Nt, Vec3f &Nb) 
{ 
    Nt = Vec3f(N.z, 0, -N.x) / sqrtf(N.x * N.x + N.z * N.z); 
    Nb = N.crossProduct(Nt); 
} 

this contruction though is not the best if the normal vector’s y-coordinate is greater than its x-coordinate. In this particular case, it is best for numerical robustness to build a vector within the plane whose x-coordinate is equal to 0. If you watch figure 3 again, you will see that this also an option:

The same logic applies here to find a possible solution for y and z. Finally we get:

void createCoordinateSystem(const Vec3f &N, Vec3f &Nt, Vec3f &Nb) 
{ 
    if (std::fabs(N.x) > std::fabs(N.y)) 
        Nt = Vec3f(N.z, 0, -N.x) / sqrtf(N.x * N.x + N.z * N.z); 
    else 
        Nt = Vec3f(0, -N.z, N.y) / sqrtf(N.y * N.y + N.z * N.z); 
    Nb = N.crossProduct(Nt); 
} 

our local coordinate system is built! now all we have to do is learn how to create samples, and then tranform these samples in this coordinate system.

create samples on the hemisphere
as u know, we can define the position of a point on the sphere (or hemisphere) using either polar coordinates or Cartesian coordinates. In polar coordinates a point or a vector (we are interested in direction here), are defined with two angles: θθ (the greek letter theta) and ϕϕ (the greep letter phi). The first angle is contained in the range [0,π/2][0,π/2] and the second angle (ϕϕ) is contained in the range [0,2π][0,2π]. Remember that the Monte Carlo integration method works if we chose the direction of the samples randomly. now this is where monte carlo theory kicks in. remember that what we want to do here is uniformly sample the hemisphere with respect to solid angle. Now we won’t explain the theory of Monte Carlo again this section. If you haven’t read the lesson on Monte Carlo yet, now is the time to do so (Mathematical Foundations of Monte Carlo Methods and Monte Carlo in Practive). Anyway, be ready: we are going to have to deal with some maths now.

Now the actual equation (which is called an estimator 估计量) that we are going to be using looks as follows (read “the monte carlo esitmator for the integral of the function f(x), with probability density function pdf(x) is”):

This is the generic form of a Monte Carlo estimator which we can use to approximate the result of an integral. What we want to integrate in this particular case is the function f(x)f(x) and as explained in the previous chapter, in our case this function can be seen as the amount of light arriving at PP the shaded point from every direction in the hemisphere oriented around the normal at P. it is maybe simpler to imagine what this function can look like if u consider the 2D case (or maybe just a slice though the hemisphere) as shown in figure 4. now as usual, the idea here is to sample that function for random directions. In the equation above, this is represented by the random variable XiXi. So the idea is that we generate a random direction DRDR in the hemisphere and evaluate the function for that direction, which in our case means “cast the ray from PP in the direction DRDR into the scene and return the color of the object that this ray intersected (if it intersected an object at all)”. That color represents the amount of light that the object the ray intersected reflects towards PP along the ray direction DRDR. Mathematically we also need to divide this result by the probability of generating the direction DRDR. This is represented in the equation above by the term called PDF (this stands for probability distribution function). The result of this function could be unique for each random direction XiXi. Which is why it is defined as a function. Hopefully, because in our case we want to generate uniformly distributed samples (or directions), all the rays will have the same PDF. In other words, the PDF in this case is constant, because all rays or all directions have the same probability of being generated (we speak of equiprobability). So let’s find what this PDF is and let’s then learn how to generate these random uniformly distributed samples or directions.

To sample a function we need to first find the PDF of that function, then compute its CDF then finally the inverse of that CDF. These steps were explained in detail in the lesson Mathematical Foundations of Monte Carlo Methods. First, keep in mind that we want to uniformly sample the hemisphere. That means that any of the sample or direction we will draw should have the same probability than any other drawn directions. The probably to draw direction A is the same than the probability to draw direction B, C, D, etc. Now, we know that the probability of drawing a sample from a given function is given by the function’s PDF (probability distribution function). We also know that all PDFs must integrate to 1 over their sample space (or domain of integration). In the case of the hemisphere, we will first define that sample space in terms of solid angle: the solid angle of a hemisphere is equal to 2π2π. And since all samples have the same chance (aka probability) of being generated, then our PDF must also be constant. If we write our PDF as p(x)p(x) and then integrate over the entire sample space (2π2π), then we get (where again we will express values in terms of solid angle to start with):

Because, as mentioned before, a PDF must integrate to 1 over the function sample space. Since p(ω)p(ω) is a constant, we can write (we replace the PDF by the constant C in the following equation):

We know (check the lesson the Mathematic of Shading in the section Mathematics of Computer Graphics) that:

From which we can find C or the PDF of our function:

So far, so good. Now it’s great to know the PDF of our function we want to sample to start with, but the problem is that we can’t really work with solid angles. Our ultimate goal is to generate random directions contained within the hemisphere so we need to express this PDF in terms of polar coordinates i.e. in terms of θθ and ϕϕ:

Now we know that the probability that a given solid angle belongs to a larger set can be written as:

This is just a generalisation of the equation we integrated above, but in the previous equation, ωω was equal to 2π+2π+, in other words ωω covered the area of the entire hemisphere. Now, if want to express the same probability density function but this time in terms of θθ and ϕϕ, we can write:

We learned in the previous lesson that a differential solid angle can be re-written in terms of polar coordinates as follows:

If we substitute this equation in the previous equation we get:

The dθdθ and dϕdϕ on each side of the equation cancel out and we get:

We just replaced p(θ)p(θ) with the value we computed above: p(ω)=1/2πp(ω)=1/2π. A little more patience, we are almost there. The function p(θ,ϕ)p(θ,ϕ) is called a join probability distribution. Now, we need to sample θθ and ϕϕ independently. In this particular case, the PDF is said to be separable. In other words, we can split the PDF in two PDFs, one for θθ and one for ϕϕ but we can also use another technique, which relies on the concept of marginal and conditional PDF. A marginal distribution is one in which the probability of a given variable (for example θθ) can be expressed without any reference to the other variables (for example ϕϕ). You can do this by integrating the given PDF with respect to one of the variables. In our case, we get (will integrate p(θ,ϕp(θ,ϕ) with respect to ϕϕ):

The conditional probability function gives the probably function of one of the variables (for example ϕϕ) when the other variable (θθ) is known. For continuous function (which the case of our PDFs in this example) this can simply be computed as the joint probability distribution p(θ,ϕ)p(θ,ϕ) over the marginal distribution p(θ)p(θ):

这里θ和ϕ是独立的，所以直接用p(θ,ϕ)除以p(θ)即可。

now that we have the PDFs for θ and ϕ, we need to compute their respective CDFs and inverse them. remember from the lesson on monte carlo that a CDF can be seen as a function that returns the probability that some random variable X will be found to have a value less or equal to x. to take a practical example, let us say "if i randomly generate a random vaue for θ, what is the probability that θ takes on a value less than or equal to say say for example π/4? we can compute this CDF as an integral of the PDF as follows:

the fundamential theorem of Calculus 微积分 says that if f(x) is continuous on [a,b] and if F(X) is any antiderivative 不定积分 of f(x), then

The antiderivative of sin⁡x is −cos⁡x, thus:

The CDF of the p(ϕ) is simpler:

Let’s finally invert these CDFs:
最后进行倒置：

另y=1-cosθ
其中θ是天顶角，如果随机一个y值，那么能得到θ的值是多少。

同样另y=ϕ/2π
给定y，也能计算出ϕ的值。

congratulation. the rest of the monte carlo theory now tells us that if we draw some random variable uniformaly distributed in the range [0,1], and repalce these random numbers in the equations above, then we get uniformaly distributed values of θ and ϕ。
如果将y在[0,1]之间进行均匀随机，那么能得到θ和ϕ的均匀随机。
which we can then use in the polar-to-Cartesian coordinates equations to compute a random uniformly distributed direction over the hemisphere. let us see how this works in practice:

Vec3f uniformSampleHemisphere(const float &r1, const float &r2) 
{ 
    // cos(theta) = r1 = y 这里说明下，因为r1是0~1之间，而根据公式1-y也是在0~1之间，所以直接将cos(theta)=r1也是可以的
    // cos^2(theta) + sin^2(theta) = 1 -> sin(theta) = srtf(1 - cos^2(theta))
    float sinTheta = sqrtf(1 - r1 * r1); 
    float phi = 2 * M_PI * r2; 
    float x = sinTheta * cosf(phi); 
    float z = sinTheta * sinf(phi); 
    return Vec3f(x, u1, z);  //这里有个错误u1，应该是r1
} 
 
... 
std::default_random_engine generator; 
std::uniform_real_distribution<float> distribution(0, 1); 
 
uint32_t N = 16; 
 
for (uint32_t n = 0; n < N; ++i) 
{ 
    float r1 = distribution(generator);  //用来得到随机的θ
    float r2 = distribution(generator);  //用来得到随机的ϕ
    Vec3f sample = uniformSampleHemisphere(r1, r2); 
    ... 
} 

there is a couple of interesting things to note about the code above. note that when we inverted the first CDF one of the equations was cosθ=1−y. the y variable is in fact the random variable uniformaly distributed between 0 and 1 that we want to use in order to randomly pick a value for θ. So 1−y or 1−ϵ if ϵ, the greek letter epsilon which represents our uniformly distributed random variable, is equal to cosθ. first cosθ is the y-coordinate of the randome direction. So we can set the y-coordinate of our random direction with 1−ϵ directly (line 9). Then note that
thus
great, it was a lot of maths, but now u know how to generate uniformly distributed random samples in the hemisphere. let us move to the next step and transform this sample to the shaded point local coordinate system.

step 3: transform the random samples to the shaded point local coordinate system
we already computed a local coordinate Cartesian system in which the normal of the shaded point is aligned with the up vector of that coordinate system. as we know from the lesson on geometry, the three vectors of that Cartesian coordinate system can be used to build a 4x4 matrix (or 3x3 matrix if we only care about directions) that will transform a direction from whatever space that direction is initially in, to the coordinate system from which we built the transformation matrix.

In this particular case, we will use the vectors of Cartesian coordinate system directly without building a matrix (if you are unsure about how this part of the code works, we recommend that you read the lesson on Geometry):

Vec3f Nt, Nb; 
 
createCoordinateSystem(hitNormal, Nt, Nb); 
 
for (uint32_t n = 0; n < N; ++i) { 
    float r1 = distribution(generator); // cos(theta) = N.Light Direction 
    float r2 = distribution(generator); 
    Vec3f sample = uniformSampleHemisphere(r1, r2); 
    Vec3f sampleWorld( 
        sample.x * Nb.x + sample.y * hitNormal.x + sample.z * Nt.x, 
        sample.x * Nb.y + sample.y * hitNormal.y + sample.z * Nt.y, 
        sample.x * Nb.z + sample.y * hitNormal.z + sample.z * Nt.z); 
    ... 
}

Keep in mind that they shading normal NN is the up vector. In other words it accounts for the y-axis in the coordinate system. If you look at the way we generated Nt and Nb in the function createCoordinateSystem, you will see that when the shading normal is equal to (0,1,0), Nt is aligned with the positive x-axis and Nb with the positive z-axis respectively.

step 4, 5 and 6: trace the indirect rays
the next step is to trace the ray, and accumulate the returned color to a temporary variable. once all the rays have been traced then we divide this temporary variable in which the result of all the rays have been accumulated (the amout of light returned along each indirect cast ray), by N, the total number of samples.

Vec3f indirectDiffuse = 0; 
for (uint32_t n = 0; n < N; ++i) 
{ 
    float r1 = distribution(generator); // cos(theta) = N.Light Direction 
    float r2 = distribution(generator); 
    Vec3f sample = uniformSampleHemisphere(r1, r2); 
    ... 
    // don't forget to apply cosine law (i.e. multiply by cos(theta) = r1).
    // we should also divide the result by the PDF (1/2*M_PI) but we can do this after
    indirectDiffuse += r1 * castRay(hitPoint + sampleWorld * options.bias, 
        sampleWorld, depth + 1, ...); 
} 
// divide by N and the constant PDF
indirectDiffuse /= N * (1 / (2 * M_PI)); 

do not forget that the light intensity returned along the direction sampleWorld also needs to be multiplied by the cosine of the angle between the shaded normal and the light direction (which is the ray direction in this case). Because our object is diffuse, we need to apply the cosine law (which we introduced in the first lesson on shading). Remember than the random variable r1 is actually equal to cos⁡θ which is the cosine of the angle between NN and the ray direction (check the code of the function uniformSampleHemisphere). So all there is to do in this case is multiply the result of castRay (the amount of light reaching PP from the direction sampleWorld) by r1 (line 8).
代码的第8行，直接乘以了r1，因为r1 is actually equal to cosθ，在0到1范围。
the generic equation we are trying to solve with this basic form of Monte Carlo integration is:

this equation is a simplied form of what we call in CG the rendering equation.渲染方程
this is a very important equation which we will introduce more formally in the next lesssons. The term L§ on the left of the equation is the color of the point P that we wish to compute. L§为要计算的P的颜色值，也是辐射率。
Technically, the value we are trying to compute is the radiance at P. On the right side of this equation, you can find the integral sign. As usual the domain of integration is the hemisphere of directions oriented about the surface normal. In this particular case, we deal with solid angles. Then we have the term L(ωi) 入射辐射率 which is the amount of light illuminating P from the direction ωi. Then we have the BRDF of the surface PP belongs to. The BRDF of a diffuse surface is ρ/π (where ρ is the surface albedo or color if you prefer). Diffuse surface are view independent. In other words the result of this function doesn’t depend on the light incoming and view outgoing direction. Though for most BRDFs, the result of the function depends on these two variables, which is why we have to pass them on to the function. Finally, we have the cosine law term, the cos⁡(ωi,N) term. The result of the BRDF should be attenuated by the cosine of the angle between the light incident direction and the surface normal.

Note that in this equation, the BRDF of the shading point and computing the cosine term are things we can easily do by accessing the scene data. Though the term L(ωi) is unknown and the way we evaluate it, is by tracing a ray in the scene.

This is a slightly more formal and compact way of looking at the problem we are trying to solve. As mentioned, we will get back to this equation in the next lessons.

Note that this is also where we should be dividing the result of castRay by the random variable PDF (which as we showed earlier, is equal to 1/(2π)). Because the PDF is a constant in this case, we don’t need to do this division for each ray (or each sample). We can do it after the loop (line 13). Though we will show later that this can be simplified even further. But don’t forget about the division by the PDF!

Step 7: Add the Contribution of the Indirect Diffuse to the Illumination of the Shaded Point
Final touch. Add the result of the indirect diffuse calculation to the total illumination of the shaded point and multiply the sum of the direct and indirect illumination by the object albedo of course.

hitColor = (directDiffuse + indirectDiffuse) * object->albedo / M_PI; 

Don’t forget that the BRDF of a diffuse surface is ρ/π where ρ stands for the suface albedo.
这里的ρ——漫反射率，或者直接用颜色代替。
So we need to divide the result by π. Though in the case of the indirect diffuse component, don’t forget that we should have divided the result of castRay by the PDF of the random variable, which as we showed earlier in this chapter is 1/(2π). Dividing indirectDiffuse by 1/(2π) is the same as multiplying this value by 2π. And since the albedo is also divided by π we can simplify the code above with:
化简公式，得到如下的代码：

// (indirectDiffuse / (1 / 2 * M_PI) + directDiffuse) * albedo / M_PI
// (2 * M_PI * indirectDiffuse + directDiffuse) * albedo / M_PI
// (2 * indirectDiffuse + directDiffuse / M_PI) * albedo
hitColor = (directDiffuse / M_PI + 2 * indirectDiffuse) * object->albedo; 

ET VOILA!

Monte Carlo Path Tracing: Results
Let’s now see some results. You can find the source code used to render the scene below in the last chapter of this lesson (as usual). In the image below you can see two renders: on the left, global illumination was turned off, and on the right it was turned on. Arrows in the image indicates part of the image in which the effect of other surfaces illuminating other objects in the scene can be clearly scene, such as for example the effect of the floor on the bottom right part of the sphere, or the effect of the green monolith on the left side of the sphere. Notice how the sphere takes on the red or green color of the floor and monolith respectively: the term color bleeding if often used to describe this phenomenon (because the color of objects bleed onto others).

the render time is already significantly slower compared to the time it takes to render direct lighting alone. Thus, in this particular example we have limited the number of bounces to 1, but you can experiment with higher number of bounces (if you are patient enough to wait for the result).

Monte Carlo Path Tracing is Slow and Noisy
蒙特卡洛方法慢而且会有造成的问题。

Now you as can see in the images below, Monte Carlo path tracing produces noisy images. This is the main problem with this technique. Why? If we go back to our poll人口统计例子 example, since each time you want to make estimate the average weight of the adult population of a country you choose randomly say 1,000 persons from that population, each time you will run the poll, you are likely to get 1,0000 totally different people. And thus, each time you compute the average of a new sample, you will get a different number (two samples of 1,000 individuals may as well have the exact same average weight but call this luck).
做两次的随机采样，都是1000的样本 ，很大的可能是得道不同的均值，但是也有小部分的概率是得到同样的均值。
So if you were to compute the average of say 10 groups of 1,000 adults, you you will get 10 different approximations. This is where the noise comes from. There is always more or less variation between what the exact result should be (which is what we want to compute with the integral) and the result of the Monte Carlo integral which is just an approximation. This variation is what produces noise in the image.
如果是采样，则或多或少的与准确的值有一定的偏差，所以我们理想的情况是使用积分得到准确的值，但是我们采用的是蒙特卡洛的近似方法。

噪声（方差）的大小和采样的样本的大小N有密切的关系。N越大，则估计越精确，噪声也就越小。根据方差的计算公式，可知，如果想让方程降低一倍，则需要让样本N，double才行。
How big is the variation depends essentially on the sample size, the number N. By increasing N, we can reduce the noise in the image. The image above shows a close up on the sphere for different values of N. As N increases, you can clearly see that the noise becomes less visible. However, obviously this has a cost, since render time increases with the number of samples used. Monte Carlo theory tells us that to reduce some level of noise (called variance) by 2, you need 4 times more samples than the number of samples used in the first place. In other words, to reduce the noise of the first image on the left by 2, you don’t need 32 samples but 64 samples, and to halve the noise in the image rendered with 128 samples, you would need 512 samples. You can find why in the lesson Mathematical Foundation of Monte Carlo methods.

Noise is one of the major problems with Monte Carlo path tracing. It is maybe less objectionable on a still image than it is for a sequence of frames. Many of the most recent CG animation feature films such as Finding Dory or Big Hero 6, use a special technique called de-noising. Once rendered, the noise can be partially removed using a special image processing technique.

How Do I Know that My Implementation Does the Right Thing? The Furnace Test