Gritty Details

Last time, I presented a "simple" algorithm that performed volume ray marching, but I didn't really explain how it was derived. We'll spend this lesson deriving the algorithm and then show some optimizations that can be done to speed up the renderer.

Volume Rendering Derivation

First, credit where it is due. This derivation is due to Dr. Jerry Tessendorf, my advisor. This is a presentation of his derivation from a course in volume rendering.

We first look at the case where there are no lights, but the medium can emit light itself.
The algorithm for ray marching in volume rendering is essentially just the numerical approximation of the following rendering equation

$L(\vec{x}_c,\hat{n}_p) \ =\ \int_0^{\infty} ds\ C_d(\vec{x}(s))\ \rho(\vec{x}(s)) \ \exp\left\{ -\int_0^s ds' \ \rho(\vec{x}(s'))\ \kappa\right\}$

where $\vec{x}_c$ is the camera position, $\hat{n}_p$ is the pixel ray direction, $\kappa$ is the scattering coefficient.
The field $C_d(\vec{x})$ is the material color of the medium, and $\rho(\vec{x})$ is the density of the medium. The path $\vec{x}(s)$ is the straight line path from the camera along the ray direction

$\vec{x}(s) \ =\ \vec{x}_c + \hat{n}_p\ s$

and $s$ is the length of the path.

Before trying to tackle this, it's useful to notice that the limits of integration $0\rightarrow\infty$ are completely overkill because the interesting, dense material is only arranged in some limited area. We can speed up the render by simply inserting a near distance $s_0$ and a maximum path length $s_{max}$, so that the integral is

$L(\vec{x}_c,\hat{n}_p) \ =\ \int_{s_0}^{s_{max}} ds\ C_d(\vec{x}(s))\ \rho(\vec{x}(s)) \ \exp\left\{ -\int_{s_0}^s ds' \ \rho(\vec{x}(s'))\ \kappa\right\}$

Numerical evaluation begins by dividing up the range $[s_0,s_{max}]$ into segments of length $\Delta s$, which we refer to as the ray march step size. We choose this step size to be small enough that we can assume the color and density fields to be constant within the segment. (The fact that this is not true at borders and other areas of high gradient is where aliasing comes in.) Using these segments, we can define a set of ray march positions $\vec{x}_i$ with the iterative definition

$\vec{x}_{i+1} \ =\ \vec{x}_i \ +\ \hat{n}_p\ \Delta s$

and the initial condition $\vec{x}_0 = \vec{x}_c +\hat{n}_p s_0$. The full integral now reduces to integrating over a collection of segments

$L(\vec{x}_c,\hat{n}_p) \ =\ \sum_{i=0}^N\ \int_{s_0+i\Delta s}^{s_0+(i+1)\Delta s} ds\ C_d(\vec{x}(s))\ \rho(\vec{x}(s)) \ \exp\left\{ - \int_{s_0}^{s} ds' \ \rho(\vec{x}(s'))\ \kappa\right\}$

The exponential term inside the integral can also be split up into intervals, but the very last interval is being integrated over. So

$\int_{s_0}^{s} ds' \ \rho(\vec{x}(s')) \ =\ \sum_{j=0}^{i-1} \int_{s_0+j\Delta s}^{s_0+(j+1)\Delta s} ds' \ \rho(\vec{x}(s')) \ +\ \int_{s_0+i\Delta s}^{s} ds' \ \rho(\vec{x}(s'))$

If we now apply the assumption that the density is constant inside a each interval, then this reduces to

$\int_{s_0}^{s} ds' \ \rho(\vec{x}(s')) \ =\ \sum_{j=0}^{i-1} \Delta s \ \rho(\vec{x}_j) \ +\ (s-s_0-i\Delta s) \ \rho(\vec{x}_i)$

Inserting this back into the integrals over the intervals, and defining

$\Delta T_i \ \equiv\ \exp\left\{ -\Delta s \ \rho(\vec{x}_i)\ \kappa \right\}$

we get the light as

$L(\vec{x}_c,\hat{n}_p) \ =\ \sum_{i=0}^N\ C_d(\vec{x}_i)\ \rho(\vec{x}_i) \ \prod_{j=0}^{i-1}\Delta T_j\ \int_{s_0+i\Delta s}^{s_0+(i+1)\Delta s} ds\ \exp\left\{ -(s-s_0-i\Delta s) \ \rho(\vec{x}_i)\ \kappa\right\}$

The remaining integral evaluates exactly to

$\int_{s_0+i\Delta s}^{s_0+(i+1)\Delta s} ds\ \exp\left\{ -(s-s_0-i\Delta s) \ \rho(\vec{x}_i)\ \kappa\right\} \ =\ \frac{1-\Delta T_i}{\rho(\vec{x}_i)\ \kappa}$

So the whole integral evaluation is

$L(\vec{x}_c,\hat{n}_p) \ =\ \sum_{i=0}^N\ C_d(\vec{x}_i) \ \left(\prod_{j=0}^{i-1}\Delta T_j\right)\ \frac{1-\Delta T_i}{\kappa}$

Since we are creating the material color field $C_d(\vec{x})$ to suit visual interest, we can choose to absorb the constant $1/\kappa$ into the color and simplify more to

$L(\vec{x}_c,\hat{n}_p) \ =\ \sum_{i=0}^N\ C_d(\vec{x}_i) \ \left(\prod_{j=0}^{i-1}\Delta T_j\right)\ \left(1-\Delta T_i\right)$

After a small bit of massaging and brain exercising, you will see that this is similar to the  algorithm presented last time. This is a full ray marching algorithm for light emitted directly from the object, without external light sources.

A fun aside is the fact that simple OpenGL "distance" fog that was/is used to hide small draw distance is basically this exact same process, but assuming that density is constant between to planes parallel to the image plane. In that case, the integral has a closed-form solution!

Adding Lights to the Equation... (All puns intended)

Adding the ability for lights (and shadowing) is not difficult if we derive the algorithm with a variable light field as well.

When there are external lights (i.e. light sources other than the material - the light(s) may actually be located within the volume), then there is a scattering term to account for light travelling from the light source into the volume, scattering and then making its way to the camera. Just a single scatter is considered. For this situation we want to allow that the volume may have two color fields, one for emission of light and one for scattering of light. So the field we have been calling $C_d(\vec{x})$ is now the emissive color, and we now label it $C_d^{emission}(\vec{x})$. The color field for scattering we will label $C_d^{scattering}(\vec{x})$

So now the rendering equation is

$L(\vec{x}_c,\hat{n}_p) \ =\ \int_{s_0}^{s_{max}} ds\ \left\{ C_d^{emission}(\vec{x}(s)) \ +\ C_d^{scatter}(\vec{x(s)})\otimes\textrm{illum}(\vec{x(s)})\right\}\ \rho(\vec{x}(s)) \ \exp\left\{ -\int_{s_0}^s ds' \ \rho(\vec{x}(s'))\ \kappa\right\}$

The new field $\textrm{illum}(\vec{x})$ accounts for the light arriving at the point $\vec{x}$ from all of the external light sources, and scattered into the direction of the camera. This illumination is a color quantity, and interacts with the material color by multiplying it, component by component. That is the meaning of the $\otimes$ symbol.

The illumination function is the sum of the light colors, weighted by the transmissivity of the volume between the point of interest and the position of the light.

$\textrm{illum}(\vec{x}) \ =\ \sum_{i=1}^{N_{lights}}\ C_L^i\ T_L^i(\vec{x})\ \beta\left( \vec{x},\ \hat{n}_L^i(\vec{x})\cdot\hat{n}_p \right)$

where $C_L^i$ is the color of light ''i'', and $T_L^i(\vec{x})$ is the transmissivity of the volume between the current position $\vec{x}$ and the position of $\vec{x}_L^i$ of light ''i''. The quantity $\beta$ is the '''scattering cross section''' which determines what fraction of the light scatters in the direction of the camera. The scattering cross section depends on the relative angle between the ray march direction and the direction of the light from the current position, i.e. on the inner product between $\hat{n}_p$ and the direction vector pointing to the light:

$\hat{n}_L^i(\vec{x}) \ =\ \frac{\vec{x}_L^i-\vec{x}}{|\vec{x}_L^i-\vec{x}| }$.

While the phase function is an important physical quantity, for out purposes it is convenient to ignore it as a constant that we can omit from the problem.

The transmissivity function is

$T_L^i(\vec{x}) \ =\ \exp\left\{ -\int_0^{|\vec{x}_L^i-\vec{x}|} ds\ \rho(\vec{y}_i(s))\ \kappa \right\}$

with a ray march path

$\vec{y}_i(s) \ =\ \vec{x} \ +\ \hat{n}_L^i\ s$

The conversion of the integral to a ray march proceeds exactly the same as for the emission only case, and produces effectively the same ray march algorithm, with only a change to the material color in the accumulation. With this, we have the exact same algorithm as presented last time. This accounts for shadowing as well. One big question now is how to speed up the light attenuation...

Deep Shadow Maps

One neat trick we can employ entails creating what is called a Deep Shadow Map, or DSM for short. The DSM is basically a grid that will store a precomputed light attenuation so that when we need it for rendering, we can simply sample the grid. This has a number of advantages. If the density is static, but the camera is moving, then we can simply compute the DSM once and re-use that lighting information for each frame. We can also use it to separate the resolution of the lighting from the underlying density. If we want fast lower-quality lighting, then we can compute a lower resolution DSM. If you want colored lighting, then the DSM can also store colored attenuation. That extension should be pretty simple...

To actually compute the DSM, we simply need to create a grid by visiting each voxel and storing the value we would have computed for lighting directly. Then afterwards, all we need to do at render-time is sample from this grid instead of performing ray marching. A really slick implementation of this will actually compute a DSM voxel's value only if it hasn't been computed yet and then reuse the computed results, but we'll stick to the simpler approach for this. Additionally, if your light rays are parallel to a face of the DSM grid, then you can vastly speed up this computation by simply using neighboring cells during integration. Perhaps this is useful if you're simulating a light which is very far away. (Sun rays going through fog near a hill for example.)

I invite you to take a look at the tutorial code again for this lesson as it should clear up most questions about all of this.

Volume Shadowing on Polygons

It's convenient now to notice that we also have a method for computing the volume shadowing effects onto a polygon. If we imagine a polygonal surface, a density above, and then a light source above that, we would imagine that the density would cast a shadow onto the polygon. This can be easily computed now. The light that reaches the polygonal surface is simply read from the DSM just as it would be for the volume rendering! This is how the lighting for the polygonal plane beneath the smoke was computed in the following video.

All for Now

I know there was a lot of actual new concepts in this lesson, but I had a lot of people asking for the theory behind this. Hopefully the derivation was clear and informative. If there are specific topics that anyone would like for me to discuss, I'm all ears. I'm going to be uploading a simple OpenCL volume tracer/fluid simulation system to GitHub soon. I'll probably give a talk on the details of GPU volume marching then. (And then maybe a tutorial on fluid simulation some time? Anybody interested?)

Here's a video of the output from the sample code for this lesson:

Again, the code is available on GitHub. Enjoy!

What the heck is volume rendering?

My first post here will explain the process of writing a simple volume renderer. Volume rendering is necessary when we need see translucent media or objects which have no well-defined surface. In practice, we're interested in rendering things like smoke, fog, fire, dust clouds, and many other media that can be characterized as a cloud of fine particles. In a later article, I'll explain the theory behind volume rendering by phrasing it as an approximate solution to the radiant transfer problem, but in this article, we'll look at a simple volume renderer which handles isotropic single scattering. We'll keep things as simple as I think they should be and then in subsequent articles keep adding complexity and richness to the system.

A few examples of some things I've written that use volume rendering...

Representation of Volume Data

As we'll see later, volume rendering forms an image by measuring the amount of scattered light that ends up on the image plane. The amount of light that scatters off of "stuff" in the scene is dictated by a lot of things, but principally by the distribution of density along the line of sight and the placement of lights in the scene.

With this in mind, one of the main things we need is a way to read how much stuff is along a the path of a ray. That stuff is usually given by a measure of density at a position in space. So we basically need a function $\rho(\mathbf{x})$ which measures density at a position in space. In empty space, this density measure would be zero, in something like smoke, the density would be larger in regions where there are more smoke particles per unit of volume.

A lot of tutorials and information that you'd find about volume rendering mention only 3D grids as if that is the only way to associate density values to points in space. It turns out that it is incredibly useful for a volume renderer to maintain a more abstract viewpoint than "grids only". We'll revisit this soon. To be fair though, a lot of the interesting stuff that is rendered from data sets (smoke simulations, MRI scans, etc.) does require the use of 3D grids.

Reading from a 3D grid is simple enough to code, and if you've ever written a bilinear interpolation code (e.g. for reading from textures), then it's not hard to see how the concept extends to 3D. I invite you to look at the sample code linked below if you have any questions about how this is implemented. Any introductory text on computer graphics will give a thorough explanation of linear sampling.

Volume Ray Marching

The volume rendering problem basically asks us "How much light arrives at a given pixel in the output image?" In the later article describing the theory of this, we'll derive the integral calculus behind the following algorithm. Suffice it to say, what is about to follow will probably look like crazy voodoo magic. If you're interested in the gory details though, just hold on to your pants and wait for the coming article.

Essentially the main idea to keep in mind is that some amount of light will leave a light source and travel outwards. As it travels around space, it will interact with material by scattering. The scattering is proportional to the product of the density at that point in space and the scattering coefficient. (We'll see later that there is more than one kind of scattering coefficient.) Some of this light will end up at the camera plane (where our pixels are). We can follow this process backwards with little effort; By backwards, I mean that we can solve the same problem by starting at pixels on the image plane and then marching back towards the scene.

We will keep track of a quantity T, which is deemed the transmissivity, which is a measure of attenuation along the view ray. It is related to the "alpha" or opacity of the final rendered image by $T=1-\alpha$.

• For each pixel...
• T = 1.0, s = 0, finalColor = black
• for(s=0, s < maxDistance, s += ds)
• rayPosition = cameraPosition + s * rayDirection
• density = scene.sampleDensity(rayPosition)
• light = scene.sampleLighting(rayPosition)
• dT =  exp(density * -ds * kappa)
• T *= dT
• finalColor += (1.0-dT) * T / kappa * someNiceColor * light
• finalColor.alpha = 1 .0- T

And of course, some explanation...

For each pixel, we will create a new camera ray which begins at the camera's location and points to a particular pixel in exactly the same way ray tracing would. We initialize the transmissivity to 1.0, indicating that there has been no attenuation yet. We also set an "output" color to black, as we will accumulate color using that color.

Our for loop will be used to marching along our view ray and sample space. It begins at s=0 (starting at the camera), steps in intervals of size ds, and ends when we've travelled a distance of maxDistance. A few things: ds is the spacing between our sampling in space and so it almost goes without saying that this will have some relationship to aliasing. There will be a speed/quality trade off here: smaller step size means your render will take longer, but will have less aliasing, and vice versa. The maxDistance shouldn't be totally pulled from a hat either. If you take a look at this procedure above, we only accumulate color if there is a positive density at that spot, so if we imagine a scene where there is some density in small region, then all of the iterations before/after that interesting region will have no contribution to the final image. You want to have the region you march in contain the interesting region that contains density, but any ray marching outside of that area will be wasted.

This is more or less the ray marching code. It updates the position that we are sampling at by walking along the view ray.

Here we're actually sampling the scene to find out what the density is at the current position in the ray march. This will be non-negative: zero density basically means there's nothing there, and increasing values of density will basically mean more "opaque".

We haven't covered this idea quite yet, but essentially this function will return a scalar value, in the range [0,1], that represents how much light arrives at this point. Since we're mutliplying by this value later, it makes sense that if this value is larger more light has arrived at this point. On the other hand, if the value is close to zero, then this point in space is in shadow.

Computing Light Attenuation

When we need to measure the amount of light in the scene, we actually need to measure density towards the light from a point in question. Another ray march is needed! Essentially what is needed is a measure of the density between the sample position and the light emitter. As we'll see next time, this is amounts to evaluating a path integral just like we did in the process above. There's a small function in the sample code that demonstrates the similarity. During the theory portion next time, how this falls together will be more obvious.

Sample Code

I've written a super basic version of this stuff for those that just want to see code. I actually think that seeing code is one of the best ways to learn, especially when something that is theoretically complicated obscures how you go about implementing it. I've added the code here on github.

The code should compile without any trouble on any remotely modern Linux system. As far as I know, it should also work in Windows, but don't hold me to that. After you compile (simply "make" on Linux), simply run the executable that's produced, and you should get

While this isn't a super-sexy example, all of the fundamental code for "Direct Volume Rendering" (DVR) is there. Next time, we'll look to speed up execution, speak a lot about the mathematical theory behind this process, and generate some more interesting volume data to render.

If this was at all helpful, or if I just left you more confused than you were when you started, please let me know! I love to hear what works and what doesn't when it comes to teaching. Till then, ciao.

Yay! I can do Inline ... $\frac \partial {\partial t}\int_\gamma \mathbf{n} \cdot \mathbf{u} = 0$ ... equations... And also single liners...

$$\int_D d\omega = \int_{\partial D} \omega$$