Introduction

We implemented a fast, efficient method for volume rendering.  This method extends part of the ATI Radeon SDK, which has some code for basic volume rendering.  We extended it to work with multiple volumes.  First the basic volume rendering method willl be described, and then we will present our extensions.

Basic Method

The basic method makes the observation that front-to-back ray casting can be equated to back-to-front polygon rasterization.  Instead of casting rays from the camera through the volume, we created a mesh of quadrilaterals inside the volume, and then draw them in the order of the furthest quadrilaterals being drawn first.  The quadrilaterals are all transparent, and they are composited using the standard OpenGL definition of compositing (see the Blue Book, glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA )). 

The colors of the vertices of the quadrilaterals are determined using an ATI extension to OpenGL 1.3.  Specifically, we used GL_EXT_texture3D because it allowed us to do the 3D texture interpolation in hardware.  This was very fast because the target video card had enough memory (64MB) to hold our 3D volume (8MB).  The basics of the texturing will be explained here, but please reference the Red Book for the full explanation.

3D Texturing

In standard 2D texturing, each axis of an arbitrary texture is mapped to [0,1].  This is done regardless of the texture resolution.  When asking for an arbitrary point in texture space [s t r q], OpenGL multiplies that point by the GL_TEXTURE matrix and then uses the [s t] coordinate to look up the correct point in texture space in the horizontal and vertical axes.  In 3D texturing, the r component is also considered to map to the depth axis of the texture.  In fact, the transform is completely generalized so that the q component may be used for arbitrary transformations.  This process is completely analogous to the perspective divide, except that q is used for the division.

We set the GL_TEXTURE matrix to map into our 3D volume, and then we used automatic texture-coordinate generation to let OpenGL figure out the correct texture coordinate for a given point in space.  See glGetTexGen in the Blue Book for the details, and read about GL_EYE_PLANE to see how to automatically generate texture-coordinates.  The graphics hardware then uses trilinear interpolation to assign a color to the vertex.

Putting it all together

Given that we can draw an arbitrary quadrilaterals and have the graphics hardware do the trilinear interpolation, we now present how the image is drawn.  Quadrilaterals are drawn in parallel planes starting furthest from the camera.  Each plane has 80x80 square quadrilaterals at evenly spaced intervals in the plane.  A total of 100 planes are drawn.  80x80x100 was chosen because we wanted to get about 10 frames per second.  Our data is 128x128x128, so obviously not all of the data is being rendered, but the trilinear interpolation makes it likely that most interesting surfaces will be visible.

Our Extension

We extended the ATI sample to work with multiple volumes of data.  To get good performance, we kept the image size in texture memory constant.  Keeping the texture at 128x128x128 meant that we had to combine the volumes in some reasonable way.  We decided to assign each volume to a different color.  The PET volumes are all red, and the CT volumes are green.

Each voxel was composited with a voxel from every volume being rendered.  This is exactly the same as doing a composite operation, but we made one change: the opacities are added.  Below is a formula that relates the current color (RGBA) of a volume C1(C1r, C1b, C1g, C1a) to the resulting color after CK(CKr, CKb, CKg, CKa) is added.

C(r) = CKr*CKa + C1r         (1)

Equation (1) is identical for the green and blue components, but it differs for the alpha component:

C(a) = CKa + C1a         (1)

This has the effect of making overlapping regions more opaque.

Results

We achieved our desired frame rate.  The images have fair quality, and the user interface is acceptable.  We added an extra clipping plane (shown as a red box) that the user can move to reveal internal structures.  Click any of the images to see a larger version.

Coronal view of the dog's chest from the front with both PET (red) and CT (green) data.  Note the heart is visible in the middle in orange because of the compositing.	Coronal view of the dog's chest from the back.  Note that the spinal column is very visibile in CT (green) and the lungs are very easy to identify in PET (red).	Right saggital view of the chest cavity.  Note that the lungs are empty because they are filled with air, which does not appear in either modality.  The spinal cord is also quite easy to see.

Time-varying Dataset

This dataset shows the effects of putting an aerosol drug into a dog's lungs.  The CT data (green) does not change, but it gives an anatomical reference.  The PET data shows the approximate location of the drug.  Each of these images was taken 5 minutes apart.  The renderings are done in the same orientation, using the clipping plane to show the internal anatomy.

This is just the CT data.	Drug introduced (0 min)	Drug spreads throughout the lungs (5 min)	Drug is absorbed into blood circulation (10 min)	Small amount remains (15 min)	(20 min)	(25 min)	(30 min)

Advantages

The best part about this method is its speed.  It's fast enough to render a 128x128x128 volume in RGBA8 mode (8 bits for each of R,G,B, and A) on a fairly cheap graphics card.  We used an ATI FireGL E1 8800 with 64MB RAM, which used to cost about $250 new.  The card has been obsoleted by further FireGL models.

This method is also useful for its relative simplicity of implementation.  The hardest programming was getting  the geometry for automatic texture coordinate generation to work.  This method is applicable to any volumetric data set, and it doesn't require too much CPU power either. 

Limitations

While this technique is very fast, it does have some limitations.  First off, there is no lighting.  Lighting would not work because the surfaces inside the volume are not found, and the quadrilaterals are drawn the same regardless of internal surfaces.  Also, the graphics hardware would need to generate surface normals based on the 3D volume, and then it would have to interpolate them, and we didn't have enough time to figure out how to do this and still maintain our desired frame rate.

Another major limitation is that this technique is not as accurate as doing real front-to-back volume rendering.  Our approximation of the volume uses a rectangular array of quadrilaterals, which use interpolation to get their colors.  Artifacts occur when they are viewed at an angle that ephasizes the distance between two planes.

Artifacts: notice the edges that appear along the surface of the CT volume when viewed at this angle.  Those edges do not actually exist in the original volume.

We assumed that the data was already registered in 3D.  We did not make any attempt to register the data.  This is appropriate because this is a technique for visualizing data, not for registering it.

Another limitation is that we used a lot of proprietary ATI extensions.  This program would not work on another graphics card without some intensive redesigning.  It is theoretically possible that any card supporting GL_EXT_texture3D could be used.

References

The OpenGL Reference Manual - The Blue Book.  Addison-Wesley Pub Co; 4th edition (March 2004).  Also available online. 

The OpenGL Programming Guide - The Red Book.  Addison-Wesley Pub Co; 4th edition (November 2003).  Also available online.