Visualization and feature extraction in isotropic Navier-Stokes turbulence

Victor M. Fernandez and Norman J. Zabusky,
Department of Mechanical and Aerospace Engineering and CAIP Center
Rutgers University, Piscataway, NJ 08855

Smitha Bhat and Deborah Silver
Department of Electrical and Computer Engineering and CAIP Center
Rutgers University, Piscataway, NJ 08855

Shi-Yi Chen
IBM T.J. Watson Research Center / Theoretical Division
and Center for Nonlinear Studies, Los Alamos National Laboratory,
Los Alamos, NM 87545

Table of contents

1. Introduction
2. Feature extraction
3. Implementation
4. Example on the turbulence field
5. Conclusions
6. Acknowledgements
7. References

Figures

Abstract

We present feature extraction and data reduction algorithms and provide an insight in the types of problems that arise in dealing with large datasets, obtained in Navier-Stokes turbulence simulations with a 512^3 mesh resolution. The developed tools are based on thresholding, object segmentation and low order ellipsoidal quantifications and are applied to the search of coherent vortex structures associated with maxima events in the turbulence field. We underline the importance of sharing tasks between the supercomputer (CM5) and the workstation (SGI Onyx), where each machine may work more efficiently at different stages of the data processing. We obtain visualizations that show the structure of the dominant coherent objects. The reduced representations employed make it possible to examine different types of fields for possible correlations. The quantification of the objects identified by the feature extraction algorithms, should contribute to the building of models that consider both coherent structures and the random background observed in Navier-Stokes turbulence.

T.o.C.

1. Introduction

The emphasis on detecting and characterizing large-scale structures in flows has been one of the most important motivations to apply digital image processing to flow visualization [1]. Much of the time of the computational aerodynamicist or fluid dynamicist will be devoted, in the future, to extracting and displaying the critical features of the computed solution [2]. Two approaches to extracting graphical information and flow visualization are currently evolving [2]. In one approach the engineer works at a graphics workstation and displays various portions of the flow field. Because the 3D dataset is large, this requires a very powerful graphics workstation. This interactivity has the advantage that one can discover information that might not be anticipated. The other approach is to program the supercomputer itself to analyze the database. In both of these approaches, algorithms must be developed to search out and display features such as shock waves, vortices, and separation lines [2]. In the case of large 3D databases, supercomputers are used not only to compute the solution, but to search the flow field itself [2]. For unsteady flows (where the data volume may not even fit in the supercomputer's memory), analysis of the flow field and extraction of the information for graphical display will have to be done ``on the fly'' as the computation is proceeding, or as a laborious process from magnetic tape [2]. Because of the needed interaction between supercomputer and workstation, high speed communication between the machines is critical [2].

Visualization is a key component in the process of obtaining quantitative and mathematical understanding of the numerical models investigated in CFD (Computational Fluid Dynamics) [3]. Visualization tools must be able to sharpen the intuition of the engineer or scientist and also contribute to building quantitative and mathematical models. The approach we take to visualization is to present data in a form where measurement (quantification) of relevant magnitudes is emphasized. This approach, called visiometrics [4], consists of several steps [5]. The basic steps are "identification, quantification and understanding/mathematization". Although this process presents a natural progression, in fact it is an iterative refinement procedure, which at the end produces the appropriate mathematical formulation [5].

Identification involves feature extraction, which has been considered and carried out in different contexts, not only to help in the process of obtaining understanding. Feature extraction can be used as a selective way to display data, which avoids the excessive, confusing visual clutter, that arises when too much information is shown [2,6]. Besides supporting the understanding process and helping to avoid visual clutter, feature extraction may be a very effective and natural way to deal with large datasets. Many existing visualization systems make performance trade-offs that assume relative small quantities of data (solution and grid fit in RAM) [7]. Unfortunately no current computer can hold some of the time dependent CFD datasets (from 5 to 162 Gigabyte), that are currently being produced [7]. The key points in dealing with these datasets are: extraction of ``scientific data'' and use of a ``persistent object database'' [7]. These objects usually correspond to ``coherent structures'' (localized objects, which persist over ``characteristic times'' [5]). Feature extraction is accompanied by large reductions in storage requirements (0.3 - 6.7% of solution size) [7]. The disadvantage is that the possibility to examine the solution domain outside the extracted domain is lost [7].

Feature extraction can be accomplished in different ways. Initially it can be just a thresholding operation and extrema tracking [8] . It can also be viewed as the process of obtaining reduced representations (ellipsoidal or skeletal representation plus vector lines released from selected release points) for the data [9][10]. These reduced representations of different variables can be used to point out causal connections between them [3,11,9,10]. The reduced representations correspond to ``identified objects'' in the data. and can be used as tools for perusal, interpretation, quantification, feature tracking [5] and also for yuxtaposition, i.e., detailed and quantitative comparison of experimental and/or computational images of similar or different functions at the same or different times [12]. Other feature extraction procedures include the use of streamlines connected through nodes or critical points, which characterize global flow topology [13,14]. These methods have been implemented as interactive tools to extract meaning from datasets [6]. The critical point approach has also been implemented as modules (TOPO [15]) in the visualization environments (FAST [16]).

In the particular case of isotropic turbulence simulations, different researchers have used a ``probe'' or ``window'' (fixed or Lagrangian) for searching and obtaining the full time history of vortex tubes [17,18]. Their object identification algorithm consists of taking a local maximum in the field, and tracing the skeleton of the vortex tube. Diagnostics include the tube's length, curvature, length to diameter ratio and circulation [19]. Algorithms for extracting ``events'' [20] have indicated that intermittent regions may be major, if not dominant contributors to global statistics in turbulence. In studies of vortex collapse and reconnection, other researchers have used a ``diagnostics box'' which surrounds regions of significant physical behavior, detected through the searching of maxima events [21,12].

In this work we use object identification algorithms to search for the dominant vortex structures in isotropic 3D turbulence. The data was produced by Dr. Shiyi Chen [22] at Los Alamos National Laboratory and has a spatial resolution of 512^3 grid points in a uniform mesh (Fig. 1). We explore the use of both a massive parallel computer (the CM5) and the Onyx workstation communicating in the frame of CMAVS/AVS networks and the use of simple feature extraction techniques (in batch mode on the CM5) and interactive postprocessing on the workstation. The examples reported below were obtained on the CM5 (1024 nodes) and the SGI Onyx at the Advanced Computing Laboratory at Los Alamos National Laboratory (ACL-LANL) and the CM5 (512 nodes) at the National Center for Supercomputing Applications (NCSA) and the SGI Onyx at the VIZLAB in Rutgers. Our work addresses the need to distribute tasks among a network of computers. We recognize that the solver, the solutions and extractors should reside and operate on supercomputers, to be able to deal with the main difficulty, which is the size of the dataset. For a $512^3$ simulation, the volume of the dataset is of the order of 0.5 Gigabyte for a scalar field (single precision). In this mode of operation, we need to be able to send different amounts of data across the network for displaying or additional postprocessing in a workstation.

Fig. 1 Turbulence 512^3 dataset. The vorticity magnitude field shown in this figure is represented by spheres corresponding to a set of thresholded points (20% of the maximum and above).

T.o.C.

2. Feature extraction

In this work we combine a feature extraction-object identification algorithm to perform a guided search for vortex structures in the turbulence field. The feature extraction is performed in three steps:

* thresholding

* object identification

* ellipsoid fitting

The objects in the dataset, or field, are defined by a scalar function, $f( \bbox{x} )$, e.g. vorticity magnitude, and a threshold, $f_{th}$. The first stage, thresholding, consist in finding and extracting the points in the dataset where $f$ is above the threshold specified by the user, including their position in space. For uniform grids, the storage required for the position of the points corresponds to a scalar value only, since it is possible to establish a map between the 3D coordinates and a 1D value. The thresholding stage usually allows a considerable reduction in the volume of information. In Fig. 2 can be observed, that the number of grid points with low vorticity magnitude values is very large. For this data set, the amount of reduction obtained by selecting the points above the threshold 20% of the maximum is from ~500 Megabyte (0.5 Gigabyte, original dataset) to ~15 Megabyte.

Object identification is performed on the thresholded grid points. The approach we take is: starting with a grid point, we find all of the points that are connected to it by checking their distances to the starting point. The new detected points are then examined to find the grid points that are connected to them. Every time a point is found to be part of an object, this point is marked to eliminate it from the subsequent search for the other objects. The algorithm is not efficient because for every grid point the connectivity search involves $O( N )$ pairwise operations, where $N$ is number of grid points (i. e., the required operations is $O( N^2 )$). Because of the large data reduction obtained by the thresholding process, this direct method still can be carried out on the workstation. In the other hand, the process could be done more efficiently on the parallel machine, where the pairwise operations can be computed and checked simultaneously.

Once the objects have been found, an object measurement (quantification) process starts, which we call the ellipsoid fitting. This is also a data reduction process that produces abstract representations of the objects: the ellipsoids. This process consists in finding the ``mass'', the centroid, the average orientation of the vector field, the tensor of second moments, the maximum and the position of the maximum inside each object $V_{ob}$ [3].

The reduced quantities $m, \bar{\bbox{x}}, \hat{\bbox{v}}, \nu_{ij}, f_{max}$ and $\bbox{x}_{max}$ are used to produce the reduced representation of the objects: the ellipsoids. the tensor of second moments $\nu_{ij}$. The position of the ellipsoids is the centroid of the objects they represent. The axes of the ellipsoids are the square root of the eigenvalues of the tensor of second moments $\nu_{ij}$, normalized, so that the ellipsoids and the objects have the same volume. The ellipsoids are oriented according to the eigenvectors of the tensor $\nu_{ij}$. The reduced representation obtained in this manner not only fits the shape of the object, but averages over the values of the scalar field in the interior, making it possible to differentiate between objects of similar shape and volume. Taking one of the reduced quantities (usually $f_{max}$), the objects are sorted and listed for further use by the user or other postprocessing programs. The use of ellipses and ellipsoids and/or the related tensor of second moments for visualization and diagnostics has been diverse in the past. The method of moments has been applied to pattern recognition problems, where characterization of images through low-order moments is equivalent to using ellipses (in 2D) [23]. Other application in image processing consisted in finding the orientation of projections of 3D objects [24]. Ellipsoids have also been used as an iconic representation of 3D tensor fields [25]. In the context of quantification, ellipses were applied as a low-order momental model for the evolution of well separated vortices [26]. Second order moments were also used as diagnostics for the evolution of vortex patches [27].

Fig. 2 Vorticity histogram of the isotropic turbulence $512^3$ dataset. The horizontal axis is the vorticity magnitude. The vertical axis is the number of grid points in each of the vorticity magnitude bins, normalized.

T.o.C.

3. Implementation

As a starting point for implementing these tools, we selected AVS [28]. AVS is a commercial package widely available, which is based on a dataflow model for visualization and control [29]. AVS is constructed based on the modularity and networking concept. Application units, called modules, are organized and made available to the user through a ``network editor''. The modules are selected by the user to form networks in the network editor working area. The networks are therefore flexible enough to meet the particular needs of the users. The modular characteristics of AVS have the advantage that the user can produce his/her own modules and insert them in networks of standard AVS modules, which makes available to the user, all of the power of the commercial product, in the user's very specific applications (the user doesn't need to rewrite the whole package again). The other advantage is that AVS has a mechanism to share tasks among different machines via ``remote modules''. In particular there is available CMAVS, which is the particular set of AVS modules for the CM5. Therefore, CMAVS/AVS provides a framework for interaction between the supercomputer and the workstation.

We tried to apply the feature extraction-object identification algorithms in two different environments. In the first one we used a remote CMAVS module running on the CM5 to send interactively selected subdomains of a large dataset, trough the network, to AVS running on a workstation, for displaying. The main motivation of this approach is to use the memory available in the supercomputer to hold the required amount of data, which can be hardly handled by the workstation. In the second implementation, the large dataset is subjected to a thresholding postprocessing operation (selective data reduction), in the supercomputer, in batch mode. The resulting reduced dataset can then be postprocessed interactively on the workstation. The reduced dataset still covers the complete domain and allows to obtain a global picture of the dataset.

3.1 Object segmentation
In order to deal with the large number of vortex structures present in the turbulence dataset, we classify them according to their relationship to maxima events, not only of vorticity magnitude but strain-rate as well. One of the tools developed with this objective is the ``object-segment'' program, which is an enhancement to the standard iso-surfacing technique,used to extract coherent structures from 3D datasets. After defining connectivity between datapoints, points above a predefined threshold are grouped into a set of objects. Currently, there are two versions of the object-segment program: the sequential version and the parallel version on the CM5 (which supports up to 256 cubed datasets).

On the CM5, a number of parallel functions are used to represent the datapoints and for operations of connectivity and membership. We demonstrate the use of this feature isolation code in Fig. 3. In this figure, the threshold value 20% of the maximum was used to detect the objects observed, however, regions are extracted based upon their connectivity. The dataset is a 256 cubed scalar field (vorticity magnitude). The output of the program consist of a list of objects sorted and colored according to the local maxima inside them, which allows the user to select "coherent" regions for further quantification. In the figure, after the objects in the field have been identified, the predominant object (colored in red) is extracted for closer examination (Fig. 3).

Fig. 3 The object-segment code separates the objects found at a given threshold. The classification of the objects is shown in the figure by coloring them according to the local vorticity magnitude maxima inside each object. The program allows the user to select single objects for further quantification.

The sequential version of the object-segment program handles both regular and curvilinear datasets. For curvilinear datasets, the periodic option has been added as for certain datasets, if an object happens to extend along a boundary, it would be split into two objects if periodicity is not enforced. For e.g, in the tokamak dataset, points on the periodic boundary can have neighbours across the boundary i.e. like a wrapping around. However in this case, the mapping function between the physical and the computational space is to be provided along with the dataset. In Fig. 4, the use of the feature isolation code on curvilinear datasets is demonstrated.

Fig. 4 Object-segmentation of a curvilinear tokamak dataset (dimensions 32 x 129 x 65). Again, a single object is selected for further quantification.

3.2 CMAVS and the interactive window in the large dataset
The large dataset produced on the supercomputer, can be accessed more easily when it still resides there. In particular, for the CM5, parallel I/O can be used if the dataset resides in the SDA (Scalable Disk Array), which provides a capacity of 25-200 Gigabytes that can be accessed at 33-254 Megabytes/second [30]. It is possible to use the CMAVS/AVS interface to access interactively the CM5 resources through the workstation. This is accomplished by using the remote modules, which run on the CM5 and are connected in a network on AVS running on the workstation. The main objective of the CMAVS module ``my_r_s_cm'' is to use the resources of the CM5 to read the data (512^3 ~ 500 Megabyte or 256^3 ~ 67 Megabyte), using parallel I/O and to hold it in its memory. After this process has been accomplished, the data reduction process is necessary to make it possible, to transfer the data through the network. Some researchers have tried to extend the processing stage on the CM5. The information passed through the network is therefore the geometries used in the rendering by the workstation [29]. The same authors have also tried to produce the rendering (the 2D pixel map) on the CM5, which may have a smaller volume of data, than the actual geometric objects forming an isosurface (for example). In our case, the data reduction process we choose consist in extracting an interactively selected cubic subdomain. This subdomain is transferred as an AVS field trough the network to the workstation, where the standard local modules can be used for displaying (e.g. isosurface) Figs. 5 and 6.

The first dial in the module my_r_s_cm (Fig. 5) allow the user to change the size of the extracted subdomain. The other three dials permit to change the position of the extracted subdomain, so that the user can browse through the data (Fig. 6). The my_r_s_cm CMAVS module can be used to view the whole dataset by adjusting the corresponding dial and using the standard CMAVS module ``downsize'', which produces a reduced dataset by decimation.

It is possible to work in this mode interactively using the whole CM5 (1024 nodes at ACL-LANL), as in fact has been done by some researchers in special circumstances (ACL Mosaic home page). Nevertheless, in practical situations, it is hard to obtain more than 128 nodes to work interactively. At the NCSA CM5, the usual amount of resources allowed for CMAVS modules are 32 nodes for 2.5min.

Another important factor for interactivity, is the amount spent in transferring data between the CM5 and the workstation. For a subdomain of $64^3$, we find that the amount of transfer time on a local ethernet network between the CM5 and the Onyx machine (at ACL-LANL) is very acceptable. The same case running on the CM5 at NCSA (Illinois) and the Vizlab Onyx machine (New Jersey), takes a few more seconds, but is still acceptable. Researchers at ACL expect to be able to send considerable larger amounts of data by using the HIPPI network.

Fig. 5 (a) CMAVS/AVS network showing the use of the remote module my_r_s_cm running on the CM5 (b) Data file browser showing list of available files in the SDA on the CM5 and input parameters of the remote module my_r_s_cm.

Fig. 6 Run of the CMAVS module: The diagnostic box is sent across the network (CM5-NCSA/Onyx-VIZLAB) with a 64^3 vorticity magnitude subdomain. The reports on the execution of the remote CMAVS module my_r_s_cm are also included.

3.3 Selective data reduction
The shape of the vorticity magnitude histogram (Fig. 2) indicates that the volume of the data can be reduced considerably by a thresholding operation, even at low thresholds. In this approach a selective data reduction process is performed on the data (in batch mode) on the CM5. The dataset is read from the SDA using a fortran subroutine. Then, in another subroutine (written in cstar), the dataset is mapped onto a 1D array, which can be used to map the 3D coordinates of the grid points to a 1D value. This value is stored as the coordinate of the grid point, from which later the 3D coordinates can be recovered again. The thresholded points (points with a scalar value greater or equal to the threshold) are marked an enumerated. The (1D) enumeration will be stored later as the coordinate value of the grid point. The array containing the enumeration is also used to simultaneously transfer the thresholded points to the reduced array by using left indexing (cstar). The threshold we select is 20% of the maximum in the dataset. The output file of this program contains the position in the 1D array and the scalar value of the thresholded points. The reduced datasets obtained in this way can have vorticity magnitude and strain rate magnitude. The thresholded points have the appearance of a set of scatter points.

The reduced dataset is the input for an AVS network that is used to search the dataset to find the relevant vortex structures in the turbulence example. This is done mainly by the three (non-standard) modules ``read_scat'' (input), ``scatter'' (feature extraction-object identification) and ``smabox'' (diagnostics box). Non-standard data types for dataflow between the modules are introduced to handle the new formats of data (scatter points and objects lists) produced by the data reduction and object identification processes. Rendering is achieved using standard AVS modules (Fig. 7).

The modules operate as follows. After read_scat reads the dataset with the thresholded points, the module scatter performs object identification (see Sec. 2) according to (a new) threshold given interactively (through a dial) by the user (Fig. 7). The first output is a list of "interesting" objects to be examined. The second output is ``ellipsoids'' representing the objects in the list, which are colored according to the local maximum inside the object. The last output of this module is an AVS geometry (``nubes'' button) representing the "filtered" or "thresholded" $512^3$ domain, which consist of spheres with size and color proportional to the magnitude of the scalar field for each of the thresholded points (Fig. 1). A dial allows to adjust the size of the spheres (Fig. 7).

The smabox module takes as input the thresholded points from read_scat and the list of interesting objects produced by the scatter module (Fig. 7). Once the user chooses from the list of interesting objects, this module takes the thresholded points available in the original reduced dataset and puts them back in a (smaller than the original) grid, called the diagnostics box, which is used to visualize the region containing the selected object, using standard AVS tools, like isosurface. This last step is not efficient because we get again the large undesired volume of data. In order to improve this process, we require the development of specific tools for the scatter point format of the reduced dataset. The position of the diagnostics box can also be changed using dials (Fig. 7) that specify the position of the center of the box. The use of this network will be illustrated in Sec. 4. The input of the scatter module can be selected using a different field than the input of the smabox module. In this way a field can be visualized (e.g. vorticity magnitude) according to the important objects found in a related field (e.g. strain rate magnitude)

3.4 Tools for closer examination of selected regions
The previous network is used for searching the field. Once an interesting region has been selected, an entire set of scalar, vector and tensor fields can be extracted from the original large dataset, covering only the interesting subdomain ($128^3$ can be handled easily by the Onyx, but for the scale of the objects in the example provided, $64^3$ or $48^3$ is better for closer examination). These smaller datasets can be examined again using the feature extraction-object identification algorithm of section 2. Now with the purpose of using the ellipsoids as pointers to local maxima, which can be also used as release points for vector line tracers. The module ``ellipsoids'' has as input a 3D scalar field. The output are fitted ellipsoids and a list of objects (Fig. 8). The list of objects is the input to the ``vtracer'' module, together with the corresponding 3D vector field. This module traces vector lines released from the interactively selected objects. The length of the lines and the integration step can also be selected interactively (Fig. 8).

Fig. 7 (a) AVS network showing the postprocessing of the thresholded points dataset. (b) Modules produced to handle the reduced dataset: read_sct, scatter and smabox, showing input parameters.

Fig. 8 (a) AVS network showing the use of the ellipsoids as pointers and release regions for vector lines. (b) Parameters of the ellipsoids and vtracer modules.

T.o.C.

4. Example on the turbulence field

Fig. 1 illustrates the output of the scatter module consisting of spheres for each of the thresholded points (20% of maximum). The simple thresholding operation allows us to obtain a representation of the $512^3$ dataset. The size and the color of the spheres is assigned according to the scalar value at each grid point (color not visible in the black and white picture). A close up of this representation is shown in Fig. 9. Fig. 10 shows the diagnostic box. After putting the thresholded points back on the 3D grid, it is possible to use the standard AVS isosurface (Fig. 10). The ellipsoids produced by the scatter module are used to guide the search for relevant vortex structures in Fig. 11. In this figure, only the objects with non-zero volume, i.e., the objects with eigenvalues larger than one grid size, have the ellipsoid representation. All of the other smaller objects are represented by lines with length proportional to the eigenvectors. The list of objects, represented by the ellipsoids, is used by the user to perform a search of coherent structures by moving the diagnostics box (Fig. 11) throughout the domain. The translation of the diagnostics box is based on the position of the identified objects, which the user can select using the classification of the objects obtained according to their local maxima (which can be vorticity magnitude, dissipation, vorticity/strain-rate alignment, etc). These reduced representations, defined by the eigenvalues and eigenvectors of the tensor of second moments of the points forming the objects, are a low order representation of the scalar field. The ellipsoids give a sense of average and orientation of the segmented regions.

The region containing the maximum vorticity magnitude is examined in Fig. 12. The vorticity magnitude isosurface at the threshold 30% of the maximum, shows the topology of this region. The ellipsoids, fitted at the threshold 45% of the maximum mark the local maxima regions inside the objects. The vector lines trace the vorticity field inside the objects. Sets of bundles are released from the three ellipsoids. The lines in the object in the tube at the center of the picture seems to be formed by two parallel tubes winding around each other, which bifurcate in the upper right corner region. The isosurface in the lower pictures corresponds to the magnitude of the strain-rate. It can be observed that the strain-rate maxima are not localized in the same position as vorticity maxima, which has been observed in other turbulence simulations [31]. Another view of the same region is presented in Fig. 13. In Fig. 14, we present vorticity and strain-rate magnitude fields for the object classified as the 12th according to the vorticity magnitude in the $512^3$ dataset. The vorticity in the two parallel vortex tubes is of opposite sign, nevertheless it doesn't seem to correspond to a case of collapse.

Fig. 9 The spheres representation is used to visualize the vorticity magnitude field inside the diagnostics box.

Fig. 10 After the diagnostics box (in this case 200^3) is produced, standard AVS modules, like isosurface in this figure, can be used.

Fig. 11 Use of the list of identified objects: The list of identified objects, represented by the ellipsoids (scatter module) is used to guide the search of relevant vortex structures using the diagnostics box (smabox module).

Fig. 12 First view of the maximum vorticity magnitude object. The isosurface in the upper figures represents the vorticity magnitude. In the lower figures, the isosurface represents the strain-rate magnitude. The lines trace the vorticity field and are released from the ellipsoids, which fit the vorticity magnitude maxima regions.

Fig. 13 Second view of the maximum vorticity magnitude object.

Fig. 14 View of the object 12, according to the classification obtained by the scatter module.

T.o.C.

5. Conclusions

The work on feature extraction we have done, gave us a first insight in the type of problems that arise when dealing with large datasets. This first experience underlined the importance of sharing tasks between the supercomputer and the workstation, where each machine may work more efficiently at different stages of the processing of the data. We have developed some tools to search the field for coherent objects, but in order to get a better understanding of the processes observed, new quantification tools are necessary. In particular we need to find measures of the identified structures that will allow us to describe them in a statistical manner. This measures should allow us to determine if the structures we observe are generic. We expect this approach will lead to a model that will consider both the coherent structures observed and the random background in the turbulence flow field.

An essential part of diagnostics corresponds to quantifications. We have not been able to perform these in detail yet. In particular it is important to determine if the vortex structures we observe are generic. In order to investigate this question, it is necessary to establish measures that will allow us to find how often this structures appear in the field. Some of these techniques have been implemented in other turbulence computations, where the radius, length and circulation of the vortex tubes have been measured [19]. Also, it is very important to follow the evolution in time of the detected structures, to determine their dynamics. Some techniques developed with these objectives have also been implemented before [17].

T.o.C.

6. Acknowledgements

Recent support was provided by DOE/ARPA. The computations were performed on the CM5 at the Advanced Computing Laboratory-Los Alamos National Laboratory and the CM5 at the National Center for Supercomputing Applications.

T.o.C.

7. References

[1] L. Hesselink. Digital image processing in flow visualization. Ann. Rev. Fluid Mech., 20:421--485, 1988.

[2] P. G. Buning and J. L. Steger. Graphics and flow visualization in computational fluid dynamics. In Proceedings of the AIAA 7th conference in CFD, AIAA paper 85-1507, June 1985.

[3] D. Silver, N. J. Zabusky, V. Fernandez, M. Gao, and R. Samtaney. Ellipsoidal Quantification of Evolving Phenomena. In N. M. Patrikalakis, editor, Visualization of Physical Phenomena. Proceedings of Computer Graphics International '91 Symposium, pages 573--588. Springer-Verlag, June 1991.

[4] F. J. Bitz and N. J. Zabusky. David and "Visiometrics": Visualizing and quantifyig evolving amorphous objects. Computers in Physics, Nov/Dec, 1990.

[5] D. Silver and N. J. Zabusky. Quantifying visualizations for reduced modeling in nonlinear science: Extracting structures from data sets. Journal of Visual Communication and Image Representation, 4(1):46--61, 1993.

[6] R. R. Dickinson. Multi-dimensonal graphics input for feature extraction of high order 3D fields. In Extracting Meaning from Complex Data: Processing, Display, Interaction, SPIE Proceedings, Vol. 1259. SPIE, February 1990.

[7] A. Globus. A software model for visualization of large unsteady 3-D CFD results. preprint, 1991.

[8] D. Silver and N. J. Zabusky. 3D visualization and quantification of evolving amorphous objects. In E. J. Farrell, editor, Extracting Meaning from Complex Data: Processing, Display, Interaction II, Proceedings of SPIE-The International Society for Optical Engineering, Vol. 1459, pages 97--108, San Jose, California, February 1991.

[9] D. Silver, M. Gao, and N. J. Zabusky. Visualizing causal effects in 4d space-time vector fields. In Proceedings of Visualization '91, pages 12--16. IEEE Computer Society, October 22-25, 1991.

[10] M. Gao. Data Extraction and Abstraction in 3D Visualization. M.S. Thesis, Rutgers University, Graduate Progam in Electrical Engineering, New Brunswick, New Jersey, 1992.

[11] N. J. Zabusky, O. N. Boratav, R. B. Pelz, M. Gao, D. Silver and S. P. Cooper. Emergence of coherent patterns of vortex stretching during reconnection: A scattering paradigm. Physical Review Letters, 67(18):2469--2472, 1991.

[12] N. J. Zabusky, D. Silver, and R. Pelz. Visiometrics, juxtaposition and modeling. Physics Today, 46(3):24--31, 1993. Also private communication by Thomas Scheidegger.

[13] A. E. Perry and M. S. Chong. A description of eddying motions and flow patterns using critical-point concepts. Ann. Rev. Fluid Mech., 19:125--155, 1987.

[14] A. E. Perry, M. S. Chong, and B. J. Cantwell. A general classification of three-dimensional flow fields. Phys. Fluids A, 2(5):765--777, 1990.

[15] A. Globus, C. Levit, and T. Lasinski. A tool for visualizing the topology of Three-Dimensional vector fields. NASA Report RNR-91-017, April, 1991.

[16] G. Bancroft, F. Merrit, T. Plessel, P. Kelaita, R. McCabe, and A. Globus. FAST: A multi-processing environment for visualization of CFD. In Proc. Visualization '90, San Francisco, 1990. IEEE Computer Society.

[17] J. Villasenor and A. Vincent. An algorithm for space recognition and time tracking of vorticity tubes in turbulence. CVGIP: Image Understanding, 55(1):27--35, 1992.

[18] M. Meneguzzi and A. Vincent. Dynamics of vortex tubes in three-dimensional turbulence. In J. T. Beale et al., editor, Vortex Flows and Related Numerical Methods, pages 279--291. Kluwer Academic Publishers, 1993.

[19] J. Jimenez, A. A. Wray, P. G. Saffman, and R. S. Rogallo. The structure of intense vorticity in isotropic turbulence. Journal of Fluid Mechanics, 255:65--90, 1993.

[20] J. G. Brasseur and Wen-Quei Lin. Structure and statistics of intermittency in homogeneous turbulent shear flow. In A. V. Johansson and P. H. Alfredsson, editors, Advances in Turbulence 3: Proceedings of the Third European Turbulence Conference, Stockholm, July 3-6 1990. Springer-Verlag.

[21] V.M. Fernandez, N.J. Zabusky, and V.M. Gryanik. Near-singular collapse and local intensification of a ``Lissajous-elliptic'' vortex ring: Non-monotonic behavior and zero-approaching local energy densities. Physics of Fluids A, 6(7):2242--2244, 1994.

[22] S.-Y. Chen and X. Shan. High-Resolution turbulent simulations using the Connection Machine-2. Computers in Physics, 6(6):643--646, 1992.

[23] M. R. Teague. Image analysis via the general theory of moments. J. Opt. Soc. Am., 70(8), 1980.

[24] D. B. Salzman. A method of general moments for orienting 2D projections of unknown 3D objects. Computer Vision, Graphics and Image Processing, 50:129--156, 1990.

[25] C. Upson and D. Kerlick. Volumetric visualization techniques. SIGGRAPH, Course notes 13: Two and Three Dimensional Visualization Workshop, August 1989.

[26] N. J. Zabusky, M. V. Melander, and J. C. McWilliams. Symmetric vortex merger in two dimensions: Causes and conditions. J. Fluid Mech., 195:303--340, 1988.

[27] J. D. Buntine and D. I. Pullin. Merger and cancellation of strained vortices. J. Fluid Mech., 205:263--295, 1989.

[28] C. Upson. The application visualization system: A computational environment for scientific visualization. IEEE Computer Graphics and Applications, 9(4):30--, 1989.

[29] M. F. Krogh and C. D. Hansen. Visualization on massively parallel cmputers using CM/AVS. preprint, 1993.

[30] Thinking Machines Corporation. CM-5 I/O System Programming Guide. 1994.

[31] G. R. Ruetsch and M. R. Maxey. The evolution of small-scale structures in homogeneous isotropic turbulence. Physics of Fluids A, 4(12):2747--2760, 1992.

T.o.C.

Victor's homepage