Structural Inferences from Massive Datasets

1 Introduction

An essential survival skill for humans is the ability to perceive objects, classify them, and predict their behav­ ior. Much of the neural machinery is devoted to these tasks. We believe that the machinery that forms the foundation of intelligence is to be found in the mech­ anisms that support, for example, the vision, the lan­ guage, and the sensorimotor faculties. Each faculty has to solve the figure-ground problem: To identify, inter­ pret , manipulate, and track salient objects.

It is important to understand how each such modality contributes to overall intelligence. Each modality draws on proprietary representations and problem-solving com­ petences. Each takes inputs from internal sources as well as from the outside world. For example, a linguistically represented problem can require us to harness the power of our visual system so as to engage the strategies of "imagistic reasoning" or visual imagination. What ap­ pears to be a high-level cognitive ability emerges from co­ operative problem-solving as peripherals work out parts If humans do have and make use of these peripheral representations, what is the form of these representa­ tions? What are the processes that manipulate them?

In this paper we focus on the development and new ap­ plications of imagistic reasoning. In particular, we pro­ pose a model to reason about spatially distributed, tem­ porally evolving 3D structures in the domain of fluid dy­ namics. Our model aims to capture the process by which human scientists interpret numerical simulation data. Our approach follows the spatial aggregation framework [Yip and Zhao, 1996]: Programs are implemented by mixing and matching a few high-level operators such as aggregate, classify, and redescribe. We report on our progress in implementing the performance model.

The research is inspired by the practical needs for ana­ lyzing large amount of data and the intellectual curiosity to understand how humans discover regularities, summa­ rize significant events, and debug conceptual models. We use scientific datasets as a testbed in order to exploit the vast amount of accumulated specialized knowledge. Sci­ entists' reasoning about fluid dynamical objects rely on the same intuitive concepts of objects: cohesiveness and continuity [Spelke et al., 1995]. An object is cohesive ifit is internally connected, externally bounded, and moves while maintaining its connectedness and boundedness. The motion of an object is not arbitrary; it traces one connected path in space-time and leaves no gaps.

Describing complex fluid phenomena in terms of such intuitive notions of cohesive objects 1 and their interac­ tion patterns is useful for qualitative understanding and almost necessary for developing new conceptual models of dynamical mechanisms.

2 Human understanding

The entire theory of incompressible fluid flow is con­ tained in the N avier-Stokes equations. However, progress in analyzing the equations analytically has been slow. Much of the understanding of flow mechanisms comes from an interplay between numerical simulation and experimental data. To see an example of this kind of understanding that fluid dynamicists arrive at after detailed investigation of simulation data, let us quote a passage in the concluding section of [Robinson, 1991]:

1 The term "coherent objects" is more commonly used in the turbulence research community.

The main conclusion is that the self­ maintaining cycle of turbulence production in the boundary layer is driven by the forma­ tion or regeneration of embedded vortical struc­ tures ... A mature vortical arch gives rise to a trailing quasi-stream wise vortex... The quasi­ streamwise vortex elements collect and lift low­ speed near-wall fluid , leaving behind a persis­ tent low-speed streak. Relatively high-speed fluid scrubbing the vortex-lifted low-speed fluid creates a shear layer which rolls up into a new vortical arch. The new arch grows outward by agglomeration and/or self-induction and circu­ lation lift, and the cycle repeat itself.

See the diagram in Figure 1 accompanying the expla­ nation.

Much can be learned from this explanation. First , the explanation is conceptual: it does not provide any new equations. Second, the explanation refers to various kinds of structures: vortices 2 , streaks, and shear layer. Third, the spatial character (e.g., arch-shaped, stream­ wise, persistent) of the structures are described and re­ lated. And finally the dynamical process is described in terms of structural morphogenesis: the growth and ag­ glomeration of structures and their effect on others (e.g., lifting, rolling up).

2 Vortex is an intuitive notion denoting some compact re­ gion of swirling fluid motion; it is related to but not the same as the vector quantiiy vorticity, the curl of the velocity, at a flow field point .

3 The Spatial Aggregation Framework

In recent years, a computational framework , spatial ag­ gregation, has been developed to unify the description of a class of imagistic problem solvers. A program written in this framework has the following properties. It takes a continuous field and optional objective functions as input, and produces high-level descriptions of structure, behavior, or control actions. It computes a multi-layer of intermediate representations, called spatial aggregates, by forming equivalence classes and adjacency relations. It employs a small set of generic operators such as aggre­ gation, classification, and localization to perform bidirec­ tional mapping between the information-rich field and successively more abstract spatial aggregates . It uses a data structure, the neighborhood graph, as a common interface to modularize computations.

Our performance model is implemented in this style. In the following, we abstract away the details of the pro­ gramming issues and focus on the incorporation of new aggregators and classifiers for manipulating 3D struc­ tures.

4 Performance Model

4.1 Objective and Strategies

The objective of the performance model is:

Given a sampled flow-field, construct a concep­ tual model to summarize the significant events implicit in the numerical data.

Numerical simulation of turbulent flow generates mas­ sive amount of data. Even a modest sampled flow-field consists of tens of gigabytes of data describing quanti­ ties such velocity, pressure, and vorticity at each sam­ pled point . From observing how some scientists work, we propose a 5-step strategy to summarize data:

1. Isolate and characterize interesting structures.
2. Determine the spatial relationships among struc­ tures.
3. Track the morphogenesis of structures backward and forward in time.
4. Describe the statistical and dynamical significance of each structural event.
5. Note correlations among structural events and as­ sign cause-and-effect relationships.

Not all of these steps have been automated. In the following we describe the pieces that have been imple­ mented.

4.2 Aggregating field points

Many categories of structures can be defined as the iso­ surface of a scalar field. An isosurface of a scalar field F is defined as the set of points (x,y,z) such that F(x,y,z) = c where c is a constant and is often called the threshold of the isosurface. Some examples are the low pressure region, the high-speed streak, and the eigenvalues of the velocity gradient.

A popular method to render isosurface is the march­ ing cube algorithm [Lorensen and Cline, 1987]. The al­ gorithm divides the field into little cubes. The inter­ section of the surface with the cube is determined by a table-lookup based on the signs of the values F(x,y,z) - c. 3 The algorithm is reasonably fast , but does not guar­ antee a topologically consistent surface when adjacent cubes share a face where the vertices on the two diago­ nals differ in sign.

Our aggregator extends the basic marching cube algo­ rithm in four ways:

• Recursive subdivision. The aggregator recursively subdivides an ambiguous cube into 8 little cubes un­ til the newly generated cubes are unambiguous. The subdivision method gives accurate results subject to the limit in sampling resolution. The values at the sub-sampled points are determined by polynomial interpolation using values from adjacent cubes.
• Multiple isosurfaces. It requires little overhead to compute surfaces corresponding to different thresh­ olds during a single march over the entire flow field .
• Connected components. The connectivity of a sur­ face patch intersecting a cube are precomputed and stored in each cube pattern. The global connected components of the surface are determined by aggre­ gating individual patches.
• Volume estimate. The volume enclosed by an iso­ surface is computed from the contributions from all the cubes on and inside the surface. Each cube pat­ tern is pre-partitioned into tetrahedra representing the volume enclosed by a surface patch.

Other structures are defined by integral curves of a vector field . For example, a vortex line, the integral curve of the instantaneous vorticity field, is the basic building block of a vortex structure. However , it is im­ portant to distinguish vorticity lines (lines everywhere parallel to the vorticity vector) and vortices (regions of nearly circular motions in the plane normal to the core of the vortex observed at the core speed). In particular, vortices are much more cohesive and less noisy than vor­ ticity lines. We adopt the following working definition for vortices:

Definition A coherent vortex is a compact bundle of ad­ jacent, high-intensity vortex lines that are geometrically similar.

Aggregating vortex lines to form coherent structures is not straightforward. Previous researchers [Moin and Kim, 1985; Robinson, 1991] have found that vortex lines are sensitive to initial conditions. Nearby vortex lines can diverge rapidly. If the initial conditions are not cho­ sen carefully, the resulting vortex lines are likely to re­ semble badly tangled spaghetti wandering over the whole flow field, making the identification of organized struc­ ture extremely difficult. This might explain why vortex lines have not been widely used for structure identifica­ tion.

3 There are 2 8 ways a cube intersects a surface. The num­ ber can be reduced to 14 distinct configurations using rota­ tion and reflection symmetries. In our implementation, we find it more convenient to use the entire table of 2 8 intersec­ tion patterns.

Our search algorithm for vortices is adaptive[Yip, 1995]: it exploits local geometry (e.g., converging or di­ verging) of the vorticity lines to decide the direction and step size of integration.

4.3 Object Cohesiveness

There is a degree of arbitrariness in defining structures in terms of isosurfaces. The shape of the object depends on the threshold chosen, but it is often unclear how to set the appropriate thresholds. For example, what is the threshold for a low pressure region? · Intuitively we prefer structures that are stable against small changes in the threshold. We compute a range of thresholds (the lower and upper bound of the thresholds are typically known) for each scalar field . Geometric properties of the re-described objects are plotted as functions of the thresholds. We choose the mean value in the largest threshold region in which the structures are stable as the desired threshold. Figure 2 shows an example of how cohesiveness of a vortex determines the choice of the threshold.

To determine the stability of the shape of a collection of surface points for threshold selection, we find it suf­ ficient to compute a gross geometric description of the points. The detailed shape representation is postponed until the surface classification stage (see next section). The algorithm is as follows . We compute a 3D bounding box that minimizes the thickness of the distribution of the points along the long axis of the box. The cloud of points is reoriented along the principal component direc­ tions. This reorientation is determined by rotating the points with the matrix of eigenvectors of the covariance matrix of the points. The reoriented points also serve as input to the detailed surface classification.

4.4 Classifying surfaces

The marching cube algorithm and the vortex line in­ tegration provide a surface representation of the fluid structures. To classify the structures, we need to re­ cover their 3D geometry. We use generalized cylinders (GC) [Binford, 1990] as our basic shape model.

A generalized cylinder consists of a spline, a cross­ section, and a sweeping rule. It is particularly useful for our problem because to a first approximation vortices are elongated slender structures which bend and turn.

Recovery of GC descriptions from images is a hard problem in computer vision because of noise and occlu­ sion ([Zerroug and Nevatia, 1996]). Our problem is dif­ ferent and is much easier. We aim to recover the GC descriptions from 3D surface points. The recovery pro­ cedure has two stages: (1) recovery of the meridians of the GC, and (2) recovery of the cross sections and the spline of the GC.

The meridians of the GC are recovered as follows:

1. Reorient the surface points S along the principal component directions.
2. Project the points onto the plane normal to the first principal component direction, i.e., the eigenvector associated with the smallest magnitude. Call the projected points Sp.
3. Compute a delaunay triangulation of Sp.
4. Remove exterior triangles of the triangulation by a sculpturing method [Boissonat, 1984]. See Figure 3.
5. Extract an oriented boundary ( anticlockwise) of the triangulation.
6. Segment the boundary at the two points corre­ sponding to the two largest local maximum curva­ tures of the boundary. Call these two segments bd1 and bd2.
7. Backproject bd1 and bd2 to obtain two meridians of the GC.

The surface points S are rotated along the last prin­ cipal component direction. The above procedure is then repeated to obtain additional meridians of the GC.

The cross-sections and the spline of the GC are recov­ ered as follows :

1. The meridians are smoothed and parametrized by an normalized arc length t from 0 to 1.
2. Align corresponding points of the meridians at fixed increments of t.
3. Fix an ellipse to these points to obtain the cross­ section.
4. The center of the ellipses are smoothed to obtain the spline of the GC.

4.5 Re-describing Surfaces

A classified surface is a compact parametric description of a large number of field points. Like a lambda abstrac­ tion, the re-describe operator encapsulates the classified surface as a primitive object for further processing. For example, the center of mass, volume, orientation, and inertia tensor of the object are computed. From these quantities we estimate the distribution of characteristic objects (e.g, stream wise slender vortices) in interesting regions of the field (e.g., near the surface of a free surface flow).

4.6 Object Persistence

Structures typically persist over time. They evolve for a while and may grow or shrink or collide with other structures. Unlike structure recovery from 2D images, the analogous correspondence problem in recovery from 3D fields is much easier to solve. For example, we do not have to deal with occlusions.

Objects smoothly deform and move in connected paths. There is a lot of redundant information - size, shape, and relative positions - that carries over from one time slice to another. We track the largest struc­ tures (currently up to 20) and discard time slices which do not contain significant changes in the spatial charac­ ter or distributions of the tracked structures. Conversely consecutive time slices that yield qualitatively different distributions are noted. These time slices are candidates for restart of numerical simulation to acquire better tem­ porally resolved data.

4. 7 Extensions

We have only automated the extraction of kinematics and statistics of a few turbulence structures (vortices and pressure regions). Much more remains to be done:

• More structures. Ejection zones, dissipation re­ gions, high stress areas, pockets, large-scale motions - these are some of the more important structures that need to be tracked.
• Spatial correlations. We need to have efficient ways to compute pairwise or even triple correlations of spatial distributions.
• Dynamics. We need to incorporate qualitative rules of interaction to determine causal relations among structures. For example, a vortex accelerates veloc­ ity {lifts) on one side and slows down on another. A bent vortex exerts a self-induced motion along its binormal and the effect is proportional to the cur­ vature of the bent [Arms and Hama, 1965] .

5 Experiment

As a testbed, we choose a direct numerical simulation 4 of free surface turbulence generated by a shear flow in a 128 3 rectangular box. The problem is important for sev­ eral reasons. First , the mechanism for scalar transport into turbulent fluids across free-shear interfaces is impor­ tant for the design and control of industrial equipments with free surfaces. Second, interest in environment and global warming raises the issue of how to accurately es­ timate mass transfer rate of C02 between atmosphere and ocean. Third, recent experiments have found evi­ dence that low-speed streaks occur near free-slip surface under sufficient high shear. Visually these streaks resem­ ble those found in wall-layer streaks. Without the com­ plication of a solid boundary, studies in the formation of these streaks near free surface may shed additional light on the phenomenon.

It has been recently proposed that a definition of vor­ tex in an incompressible flow in terms of the second eigenvalues of the symmetric tensor S 2 and 0 2 where S and n are respectively the symmetric and antisym­ metric parts of the velocity gradient tensor \7u [Jeong and Hussain, 1995]. The symmetric part determines the contribution to local pressure minimum due to rotational motion alone.

This proposed definition is attractive computationally because it refers to the isosurfaces of a scalar field , the second eigenvalues of the symmetric tensor. We compare this definition with our working definition of vortex in terms of compact bundles of vorticity lines.

We have examined dozens of structures according to both definitions. We found a surprisingly good agree­ ment between two definitions. Figure 4 is a typical result. 5 The figure shows a superposition of the vor­ tex lines recovered by laborious integration and the sur­ face patches quickly recovered by the extended march­ ing cube algorithm. Because of results like these, we now routinely use the scalar isosurfaces as surrogates for vortices to develop our matching and correlation algo­ rithms.

4 A direct simulation proceeds from the N avier-Stokes equations without extraneous modeling assumptions.

5 We are exploring why there might be such a coincidence.

6 Related Work

Our research shares some of the objectives with the vi­ sualization work in Rutgers [Silver and Zabusky, 1993]. The Rutgers group is interested in better interactive tools for helping scientists visualizing the datasets. In particular , they use ellipsoids to fit isosurfaces of vor­ ticity, the modulus of the vorticity vector. Our effort emphasizes more on building a computational model of the process by which scientists understand those data. We ask questions about how scientists think visually and how this visual thinking is connected to the rather prim­ itive knowledge of object cohesiveness, persistence, and continuity. We also use a more general shape model and explore the coincidence of different definitions of a co­ herent vortex.

Commonsense reasoning about fluids has been pro­ posed as central problem in naive physics [Hayes, 1985]. The problem is hard because fluids do not conve­ niently divide into discrete pieces that can be easily combined. Collins and Forbus have partially imple­ mented Pat Hayes' ideas for representing fluids using both the contained-liquid and piece-of-stuff ontologies [Collins and Forbus, 1987] . The work described here is closer to the scientific end of the formalization spectrum.

Our work builds on previous work by several groups on imagistic reasoning and spatial aggregation [Yip and Zhao, 1996]. Zhao[Zhao, 1994] extracts structures oftra­ jectories from 3D phase space to guide control actions. Junker[Junker and Braunschweig, 1995] reconstructs the history of wavefront propagation from 2D seismic data.

Nishida[Nishida et al., 1991] automates the qualitative analysis of 3D phase space trajectories. Our work ex­ tends the frontier to tackle full 3+1D spatial-temporal field. Unlike previous work in domains where one has a rather complete knowledge of what the interesting struc­ tures are we need to rely on informal concepts of object cohesive~ess and coincidence for structure discovery.

7 Discussion/Conclusion

We have described a framework for understanding flow fields and progress in implementation. The specific tech­ nical contributions are threefold:

• We show how aggregation algorithms (based on marching cube and adaptive line tracing) and classi­ fication algorithms (based on recovery of generalized cylinders) can be used to extract and characterize interesting 3D fluid structures.
• In collaboration with the programs, we found a sub­ stantial agreement between a recently proposed def­ inition for vortex and the more intuitive definition of vortex as compact bundles.
• We extend the library of spatial aggregation rou­ tines by incorporating new aggregators and new classifiers.

Many of these ideas are generally useful for the re­ covery of spatial-temporal patterns in other physical do­ mams.

We believe the current work also sheds light on the nature of human intelligence. Even in sophisticated sci­ entific applications, some scientists appear to think visu­ ally and rely on commonsense intuitions of object cohe­ siveness, persistence, and continuity - not unlike tho~e shown in recent cognitive psychology literature. This shared core of knowledge attests to the importance of recycling the peripheral machinery (such as vision) f~r solving otherwise symbolically-posed problems. .Trad~­ tional AI has searched for the secrets of central mtelh­ gence for years. Perhaps the central is the peripherals.