2b:null 2c:["$","title",null,{"children":"Numerical Integration/Quadrature | Sumant Pattanaik | Observable"}] 31:T547,md`# Numerical Integration/Quadrature The goal here is to numerically compute the result of an integration of the type ${tex.block`I=\int_{x_{min}}^{x_{max}} g(x) dx`}

To start with, we will assume that our domain is 1D.

The simplest quadrature method is mid-point quadrature method. In this method the result of the integration is approximated as the area of the rectangle whose width is ${tex`(x_{max}-x_{min})`} and height is the function evaluated at the midpoint. ${tex.block`I \approx (x_{max}-x_{min})g(x_{sample})`} The figure above shows the approximation of the integral where ${tex`x_{sample}={{x_{min}+x_{max}}\over 2}`}. Monte Carlo quadrature is very much similar, where the primary estimate of the integral is estimated by the area of the rectangle whose height is evaluated by computing the function at a random point uniformly sampled from the domain, i.e. ${tex`x_{sample}={x_{min}+\xi(x_{max}-x_{min})}`} where ${tex`\xi`} is a uniform random number between 0 and 1. The estimate can be improved by averaging the primary estimate over N samples, i.e. ${tex.block`I \approx {1 \over N}\sum_{i=1}^N (x_{max}-x_{min})g(x_{sample_i})`} `32:T424,md` Instead of choosing one sample inside the domain one may choose multiple numbers of samples and take the average of those function values evaluated at those sample points. i.e. ${tex.block`I =\sum_{i=1}^n (x_{i+1}-x_i)g_i = \sum_{i=1}^n (x_{i+1}-x_i){1 \over M}\sum_{j=1}^{M} g_{i,j}`} where ${tex`g_{i,j}`} is the function evaluated at the j-th random sample inside the interval ${tex`[x_i, x_{i+1}].`} The simplest approach to discretization is uniform discretization. Uniform discretization is easy in most domains. Discrete pieces of the domain are created from uniformly subdivided pieces in each dimension. Nonuniform discretization may be a preferable choice for various reasons. Once we know how many discrete pieces we want to break our domain to, then there may be a unique method of uniform discretization. Say, for example, there is a unique way to discretize a rectangular domain. However, even for such simple domains, once we decide to use nonuniform discretization there can be innumerable ways to create the same number of discrete pieces. So we must specify a criterion or basis for such discretization. In the following discussion, we will use the area-preserving criterion for the nonuniform discretization of 1D domain discretization. We will show its application afterward.`35:T691,md` Given below is an algorithm for achieving this goal. The algorithm uses the concept of cumulative function. The area under the 1D curve shown in the figures above can be written as ${tex.block` Area = \int_{x_{min}}^{x_{max}} x dx `} where [${tex`x_{min}, x_{max}`}] are the limits of the domain of x over which the function is defined. We can also write the expression for the area under the curve over a part of the domain [${tex`x_{min}, x`}] as ${tex.block`PartArea(x) = \int_{x_{min}}^{x} x' dx' `} This partial area is a cumulative function of x. When the value of ${tex`x=x_{max}`}, it gives the total area. We will now define a function by taking the ratio of this partial area to the total area. We will denote this function as F(x). A nice property of this function is that independent of f(x) this function evaluates to zero at ${tex`x=x_{min}`} and evaluates to 1 at ${tex`x=x_{max}`}. Furthermore, for a non-negative function f(x) this function is a monotonically increasing function. F(x) simply represents a fraction of the total area. So any fraction u of the total area is simply the integration shown above within the limit ${tex`[x_{min}, x]`} for some value of x. In principle, we can find that unknown x by inverting the function i.e. ${tex.block`x = F^{-1}(u)`} This is a root-finding problem. We will now make use of this principle to discretize the domain. In the example shown below, we wish to discretize the domain into 4 equal subdivisions. That means we need the value of x for the fraction u equal to 0.25, 0.5, and 0.75. The following figure illustrates the results of such computations. ${await FileAttachment("sq3.gif").image()} `36:T577,md`### Computing ${tex`F^{-1}(u)`} Some functions are analytically invertible. For example, for the constant and linear functions shown in the figures (left and middle plot) F's are respectively linear and quadratic functions. The equation for the roots of these functions are well known. This approach must be followed for analytically invertible functions. However, in general, it is not possible to find an analytical inversion. There exists a lot of root-finding methods in the literature. But, for monotonic function F, a simple table searching technique gives acceptable values at a reasonable cost. The algorithm in pseudo-code is given below. ~~~~ { // Find root for a given fraction u. 0<=u<=1. //Create an array of x by uniformly subdividing the interval [xmin, xmax]. for i = 1 to N+1{ x[i]=xmin+(i-1)*(xmax-xmin)/N compute Fx[i] for each x[i] // Note: 0 <= Fx[i] <= 1 } // Use binary search to search for u in Fx. minIndex = 1; maxIndex = N+1; found = 0; index = 0; while not(found) && (maxIndex-minIndex > 1){ mid = maxIndex; if (Fx[mid] == u){ found = 1; index = mid; }else if (Fx[mid] > u) maxIndex = mid - 1; else minIndex = mid + 1; } if (found) rootX = x[index]; else{ t = (u - Fx[minIndex])/(Fx[maxIndex] - Fx[minIndex]); rootX = x[minIndex]*(1-t) + x[maxIndex]*t; } return rootX; } ~~~~ `37:T4ba,md` We discussed uniform and nonuniform discretization at the beginning of this section. We can use either of these two approaches of discretization for dividing the domain of integration. The choice of discretization step can be kept independent of how the samples are chosen. When uniform discretization is used, the length of each discrete piece is the same. That means ${tex`x_i`} is a constant. Thus under uniform discretization, the integration is simply the average of the sample points. i.e. ${tex.block`I={{x_{max}-x_{min}}\over N} \sum_{i=1}^N g_i`} As emphasized earlier, the finer the discretization lower the error. The reason is easy to understand. We are making an assumption that the function is constant inside the sub-domain. Unless the function is originally constant, this assumption can only be true when the sub-domains are infinitesimally small. However, the finer the discretization also means that we have to carry out a larger number of additions (summation) and hence will take a larger amount of computation time. In the example shown below, we will see that careful function-specific non-uniform discretization can improve computational accuracy compared to uniform subdivisions. `38:T455,md`${await FileAttachment("sq4.gif").image()} The illustration above gives us an intuitive feeling that non-uniform discretization is likely to give us the best result. That means we will use the following equation to estimate the result of our integration. ${tex.block`I = {{x_{max}-x_{min}}\over N}\sum_{i=1}^N {{g(x_i)}\over {p(x_i)}}`} We have learned how to carry out area-preserving discretization of a domain based on some function defined over the domain. So we need a function defined over the domain [${tex`x_{min}, x_{max}`}]. Then the question arises, what function to use. The best function to use is some constant times g(), the function that is being integrated. Unfortunately, area-preserving discretization requires the total area under g(). Integration of a function is actually the area under the function. If we do not know the result of integration, then we can not know the area under the function. So using g() as the function for area-preserving discretization is ruled out. The best we can do is to use a function that has a shape somewhat similar to the function being integrated. ` 39:T505,md`The domain discretization approach generalizes well to multi-dimension. As earlier, our goal is to integrate a high-dimensional function. Instead of a 1D domain we can have a kD domain. ${tex.block`I = \int_{x=x_{min}}^{x_{max}}g(x) dx`} where x represents a multidimensional point. Instead of an area under the function, we have a hypervolume. The computation approach remains unchanged. We discretize the domain into a number of pieces and assume the function defined over each sub-domain (stratum) to be constant. So we have a box (high dimensional) over the sub-domain whose volume is the area of the base times constant height. The summation of the volume of all these boxes defined over all the strata gives an approximation to the total integration. ${tex.block`I = {{x_{max}-x_{min}}\over N}\sum_{i=1}^N g(x_i)`} As earlier, we can use uniform or nonuniform area discretization (stratification), and ${tex`g(x_i)`} can be the value of a single sample or computed from multiple samples within the domain, and the samples could be chosen at predefined or random locations of the subdomain. The discretization approach used in the 1D domain can be used in the kD domain as well. The process may not as simple though. We demonstrate a couple of 2D integration examples below.`3a:T71e,md` First we consider the problem of computing reflected radiance from a Lambertian surface illuminated by a rectangular light source of directionally uniform radiance. The figure below illustrates the scenario. ${await FileAttachment("sq5.gif").image()} The equation for the reflected light from a point on a differential area dA of the reflector is: ${tex.block`L_{out}={K\over \pi}\int_{A_{Light Source}}L(d_i){\cos\theta}{{\cos\theta_idA_i}\over{r_i^2}} `} where K is the material color of the Lambertian surface. The domain of integration is the rectangular area of the emitter and we must discretize this area to approximate the integration by the following summation: ${tex.block`L_{out} = {K\over \pi}\sum_{k=1}^N \Delta A_{Light Source,k}{L(d_{i,k})\cos\theta_k {{\cos\theta_{i,k}}\over{r_{i,k}^2}}} `} The simplest, of the discretization and the preferred discretization is to discretize according to a constant function defined over the rectangular area. So the delta areas are of equal size, i.e. ${tex`\Delta A_{Light Source,k} = {A_{Light Source} \over N}`}. ${tex.block`L_{out} = {K\over \pi}{A_{Light Source} \over N}\sum_{k=1}^N {L(d_{i,k})\cos\theta_k {{\cos\theta_{i,k}}\over{r_{i,k}^2}}} `} where ${tex`A_{Light Source} = Width \times Height`}. For a light source that is uniform along the outgoing directions ${tex`L(d_{i,k}) = constant = L(d_i)`}, and uniform over the area of the light source then the exitance of the light source ${tex.block`E_{Light Source} = \pi L(d_{i})`} and the flux of the light source ${tex.block`\Phi_{Light Source} = E_{Light Source}*A_{Light Source} = \pi L(d_{i})A_{Light Source}`} the estimate may be written as: ${tex.block`L_{out} = {K\over \pi}{\Phi_{Light Source} \over {\pi N}}\sum_{k=1}^N {L(d_{i,k})\cos\theta_k {{\cos\theta_{i,k}}\over{r_{i,k}^2}}} `} `3b:T8ed,md`We will take the help of the figure below for deriving this expression. ${await FileAttachment("sq6.gif").image()} From the definition of solid angle, a solid angle ${tex`d\omega`} is the area over the unit hemisphere. So ${tex.block` d\omega=d Height \times d Width= d\theta \times d\phi \sin\theta `} We now rewrite the hemispherical integration given earlier as ${tex.block`L_{out}={K\over \pi}\int_{Hemisphere}L(d_i) \cos\theta d\omega`} ${tex.block` = {K\over \pi}\int_{\theta={0}}^{\pi/2}\int_{\phi={0}}^{2\pi} L(d_{i}) \cos\theta d\phi \sin\theta d\theta `} Written in this fashion we convert the problem of integration over the hemispherical domain changes to problem of integration over ${tex`\theta`} in the domain [0, ${tex`\pi \over 2`}] and integration over ${tex`\phi`} in the domain [0, ${tex`2\pi`}]. We can discretize ${tex`\theta`} and ${tex`\phi`} in any suitable fashion we wish and compute the approximation as a double summation. ${tex.block`L_{out}= {K\over \pi}\sum_{j=1}^M\sum_{k=1}^N\int_{\theta={\theta_j}}^{\theta_{j+1}}\int_{\phi={\phi_k}}^{\phi_{k+1}} L(d_{i}) \cos\theta d\phi \sin\theta d\theta `} For example, - if we divide ${tex`\theta`} and ${tex`\phi`} uniformly then we can estimate ${tex`L_{out}`} as ${tex.block`L_{out} \approx {K\over \pi}{{{\pi\over 2}} \over M}\sum_{j=1}^M{{2\pi}\over N}\sum_{k=1}^N L(d_{i,jk}) \cos\theta_j \sin\theta_j`} ${tex.block`= {\pi K \over {M*N}} \sum_{j=1}^M \sum_{k=1}^N L(d_{i,jk}) \cos\theta_j \sin\theta_j`} - if we divide ${tex`\phi`} uniformly, but divide ${tex`\theta`} by equal area division of ${tex`\sin\theta`} then it will give equal solid angle division of hemisphere and we can estimate ${tex`L_{out}`} as ${tex.block`L_{out} \approx {K\over \pi}{{2\pi} \over {MN}}\sum_{j=1}^M\sum_{k=1}^N L(d_{i,jk}) \cos\theta_j`} ${tex.block`= {2K \over {MN}} \sum_{j=1}^M \sum_{k=1}^N L(d_{i,jk}) \cos\theta_j `} - if we divide ${tex`\phi`} uniformly, but divide ${tex`\theta`} by equal area division of ${tex`\cos\theta \sin\theta`} then it will give equal division of projected hemisphere and we can estimate ${tex`L_{out}`} as ${tex.block`L_{out} \approx {K\over \pi}{{\pi} \over {MN}}\sum_{j=1}^M\sum_{k=1}^N L(d_{i,jk})`} ${tex.block`= {K \over {MN}} \sum_{j=1}^M \sum_{k=1}^N L(d_{i,jk})`} `3c:T4c3,md` We use the term random to indicate something that is unpredictable. A process is said to be random if its outcomes at anytime is unpredictable, and if its outcome is independent of its prior outcome(s). We are assuming that we know the range of values the outcome can take. The value can be one from a discrete set of values or can be one from a continuous range. Let us take an example of rolling a dice. The outcome of rolling a dice can be any one of the six faces (numbered 1 to 6). So the set of outcomes in this case is finite and discrete. We can not predict the outcome of a single roll. But for a fair dice we can say for sure that the probability of getting any of the six faces is constant, and is equal to 1/6. That means if we roll the dice for a large number of times N, then the each face will appear about N/6 times. So the probability of the face is (N/6)/N = 1/6. The probability does not have to be constant. For an unfair dice the probability of each number will be different. We can estimate the probability of each face by rolling the dice a large number of times and counting the number of times the face of interest appears. The ratio of this number to the total number is the probability.`3d:T9a6,md`The lighting equation for a triangle light source remains the same ${tex.block`L_{out} \approx {K\over \pi}\sum_{k=1}^N \Delta A_{triangle,k}{L(d_{i,k})\cos\theta_k {{\cos\theta_{i,k}}\over{r_{i,k}^2}}} `} and using equal area subdivision we can express the esimation as ${tex.block`L_{out} \approx {K\over \pi}{A_{triangle} \over N}\sum_{k=1}^N {L(d_{i,k})\cos\theta_k {{\cos\theta_{i,k}}\over{r_{i,k}^2}}} `} The problem is how to carry out an equal area subdivision of the triangle, and how to get a point on that area to compute the direction for evaluating the equation. However, computing the numerical computation is different from the computation for a rectangular light source. The difference being the 2D domain of the triangle is not easily separated into two 1D domains. So equal-area division for stratification or uniform sampling is not as simple as that for the triangle. Using the u,v (barycentric coordinate) parameterization of the triangle we can represent any point on the triangle as ${tex.block`P(u,v) = P_0 + (P_1-P_0)u + (P_2-P_0)v`} where ${tex`0 \le u,v \le 1`} and ${tex`0 \le u+v \le 1`}. In the illustration below

the area of the triangle (on the left of the illustration) is ${tex`{1 over 2}\|P0P1\| \|P1P1\|sin\theta`} where ${tex`\theta`} is the angle of the triangle at ${tex`P_1`} and is related to uv triangle (shown on the right) by a scale factor ${tex`\|P0P1\| \|P1P1\| sin\theta`}. So the area of the triangle may be expressed as infinite summation of the differential area dudv (shown as shaded area on the uv triangle): ${tex.block`A_{triangle} \propto \int_{u=0}^1 \int_{v=0}^{1-u} dv du`} In this equation the limit of integration over the parameter v is dependent on u value and hence creates a problem in subdvision and hence sampling. To circumvent this problem we express the area as the infinite summation of the differential rectangle of width du. ${tex.block`A_{triangle} \propto 2 * \int_{u=0}^1 (u - {u^2 \over 2}) du`}. We can now use equiarea subvision principle of 1D function. ${tex.block`pdf(u) = 2 \int_{u=0}^1 (1-u)du`} ${tex.block`cdf(u) = 2(u-{u^2 \over 2})\mid_0^u = 2u-u^2`} So for equal area subdivision of a triangle into N parts the u extent of k-th subdivision computed by solving ${tex`2u-u^2 = {k \over N}`} is ${tex`\left[1 - \sqrt {1-{k \over N}},1 - \sqrt {1-{{k+1} \over N}}\right]`}. `3e:T9ff,plotPDFandSamples = (pdf, cdf, acdf, parameters) => { if (!parameters) parameters = {}; if (parameters.importance == undefined) parameters.importance = false; if (parameters.stratified == undefined) parameters.stratified = true; if (parameters.polar == undefined) parameters.polar = true; if (parameters.xLabelText == undefined) parameters.xLabelText = parameters.polar ? "theta" : "x"; if (parameters.yLabelText == undefined) parameters.yLabelText = parameters.polar ? "f(theta)" : "z"; if (parameters.N == undefined) parameters.N = 15; if (parameters.random == undefined) parameters.random = false; const thetaMax = Math.PI / 2; const increment = thetaMax / parameters.N; const thetas = d3.range(0, thetaMax + increment / 2, increment); const sampleThetas = thetas.map((theta, i) => { let x = parameters.stratified ? i / parameters.N : parameters.random ? Math.random() : (i + 0.5) / parameters.N; if (!parameters.importance) x = 1 - x; return acdf(x); }); const cumP = thetas.map((theta, i) => { //let x = i / parameters.N; return cdf(theta); }); //return sampleThetas; return Plot.plot({ height: 200, //width: 400, x: { label: parameters.xLabelText, domain: parameters.polar ? [0, 90] : [0, 1] }, y: { nice: true, line: true, label: parameters.yLabelText }, marks: [ Plot.line( thetas.map((theta) => parameters.polar ? [rad2deg(theta), pdf(theta)] : [Math.cos(theta), Math.sin(theta)] ) ), Plot.line( thetas.map((theta, i) => parameters.polar ? [rad2deg(theta), cumP[i]] : [cumP[i] * Math.cos(theta), cumP[i] * Math.sin(theta)] ), { stroke: "red", strokeDasharray: "3, 3" } ), Plot.ruleY([0]), Plot.dot( sampleThetas.map((theta) => { if (parameters.stratified) theta = parameters.importance ? acdf( (parameters.random ? Math.random() : 0.5) / parameters.N + Math.sin(theta) * Math.sin(theta) ) : acdf( Math.cos(theta) - (parameters.random ? Math.random() : 0.5) / parameters.N ); return parameters.polar ? [rad2deg(theta), 0] : [Math.cos(theta), Math.sin(theta)]; }), { fill: "black", opacity: 0.25 } ) ] }); }3f:T5f4,sampleTriangle = (inputs) => { const P0 = [1 / 4.5, 1 / 8], P1 = [6 / 9, 1 / 8], P2 = [1 / 2, 6 / 9]; const P = (u, v) => [ P0[0] + (P1[0] - P0[0]) * u + (P2[0] - P0[0]) * v, P0[1] + (P1[1] - P0[1]) * u + (P2[1] - P0[1]) * v ]; const points = [P0, P1, P2, P0]; const N = inputs.N; const uSamples = d3.range(0, N, 1).map((k) => 1 - Math.sqrt(1 - k / N)); const stratifiedSamples = uSamples .map((u_k) => { return d3.range(0, N, 1).map((j) => { const xi_u = inputs.uRandom ? Math.random() : 0.5; const u = 1 - Math.sqrt((1 - u_k) * (1 - u_k) - xi_u / N); const xi_v = inputs.vRandom ? Math.random() : 0.5; return P(u, (1 - u) * (j / N + xi_v / N)); }); }) .flat(); const ulinkEnds = uSamples.map((u) => [P(u, 0), P(u, 1 - u)]); const vlinkEnds = d3.range(0, N, 1).map((i) => [P(0, i / N), P(1, 0)]); //return stratifiedSamples; return Plot.plot({ x: { ticks: 0 }, y: { ticks: 0 }, marks: [ Plot.line(points), Plot.link(ulinkEnds, { x1: (d) => d[0][0], y1: (d) => d[0][1], x2: (d) => d[1][0], y2: (d) => d[1][1], stroke: "steelblue", strokeDasharray: "1, 3" }), Plot.link(vlinkEnds, { x1: (d) => d[0][0], y1: (d) => d[0][1], x2: (d) => d[1][0], y2: (d) => d[1][1], stroke: "steelblue", strokeDasharray: "1, 3" }), Plot.dot(stratifiedSamples, { fill: "black" }) ] }); }40:[{"id":"b0378baf67a4d16eb9c3472be79d3cda95090d458a8686d7b7a2c25a54f9492124ce2b6a4713b25ff1cfed577175bf026bb6f4e0099124764af772d0b1936e09","href":"https://static.observableusercontent.com/files/b0378baf67a4d16eb9c3472be79d3cda95090d458a8686d7b7a2c25a54f9492124ce2b6a4713b25ff1cfed577175bf026bb6f4e0099124764af772d0b1936e09","name":"sq1.gif","type":"image/gif","lastModified":1647553831633.609,"size":4119},{"id":"95f9697e3c09d571c89d356c04a6e8922c9eddcd7e9a1960e10cb1a3935f071ac5e5a5314831d290b33730dd4d736e0cc033efc7c70d44feb3bb29d2c5c29b80","href":"https://static.observableusercontent.com/files/95f9697e3c09d571c89d356c04a6e8922c9eddcd7e9a1960e10cb1a3935f071ac5e5a5314831d290b33730dd4d736e0cc033efc7c70d44feb3bb29d2c5c29b80","name":"sq6.gif","type":"image/gif","lastModified":1647553831680.789,"size":10992},{"id":"d8bfb89775318b6bf85d09576f4f9b9a9e72cd13e32ea63ff6234cf7f6c852742c2a70656a77cfac84a4b5451dd1ff298460aa0a08b4cbbf4e01d6eee0b39653","href":"https://static.observableusercontent.com/files/d8bfb89775318b6bf85d09576f4f9b9a9e72cd13e32ea63ff6234cf7f6c852742c2a70656a77cfac84a4b5451dd1ff298460aa0a08b4cbbf4e01d6eee0b39653","name":"sq2.gif","type":"image/gif","lastModified":1647553831634.218,"size":4959},{"id":"d00c9c39bdeca856fcdad98fdc244b53cd9a56927a6a2ddb63a8fccda0a79a48f7d7966f6b0c2a89e21a2db025f9106459b5eea30703f6e97fb4812d3d2b5c5a","href":"https://static.observableusercontent.com/files/d00c9c39bdeca856fcdad98fdc244b53cd9a56927a6a2ddb63a8fccda0a79a48f7d7966f6b0c2a89e21a2db025f9106459b5eea30703f6e97fb4812d3d2b5c5a","name":"sq3.gif","type":"image/gif","lastModified":1647553831657.459,"size":5073},{"id":"b4619cc7f119eb6b4874c6594f32c16a2328d09f6a78da180e134e50510c6aa28562dba6a25b49df09c5dc059bd83bd7e26ea3faf68fb0eed8ddd00c4007f037","href":"https://static.observableusercontent.com/files/b4619cc7f119eb6b4874c6594f32c16a2328d09f6a78da180e134e50510c6aa28562dba6a25b49df09c5dc059bd83bd7e26ea3faf68fb0eed8ddd00c4007f037","name":"triangleArea@2.png","type":"image/png","lastModified":1647887161239.414,"size":59031},{"id":"0c00ce247ba328c857422d81dc585f724835e0e9df40ccaad69b0408a78e83f534847b35e75593d89cfc8710e7bacae41f06dfc90798c065453b5f498017da1b","href":"https://static.observableusercontent.com/files/0c00ce247ba328c857422d81dc585f724835e0e9df40ccaad69b0408a78e83f534847b35e75593d89cfc8710e7bacae41f06dfc90798c065453b5f498017da1b","name":"quadratureMidPointRule@2.png","type":"image/png","lastModified":1647915299319.508,"size":43722},{"id":"2ce66c925e6ff4b210678d2e0bcc411b8b50b70d1ff260e60d968bfe72808bc96f1edd01ec3b0db87c93119157b18cf899f4186cc9b9b952fedc277c30fea725","href":"https://static.observableusercontent.com/files/2ce66c925e6ff4b210678d2e0bcc411b8b50b70d1ff260e60d968bfe72808bc96f1edd01ec3b0db87c93119157b18cf899f4186cc9b9b952fedc277c30fea725","name":"sq5.gif","type":"image/gif","lastModified":1647553831662.085,"size":5975},{"id":"d67bd047824b327cd3895d31df0c71349f08706e1edcea70822b0d490933e40b1297681e272fb24f6fc3a900563eff6e36316749bddf2b99ce89c480e2005eb5","href":"https://static.observableusercontent.com/files/d67bd047824b327cd3895d31df0c71349f08706e1edcea70822b0d490933e40b1297681e272fb24f6fc3a900563eff6e36316749bddf2b99ce89c480e2005eb5","name":"sq4.gif","type":"image/gif","lastModified":1647553831681.804,"size":22387},{"id":"0570eee24d010f60f61eefb1e819b271fd55996d4882341e66dc9e96bdbc2445fb2096a9f22fcf992079ef5063c4f6e8f9cf8aacdb980a5347b93cd6186560f4","href":"https://static.observableusercontent.com/files/0570eee24d010f60f61eefb1e819b271fd55996d4882341e66dc9e96bdbc2445fb2096a9f22fcf992079ef5063c4f6e8f9cf8aacdb980a5347b93cd6186560f4","name":"integration@1.png","type":"image/png","lastModified":1647913656062.909,"size":34145},{"id":"79400ebd33fddf6e541dd555c6741b321ff6a8eccdc7577f87d12e878c485224fa643f49c0e66fd5111fb4c81db60ccf1e3c61ebac1567d8c89189d5fa66f703","href":"https://static.observableusercontent.com/files/79400ebd33fddf6e541dd555c6741b321ff6a8eccdc7577f87d12e878c485224fa643f49c0e66fd5111fb4c81db60ccf1e3c61ebac1567d8c89189d5fa66f703","name":"stratified.png","type":"image/png","lastModified":1647988764055.782,"size":129855},{"id":"15a43f823dba3d87cc61246d111ec33cb1b8c28a1d66601c6c8ae31f148c65acd9e21c70b6493194885b4707d7e61606547044e1eaa56f096d4d9e0b1789a952","href":"https://static.observableusercontent.com/files/15a43f823dba3d87cc61246d111ec33cb1b8c28a1d66601c6c8ae31f148c65acd9e21c70b6493194885b4707d7e61606547044e1eaa56f096d4d9e0b1789a952","name":"uniformVsNonUniform.png","type":"image/png","lastModified":1647990128677.833,"size":160365}] 2f:{"id":"e88e364471d335f0","slug":"numerical-integration-quadrature","title":"Numerical Integration/Quadrature","thumbnail":"c2ef8d1906d7af56af90881e434e814fb6b1e395d16feef2667dc85a3a7a35f2","access_level":"public","create_time":"$D2022-03-17T21:41:42.252Z","update_time":"$D2022-03-23T15:21:09.056Z","view_time":null,"first_public_time":"$D2022-03-19T12:58:01.386Z","fork":false,"stars":1,"starred":false,"owner":{"type":"team","id":"91415db2a19585e8","login":"spattana","name":"Sumant Pattanaik","avatar_url":"https://avatars-next.observableusercontent.com/5c3b1a0caf2ce230772f04126acb360531061440f39e4326260bd62e4a5e16f0","tier":"starter_2024"},"version":1048,"publish_version":1048,"role":null,"authored":false,"authors":[{"type":"individual","id":"a51e58c27848cfd7","login":"spattana","name":"Sumant Pattanaik","avatar_url":"https://avatars-next.observableusercontent.com/5c3b1a0caf2ce230772f04126acb360531061440f39e4326260bd62e4a5e16f0","tier":null}],"members":[{"type":"individual","id":"a51e58c27848cfd7","login":"spattana","name":"Sumant Pattanaik","avatar_url":"https://avatars-next.observableusercontent.com/5c3b1a0caf2ce230772f04126acb360531061440f39e4326260bd62e4a5e16f0","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":{"id":"a86a4a88340c0df9","slug":"ray-tracing","title":"Ray Tracing","owner":{"id":"91415db2a19585e8","login":"spattana","name":"Sumant Pattanaik"}},"type":"notebook","body":{"title":"Numerical Integration/Quadrature","cells":[{"id":30,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":926,"value":"md` The estimate/approximation of the above quadrature method is improved by dividing the domain into a number of pieces, (subdomains, or strata), and computing the integration as a sum of approximation/estimate over the smaller domains.\n

\n${tex.block`I = \\int_{x_1=x_{min}}^{x_2} g(x) dx + \\int_{x_2}^{x_3} g(x) dx +...+\\int_{x_n}^{x_{n+1}=x_{max}}g(x) dx`} or equivalently ${tex.block`I = \\sum_{i=1}^n \\int_{x_i}^{x_{i+1} }g(x) dx \\approx \\sum_{i=1}^n (x_{i+1}-x_i)g_i`}\nwhere ${tex`(x_{i+1}-x_i)g_i`} is the estimate of the integration in the smaller domiain. \nIf we make a fine enough subdivision of the domain then the function defined over the smaller domain approaches a constant function and hence the sum of the estimate approaches the exact value of the integration. \n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":92,"value":"md` The sample point may be chosen from a predefined location (say the midpoint of the domain). Then, such samplings are called deterministic sampling and the quadrature is deterministic quadrature. The sample point may be chosen at random. Then such sampling is called Monte Carlo sampling or random sampling and the quadrature is called Monte Carlo quadrature.\nOne way to choose a random sample is to make use of a random number generator as shown below.\n${tex.block`\ng_i= g\\left(x_i+(x_{i+1}-x_i)*random()\\right)\n`}\nwhere random() is a function that returns a random floating-point number within the interval [0,1].\n\nBecause we evaluate the function only at the sample points, sampling may be said to be the main operation in numerical integration. `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":111,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":122,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":100,"value":"md` Pdf's that are functionally similar to the integrand g(), provide better approximations to the actual solution. Such functions are then said to be *importance functions*, and sampling such functions is known as *importance sampling*. This importance sampling is carried out by computing the cumulative density function (cdf) of the pdf, followed by applying the inverse cumulative function to uniformly generated random numbers.\n${tex.block`cdf(x) =\\int_{x_{min}}^x p(x') dx'`}\nSample ${tex`x_i`}s are generated from the uniform random numbers as ${tex.block`x_i=cdf^{-1}\\left(random()\\right)`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":140,"value":"md`## Discretization for Numerical Quadrature`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":126,"value":"md` In the rest of the writeup we will assume that we carry out the integral evaluation by stratifying the domain and using a single sample inside the subdomain (or strata). We will next discuss the discretization process appropriate for quadrature.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":915,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":910,"value":"md`We show here 3 functions defined over a 1D domain [0,1]. \n${await FileAttachment(\"sq1.gif\").image()}\nThe first one is a constant function, the second one is a linear function and the third one is a quadrant of a circle (see illustration). Our goal is to find a suitable subdivision of our 1D domain so that the areas under the curve of those small pieces are all equal.\n\n${await FileAttachment(\"sq2.gif\").image()}\n\nUniform subdivision of the domain creates uniform areas for the first case (left plot in the illustration) but creates nonuniform areas in the remaining two cases. A non-uniform subdivision shown below is required to achieve the goal.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":904,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":896,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":894,"value":"md`### Computing F(x)\nF(x), the fractional area function, can be computed by integrating f(x). In all the examples shown in the above figures, we can compute F(x) analytically. For some functions, it may not be possible to do so, and even some of the f(u) may be available only in tabular form. Complicated functions may be tabulated. Tabulated F(x) values can be directly computed by simple summation. The code snippet for tabulated computation is given below.\n~~~~\n{\n // Let fx be an array of f(x) values at N+1 different x uniformly sampled in the interval [xmin, xmax].\n Fx[1] = fx[1];\n for i = 2 to N+1\n Fx[i]= Fx[i-1]+fx[i]\n}\n~~~~`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":143,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":148,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":157,"value":"md`Let us say for example that \n${tex`g(x) = h(x)*\\cos(\\pi/2*x)`} where h() is a some complicated function. In spite of the shape of h(), its product with cos() makes g() closer to cos(). We can thus very well use cos() as the function for our area preserving discretization, and intuitively for the same number of subdivisions the discretization based on cos() will give more accurate result compared to the one based on uniform subdivision.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":162,"value":"md`### Discretization in Higher Dimension`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":160,"value":"$39","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":212,"value":"md`### Direct Lighting Due to Area Light Source`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":153,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":215,"value":"md`### Direct Lighting due to Environment Light Source`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":177,"value":"md` Next we consider the problem of lighting due to environment light source. Let us assume that the reflector is Lambertian as earlier. The equation for such lighting is:\n\n${tex.block`L_{out}={K\\over \\pi}\\int_{Hemisphere}L(d_i) \\cos\\theta d\\omega`} \n\nHere the integration domain is a solid angle due to the hemisphere of directions around the reflection point. The numerical approximation of this integration will be \n\n${tex.block`L_{out}= {K\\over \\pi}{1 \\over N} \\sum_{k=1}^N \\int_{\\Delta\\omega_{i,k}}L(d_i) \\cos\\theta d\\omega \\approx {K \\over \\pi}{1 \\over N} \\sum_{k=1}^N \\Delta{\\omega_{i,k}} L(d_{i, k}) \\cos\\theta_k`}\n\nwhere ${tex`\\Delta{\\omega_{i,k}}`}'s are the discrete solid angle pieces or equivalently, are small pieces of area over the unit hemisphere. \nSo, how do we uniformly discretize the surface of a hemisphere? We do not know a simple way to do so.\nWe will give here an answer where the integration over the hemisphere is expressed as a double integration over two 1D domains. `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":248,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":25,"value":"md`### Random Process`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":23,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":302,"value":"md`### Probability Density Function (pdf)`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":17,"value":"md`Next we will take the example when the outcome range is continuous. One such example is the duration of wait at a traffic stop. Traffic stops are electronically set to turn from red to green in a fixed amount of time. So your wait at the traffic stop could be anything from 0 to setTime, depending on when you arrive. Though each outcome is unpredictable, it is possible to associate a probability distribution of the wait time. The probability distribution is specified in terms a density function f(t) (called probability density function) that gives us the probability of wait time in the infinitesimal time interval [t, t + dt] as f(t)dt. The total probability is f(t)dt=1.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":306,"value":"md`### Cumulative Probability Function (CDF)`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":14,"value":"md`\nIn the remaining part of this section we will discuss random sampling, or the simulation of the outcome x of a random process where x is in the continuous range [${tex`x_{min}, x_{max}`}] and whose outcome distribution is defined by a probability density function f(x).\n${tex.block`\\int_{x_{min}}^{x_{max}} f(x) dx = 1`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":327,"value":"md`### Random Sampling a PDF`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":329,"value":"md`We will define a cumulative probability function ${tex`F(x)`} as\n${tex.block`F(x) = \\int_{x_{min}}^{x} f(x') dx' = 1`}\n${tex`F(x)`} is the probability that outcome has a value less than equal to x. As it is a probability, its value is a fraction between 0 and 1. \nGiven a fraction 'a' we can find x as ${tex`F^{-1}(a)`} such that the probability of outcome less than x is 'a'. If we uniformly choose 'a' as a random number in [0, 1] then the distribution of ${tex`F^{-1}(a)`} will be equal to f(x). This makes a method for random sampling a density function, and the method is known as 'importance sampling'. `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":314,"value":"md` Most programming languages provide a pseudo uniform random number generator as a function through their math library. The function is a pseudo-random number generator because the random numbers are computed algorithmically. For the same starting seed, the algorithm generates an identical pattern of numbers. Using such a random number generator we give below a sampling algorithm.\n~~~~\nfunction sampling(pdf, xmin, xmax){\n F = cdf;\n a = random();\n return Finverse(a);\n}\n~~~~\nWe have discussed methods of computing cumulative function ${tex`F`}, and of computing ${tex`F^{-1}(a)`} earlier in the context of domain discretization. The same methods are used here. Wherever possible F and ${tex`F^{-1}(a)`} are evaluated analytically. Most often tabulation and binary search on the table are used for ${tex`F`} and ${tex`F^{-1}(a)`} evaluation.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":621,"value":"viewof inputs1 = Inputs.form({\n stratified: Inputs.toggle({ value: true, label: \"Stratified\" }),\n polar: Inputs.toggle({ value: true, label: \"Polar Plot\" }),\n random: Inputs.toggle({ value: true, label: \"Random\" }),\n N: Inputs.range([1, 100], {\n step: 1,\n value: 15,\n label: \"Number of Samples\"\n })\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":350,"value":"md`Example: Sampling of ${tex`\\sin\\theta`} for ${tex`\\theta \\in [0, \\pi/2]`} \n${plotPDFandSamples(\n (x) => Math.sin(x),\n (x) => 1 - Math.cos(x),\n (x) => Math.acos(x),\n {\n importance: false,\n polar: inputs1.polar,\n stratified: inputs1.stratified,\n random: inputs1.random,\n N: inputs1.N\n }\n)}\nConfirm ${tex`\\sin\\theta`} is a pdf. \n${tex.block`\\int_0^{\\pi/2} \\sin\\theta d\\theta = -\\cos\\theta \\vert_0^{\\pi/2} = -cos{\\pi\\over 2} + cos0 = 1`}\nSo ${tex`pdf(\\theta) = \\sin\\theta`}. \nCompute cdf ${tex`F`}:\n${tex.block`F(\\theta) = \\int_0^\\theta \\sin\\theta' d\\theta' = - \\cos\\theta \\vert_0^\\theta = -\\cos\\theta + \\cos0 = 1-\\cos\\theta`}\nLet ${tex`\\zeta`} be a random number. Then ${tex`1-\\cos\\theta = \\zeta`}. \nA random sample from pdf ${tex`\\sin\\theta`} is:\n${tex`\\theta = \\cos^{-1}(1-\\zeta)`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":619,"value":"viewof inputs2 = Inputs.form({\n stratified: Inputs.toggle({ value: true, label: \"Stratified\" }),\n polar: Inputs.toggle({ value: true, label: \"Polar Plot\" }),\n random: Inputs.toggle({ value: false, label: \"Random\" }),\n N: Inputs.range([1, 100], {\n value: 15,\n step: 1,\n label: \"Number of Samples\",\n width: 200\n })\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":443,"value":"md`Example: Sampling of ${tex`\\sin\\theta\\cos\\theta`} for ${tex`\\theta \\in [0, \\pi/2]`} \n${plotPDFandSamples(\n (x) => 2 * Math.sin(x) * Math.cos(x),\n (x) => Math.sin(x) * Math.sin(x),\n (x) => Math.asin(Math.sqrt(x)),\n {\n importance: true,\n N: inputs2.N,\n polar: inputs2.polar,\n stratified: inputs2.stratified,\n random: inputs2.random\n }\n)}\nConfirm ${tex`\\sin\\theta\\cos\\theta`} is a pdf. If not make it a pdf. \n${tex.block`\\int_0^{\\pi/2} \\sin\\theta \\cos\\theta d\\theta = {1 \\over 2}\\sin^2\\theta \\vert_0^{\\pi/2} = {1 \\over 2}(1-0) = {1 \\over 2}`}\nSo ${tex`pdf(\\theta) = 2\\sin\\theta \\cos\\theta`}. \nCompute cdf ${tex`F`}:\n${tex.block`F(\\theta) = \\int_0^\\theta 2\\sin\\theta' \\cos\\theta' d\\theta' = \\sin^2\\theta' \\vert_0^\\theta = \\sin^2\\theta`}\nLet ${tex`\\zeta`} be a random number. Then ${tex`\\sin^2\\theta = \\zeta`}. \nA random sample from pdf ${tex`2\\sin\\theta\\cos\\theta`} is:\n${tex`\\theta = \\sin^{-1}\\sqrt \\zeta`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":686,"value":"md`#### Direct Lighting due to Triangle Light Source`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":699,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":807,"value":"md`Similarly, for uniform sampling of the triangle the samples must be drawn at \n${tex`u = 1 - \\sqrt {1-\\xi_u}`} or equivalently at ${tex`u = 1 - \\sqrt {\\xi_u}`}, \nand v may be drawn uniformly between 0 and 1-u, or \n${tex`v = (1-u)\\xi_v=\\sqrt {\\xi_u}\\xi_v`} \nwhere ${tex`\\xi_u`}, ${tex`\\xi_v`} are uniform random numbers between 0 and 1. \n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":767,"value":"md`For stratified sampling, the pdf(u) in the rectangular k-strata of N equal subdivision, will be\n${tex.block`pdf(u) = 2N \\int_{u=u_k}^{u_{k+1}} (1-u)du`}\n${tex.block`cdf(u) = N(2u-u^2)\\mid_{u_k}^u= N \\left((2u-u^2)-(2u_k-u_k^2)\\right)`}\nGiven this CDF, the sample u drawn uniformly in the strata can be computed by setting\n${tex.block`N \\left((2u-u^2)-(2u_k-u_k^2)\\right)=\\xi_u`}\n${tex.block`\\Rightarrow 2u-u^2={\\xi \\over N} + (2u_k-u_k^2)`}\n${tex.block`\\Rightarrow u^2-2u+1=1-{\\xi_u \\over N} - (2u_k-u_k^2)`}\n${tex.block`\\Rightarrow (u-1)^2=(1-u_k)^2 - {\\xi_u \\over N}`}\n${tex.block`\\Rightarrow u=1 - \\sqrt{(1-u_k)^2-{\\xi_u \\over N}}`}\nand \n${tex.block`v = \\left(\\sqrt{(1-u_k)^2-{\\xi_u \\over N}}\\right)\\xi_v`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":858,"value":"viewof inputs3 = Inputs.form({\n uRandom: Inputs.toggle({ value: false, label: \"Random along u\" }),\n vRandom: Inputs.toggle({ value: false, label: \"Random along v\" }),\n N: Inputs.range([1, 10], {\n step: 1,\n value: 5,\n label: \"Number of Samples along u and V\"\n })\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":825,"value":"sampleTriangle(inputs3)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":477,"value":"md`### Used functions`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":479,"value":"rad2deg = (theta) => (theta * 180) / Math.PI","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":390,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":787,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":586,"value":"md`### Imports`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":585,"value":"import {columns} from \"@bcardiff/observable-columns\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W40","stdlib":"1","titling":"auto"}} 30:"$2f:body"