Abstract:
In order to solve the need for an accurate two-dimensional representation of a three-dimensional image, the invention provides an improved method and apparatus for processing voxel data which avoid the problem of color artifacts coming from zero opacity regions. The three-dimensional image is stored in a computer memory as a three-dimensional data array having object data values associated with sample points in a three-dimensional space. A central processing unit (CPU) processes the object data values to determine voxel colors C and voxel opacities α for the plurality of sample points and then opacity weights the voxel colors to produce a set of opacity weighted colors  {tilde over (C)} . Next, the voxel opacities α and the opacity weighted colors  {tilde over (C)}  from each sample point are composited in a compositor to form the two-dimensional representation. Finally, a display unit visually displays the two-dimensional representation. In this manner, artifacts resulting from zero opacity areas contributing to the overall color are eliminated because of the opacity weighting performed prior to compositing.

Description:
BACKGROUND OF THE INVENTION 
     The present invention relates to computer image generation and more particularly to methods for generating two-dimensional displays of an image based on three-dimensional image data. 
     Many imaging systems produce data representing physical entities or properties in three-dimensional space. Examples include CAT scans, MRI scans, ultrasound images and seismic test images. The process of converting the three-dimensional arrays of data into two dimensional image displays is referred to as volume-rendering. The term voxel processing is also used, because the 3-D space is ordinarily broken down into small volumes, referred to as voxels. Each voxel has optical properties, ordinarily expressed in terms of opacity (or its inverse, transparency) and color (or grey-scale in black and white imaging systems). These properties contribute to the overall image when it is to be computed by mathematical projection and sampling and converted into pixels to be displayed on a two dimensional display surface. 
     Volume rendering employs a class of algorithms for visualizing or representing three-dimensional grids of data. Such algorithms are based on the physics of light interaction with modeled particles in a volume. The beginnings of a particle model for graphics rendering were developed by Blinn (Bli82) and extended by Kajiya et al. (KH84) to nonhomogeneous densities of particles. The common algebraic technique for combining ray samples is compositing, derived separately by Blinn (Bli82) and Porter et al. (PD84). Modem direct volume rendering algorithms are largely based on Levoy (Lev88), Drebin et al. (DCHB88), Sabella (Sab88), and Upson et al. (UK88). Up-to-date surveys may be found in Max (Max95) and Kaufman (Kau91). 
     Levoy in (Lev88) proposed a volume rendering pipeline that classifies and shades voxels before interpolation. This approach can lead to two computational advantages compared to classifying sample points interpolated from voxel data. First, classification can be treated as a preprocess and need not be repeated per view. Second, the number of voxels in a data set can be significantly less than the number of sample points when computing a view, resulting in lower classification cost. Commonly-assigned U.S. Pat. No. 5,557,711 discloses a multiprocessor and parallel method for volume rendering by preprocessing transparency and luminescence signals in subsets or regions of an object space and then compositing the signals from such regions for the object space as a whole. Unfortunately, classification first followed by interpolation as described in (Lev88, Lev89, Lev90) erroneously interpolates shaded color values independently from opacity. 
     The wide use of volume rendering software from the seminal Levoy references, and the inherent complexity of the underlying model, have meant that implementations of volume rendering can include shading artifacts due to improper classification, shading, and interpolation. Drebin et al. (DCH88) specify opacity weighting of the colors. Others have noted the difficulties associated with opacity and color, notably Wilhelms (WG91) and Blinn (Bli94), but not as it relates to resampling of shaded color values in volume ray tracing. Wilhelms notes that different answers will result if one interpolates colors versus data, but the former may be faster. Blinn discusses the issues in filtering an image that has opacities at each pixel. He shows an example of overlaying and downsampling, and the different answer that results if downsampling is done first without opacity weighting the colors. He calls the opacity weighted colors “associated colors.” Wilhelms, Van Gelder et al. (Wil91,VGK96) discuss the issues in interpolating colors instead of materials, and resulting artifacts, but do not clarify that additional artifacts result if opacity weighted colors are not properly interpolated. 
     Accordingly, a need remains for a better way to process voxel data to interpolate color values in proper relation to opacity values, so as to avoid introducing artifacts. 
     SUMMARY OF THE INVENTION 
     In order to solve the need for an accurate two-dimensional representation of a three-dimensional image, the invention provides an improved method and apparatus for processing voxel data which avoids the problem of color artifacts coming from zero opacity regions. 
     Particularly, according to the invention, a computer system has a processing subsystem, a memory and a display for displaying an image in a two-dimensional array of pixels. The three-dimensional image is stored in the computer memory as a three-dimensional data array comprising object data values associated with a plurality of sample points in a three-dimensional space. A method for generating a pixel color value to be displayed in a subset of the array of pixels as a part of the two-dimensional representation comprises the steps of, first, processing the object data values to determine voxel colors C and voxel opacities α for the plurality of sample points. Second, the voxel colors are opacity weighted to produce a set of opacity weighted colors  {tilde over (C)} . Third, the voxel opacities α and the opacity weighted colors  {tilde over (C)}  from each sample point are composited to form the two-dimensional view. In this manner, artifacts resulting from zero opacity areas contributing to the overall color are eliminated because of the opacity weighting performed prior to compositing. 
     An apparatus for mathematically calculating and visually displaying a two-dimensional view of a three-dimensional image is also provided. In this apparatus, a two-dimensional array of pixels for representing the three-dimensional image is created by generating a pixel color value to be displayed in a subset of the array of pixels. The apparatus comprises a memory for storing the three-dimensional image as a three-dimensional data array comprising object data values associated with a plurality of sample points in a three-dimensional space. A central processing unit (CPU) is provided for processing the object data values to determine voxel colors C and voxel opacities a for the plurality of sample points. The processing unit must also be capable of mathematically calculating a set of opacity weighted colors  {tilde over (C)}  by opacity weighting the voxel colors. A compositor composites the voxel opacities α and the opacity weighted colors  {tilde over (C)}  from each sample point to form the two-dimensional view, and a display unit visually displays the two-dimensional view. 
     The foregoing and other objects, features and advantages of the invention will become more readily apparent from the following detailed description of a preferred embodiment of the invention which proceeds with reference to the accompanying drawings. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1A is a flowchart of a volume rendering pipeline analogous to a prior art system but using opacity weighted colors for interpolation implemented by modified and added steps shown in bold lines. 
     FIG. 1B is a block diagram of a voxel processing subsystem for implementing the present invention. 
     FIGS. 2A and 2B are ray trace diagrams of a ray leaving a fully opaque material into empty space showing sample points and compositing results, differing values underlined for clarity; FIG. 2A showing the result with independent interpolation of color and opacity, as is commonly done in the prior art, and FIG. 2B showing the correct results with opacity weighted interpolation. 
     FIGS. 3A and 3B are black/white/grey-scale renderings of a torus, the view in FIG. 3A including artifacts of separate interpolation of colors according to the prior art and the view of FIG. 3B showing the improved image resulting from using opacity-weighted color interpolation according to the invention. 
     FIGS. 4A,  4 B and  4 C are black/white/grey-scale reproductions of color renderings of sections of human vertebrae collected via spiral CT and rendered using raycasting with preclassification. 
     FIG. 5 is block diagram of voxel processing subsystem similar to FIG. 1B, like numerals representing similar elements, showing an example of an architecture which can be used to implement the opacity weighting method of the present invention in a full color system. 
    
    
     DETAILED DESCRIPTION 
     I. Volume Rendering Pipeline with Opacity Weighted Colors 
     FIG. 1A shows a volume rendering pipeline  10  analogous to the prior art system described by Levoy (Lev88) but with modified and added blocks (shown in bold) to implement opacity weighting of colors according to the present invention. Opacity weighting must occur prior to interpolation for the ray resampling, and hence cannot be corrected for in the compositing stage following interpolation. The volume rendering pipeline  10  is arranged to avoid this difficulty, as well as the problem further outlined below. 
     The volume rendering pipeline  10  begins, as in the prior art, with a set  12  of acquired values f(x) for an object space, and the data is processed in a data preparation step  13  to produce a set  14  of prepared values f(x). The prepared values are then split by shading and classification steps  15 ,  17  into a set  16  of voxel colors C(x) and a set  18  of opacities α. In the prior art, these values would be processed separately through the two processing branches before being combined in a final compositing step. According to the invention, however, the voxel opacities are fed back to the color processing branch to be used in an opacity weighting step  19 , to produce a set  20  of opacity weighted colors  {tilde over (C)}  (x). Then, while a conventional ray-tracing step  21  proceeds on the opacity branch of the pipeline to produce a set  23  of sample opacities α, a ray tracing/resampling step  22  proceeds on the color branch, using the opacity weighted colors  {tilde over (C)}  (x). This step produces a set  24  of sample opacity weighted colors  {tilde over (C)}  (u). A final compositing step  25  composites the sample opacities α and sample opacity weighted colors  {tilde over (C)}  (u) to produce a set  26  of image pixels. 
     FIG. 1B shows an architecture of a voxel processing circuit  30  in which the foregoing pipeline can be implemented. The circuit  30  is designed to be used in a conventional workstation or personal computer, plugged into the main bus to communicate with the computer&#39;s central processing unit (CPU), random access and disk memories, and display subsystem and monitor. Accordingly, the circuit  30  has a bus interface  32  which communicates to the computer&#39;s main bus to receive voxel block data (object space address, opacity and unweighted color data for each sample point in object space, which can be in a compressed format) and ray data (vectors defining rays through the object space). The voxel block data is received by a decompression engine  34 , decompressed, and then forwarded to a raw voxel cache or memory  36 . The voxel data is then passed to a classifier/shader  38  which, using a classification RAM  40  as a scratch memory, processes the raw voxel data for each sample point into an opacity α and a color value C, which is a single quantity in a black/white/grey-scale system, or three values (C r , C g , C b ) in a full color system (see FIG.  5 ). The classified voxel data are stored in a classified voxel cache  42  which includes a cache for the opacity data and, for a full color system, three caches for the three sets of color data (see FIG.  5 ). These data are then held pending selection of a view of the data by the user of the computer system. 
     Viewing software running on the CPU generates rays defining a direction for viewing through the object space responsive to a user inputting a selected viewing direction. These rays are transmitted via the bus interface  32  to a ray sequencer  44  which controls operation of an interpolator  46 . The interpolator  46  comprises an interpolation engine for processing the opacity data (α) and, for a full color system, three interpolation engines for processing the three sets of color data (C r , C g , C b ) (see FIG.  5 ). Each interpolation engine operates to resample the voxel data along the rays, interpolating between the sample points as needed to determine the most likely opacity and color values at each sample point along the rays. Associated with each interpolation engine is a compositor  48  which composites the resampled data along the rays, preferably in a back-to-front manner. The compositors produce a set of pixel data defining a two-dimensional view of the object space as seen along the rays. This pixel data is stored in an accumulation buffer  50  for output to the frame buffer of the computer display subsystem. 
     Section II discusses the root of the problem and Section III describes the present invention in further detail. Sections IV and V show comparison results and conclusions. We demonstrate the problem encountered by separate interpolation of color and opacity in the context of ray cast volume rendering. The suitability of the interpolation method presented is not limited to ray-casting, however, and applies to other techniques such as splatting (Wes90) and three-dimensional texture mapping (VGK96) as well. 
     II. Problem Example 
     We follow Blinrn&#39;s (Bli94) notation for an opacity-weighted color, which he calls an associated color, and denotes as  {tilde over (C)} =αC, where C is simply the nonweighted color. The compositing equations for back-to-front compositing are: 
     
       
           {tilde over (C)}new =(1−α front )  {tilde over (C)}back + {tilde over (C)}front   (1) 
       
     
     
       
         α new =(1−α front )α back +α front   (2) 
       
     
     If we were to use nonweighted colors, the corresponding equation would be (opacity is calculated the same as eq. 2): 
     
       
           {tilde over (C)}new =(1−α front )C back α back +C front α front   (3) 
       
     
     To illustrate the difficulties encountered in performing separate interpolation of opacity and color values, compare FIGS. 2A and 2B, showing sample points and compositing results for the prior art and the present invention in the case of a ray leaving a fully opaque material into empty space. FIG. 2A shows the result with independent interpolation of color, and opacity, as is conventionally done. This approach incorrectly allows colors associated with an empty region of space to contribute to the ray, thereby resulting in unwanted artifacts. FIG. 2B shows the correct results using opacity weighted interpolation. Differing values are underlined for clarity. 
     First, consider FIG. 2A, which assumes the back-to-front compositing formula given as in equation 3. In FIG. 2A, the ray is shown exiting a fully opaque region, into a region of empty space. Voxels are classified and shaded before interpolation. The difficulty arises in the middle sample point which is positioned between opaque and transparent samples. The conventional method yields a color value of C=0.5 since it averages the material and empty space colors. Since the opacity value is 0.5, this sample does contribute color to the ray. However, the color that it contributes is incorrect since it has been partially affected (50%) by the color that the empty space takes on. The color assignment to regions with zero opacity should not contribute to the final image. This artifact arises because the interpolation that is conventionally used is not a material-based interpolation. The resultant ray color (or grey-scale value)  {tilde over (C)}r , as affected by zero opacity regions, is shown in the lower left hand corner of FIG.  2 A. 
     FIG. 2B shows the interpolation of opacity weighted colors, as discussed further in the next section. Note that the resultant ray color is not affected by the color assignment in empty space. This interpolation method essentially treats opacity as a material density and weights the color contributions of each voxel by this density. From a particle model perspective (KH84), the greater number of particles from a high density material influences color more than a lower density material&#39;s contribution. The resultant ray capacity  {tilde over (C)}r , unaffected by the zero opacity regions, is shown in the lower left hand corner. 
     III. Opacity Weighted Interpolation of Colors 
     Proper calculation of resampled colors is done by weighting each source voxel color value, (C r , C g , C b ) by the opacity α. These weighted colors are then interpolated, as are the opacity values. A ray sample opacity and color are calculated using a weighting factor w as:              α   =       ∑   i            w   i          α   i                 (   4   )                 C   ~     =         ∑   i            w   i          α   i          C   i         =       ∑   i            w   i            c   ~     i                   (   5   )                                
     The compositing is then performed as above (equations 1, 2), using these opacity weighted colors at each sample point. At each source voxel, we choose to store opacity α and unweighted color C, and perform the opacity weighting during interpolation. For the case where α is 0, the products with (C r , C g , C b ) are zero, therefore, the colors can be set to 0. Since the products w i α i  are used for each channel,they can be computed once as ω i =w i α i , resulting in:              α   =       ∑   i          ω   i               (   6   )                 C   ~     =       ∑   i            ω   i          C   i                 (   7   )                                
     This new method uses efficient calculation, with minimal storage requirements, and computes correct ray colors. In fact, compared to separate interpolation of color and opacity, the new method involves one less multiplication per sample point along a ray, per color channel, due to the simpler compositing equation (1) versus (3). 
     There are three alternative approaches to avoiding the interpolation artifacts described in Section II. These approaches are summarized and compared in the following table, Table 1. 
     
       
         
               
             
               
               
               
               
             
               
               
               
               
             
           
               
                 TABLE 1 
               
             
             
               
                   
               
               
                 Volume rendering approaches that avoid interpolation artifact. 
               
               
                 Volume Rendering Approaches 
               
             
          
           
               
                   
                 Method 
                 Advantages 
                 Disadvantages 
               
               
                   
                   
               
             
          
           
               
                 (1) 
                 Interpolate then 
                 Sharper surface detail 
                 Classification: 
               
               
                   
                 classify 
                 as compared to (2), (3) 
                 *view dependent 
               
               
                   
                   
                   
                 *needed per view 
               
               
                 (2) 
                 Store opacity 
                 Shorter render 
                 Requires more 
               
               
                   
                 weighted colors 
                 time compared to (1) 
                 bits/voxel than (3) 
               
               
                 (3) 
                 Store Unweighted 
                 Short render time 
                 Softer imagery 
               
               
                   
                 Colors, α weighted 
                 less bits/voxel 
                 compared to (1) 
               
               
                   
                 interpolation 
               
               
                   
               
             
          
         
       
     
     The first alternative is to store raw data values themselves, interpolate these down to ray samples, and then classify the interpolated results. This is sometimes referred to as post classification. Opacity weighting problems are avoided by this approach. In addition, image quality is often superior since sub-voxel surfaces are reconstructed by this approach. There are, however, certain drawbacks. First, classification must be performed for every frame being computed, unlike the preclassification case where classification is only performed when the classification function changes. Computational cost is also typically higher since there are often more ray sample points than voxels involved in a rendering. Second, since classification acts on sample points whose positions are view dependent, the classification is itself view dependent. This approach can introduce view-dependent artifacts. 
     The second alternative is to store classified-shaded-color values preweighted by opacity, along with opacity at each voxel, as was done in Drebin et al. (DCH88). Interpolation is then performed on preweighted colors. Preweighted colors, however, require more bits of precision than unweighted colors for the same image quality. This is because opacity weighting returns a fraction of the shaded values, which in fixed point will reduce the accuracy of subsequent calculations. In the proposed alternative a higher precision representation can be used during weighting which results in greater accuracy. 
     The third alternative is opacity weighted color interpolation as described herein. Color and opacity values are stored independently. Opacity weighting is performed during interpolation for memory savings (bit precision), computational advantages of preclassification, and comparable speed to preweighting colors. Through examination of source code and inquiry we have determined that Lacroute et al. (LL94) perform weighting when loading scanlines for a similar solution. Three-dimensional texture mapping using hardware acceleration (VGK96) does not opacity weight colors, although it could be done by preweighting. Unfortunately, that may result in difficulties due to the fixed point texture storage formats. The next section demonstrates the artifacts incurred by interpolation of unweighted colors. 
     IV. Comparative Results 
     As illustrated in FIG. 2A, the interpolation of unweighted colors miscalculates the shading of new sample points along the rays. In volume rendering, such errors can result in color shifting, darkening, greater aliasing, and contribution of color from null regions. We show two examples of artifacts resulting from this error. 
     FIGS. 3A and 3B show two volume renderings of a volumetric torus. The data set is a sampled analytic torus function that has been low-pass filtered to 8 bit voxels for antialiasing. In this example, color consists of grey-scale shading. This volume is then classified before rendering using a standard isosurface classifier (Lev89). The rendering process itself uses a step size of one-half of the voxel spacing. The rendering of FIG. 3A, using separate interpolation of colors and opacity, causes severe banding artifacts since the colors assigned to empty space contribute to the ray with this scheme. The rendering using opacity weighted color interpolation, FIG. 3B, does not show these artifacts. Although voxelization of the torus is still somewhat apparent, this is due to the limited voxel resolution (256×256×128) as opposed to any interpolation artifacts. 
     The second example, shown in FIGS. 4A,  4 B and  4 C, compares interpolation strategies on a rendering of a computed tomography data set. This 8 bit data set (100×96×249 voxels) includes several human vertebrae and was collected with spiral CT (data is courtesy of Dr. Ramani Pichumani, Stanford University). Classification was also performed as a preprocess. Empty space (opacity=0.0) surrounding the thresholded vertebrae was colored red. This coloring readily labels surface regions where inappropriate colors are being introduced. The data was rendered using volumetric ray-casting with a step size of one half of the voxel spacing. 
     FIG. 4A shows the results of separate interpolation of color and opacity using the prior art method. In the color original, this image has a pinkish cast (shown as a grey tone in FIG. 4A compared to FIG. 4B) imparted erroneously by the prior art method from the transparent surrounding space, which was colored red to demonstrate the anomalous behavior of the prior art method. FIG. 4B shows opacity weighted interpolation of colors according to the present invention, using the same data set as for FIG. 4A, including the red coloration in the transparent space surrounding the vertebra. FIG. 4C shows a normalized difference image between FIGS. 4A and 4B. This view shows the extent of image data erroneously imparted into FIG. 4A from the transparent but red colored surrounding space. Note that the worst artifacts appear (red or grey areas in the image of FIG. 4A) when rays are able to graze along the surface for some distance, leading to a result that is dependent on surface orientation relative to the viewer. Animation also shows that the artifact causes banding and is view dependent. 
     V. Conclusions 
     The foregoing comparisons show what appears to be a common artifact in volume rendering software and algorithms. The origins of the artifact arise from an interpolation method that has been widely described in the literature and used, but not previously recognized to be improper. This artifact occurs when classifying and shading before resampling, and the present invention provides a solution for this case: opacity weighted color interpolation. Prior solutions have either interpolated ray values, and then classified and shaded, or stored entire opacity weighted volumes, which requires more storage precision. We have found that the error magnitudes are dependent on data set, view, classification, and shading, which may account for the oversight. Also, because of the number of steps in volume rendering, other artifacts may dominate, such as reconstruction error, gradient calculation, and shading, thereby reducing the perceived importance of opacity weighting. While the effects of misuse of opacity and color in volume rendering can be subtle, this invention clarifies an important nuance in a widely used algorithm, and will improve the quality and understanding of future implementations. 
     Having described the principles of the invention in a preferred embodiment thereof, it should be apparent that the invention can be modified in arrangement and detail without departing from such principles. We therefore claim all modifications and variations coming within the spirit and scope of the following claims. 
     VI. References 
     (Bli82) J. Blinn, Light reflection functions for simulations of clouds and dusty surfaces. In Proc. of SIGGRAPH, pages 21-29, July 1982. 
     (Bli89) J. Blinn, Return of the Jaggy. IEEE Computer Graphics and Applications, pages 82-89, March 1989. 
     (Bli94) J. Blinn, Compositing, part I: Theory. IEEE Computer Graphics and Applications, pages 83-87, September 1994. 
     (DCH88) R. A. Drebin, L. Carpenter, and P. Hanrahan, Volume Rendering In Proc. of SIGGRAPH, pages 65-74, August 1988. 
     (Kau91) A. Kaufman, editor, Volume Visualization. IEEE Computer Society Press, Los Alamitos, Calif., 1991. 
     (KH84) J. T. Kajiya and B. Von Herzen, Ray tracing volume densities, In Proc. of SIGGRAPH, pages 165-174, July 1984. 
     (Lev88) M. Levoy, Display of surfaces from volume data, IEEE Computer Graphics and Applications, 8(5):29-37, May 1988. 
     (Lev89) M. Levoy. Display of Surfaces From Volume Data. Ph.D. thesis, University of North Carolina at Chapel Hill, 1989. 
     (Lev90) M. Levoy, Efficient ray tracing of volume data, ACM Transactions on Graphics, 9(3):245-261, July 1990. 
     (LL94) P. Lacroute and M. Levoy, Fast volume rendering using a shear-warp factorization of the viewing transformation, In Proc. of SIGGRAPH, pages 451-458, Orlando, Fla., July 1994. 
     (Max95) N. Max, Optical models for direct volume rendering, IEEE Transactions on Visualization and Computer Graphics, 1(2):99-108, June 1995. 
     (PD84) T. Porter and T. Duff., Compositing digital images, In Proc. of SIGGRAPH, pages 253-259, August 1984. 
     (Sab88) P. Sabella, A rendering algorithm for visualizing 3d scalar fields, In Proc. of SIGGRAPH, pages 51-58, August 1988. 
     (UK88) C. Upson and M. Keeler, V-buffer: Visible volume rendering, In Proc. of SIGGRAPH, volume 22, pages 59-64, July 1988. 
     (VGK96) A. Van Gelder and K. Kim, Direct volume rendering with shading via 3D textures, In Symposium on Volume Visualization, San Francisco, Calif., Oct. 1996. 
     (Wes90) L. Westover: Footprint Evaluation for Volume Rendering, in Proc. of SIGGRAPH, volume 24, no. 4, pages 367-376, July 1990. 
     (WG91) J. Wilhelms and A. Van Gelder, A coherent projection approach for direct volume rendering, In Proc. of SIGGRAPH, pages 275-284, 1991. 
     (Wil91) J. Wilhelms, Decisions in Volume Rendering, In SIGGRAPH 91 Course Notes 8-State of the Art in Volume Visualization, pages I.1-I.11, 1991.