Jason Freund, Silicon Graphics, Inc.

Kenneth Sloan ², University of Alabama at Birmingham

One of the challenges of volume rendering is speed. It would be useful to have real-time control to help focus quickly on a desired aspect. While surface representations allow nearly instantaneous positioning of a model at a desired viewpoint, the next level of expectation is to be able to maintain this speed, but produce volume renderings of the entire dataset. Unfortunately, volume data is immense compared to geometric data. A typical geometric representation might require a few hundred objects; a similar volume database may contain millions of voxels. Although volume data is immense, there are many special coherency factors inherent to volume data which make these elevated rendering goals tempting.

Much work has been done to accelerate volume rendering in the image-order (ray-casting) and object-order (back-projection) domains. Early accelerated volume rendering techniques used hierarchical data structures such as K-d trees [12] and octrees [10,11] to quickly traverse empty regions of the volume data, thus reducing the number of samples needed to construct an image.

Later, several papers used hierarchical datastructures to accelerate rendering of homogeneous as well as empty regions [1,5,9,15]. Homogeneity-accelerated ray-casting techniques accelerate a ray through portions of the volume data in exchange for error. The use of a hierarchical datastructure, has several drawbacks. Each method must choose a single, global parameter for their octrees to determine where the error should be distributed. They differ in the choice of this parameter. Arvo [1] uses density, meaning that regions of low density are automatically assumed to be unimportant. Van Gelder, et al [15] use the average deviation from the average intensity value in a node. Laur and Hanrahan [9] store average intensity values in octnodes to control the sampling rate. Danskin and Hanrahan [5] use density ranges to gain a speedup from regions which are relatively uniform. None of the techniques give the user the ability to change an importance classification independently of the global

To circumvent the problems of hierarchical datastructures, later techniques explored other forms of encoding such as a look-aside buffer, flat pyramid or proximity cloud to encode empty space [4,17,16], amorphous "shell" encodings to allow interactive thresholding [13], or run-length encodings to precisely skip empty regions [8]. The success of these techniques over hierarchical methods stems from their encoding scheme where encodings or "skips" are indexed with the same indices used for the volume data. Most recently, the greatest time savings have been made with Lacroutes [8] shear-warp back-projection algorithm which precisely skips empty regions of space.

This paper extends the successful idea of directly encoding empty regions to include homogeneous regions. In Section 3, we show a method for constructing homogeneous regions and a datatype for efficiently encoding them. In Section 4, we show how the direct encoding of homogeneous regions can be used to accelerate volume ray-casting beyond speeds realized by hierarchical schemes or schemes which directly encode only empty regions. In Section 5 we present a new algorithm which extends Lacroutes fast shear-warp back-projection algorithm [8] to use our encoded homogeneous regions.

Proximity clouds [4,16] use a 3-D distance transform [2,3,7,17] to accelerate rays through empty space. A "proximity'' or "skip" field is calculated by a distance transform algorithm which computes the distance from each empty voxel, to the nearest opaque voxel, making each voxel the center of its own empty region. This "skip" field is the distance that may be safely skipped along any ray which samples that voxel. Based on a user-defined transparency threshold, the distance transform calculation is an inexpensive linear time operation, calculated once for the entire volume dataset.

There are several ways to compute an empty space distance transform, and all of them require a two-step initialization phase to compute a distance map. The distance map generated empty space distance transforms record, for each voxel, the distance to the nearest opaque voxel. First, voxels are classified as being transparent or non-transparent. Those which are transparent are called background, and are initialized to the maximum distance in a distance map. The rest are called foreground and are initialized to zero (or the minimum distance). Next, interior voxels are removed by overlaying a simple six-connected neighbor mask at each voxel. If the center voxel and all six neighbors are foreground, the center voxel is changed to background. The result of this initialization phase is a distance map where voxels on the surface of non-transparent regions are initialized to the minimum distance.

After the distance map has been initialized, there are several ways to compute the distance transform. Our implementation uses a fast distance transform calculation presented in [3] which propagates a wavefront of Euclidean distances represented by x and y components, through the volume. A wavefront (represented as a queue) of distances is grown from each foreground voxel until voxels in the distance queue reach the maximum distance or a border. After the queue is exhausted, the Pythagorean theorem is used to compute integer approximations to the true Euclidean distance from the x and y components.

With a slight modification to the initialization phase, any distance transform algorithm can be made to generate distances for homogeneous as well as empty regions. Instead of using an isovalue determination to label a voxel as transparent or non-transparent, each voxel is tagged or colored according to what kind of region it lies in. There are many possible ways of identifying regions. Regions could be identified by thresholding, by using a contouring scheme, using local information such as a gradient operator, or by using any number of graph, split-and-merge, or region growing techniques. Because the choice of homogeneous region identification algorithms is so varied and application-dependent, this choice is left open.

We separate the question of region boundary definition from the volume data. Instead of using a parameter of the volume such as density, opacity, or homogeneity to control error bracketing and sample spacing, our direct encoding scheme allows complete control over how regions are identified. This allows the user to have flexibility in changing parameters. The disadvantage is that the user is responsible for choosing a suitable technique to identify region boundaries.

After voxels are identified according to their region, interior voxels must be removed using a similar kind of neighbor-checking filter to that used by empty space encoding distance transforms. If a voxel has been identified as belonging to the same region as its six neighbors, it is set to background in the corresponding distance map. At the end of this initialization phase, only voxels corresponding to the surfaces of identified regions will be set to the foreground (see Plate 1a). After the homogeneous region initialization phase is complete, any of the empty space distance transform computation algorithms mentioned in the beginning of Section 3.1 can finish the job. The resulting distance map will record the distance from each voxel to the nearest boundary voxel whether that voxel happens to be inside an empty or homogeneous region (see Plate 1b).

Once the volume is pre-processed to find homogeneous regions, skip fields must be stored in an efficient manner. We tested several skip encoding methods for speed and storage efficiency which are all similar in their purpose, but differ in implementation details which lead to different optimizations in the volume traversal portion of rendering algorithms. The most straightforward way to encode skip information is to create a look-aside volume [17] of skip fields which correspond to each voxel. Although this encoding scheme requires a secondary volume of byte-sized skip fields, it does not require extra overhead for indexing because the same indices into the voxel database will index the look-aside volume.

Proximity clouds [4,16] do not require auxiliary storage for encoding skip fields. Proximity clouds store skip fields in place of voxel intensities wherever a voxel contains air. Non-air voxels simply store the original voxel intensity. This duality leads to a minor dilemma. It is not possible for a voxel to contain an intensity less than the maximum possible skip size. Fortunately, this should not cause concern because low voxel densities are usually classified as transparent.

Unlike previous hierarchical encoding schemes mentioned in Section 2, this hybrid datatype requires no auxiliary skip storage, no additional indexing, and no additional memory accesses to fetch skip fields. Datasets consisting of 8 bit samples incur a two fold memory penalty because an extra byte for each voxel

must be used to store the skip field. The maximum skip size of 15 provided by these upper four bits is more than adequate for recording skip fields because our experiments show the tradeoff in speed for skip size falls logarithmically as the speed bottleneck shifts from traversing and sampling the volume to other parts of the rendering algorithm.

Using a ray-caster, the different encoding techniques were compared for their inherent efficiency on three volume datasets. The 3dhead and 3dknee [14] are 256x256x109 MRI datasets of 8 bit intensity values. The hipip [14] dataset is a 64x64x64 electron density map of a high-energy protein molecule. Table 1 shows a comparison of the inherent efficiency of the skip encoding datatypes when tested with our three datasets. Because proximity cloud and look-aside datatypes can only encode empty regions, Table 1 compares our hybrid encoding scheme when only transparent (and not homogeneous) regions have been encoded in the datasets.

Skip fields can easily be used to accelerate ray-casting. Skip fields indicate how large a jump any ray sampling that voxel may take in any direction. The skip is also used as a weight for sampling that voxel so the sample will represent the linear span of voxels the ray traverses. The upper four skip size bits of a hybrid voxel conveniently double as an index into pre-computed higher-weighted portions of the anatomy classification table, giving us sampling weight multiplications for the price of the anatomy classification lookup.

Our homogeneous region ray-caster was tested by rendering a clean image of each of the datasets mentioned in Section 3.3. For each dataset, we classified anatomy by assigning colors, opacities, and a transparency threshold. The wavefront [3] algorithm was used to generate Euclidean skip fields which were then encoded using our "hybrid" datatype. When zero error is required, our method can render a perfect image using the same amount of pre-processing time as other ray-casting algorithms which directly encode empty regions. The time used by our algorithm to render an error-free image is comparable to [16] (refer back to Table 1).

If some image error can be tolerated, then image quality can be traded for a rendering speedup. We tested the amount of time which could be saved during ray-casting against the amount of image error generated by encoding various sized homogeneous regions. Our experiment re-processed each of the three volume datasets 15 times for a range of uniformly relaxed to incrementally more stringent gradient threshold values. A simple gradient threshold was used to identify region boundaries, and skip fields were created using our homogeneous region distance transform, and stored using our hybrid datatype. Danskin and Hanrahan [5], from whom our error metric was taken, noted that a 20% error in image quality was usually the maximum useful limit. The graph in Figure 2 shows the tradeoff between speed and accuracy for our homogeneous region accelerated ray-caster on each of the datasets. All experiments generate 256² pixel images.

Figure 2

The first phase of pre-processing for our homogeneous region scanning algorithm computes and stores skip fields for homogeneous regions, as in Section 3.2 (see Plates 1a,b). Next, the fewest number non-overlapping homogeneous regions, covering the most amount of space, are extracted from the volume data, (see Plate 1c). This is not a precise problem description because we do not need to compute an optimal solution, nor can we do so in a reasonable amount of time (the perfect solution is probably NP-hard). Our imperfect solution favors extracting larger regions rather than smaller ones because fewer large regions can fit inside a volume. It is easy, however, to think of a counterexample where it would not be preferable to extract the largest regions first.

A brute force approach to extracting non-overlapping homogeneous regions makes a sweep through the volume once for each possible skip size (in order from largest to smallest), pulling out regions of that skip size which do not overlap any other region. Certain information associated with extracted regions is appended to a list of regions datastructure. Each voxel which is part of an extracted homogeneous region must be tagged so that other regions will be able to check against that voxel for possible overlap.

The brute force algorithm is effective at finding sets of large, non-overlapping homogeneous regions. Unfortunately, it is very slow because it traverses the entire volume once for each possible skip size, and it must visit every voxel within the shells neighborhood to determine if one region overlaps with another. On a 256 x 256 x 100 MRI volume dataset, the brute force homogeneous region extraction algorithm typically runs from 15 seconds to several minutes on a Pentium Pro 200 computer. Time varies a great deal, depending on both the number and the size of non-overlapping homogeneous regions found [6].

A faster version of the greedy algorithm can be developed based on the observation that the next largest non-overlapping region can be found in logarithmic time [6]. The fast homogeneous region extraction algorithm uses an octree to direct queries to the next largest homogeneous region. After creating the initial octree, this algorithm runs approximately six to eight times faster than the brute force method.

The list of regions datastructure simply stores certain attributes of extracted regions that will later be required for rendering. The color of the region is pre-computed as the average color of all voxels contained in the region. Because it is possible for a region to contain different shades of anatomy, we only need to keep track of the average color (not intensity) value over the entire region, computed as:

where n is the number of samples enclosed in the region and c_i is the pre-multiplied color and opacity of a sample in the region. The only other information about the extracted region that is required is the location of its center and its skip size.

Our accelerated homogeneous region back-projection algorithm operates in two phases. The first phase is identical to Lacroutes fast shear-warp back-projection algorithm [8]. Run-length encoded voxels are projected onto an intermediate image which is then 2-D warped to become the final image. In the second phase, footprints of the extracted homogeneous regions are blended into the final image.

The image produced by warping the intermediate image will be a perfect rendering of the non-transparent, non-occluded voxels, minus extracted homogeneous regions. The last step of the rendering process adds extracted homogeneous regions back into the final image. Our algorithm saves time by drawing 2-D images of the 3-D extracted homogeneous regions. Extracted regions are blended into the final image from the list of regions structure generated during pre-processing. This is done by first projecting each region on to the image plane in order to find out its image space coordinates. Then the average region color is blended into the image according to weights for each position in a re-sampling mask. This weight indicates how deep the region is at a particular point in the mask.

We tested the amount of time which could be saved by our object-order renderer against the amount of image error generated by encoding various sized homogeneous regions. Again, a simple gradient threshold was used to identify region boundaries, and skip fields were created using our homogeneous region distance transform. Our experiment re-processed each of the three volume datasets 15 times for a range of uniformly relaxed to incrementally more stringent gradient threshold values. Our homogeneous region accelerated object-order renderer is capable of pre-processing volume data and rendering perfect images in the same amount of time as Lacroutes algorithm [8].

If, however, image error is acceptable, we were able to cut rendering time nearly in half with less than 4% image error (Figure 4 and Plate 2d). Speed and image quality might be further improved by using more sophisticated region boundary identification schemes, or by focusing on detail in specific regions rather than using a global, uniform gradient threshold as in our experiments. The group of dots in the upper-left corner of the hipip graph in Figure 4 represents attempts to relax region thresholds in order to reduce rendering time. Because voxel intensities from this dataset represent smooth transitions from high to low potential, the dataset demonstrates the worst-case behavior of homogeneous region accelerated renderers.

Shading is a slight problem in the object-order rendering algorithm. Tissue highlights generated by a shading function can get "averaged out" in extracted homogeneous regions. A partial solution to this is to simply use the average pre-shaded and pre-multiplied color and opacity value when extracting homogeneous regions. Fortunately, shading should not pose much of a problem because most shading algorithms affect only tissue surfaces whose true shaded color values can stay intact with the RLE volume.

Error is introduced at three stages during pre-processing and rendering:

• During the two re-sampling steps as the volume is transformed into sheared-object space,
• By approximating the color of a region from its average color,

• By approximating the thickness of a region from any angle using a single pre-computed mask.

The first error source is inherent to any object-order rendering technique which uses the shear-warp factorization. Lacroute reported that error due to resampling is negligible [8]. The second source of error depends on the users region definitions. If homogeneity definitions are strict, then relatively small amounts of error will be incurred. The last source of error is dependent on several factors such as the viewing angle, the distance metric chosen to encode homogeneous regions during pre-processing, and the weights used to compute the depth of each point in the re-sampling mask. We chose the Euclidean distance metric so that regions would look roughly symmetrical from any angle, however any distance metric should work as long as weights are computed properly. If error generated at this stage is of concern, then new region re-sampling masks could be re-computed each time the viewing angle changes.

The homogeneous region accelerated back-projection algorithm only works for parallel projections right now. This is not a limitation of sheared-object space. Lacroute [8] developed a transformation to sheared-object space which allows for perspective projections. In order to allow for perspective projections in our algorithm, it would be necessary to compensate for the scaling effect on our omnidirectional homogeneous regions. Another limitation of our object-order renderer is that it currently does not support the use of the "over" operator. Instead, samples must be collected by blending rgba values because homogeneous regions are composited after scanning. The over operator might be incorporated into our back-projection algorithm if extracted homogeneous regions were to be sorted by depth for each principal viewing axis. This way, scanning could be interrupted each time a region is encountered so the region could be drawn into the image at the correct time.

If no image error can be tolerated, our homogeneous region rendering algorithms can easily be made to fall back on a conventional empty space distance transform to produce a perfect rendering in the same amount of time as previous algorithms. If a small amount of image error is tolerable, and homogeneous regions exist in the volume, a significant additional speedup can be gained with little or no extra cost in pre-processing time. In Section 3.2, standard distance transforms which had previously been used to encode empty space distances, were used to compute distance transforms for homogeneous regions. Direct encoding of homogeneous regions provides flexibility over hierarchical based encoding techniques, giving the user control over importance for each anatomy. In Section 3.3, we showed that during ray-casting, our hybrid datatype is just as efficient as other direct encoding datatypes, while allowing for both empty and homogeneous regions to be directly encoded