Riemann Sum

A simple (albeit inaccurate) method for numerical approximations of integrals is the Riemann sum. The idea behind the Riemann sum is to take an integral, for example of the form

$$\int_a^b f(x) dx$$ (25)

and break it up into a finite amount of pieces that we sum together. The sum is

$$\sum_{i=1}^n f(x_i) \Delta{}x$$ (26)

where $\Delta{}x = (b-a)/n$ and the sample $x_i$s are chosen such that $(i-1)\Delta{}x \leq x_i \leq i\Delta{}x$. As $\Delta{}s$ approaches zero (and $n$ approaches infinity), Equation 1.8 converges to Equation 1.7. Consequently, using more terms in the sum results in a more accurate numerical approximation.

Max [63] gives an overview of how we may apply the Riemann sum to the volume rendering integral. The volume rendering integral (Equation 1.6) actually has several integrals in it. The first is the integral attenuating the incoming light. We approximate $\exp\left(-\int_0^D \tau (t) dt\right)$ as follows.

$$\exp\left(\sum_{i=1}^n \tau (t_i) \Delta{}t\right) = \prod_{i=0}^n \exp\left(\tau (t_i) \Delta{}t\right) = \prod_{i=0}^n \zeta_i$$ (27)

Note that in Equation 1.9 I have substituted $\zeta_i$ for $\exp\left(\tau (t_i) \Delta{}t\right)$. Assuming that $\Delta{}t$ is fixed, we can precompute values of $\zeta_i$ based on the associated values for $\tau (t_i)$.

We approximate the second outer integral of Equation 1.6 in much the same way. $\int_0^D L(s) \tau (s) \exp\left(\int_s^D \tau (t) dt\right)$ becomes

$$\sum_{i=1}^n g_i \prod_{j=i+1}^n \zeta_j$$ (28)

Putting Equation 1.9 and Equation 1.10 together, we get the following approximation for the volume rendering integral.

$$I(D) \approx {I_0 \prod_{i=0}^n \zeta_i + \sum_{i=1}^n g_i \prod_{j=i+1}^n \zeta_j }$$

$$= { g_n + \zeta_n(g_{n-1} + \zeta_{n-1}(g_{n-2} + \zeta_{n-2}(g_{n-3} + \cdots \zeta_2(g_1 + \zeta_1 I_0) \cdots ))) }$$ (29)

Equation 1.11 yields simple front-to-back or back-to-front methods for computing it.

Until now, I have been intentionally vague on the form of $g_i$. Max [63] defines it as the traditional Riemann sum form dictates it: $g_i = L(s_i)\tau (s_i)\Delta{}s$. However, in Section 1.3.3, I show that the output intensity of a ray segment of length $\Delta{}s$ is

$$L(s_i) \left(1 - e^{-\tau (s_i) \Delta{}s}\right) = L(s_i) (1 - \zeta_i)$$ (30)

Using Equation 1.12 for $g_i$ results in a more accurate approximation of the integral. Nevertheless, we may precompute Equation 1.12 just as we can precompute $\zeta_i$. Furthermore, consider what happens when we substitute Equation 1.12 for the $g_i$s in Equation 1.11.

$$I(D) \approx { L_n (1-\zeta_n) + \zeta_n(L_{n-1} (1-\zeta_{n-1}) + \cdots \zeta_2(L_1 (1-\zeta_1) + \zeta_1 I_0) \cdots ) }$$ (31)

We can use graphics hardware to perform the basic operation $L_i (1-\zeta_i) + \zeta_i(\cdots)$ used in Equation 1.13 as described in Section 1.2.4.

Volume rendering systems use the Riemann sum method most often when sampling of the data set along viewing rays is convenient and does not result in a large loss of information. Such situations occur when rendering rectilinear grids using methods such as those described in Section 1.1.3.