Linear Interpolation of Luminance and Intensity

In Section 1.3.5, I solve the volume rendering integral for linearly varying luminance and attenuation. I parameterize the integral with the front and back values of the luminance and attenuation, as shown in the following two equations.

$$L(s) { = {L}_{\,\mathrm{b}}\frac{D-s}{D} + {L}_{\,\mathrm{f}}\frac{s}{D} }$$ (35)

$$\tau (s) { = {\tau }_{\,\mathrm{b}}\frac{D-s}{D} + {\tau }_{\,\mathrm{f}}\frac{s}{D} }$$ (36)

The solution to the volume rendering integral using the linear terms demonstrated in Equations 1.17 and 1.18 is complicated. If we wish to compute the volume rendering integral using only finite, real values, we need at least three forms of the solution. The first form has real and finite terms when ${\tau }_{\,\mathrm{b}} > {\tau }_{\,\mathrm{f}}$.

$$\begin{multline} I(D) = { I_0 e^{-D\frac{{\tau }_{\,\mathrm{b}}+{\tau }_{\,\mat... ...tau }_{\,\mathrm{b}} -{\tau }_{\,\mathrm{f}})}}}\right)\right] } \end{multline}$$

The second form has real and finite terms when ${\tau }_{\,\mathrm{b}} < {\tau }_{\,\mathrm{f}}$.

$$\begin{multline} I(D) = { I_0 e^{-D\frac{{\tau }_{\,\mathrm{f}}+{\tau }_{\,\mat... ...tau }_{\,\mathrm{f}} -{\tau }_{\,\mathrm{b}})}}}\right)\right] } \end{multline}$$

The third form is valid when ${\tau }_{\,\mathrm{b}} = {\tau }_{\,\mathrm{f}} = \tau$.

$$I(D) = { I_0 e^{-\tau D} + {L}_{\,\mathrm{b}}{\left(\frac{_{1}}{^... ...\left(1 + \frac{_{1}}{^{\tau D}} e^{-\tau D} - \frac{_{1}}{^{\tau D}}\right)} }$$ (37)

Computing Equations 1.19 through 1.21 is not straightforward. First is the problem of computing the functions $\mathrm{erf}$ and $\mathrm{erfi}$ with high numerical accuracy. Second is the problem where the value in the brackets can become quite small while the value it multiplies with becomes exceptionally large, which leads to numerical instability. Williams, Max, and Stein [105] solve these problems, and I review their solutions here.

We first consider Equation 1.19. In particular, consider the last term. For convenience, let us use the groupings ${\delta}_{\,\mathrm{b}} \equiv {\tau }_{\,\mathrm{b}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{b}} -{\tau }_{\,\mathrm{f}})}}}$ and ${\delta}_{\,\mathrm{f}} \equiv {\tau }_{\,\mathrm{f}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{b}} -{\tau }_{\,\mathrm{f}})}}}$. The third term of Equation 1.19 becomes

$$({L}_{\,\mathrm{b}} - {L}_{\,\mathrm{f}}) \frac{_{1}}{^{\sqrt{D({... ...{\,\mathrm{b}}\right) - \mathrm{erf}\left({\delta}_{\,\mathrm{f}}\right)\right]$$ (38)

In particular, note that if ${\delta}_{\,\mathrm{b}}$ and ${\delta}_{\,\mathrm{f}}$ are moderately large, then the value inside the brackets is close to zero and the exponent outside the brackets is large. The disparity between the two factors leads to numerical errors.

To compute the $\mathrm{erf}$ functions, Williams uses an approximation given by Press and colleagues [75]. (Press actually gives his approximation for the complementary error function, defined as $\mathrm{erfc}(x) = 1 - \mathrm{erf}(x)$, but the conversion is trivial.) Press demonstrates approximating $\mathrm{erf}$ as

$$\mathrm{erf}(x) \cong 1 - u(x) e^{-x^2 + p(u(x))}$$ (39)

where $u(x) = \frac{1}{1+0.5x}$ and $p(z)$ is a ninth degree Chebyshev polynomial selected to give an accurate fit to the tail of the error function. Specifically, $p(z)$ is

$$p(z) = -1.26551223 + z*(1.00002368 + z*(0.37409196$$

$$\qquad + z*(0.09678418 + z*(-0.18628806 + z*(0.27886807$$

$$\qquad \qquad + z*(-1.13520398 + z*(1.48851587 + z*(-0.82215223$$

$$\qquad \qquad \qquad + z*0.17087277))))))));$$ (40)

Substituting Equation 1.23 into Equation 1.22, we get more flexibility.

$$({L}_{\,\mathrm{b}} - {L}_{\,\mathrm{f}}) \frac{_{1}}{^{\sqrt{D({... ...mathrm{f}})e^{-{\delta}_{\,\mathrm{f}}^2+p(u({\delta}_{\,\mathrm{f}}))} \right]$$

Now we cancel the 1s and move the exponent inside the brackets.

$$({L}_{\,\mathrm{b}} - {L}_{\,\mathrm{f}}) \frac{_{1}}{^{\sqrt{D({... ...,\mathrm{f}}^2-{\delta}_{\,\mathrm{b}}^2+p(u({\delta}_{\,\mathrm{b}}))} \right]$$ (41)

Equation 1.25 does not have the disparate factors as before and therefore is generally far more accurate when evaluated.

We next consider Equation 1.20. I rewrite the last term using the groupings ${\delta}_{\,\mathrm{b}}' \equiv {\tau }_{\,\mathrm{b}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{f}} -{\tau }_{\,\mathrm{b}})}}}$ and ${\delta}_{\,\mathrm{f}}' \equiv {\tau }_{\,\mathrm{f}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{f}} -{\tau }_{\,\mathrm{b}})}}}$.

$$({L}_{\,\mathrm{b}} - {L}_{\,\mathrm{f}}) \frac{_{1}}{^{\sqrt{D({... ...\mathrm{f}}'\right) - \mathrm{erfi}\left({\delta}_{\,\mathrm{b}}'\right)\right]$$ (42)

Rather than compute $\mathrm{erfi}$ directly, Williams instead computes Dawson's integral, defined as $D(x) \equiv e^{x^2} \int_0^x e^{y^2} dy$. Dawson's integral relates to $\mathrm{erfi}$ as $D(x) = \frac{1}{2}\sqrt{\pi}e^{-x^2}\mathrm{erfi}(x)$. Rybicki [82] gives a good numerical method for computing Dawson's integral. Press and colleagues [76] also summarize the method.

Substituting $\mathrm{erfi}(x) = \frac{2}{\sqrt{\pi}}e^{x^2} D(x)$ into Equation 1.26, we again get more flexibility.

$$({L}_{\,\mathrm{b}} - {L}_{\,\mathrm{f}}) \frac{_{1}}{^{\sqrt{D({... ...^{\sqrt{\pi}}}e^{{\delta}_{\,\mathrm{b}}'^2} D({\delta}_{\,\mathrm{b}}')\right]$$

$$({L}_{\,\mathrm{b}} - {L}_{\,\mathrm{f}}) \frac{_{\sqrt{2}}}{^{\s... ...\,\mathrm{b}}'^2-{\delta}_{\,\mathrm{f}}'^2} D({\delta}_{\,\mathrm{b}}')\right]$$ (43)

Equation 1.21 contains no form of the error function and is numerically stable enough to compute directly with the exception of when the quantity $\tau D$ approaches zero. However, it is straightforward to show that $\lim_{\tau D = 0} I(D) = I_0$. Therefore, this is just a simple special case.

So, in summary of Williams and colleagues' work [105], we can accurately compute the volume rendering integral with linearly interpolated attenuation and luminance with the following equation.

$$I(D) \cong \begin{cases}\begin{gathered}{ I_0 e^{-D\frac{{\tau ... ...} = {\tau }_{\,\mathrm{f}} = \tau ) \cap (\tau D = 0) \end{cases}$$ (44)

where ${\delta}_{\,\mathrm{b}} \equiv {\tau }_{\,\mathrm{b}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{b}} -{\tau }_{\,\mathrm{f}})}}}$, ${\delta}_{\,\mathrm{f}} \equiv {\tau }_{\,\mathrm{f}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{b}} -{\tau }_{\,\mathrm{f}})}}}$, ${\delta}_{\,\mathrm{b}}' \equiv {\tau }_{\,\mathrm{b}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{f}} -{\tau }_{\,\mathrm{b}})}}}$, ${\delta}_{\,\mathrm{f}}' \equiv {\tau }_{\,\mathrm{f}} \frac{_{\sqrt{D}}}{^{\sqrt{2({\tau }_{\,\mathrm{f}} -{\tau }_{\,\mathrm{b}})}}}$, $u(x) = \frac{1}{1+0.5x}$, $p(x)$ is defined in Equation 1.24, and $D(x)$ is Dawson's integral.

The listing for a fragment program that calculates Equation 1.28 resides in Appendix 1.3.