Abstract: High-speed and effective algorithm of volume rendering is proposed. New methods of object representation by using unified freeforms and arbitrary algebraic surfaces are suggested as well. Proposed algorithm allows to lower complexity of computations, memory and bandwidth requirements in comparison with known methods. The prototyping of hardware system allows operating in real-time with high-quality images. Proposed systems are targeted for scientific visualization, medical applications such as virtual bronchoscopy, endoscopy, and tomography, fly through a virtual patient; virtual reality-based systems and high fidelity conventional image generators. Primary goal of the project is the development of algorithm and 3D data presentations for fast volume-rendering, evaluation of developed algorithms in order to find optimal hardware architecture of real-time image generator. 1. Introduction Virtual reality systems, where viewer merges into the model world, require high quality visualization, for example, on the level provided by Onyx2 Reality Monster by Silicon Graphics, priced around 100 thousand dollars. In addition there is a number of application area specific problems, for example, dynamic distortion correction, which is required for image projection to non-planar screens taking into account viewer and projection device *movements. These problems are solved in high-end real-time visualization systems ESIG and HARMONY of Evans & Sutherland). More recently, voxel-based methods have been developed which create a 3D view directly from the volume data. For the first time on PC-class computers, the vg500 of Mitsubishi Electric allows high quality, real time volume rendering (30 frames/sec.). The vg500 computes 500 million voxels per second with associated 3D gradient estimations. The surface-based methods first create on intermediate surface representation of the object to be shown. Both surface-and voxel-based methods have their merits; the decision which one should be used for a particular application depends on the visualization goals. We are present a new project on volume modeling and rendering that combines volume representations by voxel data and by real continuous functions. Image generators traditionally use polygons as database primitives [1]. It is difficult to use polygons to make convincing models of clouds, smoke, trees and anything else for which no simple surface representation exists. Research system that uses voxels as the database primitive has the complementary behavior. An image generator architecture that could process both polygons and voxels in an elegant way would have the advantages of both approaches. We suggest expanding the notion of primitives and making it possible to process them by easy and effective method without approximation by polygons. The analysis of possible directions of evolution of a real-time visualization systems shows that the easiest way to improve picture quality, i.e. to increase number of polygons rendered per frame, is not the most effective one. Going this way, the qualitative changes can be hardly achieved. For instance, introduction of texture was such a change. Even displaying much more polygons, an image without textures will be poorer than an image with fewer polygons but with color textures. Moreover, since Phong shading renders the higher quality with a much smaller number of polygons, one needs much less computing power for geometry processing; consequently, the overall cost of the Phong-based system is significantly less than that of a comparable Gouraud-based system. Thus, changing polygon "quality" holds more promise for evolution of visualization systems rather than increasing polygon quantity. Nevertheless imaging realistic mountains and terrain in general requires big number of polygons [2]. The same is with sculptural surfaces. Even if NURBS (non-uniform rational B-splines) or Bezier patches are applied, a comparatively big number of them is required. All this proves need for processing a lot of primitives in visualization systems. Existing algorithms work either with polygonal or with polynomial representations of geometry. Surfaces of high order are approximated by polygons or patches and these approximations again require lots of primitives. On the other side using polynoms of high order is not convenient due to complexity of computations and accuracy limitations. In this article a technique for rendering free form surfaces and volumes is proposed, which provides quality impossible with polygonal representation even with very big number of polygons. From a practical point of view an versatile system providing solution of wide range of tasks and minimal replacement of modules for application area specific tasks is the most efficient. Such versatility is one of requirements to system of new generation - Voxel Volumes (VxVl), which is described below. 2. Volume oriented visualisation There are five basic techniques of volume visualization 2.1. Polygonal rasterization Commercially available hardware for mapping 3D texture to polygons allows one to visualize (though with limited efficiency and quality) 3D volume data. Image of 3D data is presented by layered billboard polygons. Silicon Graphics has optimized OpenGL Volumizer for all current Silicon Graphics workstations so application developers will be able to write their application once and achieve maximum performance on all systems. OpenGL Volumizer is built upon the concept of the voxel as the fundamental small-scale volumetric primitive, and the tetrahedron- a four-sided three-dimentional solid object- as the fundamental large-scale volumetric primitive. 2.2. Object-space rendering or forward viewing methods In Maspar MP-1 [3] è CM-2 systems [4] a technique of multiple transformations was implemented, the earliest implementation of splatting algorithm was done on network of several Sun workstations [5], analogous algorithm was implemented in Pixel Planes 5 [6],[7] where particular screen area is assigned to a rendering processor. Efficient shift algorithm was implemented for shared memory architecture in Challenge system of Silicon Graphics [8], in which synchronization time is minimized by dynamic overload balancing and such task subdivision that decreases number of events requiring synchronization. The algorithm is based on support interprocessor connection with low wait time and on hardware caching for wait time minimization. Forward viewing rendering is implemented in Voxel Processor, Voxel Scope II [9], TMC CM-5 [10] etc. 2.3. Backward viewing methods Each node of multiprocessor computer DASH [11] is assigned to particular part of screen and implements ray casting algorithm, loading non-local data is done by efficient caching subsystem. Ray tracer implementation on Intel iPSC/2 [12] used static scheme of data distribution. Ray tracer also was implemented on MPP-system Fujitsu AP1000 MIMD [13]. Three coherent algorithms were designed for Cray T3D [14]. Some architectures with shared distributed memory implement CellFlow technique [15], [16]. 2.4. Hybrid methods Method of inverse incremental transform to object space was implemented on IBM Power Visualization System (PVS) [17]. 2.5. FVR-Fourier Volume Rendering This method uses not volumetric object description but their frequency domain representation. It was developed in Hewlett-Packard [18]. 3. Voxel-Volumes We developed new generation visualization system for volume rendering. The system provides high level of versatility and high image quality. It is based on: a) New surface representation by free-form surfaces without using polynoms or polygonal/patch approximation and representation volume areas by arrays b) Efficient algorithms of rasterization 4. New surface representation Proposed surface representation uses quadrics-surfaces of second order as basic primitives. Quadrics are foundation for free-form surfaces. 4.1. Implicit shape driving functions. It is possible to describe complex geometry forms by specifying surface deviation function (of second order) in addition to surface base function of second order. Generally a function F(x,y,z) specifies surface of second order – quadric: F(x,y,z) =Ax2 + By2 + Cz2 + Dxy + Exz + Fyz + Gx + Hy + Iz + K. Free form surface is a composition of surface base function and surface deviation (deformation) function: F’(x,y,z) = F(x,y,z) + R(x,y,z), where deviation R(x,y,z) is defined:
Fig. 1.
Surfaces produced by composition of base and deviation function are smooth. Often only few such compositions are needed to construct complex surfaces (Fig. 1. Scene defined by 45 quadrics). 4.2. Scalar shape driving functions The most novel real-time visual systems have terrain-skinning processors, which guarantee continuous level of detail (LOD), change with respect to surface roughness and viewpoint distance. Many terrain-skinning systems use regular terrain elevation grid with square cells. It leads to algorithmic simplicity of computations, database uniformity, and strict definition of relationship between adjacent LOD, which results in database generation simplification. Terrain skins constructed from regular grids are efficient to generate, prioritize, query (for culture conforming) and store. However, features with high spatial frequency context (ridgelines and canyons) require large numbers of polygons to meet a specified level of terrain accuracy. Non-regular terrain algorithms are more complex but have the potential to reduce the number of polygons that the system must process. However, for many terrains with limited elevation gradient, the average expected reduction in the image generator load is small for irregular grid as compared to regular grid. Each irregular grid node consumes much more computational time and memory than that of regular. Systems, which employ irregular grids, are very limited in the number of LODs that can be handled. And it is very difficult to solve problems concerned with terrain deformation when using irregular grid (explosions, pits in the ground).
Fig. 2.
The dream of every simulator designer is to render a terrain as easily as a texture. The algorithm considered below does this (Fig. 2. Mountain landscape without feature texture generated for altitude/height map 64x64). A terrain model is coded as differential height map, i.e. the carrier surface is defined by algebraic means and only deviation from this basic surface is stored in the each node. Such a modeling method simplifies creation of smooth detail levels and shading. The data of height grid is not subject to geometry transformations as the triangle vertices do. The geometry transformations are only required for the carrier surface. During the recursive voxel subdivision on each level, we project the centers of the voxels onto some plane. The computed coordinates, as well as in the case of ordinary RGB texture map, will define address in the so-called “altitude map” or “shape texture” [19]. We calculate the altitude corresponding to this address and a level of details, and use it to modify coefficients of the plane or quadric equation. As a result we will obtain a smooth surface of arbitrary shape modulated with the values from the altitude map. But the problems solved by this algorithm require much more complicated methods within the traditional approach. Indeed, the common way to present terrain with polygons requires an abundance of polygons. Besides, the number of additional problems arises such as high depth complexity, hidden polygons removal, priorities, switching between levels of detail, clipping polygons by the pyramid of vision, etc. Such problems do not appear in the proposed method. The Geometry Processor works with the single plane. The corresponding traversing of the tree and the set of masks provide the right priority order. The backside of a terrain is rejected automatically. The clipping terrain by the pyramid of vision becomes unnecessary since sampling of just required altitudes from the altitude map is provided automatically by the rendering algorithm. To switch between levels of detail the same procedure is used as for the ordinary texture. 5. Definition of volume areas by scalar arrays To visualize semi-transparent structures with internal density distribution, 3D textures, boundary and internal voxels are used during rendering. It is possible to render such structures when they are bounded by non-trivial surfaces. Light penetration through materials and semi-transparency is modeled using color blending. 5.1. Computation of pixel color. In ray casting approach secondary (scattered and refracted) rays are not traced to increase computation speed. In this aspect ray casting differs from ray tracing. We also do not trace secondary rays when modeling light penetration through semi-transparent environments. Reflected rays and attenuation is taken into account in our lighting model. Computation of pixel color is done by formula [20]:
,
is final pixel color, 
is either r, g or b (i.e. red, green or blue), 
is intensity computed in n-th voxel by Phong model, 
is the transparence (opacity) of n-th voxel. 
is reflected light at the first point on a ray, 
is the background color and 
. Accumulation threshold test: if on k-th step summed transparency 
is less than 
, then all voxels after k can be discarded.
5.2. CSG set operations for volumes
In case when there are several semi-transparent objects composed as union an object priority for overlapping is introduced. Union of several objects when at least one of them is semi-transparent differs from union of opaque objects (Fig. 3. The objects with different priorities overlapping). If the first object of two has greater priority then for union operation it is drawn over the second object. When more than one object is found in voxel then voxel is assigned attributes (color, transparency) of the object with higher priority.
Fig. 3.
5.3. Considering perspective for transparency computation Described algorithm renders only parts of objects which are inside cube, in [21] an algorithm was implemented, that applies perspective transformation, which generalizes described algorithm for pyramidal volumes. It allows one to synthesize images with perspective. To take into account perspective it is necessary to make correction to algorithm of pixel color computation. Voxel sizes after geometry transformation of primitives depend on Z coordinate. The more Z is the more voxel is, so when color is computed voxel transparency is recomputed taking into account change of voxel depth. 5.4. 3D Texture When surfaces are semi-transparent free-form surfaces are used to bound 3D texture (Fig.4. Semi-transparent lid defined by 4 quadrics and filled by 3D-texure). In general 3D texture in this case can be considered as separate object returning in all space positive result when requested about intersection. Mapping texture to object can be seen as intersection of the object with texture object, color and transparencies are taken not from object attributes but from texture map.
Fig.4.
Texture is generated at visualization time using three-dimensional density map. Transparency of elements is computed depending on density value and normal to iso-surface (Fig.4.Cylinder with 2 oval shells, normal is needed in lighting model) is computed as density gradient [20]. For applications with high frame rate (real-time visualization systems) efficient antialiasing is required. To accomplish this 3D variant of mipmapping technique should be used as well as quadrolinear interpolation for texture filtering. Important property of algorithm should be an ability of fast rendering anisotropic and pseudo-anisotropic objects for which it is not necessary to scan volume down to the last recursion level, but anisotropic areas should be skipped along z-coordinate, colors and transparency are computed immediately - visualization of atmospheric effects – 3D-clouds, smoke, fog (Fig. 5).
Fig. 5.
6. New rasterization technique Along with possibility to rasterize 2D space, the main feature of the Voxel Volumes visualization system is rasterizing made by 3D-space quadtree subdivision of pyramids of different levels, which constitute the whole pyramid of vision. Then a pyramid of the lowest level is binary-tree subdivided into voxels of the lowest level - Recursive Multilevel Ray Casting (RMLRC). In the latter case extent space regions (in depth) can be masked out. The technique of RMLRC allows determining an intersection of a ray (pyramid) of any level with a surface effectively. It is also suitable for fast culling of a spatial region outside an object. The core of this approach is effective search of volume elements (hereinafter voxels), involved in current frame generation, fused with direct projection. In other words, we traverse scene space (unit cube) by quad-tree, leafs of which (slices) is a roots of the binary subtrees. Analyzing mutual location of i-level cell Vi and object, we can differentiate the following cases: Vi is outside the object. All screen space corresponding to Vi is omitted from further consideration. Vi is inside the object, that is, all Vi voxels belong to the object. Vi is called “intersected” and is subdivided into cells of (i+1) level. New cells should be analyzed further. Identification of “inner” cells on the each step of subdivision can be used for masking or filling 3D texture.
Fig. 6.
6.1. Transformation of coefficients to local coordinates of cell Vi during subdivision of Vi-1 Quadric object is base for building all other objects. Its surface is zero level surface of quadric three-dimensional function, i.e. function of its surface is defined in non-explicit form using ten quotients (A, B, C, … K): Q(x,y,z)= Ax2 + By2 + Cz2 + Dxy + Exz + Fyz + Gx + Hy + Iz + K = 0, where x, y è z are spatial variables. It is possible to write this equation in matrix form:
The space in which the algorithm subdivides cubic volume is called further work or model space. M will denote this space further. Now introduce coordinate system immediately related to real space metrics (viewer or camera coordinate system). Four planes and depth of view field define viewing frustum in this space. Frustum volume is called pyramidal volume and denoted by P.
If geometry transform is represented by matrix C then new, transformed, matrix of quotients Q’ = CT * Q * C is applied in space M according to matrix Q in space P.
In particular on each subdivision step along quadtree, when scaling by 2 and shift by ±1 along two coordinate axes is performed, recursive transformation of quadric quotients looks like:
A’=A/4, B’=B/4, C’=C, D’=D/4, E’=E/2, F’=F/2,G’=G/2+i*A/2+j*D/4,H’=H/2+i*D/4+j*B/2
I’=I+i*E/2+j*F/2,K’=K+i*G/4+j*H/4,K’’=K+i*G’/2+j*H’/2 Quotients without ’ are taken from previous recursion step; variables i, j = ±1 are defined in algorithm depending on recursion direction (i - along x-axis, j – along y-axis). Obtained quotients are used in intersection test. 6.2. Criteria of zero existence at the d -area of the coordinate origin Consider an intersection of some object defined as f(x,y,z) ³ 0 with a subdivision cell V(P, d )={Xi|Xi-P|£ d }. It is evident that an accurate solution of this problem (set of inequalities) is possible only for sufficiently simple functions f(x,y,z). Besides, there is no ultimate need in the accurate solution. As though the certain cases (a cell is totally overlapped or totally clean) are solved by relatively facile formulas, it will be enough to determine whether a point of a contour belongs to the area. Definition: a set of points V(P, d ) = {X:dist(X,P) £ d } is called an area of point P in the metric space. All further consideration will be made for 3D space and the metrics dist(X,Y) = max{|X|,|Y|}. Let f(x,y,z) be an analytical function. Then the function f(x,y,z) can be expanded in Taylor series at V((0,0,0), d ).
,
rewritten equation 1
,where
,
condition:
(1)
and 
if 
then 
,
because F(e =0) = F0 = | f(0,0,0) | > 0 (2)
F(e ) decreases when
, besides Fd > 0. Since F(e ) > 0 for 0 < e < d d , we obtain sufficient condition of not having zeros in V((0,0,0),e ) for f(x,y,z). Since F(e ) > 0 when e = 0 it is sufficient to test sign of F(e ,e =d ). If following condition is true then f(x,y,z) has zeros in V((0,0,0), d ):
Û
(3)
In case of unite area, V((0,0,0), d =1), we obtain condition: “module of free member compares to sum of modules of other coefficients”.
6.3. Determination of Vi cell position
Intersection test for base quadric is performed as follows:If K+R<0 and (|A|+|B|+|C|+|D|+|E|+|F| +|G|+|H|+|I| < -(K+R)) then slab is outside;
or C*M=P, where C is the transformation matrix, (M) are homogenous coordinates of point in M space, and P are corresponding coordinates in P space.
Many military simulator applications require the generated image to be free of distortion regardless of where in some allowable volume the pilot’s eyepoint lies. Since the distortion mapping is in general different for each eyepoint, some means of locating the eyepoint is necessary. A head-tracking device fitted to the pilot’s helmet is the usual solution. A more difficult problem lies in determining what the distortion mapping looks like from each viewpoint and inverse mapping. To perform correction for changing eyepoint in a spherical dome, extra computations are required.
This approach takes into account dynamic correction of projection system distortion.
Viewing frustums, which form distorted spherical object space, intersects orbits around viewer (360 degrees in all directions).
Simulators require efficient processing large databases describing modeled environment in real time. Algorithms that solve correctly tasks in this application area can be defined as algorithms of visualization of open database. Term open database means that whole data base does not fit in view field. New and promising arm systems make principally new requirements to computer image generators (CIGs) used in simulators for environment synthesis [1]. Most complicated tasks are low height flights following terrain, officer staff training, helicopter pilot training etc. Design of new CIGs require solution of such problems as:
- Improving realism of synthesized images due to increasing system performance and due to implementing new visual effects.
- Imaging large terrain areas (400õ400 km and more).
- Projecting images on spherical screen of large size is challenging and complex task. Its solution is required primarily for army simulator applications. Complexity and volume of computations increase dramatically for moving viewer and projection device.
- Number of surfaces for rendering non-planar objects is greatly decreased. Defining volumes by free-form surfaces allows one to compress data base size 100 and more times as compared with polygonal representation.
- Distortion correction is greatly facilitated. Let us note that there is a majority of non-planar surfaces in visualization system implementing distortion correction, since after distortion of planar surfaces on screen they become non-planar. Yet non-planar surfaces can hardly become planar ones after distortion on the screen.
- Geometry processor load is decreased greatly as well as data flow from geometry to video processor.
- Voxel-based rendering solves almost all problems pertaining to photorealistic imaging of terrain by elevation grid without preliminary triangulation [24-28].
- Time of landscape processing and rendering does not depend on elevation map resolution.
- Simplicity of surface animation and deformation.
- Possibility to visualize volumes.
- Wide range of possible application areas.
11. References
- Michael A.Cosman, Allen E.Mathisen, John A.Robinson. Evans & Sutherland. "A New Visual System To Support Advanced Requirements". Presented at the IMAGE V Conference Phoenics, 19-22 June 1990.
- Michael A.Cosman. Evans & Sutherland. "Mission Rehearsal Modeling". Presented at the IMAGE VI Conference Scottsdale, Arizona 14-17 1992.
- Vezina,G., Fletcher, P., and Robertson, P. K., "Volume Rendering on the MasPar MP-1", Proceeding of 1992 Workshop on Volume Visualization, Boston, MA, 3-8, October 1992.
- Schroeder, P. and Salem, J.B.,"Fast Rotation of Volume Data on Data Parallel Architecture", Proceeding of Visualization'91, San Diego, CA, 50-57, October 1991.
- Westover, L.,"Footprint Evaluation for Volume Rendering", Computer Graphics, 24(4): 37-376, August 1990.
- Fuchs, H., J.Poulton, J.Eyles, T.Greer, J.Goldfeather, D.Ellsworth, S.Molnar, G.Turc, B.Tebbs, and L. Israel, "Pixel-Planes 5": "A Heterogeneous Multiprocessor Graphics System Using Processor-Enhanced Memories, "Computer Graphics, Vol. 23, No. 3, pp. 79-88, 1989.
- Fuchs H., "Direct visualisation of Volume Data. Volume-based Anatomical Application", IEEE CG & Applications, 12(4):63-71, July 1992.
- Lacroute P.,"Real-Time Volume Rendering on Shared Memory Multiprocessors Using the Shear-Warp Factorization", Proceedings 1995 Parallel Rendering Symposium, Atlanta, Georgia, October 1995, pp.15-22.
- Goldwasser S.M., Reynolds R.A., "Real-time display and manipulation of 3D medical objects: The voxel processor architecture". CVG&IP 39, 1-27, 1987.
- Amin A.B., A. Grama, and V.Singh,"Fast Volume Rendering Using an Efficient, Scalable Parallel Formulation of the Shear-Warp Algorithm", Proceedings 1995 Parallel Rendering Symposium, Atlanta, Georgia, October 1995, pp.7-14.
- Neih, J. and Levoy, M.,"Volume Rendering on Scalable Shared Memory Architecture", Proceeding of 1992 Workshop on Volume Visualization, Boston, MA, 17-24, October 1992.
- Montani, C., Perego, R., and Scopingo, R.,"Parallel Volume Visualization on a Hypercube Architecture", Proceedings of 1992 Workshop on Volume Visualization, MA, 9-16, October 1992.
- Corrie, B. and Mackerras, P.,"Parallel Volume Rendering and Data Coherence on the Fujitsu AP1000", Technical Report TR-CS-92-11, Department of Computer Science, Australian National University, Canberra, ACT, Australia, 1992.
- Law, A. and Yagel, R.,"Exploiting Spatial, Ray, and Frame Coherency for Efficient Parallel Volume Rendering", GraphiCon'96, Saint-Petersburg, Russia, July 1996.
- Law, A., Yagel, R., and Jayasimha, D. N.,"VoxelFlow: A Parallel Volume Rendering Method for Scientific Visualization", ISCA Journal of Computers and Their Applications, 1996.
- Law, A., and Yagel R.,"CellFlow: A Parallel Rendering Scheme for Distributed Memory Architectures", Proceedings of International Conference on Parallel and Distributed Processing Techniques and Applications (PDPTA'95), November 1995, pp.1-10.
- Yagel, R., and R. Machiraju, "Data Parallel Volume Rendering Algorithms", The Visual Computer, 11(6):319-338,1995.
- Malzbender, T.,"Fourier Volume Rendering", ACM Trans. On Gr.,Vol. 12, N3, July 1993.
- Sergei I. Vyatkin, Boris S. Dolgovesov, Valerie V. Ovechkin, Sergei E. Chizhik, Nail R. Kaipov “Photorealistic imaging of digital terrains, freeforms and thematic textures in realtime visualization system Voxel-Volumes”, GraphiCon ’97, Moscow.
- G.Knittel “Voxel Engine for Real-time Visualization and Examination” Eurographics ’93, Volume 12, Numb. 3, pp 39-47.
- S.I. Vyatkin, B.S. Dolgovesov, S.E. Chizhik, “Synthesis of virtual environment with object-space recursive subdivision” - Graphicon’98 Proceedings, S. Klimenko et al. (Eds). 1998.
- R.I. Velickohatny, S.I. Vyatkin, O.Yu. Ghimautdinov, B.S. Dolgovesov, N.R. Kaipov, A.V. Romanovsky, S.E. Chizhik, “ARIUS” - Family of real-time 3D graphics systems for PC platforms, Graphicon’97 Proceedings, S. Klimenko et al. (Eds). 132-135, 1997.
- Savchenko V.V., Pasko A.A., Vyatkin S.I., Dolgovesov B.S. New approach in geometric modeling: distributed and hardware implementation perspectives.// International Conference on Computers and Devices for Communication CODEC-98. Calcutta, India, 1998, pp. 285-290.
- D. Cohen and A. Shaked, "Photo-realistic Imaging of Digital Terrains", In Proceedings of Eurographics '93, 363-373, UK, 1993.
- D.W. Paglieroni and S.M.Petersen, "Parametric Heights Fiels Ray Tracing", In Proceedings of Graphics Interface '92, 192-200, 1992.
- G.Vezina and P.K.Robertson, "Terrain Perspectives on a Massively Parallel SIMD Computer", In Proceedings of CG International '91, 163-188, 1991.
- P.Pitot, Y.Duthen and R.Caubet, "A parallel architecture for ray-casting", In Computer Graphics '89, 463-472, 1989.
- G.Agranov and C.Gotsman, "Algorithms for Rendering Realistic Terrain Image Sequences and Their Parallel Implementation", In Proceedings of Graphicon '95, 153-161, 1995.
Computer Graphics & Geometry