Visualizing Flow Over Curvilinear Grid Surfaces Using Line Integral Convolution Lisa K. Forssell Computer Sciences Corporation, NASA Ames Research Center and Stanford University lisaf@cs.stanford.edu Abstract Line Integral Convolution (LIC), introduced by Cabral and Leedom in Siggraph '93, is a powerful technique for imaging and animating vector fields. We extend the LIC paradigm in three ways: 1. The existing technique is limited to vector fields over a regular Cartesian grid. We extend it to vector fields over parametric surfaces, specifically those found in curvilinear grids, used in computational fluid dynamics simulations. 2. Periodic motion filters can be used to animate the flow visualization. When the flow lies on a parametric surface, however, the motion appears misleading. We explain why this problem arises and show how to adjust the LIC algorithm to handle it. 3. We introduce a technique to visualize vector magnitude as well as vector direction. Cabral and Leedom have suggested a method for variable-speed animation, which is based on varying the frequency of the filter function. We develop a different technique based on kernel phase shifts which we have found to show substantially better results. Our implementation of these algorithms utilizes texture- mapping hardware to run in real time, which allows them to be included in interactive applications. 1. Introduction Providing an effective visualization of a vector field is a challenging problem. Large vector fields, vector fields with wide dynamic ranges in magnitude, and vector fields representing turbulent flows can be difficult to visualize effectively using common techniques such as drawing arrows or other icons at each data point, or drawing streamlines[2]. Drawing arrows of length proportional to vector magnitude at every data point can produce cluttered and confusing images. In areas of turbulence, arrows and streamlines can be difficult to interpret. Various techniques have been developed which attempt to address some of these problems. Max, Becker, and Crawfis[13], Ma and Smith[12], and Max, Crawfis, and Williams[14] have implemented systems which advect clouds, smoke, and flow volumes. These techniques show the flow on a coarse level but do not highlight finer details. Hin and Post[11] and van Wijk[16] have visualized flows with particle-based techniques, which show local aspects of the flow. Bryson and Levit [4] have used an immersive virtual environment for the exploration of flows. Helman and Hesselink[10] have generated representations of the vector field topology, which use glyphs to show critical points in the flow. In this paper we discuss a new technique for visualizing vector fields which provides an attractive alternative to existing techniques. Our technique makes use of Line Integral Convolution (LIC)[5], which is a powerful technique for imaging and animating vector fields. The image of a vector field produced with LIC is a dense display of information, and flow features on the surface are clearly evident. The LIC algorithm as presented by Cabral and Leedom in [5] is applicable only to vector fields over regular 2- dimensional Cartesian grids. However, the grids used in computational fluid dynamics simulations are often curvilinear. In this paper we show how to extend the LIC algorithm to visualize vector fields over parametric surfaces. Thus, for example, our extended algorithm allows us to visualize the flow over the surface of an aircraft or turbine. In the original work on LIC, a technique for animation of vector field visualizations is presented. Our work extends this animation technique to apply to the parametric surfaces found in curvilinear grids as well. Lastly, we present a new technique for displaying vector magnitude which can be applied to both 2- dimensional regular grids and parametric surfaces. Our method varies the speed of the flow animation to give an intuitive representation of vector magnitude. In the next section we discuss the basic LIC algorithm. In section 3 we describe our extension to curvilinear surfaces. In section 4, we discuss the implementation of animation for curvilinear grid surfaces. In section 5, we introduce our technique for displaying vector magnitude. In section 6, we describe our implementation of all the algorithms in the paper. We conclude with a brief discussion of directions for further applications of the LIC algorithm in vector field visualization. 2. Background The Line Integral Convolution (LIC) algorithm takes as input a vector field lying on a Cartesian grid and a texture bitmap of the same dimensions as the grid, and outputs an image wherein the texture has been "locally blurred" according to the vector field. There is a one-to-one correspondence between grid cells in the vector field, and pixels in the input and output image. Each pixel in the output image is determined by the one-dimensional convolution of a filter kernel and the texture pixels along the local streamline indicated by the vector field, according to the following formula: Cin(i,j)= sum ( Cin(p) * h(p)) p c t where t = the set of grid cells along the streamline within a set distance 0.5 l from the point (i, j), shown as the shaded cells in Figure 1. l = the length of the convolution kernel Cin(p) = input texture pixel at grid cell p h(p) = integral[a,b](k(w) dw) where a = the arclength of the streamline from the point (i, j) to where the streamline enters cell p b = the arclength of the streamline from the point (i,j) to where the streamline exits cell p k(w) = the convolution filter function [Figure 1] Thus each pixel of the output image is a weighted aver- age of all the pixels corresponding to grid cells along the streamline which passes through that pixel's cell. Section 4 of [5] provides the complete details of the algorithm. When this algorithm is applied at every pixel, the resulting image appears as if the texture were "smeared" in the direction of the vector field. 3. Curvilinear Grid Surfaces Because of the one-to-one correspondence between grid cells and pixels in the input/output images, the algorithm described above requires that the vector field lie on a regular, Cartesian grid. Here we show how to use the algorithm on 2-dimensional slices of structured curvilinear grids, which describe parametric surfaces. [Figure 2] We denote the curvilinear space coordinates of a point i as (xi, eta, zeta) and the physical space coordinates as (x, y, z). The vector which describes the velocity of the flow at each point is [dx/dt, dy/dt, dz/dt]T. We transform the vector field to the coordinate system of the curvilinear grid, hereafter called "computational space." The transformation from physical space to computational space is performed by multiplying the physical-space velocity vectors by the inverse Jacobian matrix s.t. | dxi | |dx dx dx |-1 |dx | |---- | |-- -- -- | |-- | | dt | |dxi deta dzeta | |dt | | | | | | | |deta | |dy dy dy | |dy | |---- | = |-- -- -- | |-- | | dt | |dxi deta dzeta | |dt | | | | | | | |dzeta| |dz dz dz | |dz | |-----| |-- -- -- | |-- | | dt | |dxi dzeta dzeta | |dt | The computational-space vectors give velocity in grid- cells per unit time. Because the data points are given for integer coordinates in computational space, this constitutes a regular Cartesian grid. We can compute a LIC-image of any 2-dimensional slice of the grid by projecting the vector field onto it. For example, if we want to examine the k=1 plane of the computational grid, which in many CFD data formats usually lies on the surface of the object about which the flow is being simulated, we drop the (z/t) term and use the 2- dimensional vector [x/t, h/t]T in the LIC algorithm. The resulting image, which is a visualization of the vector field in computational space (see Figure 3) is then mapped onto the surface in physical space using a standard inverse mapping algorithm, such as that described in [9]. The inverse mapping converts the vector field representation back into physical space (Figure 3). The final result is a visualization of the flow which is dense, easily interpreted, and effectively handles the complicated areas of the flow. 4. Animation While the image described above and shown in Figure 3 correctly shows the streamline direction of the vector field, the visualization is ambiguous in regards to whether the flow is moving forward or backward along the lines indicated. To disambiguate the direction of flow, animation is useful. Also, animating a flow visualization is physically meaningful. [Figure 3] As Cabral and Leedom discuss in [5], periodic motion filters [6] can be used together with LIC to create the impression of motion, such that a flow appears to be moving in the direction of the vector field. A small number n of LIC images are computed, where in frame i the filter kernel is phase shifted by is/n, where s is the period of the filter function. When played back, these images cause the appearance of ripples moving in the direction of the vector field. Because the filter kernel is periodic, the n frames can be cycled through continually for smooth motion. On a parametric surface, the images are `played' by texture mapping each in turn onto the surface. However, additional steps must be taken to ensure that the animation does not introduce misleading information into the visualization. The conversion from computational space to physical space maps square grid cells into quadrilaterals of varying dimensions. Therefore, the length of the convolution filter, which is measured in computational space units, is mapped to varying lengths in physical space. The length of the periodic filter determines the size and speed of the "ripples" in the animation. The speed is given by the amount of phase shift in physical space per unit time. Thus, if the period of one filter function is longer in physical space than another, that ripple appears to move faster than a shorter filter. As a result of the warping that occurs in the mapping from computational to physical space, the animation appears uneven and erratic. In areas where the grid is sparse, the flow appears as little ripples moving fast, because the convolution kernel has been compressed, and in areas where the grid is dense, the flow appears as large ripples moving slowly. Since there is no correlation between apparent speed and actual speed of the flow, this motion is highly misleading. The situation can be corrected by varying the length of the convolution filter while computing the LIC image. The length of the convolution filter must vary inversely with the grid density in the direction of the flow. Where the grid is sparse in physical space, we want to use a narrow convolution filter in computational space, as it will be stretched out when mapped. Likewise, where the grid is dense in physical space, we want to use a wide filter in computational space, as it will be compressed when mapped. See figure 4, next page. We compute frames where the length of the convolution kernel used in the LIC algorithm at each grid cell p is given by l(p) = a + b/r(p) where a is the minimum length of the kernel, measured in grid cells b controls the range of possible kernel lengths, and r is the grid density at grid cell p in the direction of the flow. a must greater than 1; if the length of the filter is 1 or less, the LIC algorithm simply returns the input texture pixel, unaffected by the vector field. b must be set to a finite length which will vary with the particular grid. r(p) for each grid cell is given by . <-- indicates that vector has been |dx dx dx |-1 |dx | normalized to unit length |-- -- -- | |-- | |dxi deta dzeta | |dt | | | | | |dy dy dy | |dy | |-- -- -- | * |-- | |dxi deta dzeta | |dt | | | | | |dz dz dz | |dz | |-- -- -- | |-- | |dxi dzeta dzeta | |dt | r(p) for the entire grid is computed by the following steps: 1) Normalize the vector field to unity in physical space. 2) Convert to computational space using the inverse Jacobian as described in section 3. 3) Take the magnitude of the computational-space vectors. [Figure 5] Figure 5 shows a single texture frame of an animation sequence computed in this way. When the 10 texture frames are played back, the flow appears smooth and even everywhere on the surface, rather than uneven and erratic. Thus we are able to use periodic motion filters even on parametric surfaces from curvilinear grids. 5. Variable Speed The next step in flow animation, whether on a regular grid or on a parametric surface, is to give a visualization of vector magnitude as well as vector direction. Thus, in a CFD flow visualization, the periodic motion should be slow where the flow has low velocity and quick where the flow has high velocity. Cabral and Leedom [5] suggest achieving this effect by varying the frequency of the filter function, while keeping its length constant. However, the limited dynamic range (experimentation shows only between 2 and 4 ripples per kernel are interpretable) and the artifacts caused by changing the shape of the filter make it difficult to use this approach for meaningful results. We have found that a better solution is to vary the amount of filter function phase shift at each grid cell in proportion to the physical-space vector magnitude. The amount of phase shift is what determines the apparent speed, given a uniform-length filter kernel. An infinitesimally small phase shift will appear not to move at all. Likewise, a 90-degree phase shift in every frame will produce a full cycle in four frames, which appears to move very quickly. (At anything greater than 180 degrees, temporal aliasing occurs.) Phase shifts ranging from 0 to 90 degrees can be mapped to the actual range of physical vector magnitude for a convincing variable-speed animation. In a frame from a variable-speed animation sequence, each pixel will be computed with a convolution kernel that has a phase shift proportional to the corresponding grid cell's physical vector magnitude. Therefore the period of the filter function is different at each pixel, and there is no fixed number of frames that can be used in a cyclic animation. Therefore, we adopt the following strategy of sampling the "real" solution and interpolating to find the pixel values which we will display. In practice, the texture frames are computed as follows: 1) First compute N LIC images, such that in image i, where qi is the amount of filter phase shift, qi= is/N. As in section 4 s is the period of the filter function. The larger N, the more accurate the visualization. The intensity of pixel p in image i is defined as T(i, p). 2) For each grid cell p , let q = y - min(y) ----------- max(y) - min(y) where y denotes physical vector magnitude at grid cell p. q is a real number in [0,1] that gives the vector magnitude in cell p relative to the magnitudes in the whole grid. 3) The intensity of pixel p in frame j of the displayed image, I(j,p), is found by interpolating linearly between the two LIC images from step (1) closest to it: Let a = q(N/4)jmodN - floor(q(N/4)j)modN at cell p. Then I(j, p) = (1-a) T(floor(q(N/4)j) mod N, p) + aT(ceil(q(N4)j) mod N, p) 6. Implementation We use the texture-mapping capabilities of a high-end workstation [8] to display surfaces with LIC-images mapped onto them in an interactive program. All LIC images are computed prior to running the interactive program, since the LIC algorithm is fairly compute- intensive (images take on the order of several seconds to minutes to compute). The hardware is capable of switching between pre-loaded textures quickly enough that we are able to run animation in real time, while the user manipulates the surface. For single-speed periodic motion, we find that 5 - 12 texture frames is sufficient for smooth animation. For variable-speed animation, we are no longer able to precompute a finite number of frames and cycle through them, because the amount of phase shift varies at every grid cell. However, steps 2 and 3 of the algorithm described in section 5 are not compute intensive, once the N LIC-images of step 1 have been precomputed. This implies that a real time implementation of variable-speed animation should be possible. Unfortunately, we have not been able to achieve this with our hardware, an SGI Reality Engine, because the time required to load a new texture into the texture cache is too long to permit good frame rates. However, we expect that within the foreseeable future texture-mapping hardware will allow fast texture definition and this feature will be possible. In the meantime, we have experimented with two alternative solutions. (1) Calculate and store a large number of texture frames using a real phase shift at every grid cell which is a linear function of the physical vector magnitude in that grid cell. The number of frames required for anything more than a few seconds of animation using this solution is so large that storage requirements quickly exceed the memory capabilities of a workstation. Therefore either (a) the short sequence of animation must be continually restarted, which causes the flow to appear to "jump" every few seconds, or (b) the animation must be stored on video or another digital playback device, and the real time interactivity possible with all other techniques described in this paper are forfeited. We show an implementation of (a). The animation is smooth in spurts of 5 seconds, and the speed of the flow is clearly varying across the surface (see video). (2) Approximate the continuous solution by choosing a minimum phase shift, f, and quantizing all phase shifts as integer multiples of f. In this solution, only M = s/f frames need to be precomputed, because the M frames will form a complete cycle. The drawback of this solution is that it is susceptible to aliasing and rasterization caused by the sampling and quantization of phase shifts. The advantage is that it can be played continually in real time on a workstation, without the jumps of solution (1). To compute the frames for this discretized solution, we follow the same steps as described for the continuous solution in section 5, but round q to the nearest multiple of 1/M. In this case there is no interpolation and M frames suffice to form a cycle of animation. In our implementation of solution (2), we see that while the speed of the flow is clearly varying, aliasing and rasterization artifacts do appear. 7. Future Work There are several promising directions for future work. First, in order for this technique to become useful for practi- cal applications, a number of extensions must be imple- mented. Foremost among these are multigrid solutions, unsteady flows, and unstructured grids. Also, we hope to extend this technique to the visualiza- tion of 3-dimensional vector fields. While the LIC algo- rithm in itself extends easily to a 3-dimensional Cartesian grid, the output image data requires additional processing before a useful image is produced. Cutting planes, isosur- faces, or volume rendering techniques will be necessary for this extension. Experimentation with input textures and convolution filters will be needed to achieve effective im- ages. Furthermore, new algorithms will be required to han- dle curvilinear grids in this situation as well. 8. Summary We have presented several extensions to Line Integral Convolution. First, we have described how to use the LIC algorithm on curvilinear grid surfaces. We have shown how to solve the problems that arise when using periodic motion filters in LIC on a curvilinear surface. Lastly, we have in- troduced a method of incorporating visualization of vector magnitude into the LIC algorithm, by showing the anima- tion at variable speeds. All algorithms are designed such that with modern graphics hardware the surfaces can be dis- played, animated, and manipulated in real time. Our visualization technique provides intuitive and accu- rate information about the vector field, and thus is a useful complement to other visualization techniques. Acknowledgments The author wishes to thank David Yip, who provided support of many forms throughout this research, the Ap- plied Research Branch staff of NASA Ames, who offered helpful discussions, Marc Levoy, who gave advice on the research and on the writing of this paper, and Brian Cabral and Casey Leedom, who made their LIC code publicly available. References [1] D. Asimov, "Notes on the Topology of Vector Fields and Flows," NASA Ames Research Center Report RNR- 93-003, February, 1993 [2] M. Bailey, C. Hansen, Introduction to Scientific Visual- ization Tools and Techniques, Course Notes, ACM SIG- GRAPH (1993). [3] G. Bancroft, F. Merritt, T. Plessel, P. Kelaita, R. Mc- Cabe, A. Globus, "FAST: A Multi-Processing Environment for Visualization of CFD," Proc.Visualization `90, IEEE Computer Society, San Francisco (1990). [4] S. Bryson, C. Levit, "The Virtual Windtunnel: An Envi- ronment for the Exploration of Three-Dimensional Un- steady Flows," Proc. Visualization `91, IEEE Computer Society, San Diego (1991). [5] B. Cabral, C. Leedom, "Imaging Vector Fields Using Line Integral Convolution," Computer Graphics Proceed- ings `93, ACM SIGGRAPH (1993) [6] W.T. Freeman, E.H. Adelson, D.J. Heeger, "Motion Without Movement," Computer Graphics Proceedings `91, ACM SIGGRAPH (1991) [7] A. Globus, C. Levit, T. Lasinski, "A Tool for Visualiz- ing the Topology of Three-Dimensional Vector Fields," Proc.Visualization `91, IEEE Computer Society, San Diego (1991). [8] Graphics Library Programming Guide, Volume II, Sili- con Graphics, Inc. (1992) [9] E. Haines, "Essential Ray Tracing Algorithms", Chapter 2 in A. Glassner,Ed., An Introduction to Ray Tracing, Aca- demic Press (1989) [10] J. Helman, L. Hesselink, "Surface Representation of Two- and Three- Dimensional Fluid Flow Topology," Proc.Visualization `90, IEEE Computer Society, San Fran- cisco (1990). [11] A. Hin, F. Post, "Visualization of Turbulent Flow with Particles," Proc. Visualization `93, IEEE Computer Society, San Jose(1993). [12] K-L. Ma, P. Smith, "Virtual Smoke: An Interactive 3D Flow Visualization Technique," Proc. Visualization `92, IEEE Computer Society, Boston (1992). [13] N. Max, B. Becker, R. Crawfis, "Flow Volumes for In- teractive Vector Field Visualization," Proc. Visualization `93, IEEE Computer Society, San Jose(1993). [14] N. Max, R. Crawfis, D. Williams, "Visualizing Wind Velocities by Advecting Cloud Textures," Proc. Visualiza- tion `92, IEEE Computer Society, Boston (1992). [15] P. P. Walatka, P. G. Buning, PLOT3D User's Manual, NASA Technical Memorandum 101067, NASA Ames Re- search Center. [16] J. van Wijk, "Rendering Surface-Particles," Proc. Vi- sualization `92, IEEE Computer Society, Boston (1992). [17] W.H. Press, et al., Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, Cam- bridge (1988).