# Visualizing High-Order Symmetric Tensor Field Structure with Differential Operators

Authors: Tim McGraw; Takamitsu Kawai; Inas Yassine; Lierong Zhu
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

The challenge of tensor field visualization is to provide simple and comprehensible representations of data which vary both directionallyand spatially. We explore the use of differential operators to extract features from tensor fields. These features can be used to generate skeleton representations of the data that accurately characterize the global field structure. Previously, vector field operators such as gradient, divergence, and curl have previously been used to visualize of flow fields. In this paper, we use generalizations of these operators to locate and classify tensor field degenerate points and to partition the field into regions of homogeneous behavior. We describe the implementation of our feature extraction and demonstrate our new techniques on synthetic
data sets of order 2, 3 and 4.


---

## Body

## 1. Introduction

Many approaches to the visualization of vector and tensor fields involve reducing the dense input data to a sparse set of features that are more easily displayed and understood. Topological approaches to vector field visualization also attempt to reduce the input data to a simpler representation of the structure of the field in terms of critical points (such as foci of sources and sinks, or the centers of vortices for, example) and a few separating streamlines (separatrices). These points and curves comprise a skeleton description of the flow field [1]. Such flow skeletons are concise and intuitive representations of vector fields. In contrast, direct visualization of the vector field in terms of glyphs or streamlines often results in distracting visual clutter. To improve the robustness of topological feature extraction, some vector field visualization methods apply feature extraction to a decomposition of the original field. One particular decomposition, the Helmholtz decomposition, allows a flow field to be separated into divergence-free (solenoidal) and curl-free (irrotational) parts. These parts may then be analyzed separately to robustly identify different types of critical points in the field.Some developments in topological tensor field visualization have proceeded by generalizing the concepts of vector field topology. Degenerate points (in 2D tensor fields) and degenerate lines (in 3D tensor fields) have commonly been defined in terms of eigenvectors of the tensors. Separatrices in the tensor case are hyperstreamlines, or integral curves of the eigenvector field. A simplified tensor field representation may also be generated by computing scalar fields of tensor invariants, and by extracting ridge lines of these scalars.In this paper, we present a novel approach to tensor field visualization which builds upon our previously developed tensor decomposition method [2]. We do not depend on the computation of eigenvalues or eigenvectors as many other approaches do. Nor do we trace deterministic or probabilistic hyperstreamlines. Instead, we generate scalar fields based on differential operators which have been generalized from existing vector field operators. With these scalar fields, we are able to locate and classify features in the field and partition the field into regions of homogeneous behavior. In this method, feature points are local extrema and segmenting curves are isocontours. This approach has the benefit of being very general with respect to tensor order. In contrast to our earlier paper [2], this work describes a visualization method which segments the data and explicitly extracts feature points. We are also now able to identify additional features (helices and saddles), whereas the earlier work only concerns sources and sinks in the field.The rest of this paper is organized as follows. In Section2, we review previous work on visualization of vector and tensor fields, with special attention being paid to vector field methods utilizing the Helmholtz decomposition. In Section 3, we describe the differential operators and their generalization to high-order tensors. In Section 4, we present the detailed formulation of the generalized Helmholtz decomposition for tensor fields. In Section 5, we describe our experiments in visualizing synthetic tensor fields and discuss the performance of the algorithm. Lastly, in Section 6, we present conclusions and discuss areas for future work.


## 2. Previous Work

The Helmholtz decomposition has been used in the solution of many problems in fluid mechanics and electromagnetism and has recently also proven to be useful in the topological analysis of vector fields. Polthier and Preuss [3, 4] used a discrete Helmholtz decomposition to robustly locate singularities in vector fields. Since they were dealing with vector fields defined on irregular meshes, they were required to carefully define divergence and curl operators. In our approach, since we consider data only on regular Cartesian grids, we use standard finite difference approximations to differential operators. Tong et al. [5] described vector fields in a multiscale framework by defining a vector field scale space in terms of the separate scale spaces of the solenoidal and irrotational parts of the field. They then smoothed the vector field by separately smoothing these constituent fields. They were able to show that this approach better preserves singularities than smoothing the input field directly. They also enhanced features of the field by separately amplifying the components of the decomposed field. Li et al. [6] used the Helmholtz decomposition to segment 2D discrete vector fields. They used Green's function method to compute the decomposition, then for each critical point, they found the region of influence using graph cuts. In their work, the critical points were defined in terms of scalar stream and potential functions.Several approaches to topological tensor field visualization have been described in previous literature. Many have considered the topology of the dominant eigenvector field [7, 8] and so defined degenerate points as locations where two or more eigenvalues are equal to each other. Zheng et al. [9, 10] described categories of feature points and numerically stable methods for extracting them. The points were then joined to form feature lines. By using discriminant analysis, they were able to find the locations of critical features with subvoxel accuracy.Approaches specific to diffusion tensor MRI have traditionally considered the topology of scalar fields of tensor invariants as defined by crease lines. Tricoche et al. [11] used this framework applied to tensor mode (which is related to the skewness of eigenvalues), and Kindlmann et al. [12, 13] used fractional anisotropy (which is related to the variance of eigenvalues) [14]. Another approach based on degenerate lines derived from probabilistic tractography has been described by Schultz et al. [15], as well as robust method for extracting crease surfaces from 3D data. In a sense, our approach is the opposite of several of these earlier approaches. Rather than extracting a scalar field and analyzing its critical points with differential operators, we apply the differential operators to the tensor field itself and then compute scalar values from the results.


## 3. Background

In the following subsections, we will define the divergence, curl, and gradient of Cartesian tensors as given by Heinbockel [16]. We will also demonstrate the usefulness of these differential operators in extracting critical points and classifying tensor field features. The features we consider are shown in Figure 1 along with the functions we use to classify them. These features are shown as a single slice of the 3D second-order tensor fields describing a helix, a saddle point, a source, a spiral, and a vortex.Second-order tensor field features and surface plots of the functions whose critical points identify them. From top to bottom, the functions are trace of helicity, divergence of the Lamb vector, trace of irrotational component, and trace of solenoidal component. The irrotational and solenoidal parts are computed from the generalized Helmholtz decomposition.
(a)
Helix(b)
Saddle(c)
Source(d)
Spiral(e)
Vortex


### 3.1. Tensor Notation

Theorder of the tensor (referred to as rank in some literature) is the number of indices required to specify each element. Tensors of order0 are represented by scalars, and tensors of order 1 are represented by vectors. An order 2 tensor can be represented as a matrix. Higher-order tensors can be represented as multiway arrays. For tensors in d-dimensional space, each index can take one of d different values. In 3 dimensions, an order-ℓ tensor then has 3ℓ components. In general, a tensor may also have two different types of indices, covariant and contravariant, usually denoted using subscripts and superscripts. For Cartesian tensors, these are equivalent, so they will denote indices using only subscripts.In writing an expression containing tensors, we will use the Einstein summation convention. This means that repeated indices are to be multiplied pairwise and summed over all possible values,(3.1)Ai1i2⋯ilBi1i2⋯il=∑i1=1d∑i2=1d⋯∑il=1dAi1i2⋯ilBi1i2⋯il.Partial derivatives will be denoted by the∂ symbol, where(3.2)∂i=∂∂xi.This notation permits tensor equations to be expressed in a very compact manner. For example, the familiar vector field divergence operator in 3D(3.3)div⁡v=∂v1∂x+∂v2∂y+∂v3∂z,
which is given as ∇·v in vector calculus notation, would be given as ∂ivi in index notation.


### 3.2. Tensor Field Divergence

In general, the divergence of an orderℓ tensor field is an order (ℓ-1) tensor field. For ℓ=4, the divergence is given in Einstein notation as(3.4)div⁡D=∂iDijkl.
This notation indicates that for all possible values of index i, the tensor components are differentiated with respect to that index and summed over. Note that when the field consists of totally symmetric tensors, the divergence tensor is also totally symmetric.The divergence operator frequently appears in conservation laws. We can understand the meaning of tensor field divergence in the context of conservation of mass. Consider the volume element shown in Figure3 and the velocity field u=[ux,0,0]T. The mass flux into the left face is jin=ρuxΔyΔz, and the flux out of the right face is jout=ρ(ux+(∂ux/∂x)Δx)ΔyΔz, where ρ is the material density. If mass is to be conserved, it is required that jout=jin, or ∂ux/∂x=0. For more general flows with nonzero y and z components, if we take the limit as the volume, V=ΔxΔyΔz, of the box approaches 0, this conservation property can be written in terms of the integral(3.5)lim⁡V→01V∫∫u⋅ndS,
which is the definition of the divergence of u.Now consider the case of a flow fieldDu resulting from the action of a second-order tensor field D on u. The flux into the left side of the volume will be given by jin=ρ(Dxxux+Dxyuy+Dxzuz)ΔyΔz, and the flux out of the right side is given by jout=ρ(Dxxux+Dxyuy+Dxzuz)+∂(Dxxux+Dxyuy+Dxzuz)/(∂xΔxΔyΔz). Conservation of mass for the one-dimensional case leads to the condition ∂/∂x(Dxxux+Dxyuy+Dxzuz)=0. By applying the product rule for differentiation, we get terms which depend on the rate of change of D, and terms which depend on the rate of change of u. The terms which depend on the rate of change of D, are characterized by the divergence of D.In index notation, applying the product rule for differentiation toDu,(3.6)div⁡(Du)=∂i(Dijuj)=uj∂iDij+Dij∂iuj.
The first term on the right hand side, uj∂iDij, depends on the divergence of D, and the second term Dij∂iuj depends on the gradient of u. This can also be rewritten in vector calculus terms as u·divD+tr(D∇u).In the context of diffusion, the divergence appears in Fick's second law(3.7)div⁡(D∇c)=0,
which is a statement of conservation of mass for the diffusion process governed by concentration gradient ∇c and diffusion tensor field D. Expanding (3.7) using the product rule in (3.6), it is easy to see that if the concentration gradient is nonzero and constant, then mass can only be conserved if the divergence of the tensor field is zero.We consider extrema (local minima and maxima) of the magnitude of the divergence of the tensor field to be one type of critical point which describes the structure of the tensor field.


### 3.3. Tensor Field Curl

The curl of an orderℓ tensor field in 3 dimensions is an order ℓ tensor field. For the order 4 case, it is defined as(3.8)curl⁡D=ɛijk(∂jDklmn),
where ɛijk is the permutation tensor(3.9)ɛijk={+1(i,j,k)is an even permutation of indices,-1(i,j,k)is an odd permutation of indices,0otherwise.
The permutation tensor is often used to define the vector cross product u×v=ɛijkujvk. This leads to the vector calculus notation curlu=∇×u. The tensor curl can be analyzed just as the divergence was in (3.6). In this case, the result is that curl(D) characterizes the part of curl(Du) which is due to the spatially varying tensor field D,(3.10)curl⁡(Dklul)=εijk∂j(Dklul)=εijk(ul∂jDkl+Dkl∂jul)=ul(εijk∂jDkl)+εijk(∂jul)Dkl.Extrema of the magnitude of the curl of the tensor field are the second type of critical point we consider.


### 3.4. Tensor Field Gradient

The gradient of an orderℓ tensor field is an order (ℓ+1) tensor field. For n=4, the gradient is given by(3.11)gradD=∂iDjklm.
For order-0 tensors, the vector calculus notation is gradϕ=∇ϕ, and for order-1 tensors, the gradient is equivalent to the Jacobian matrix. We do not analyze the tensor field directly in terms of gradient but use it within the generalized Helmholtz decomposition. Although the gradient may contain useful information, the fact that the resultant tensor field is one order higher than the field being studied makes it costly to compute and store.


### 3.5. The Lamb Tensor

The Lamb vector is commonly used to analyze turbulent flows [17]. It is given by v×∇×v where v is the velocity. The Lamb vector for the flow around a cylinder is shown in Figure 4.The Lamb vector can be generalized by considering the action of a tensor field on a vector field and expanding the definition of the Lamb vector using the product rule. The part of that expansion which is due to the changing tensor field is given by(3.12)Lpnl=εpmiDmnεij∂jDkl.
Sample Lamb tensor plots are shown in the top row of Figure 2. It can be seen from (3.12) that in general the Lamb tensor of an order ℓ tensor field is of order 2ℓ-1. In order to reduce the storage requirements for L, we perform order reduction to obtain a lower-order tensor field which approximates L. Analysis of L is done in terms of its divergence, which again results in reduction of storage requirements. We have found that reduction of L to its order-one (vector) approximation is useful for identifying saddle points in tensor fields.Lamb tensor (top row) and Helicity tensor (bottom row) of synthetic second-order tensor fields. The Lamb vector characterizes shear in the tensor field. Positive lobes are plotted in blue, and negative lobes are red. Note that the shear is primarily in an outward radial direction for all fields except for the saddle where the shear is directed inward. This makes the saddle feature easy to identify using the divergence of the Lamb tensor. The helicity tensor has both positive and negative lobes (and trace = 0) for all features except the helix.
(a)(b)(c)(d)(e)(f)
Helix(g)
Saddle(h)
Source(i)
Spiral(j)
VortexFigure 3
Flux in equals flux out for the control volume.A slice of a 3D velocity field for fluid flow around a cylinder (a) and streamlines of the corresponding Lamb vector field (b).
(a)(b)The order reduction process was described by Ozarslan et al. in the context of even-order diffusion tensors, but the same idea holds here. The homogeneous polynomial which is represented by a tensor can be decomposed into the spherical harmonic basis. The spherical harmonic coefficients for both tensors can be written as a linear combination of tensor components, resulting in a linear system which may be solved by least squares. The reduction of a symmetric third-, fifth- and seventh-order tensor to a vector is given by (3.13), (3.14), and (3.15) respectively,(3.13)vx=610(Dxxx+Dxyy+Dxzz),vy=610(Dxxy+Dyyy+Dyzz),vz=610(Dxxz+Dyyz+Dzzz),(3.14)vx=135(15Dxxxxx+3Dxxxyy+3Dxxxzz+3Dxyyyy+3Dxzzzz+Dxyyzz),vy=135(15Dyyyyy+3Dxxyyy+3Dyyyzz+3Dxxxxy+3Dyzzzz+Dxxyzz),vz=135(15Dzzzzz+3Dxxzzz+3Dyyzzz+3Dxxxxz+3Dyyyyz+Dxxyyz),(3.15)vx=1105(35Dxxxxxxx+5Dxxxxxyy+5Dxxxxxzz+3Dxxxyyyy+Dxxxyyzz+3Dxxxzzzz+5Dxyyyyyy+Dxyyyyzz+Dxyyzzzz+5Dxzzzzzz),vy=1105(5Dxxxxxxy+3Dxxxxyyy+Dxxxxyzz+5Dxxyyyyy+Dxxyyyzz+Dxxyzzzz+35Dyyyyyyy+5Dyyyyyzz+3Dyyyzzzz+5Dyzzzzzz),vz=1105(5Dxxxxxxz+Dxxxxyyz+3Dxxxxzzz+Dxxyyyyz+Dxxyyzzz+5Dxxzzzzz+5Dyyyyyyz+3Dyyyyzzz+5Dyyzzzzz+35Dzzzzzzz).
We use these relations to reduce the order of the Lamb tensor of order 2, 3, and 4 fields, respectively.Since the Lamb vector is not guaranteed to be symmetric, we apply the order reduction to the symmetric part of the tensor. The symmetric part of a tensor is given by(3.16)SymD=1l!∑j=1l!D(i1i2⋯il)j,
where (i1i2⋯iℓ)j denotes the jth permutation of the indices.In the appendix, we present an analytical demonstration of the ability of the Lamb tensor to differentiate saddles from other tensor field features.


### 3.6. Tensor Field Helicity

Helicity is motion along a helical, or corkscrew path. The equation for helicity is similar to the Lamb vector but involves the inner product of the vector field and its curl [18]. As a result, the helicity of an order ℓ tensor field is an order 2ℓ-2 tensor field. For a second-order tensor field, the helicity tensor is given by(3.17)Hml=Dmiεijk∂jDkl.Sample helicity plots are shown in the bottom row of Figure2. We use the trace of the helicity tensor to quantify the magnitude of helicity. Extrema of tr(H) indicate the location of centers of helices. In the appendix, we present an analytical justification of this choice.


## 3.1. Tensor Notation

Theorder of the tensor (referred to as rank in some literature) is the number of indices required to specify each element. Tensors of order0 are represented by scalars, and tensors of order 1 are represented by vectors. An order 2 tensor can be represented as a matrix. Higher-order tensors can be represented as multiway arrays. For tensors in d-dimensional space, each index can take one of d different values. In 3 dimensions, an order-ℓ tensor then has 3ℓ components. In general, a tensor may also have two different types of indices, covariant and contravariant, usually denoted using subscripts and superscripts. For Cartesian tensors, these are equivalent, so they will denote indices using only subscripts.In writing an expression containing tensors, we will use the Einstein summation convention. This means that repeated indices are to be multiplied pairwise and summed over all possible values,(3.1)Ai1i2⋯ilBi1i2⋯il=∑i1=1d∑i2=1d⋯∑il=1dAi1i2⋯ilBi1i2⋯il.Partial derivatives will be denoted by the∂ symbol, where(3.2)∂i=∂∂xi.This notation permits tensor equations to be expressed in a very compact manner. For example, the familiar vector field divergence operator in 3D(3.3)div⁡v=∂v1∂x+∂v2∂y+∂v3∂z,
which is given as ∇·v in vector calculus notation, would be given as ∂ivi in index notation.


## 3.2. Tensor Field Divergence

In general, the divergence of an orderℓ tensor field is an order (ℓ-1) tensor field. For ℓ=4, the divergence is given in Einstein notation as(3.4)div⁡D=∂iDijkl.
This notation indicates that for all possible values of index i, the tensor components are differentiated with respect to that index and summed over. Note that when the field consists of totally symmetric tensors, the divergence tensor is also totally symmetric.The divergence operator frequently appears in conservation laws. We can understand the meaning of tensor field divergence in the context of conservation of mass. Consider the volume element shown in Figure3 and the velocity field u=[ux,0,0]T. The mass flux into the left face is jin=ρuxΔyΔz, and the flux out of the right face is jout=ρ(ux+(∂ux/∂x)Δx)ΔyΔz, where ρ is the material density. If mass is to be conserved, it is required that jout=jin, or ∂ux/∂x=0. For more general flows with nonzero y and z components, if we take the limit as the volume, V=ΔxΔyΔz, of the box approaches 0, this conservation property can be written in terms of the integral(3.5)lim⁡V→01V∫∫u⋅ndS,
which is the definition of the divergence of u.Now consider the case of a flow fieldDu resulting from the action of a second-order tensor field D on u. The flux into the left side of the volume will be given by jin=ρ(Dxxux+Dxyuy+Dxzuz)ΔyΔz, and the flux out of the right side is given by jout=ρ(Dxxux+Dxyuy+Dxzuz)+∂(Dxxux+Dxyuy+Dxzuz)/(∂xΔxΔyΔz). Conservation of mass for the one-dimensional case leads to the condition ∂/∂x(Dxxux+Dxyuy+Dxzuz)=0. By applying the product rule for differentiation, we get terms which depend on the rate of change of D, and terms which depend on the rate of change of u. The terms which depend on the rate of change of D, are characterized by the divergence of D.In index notation, applying the product rule for differentiation toDu,(3.6)div⁡(Du)=∂i(Dijuj)=uj∂iDij+Dij∂iuj.
The first term on the right hand side, uj∂iDij, depends on the divergence of D, and the second term Dij∂iuj depends on the gradient of u. This can also be rewritten in vector calculus terms as u·divD+tr(D∇u).In the context of diffusion, the divergence appears in Fick's second law(3.7)div⁡(D∇c)=0,
which is a statement of conservation of mass for the diffusion process governed by concentration gradient ∇c and diffusion tensor field D. Expanding (3.7) using the product rule in (3.6), it is easy to see that if the concentration gradient is nonzero and constant, then mass can only be conserved if the divergence of the tensor field is zero.We consider extrema (local minima and maxima) of the magnitude of the divergence of the tensor field to be one type of critical point which describes the structure of the tensor field.


## 3.3. Tensor Field Curl

The curl of an orderℓ tensor field in 3 dimensions is an order ℓ tensor field. For the order 4 case, it is defined as(3.8)curl⁡D=ɛijk(∂jDklmn),
where ɛijk is the permutation tensor(3.9)ɛijk={+1(i,j,k)is an even permutation of indices,-1(i,j,k)is an odd permutation of indices,0otherwise.
The permutation tensor is often used to define the vector cross product u×v=ɛijkujvk. This leads to the vector calculus notation curlu=∇×u. The tensor curl can be analyzed just as the divergence was in (3.6). In this case, the result is that curl(D) characterizes the part of curl(Du) which is due to the spatially varying tensor field D,(3.10)curl⁡(Dklul)=εijk∂j(Dklul)=εijk(ul∂jDkl+Dkl∂jul)=ul(εijk∂jDkl)+εijk(∂jul)Dkl.Extrema of the magnitude of the curl of the tensor field are the second type of critical point we consider.


## 3.4. Tensor Field Gradient

The gradient of an orderℓ tensor field is an order (ℓ+1) tensor field. For n=4, the gradient is given by(3.11)gradD=∂iDjklm.
For order-0 tensors, the vector calculus notation is gradϕ=∇ϕ, and for order-1 tensors, the gradient is equivalent to the Jacobian matrix. We do not analyze the tensor field directly in terms of gradient but use it within the generalized Helmholtz decomposition. Although the gradient may contain useful information, the fact that the resultant tensor field is one order higher than the field being studied makes it costly to compute and store.


## 3.5. The Lamb Tensor

The Lamb vector is commonly used to analyze turbulent flows [17]. It is given by v×∇×v where v is the velocity. The Lamb vector for the flow around a cylinder is shown in Figure 4.The Lamb vector can be generalized by considering the action of a tensor field on a vector field and expanding the definition of the Lamb vector using the product rule. The part of that expansion which is due to the changing tensor field is given by(3.12)Lpnl=εpmiDmnεij∂jDkl.
Sample Lamb tensor plots are shown in the top row of Figure 2. It can be seen from (3.12) that in general the Lamb tensor of an order ℓ tensor field is of order 2ℓ-1. In order to reduce the storage requirements for L, we perform order reduction to obtain a lower-order tensor field which approximates L. Analysis of L is done in terms of its divergence, which again results in reduction of storage requirements. We have found that reduction of L to its order-one (vector) approximation is useful for identifying saddle points in tensor fields.Lamb tensor (top row) and Helicity tensor (bottom row) of synthetic second-order tensor fields. The Lamb vector characterizes shear in the tensor field. Positive lobes are plotted in blue, and negative lobes are red. Note that the shear is primarily in an outward radial direction for all fields except for the saddle where the shear is directed inward. This makes the saddle feature easy to identify using the divergence of the Lamb tensor. The helicity tensor has both positive and negative lobes (and trace = 0) for all features except the helix.
(a)(b)(c)(d)(e)(f)
Helix(g)
Saddle(h)
Source(i)
Spiral(j)
VortexFigure 3
Flux in equals flux out for the control volume.A slice of a 3D velocity field for fluid flow around a cylinder (a) and streamlines of the corresponding Lamb vector field (b).
(a)(b)The order reduction process was described by Ozarslan et al. in the context of even-order diffusion tensors, but the same idea holds here. The homogeneous polynomial which is represented by a tensor can be decomposed into the spherical harmonic basis. The spherical harmonic coefficients for both tensors can be written as a linear combination of tensor components, resulting in a linear system which may be solved by least squares. The reduction of a symmetric third-, fifth- and seventh-order tensor to a vector is given by (3.13), (3.14), and (3.15) respectively,(3.13)vx=610(Dxxx+Dxyy+Dxzz),vy=610(Dxxy+Dyyy+Dyzz),vz=610(Dxxz+Dyyz+Dzzz),(3.14)vx=135(15Dxxxxx+3Dxxxyy+3Dxxxzz+3Dxyyyy+3Dxzzzz+Dxyyzz),vy=135(15Dyyyyy+3Dxxyyy+3Dyyyzz+3Dxxxxy+3Dyzzzz+Dxxyzz),vz=135(15Dzzzzz+3Dxxzzz+3Dyyzzz+3Dxxxxz+3Dyyyyz+Dxxyyz),(3.15)vx=1105(35Dxxxxxxx+5Dxxxxxyy+5Dxxxxxzz+3Dxxxyyyy+Dxxxyyzz+3Dxxxzzzz+5Dxyyyyyy+Dxyyyyzz+Dxyyzzzz+5Dxzzzzzz),vy=1105(5Dxxxxxxy+3Dxxxxyyy+Dxxxxyzz+5Dxxyyyyy+Dxxyyyzz+Dxxyzzzz+35Dyyyyyyy+5Dyyyyyzz+3Dyyyzzzz+5Dyzzzzzz),vz=1105(5Dxxxxxxz+Dxxxxyyz+3Dxxxxzzz+Dxxyyyyz+Dxxyyzzz+5Dxxzzzzz+5Dyyyyyyz+3Dyyyyzzz+5Dyyzzzzz+35Dzzzzzzz).
We use these relations to reduce the order of the Lamb tensor of order 2, 3, and 4 fields, respectively.Since the Lamb vector is not guaranteed to be symmetric, we apply the order reduction to the symmetric part of the tensor. The symmetric part of a tensor is given by(3.16)SymD=1l!∑j=1l!D(i1i2⋯il)j,
where (i1i2⋯iℓ)j denotes the jth permutation of the indices.In the appendix, we present an analytical demonstration of the ability of the Lamb tensor to differentiate saddles from other tensor field features.


## 3.6. Tensor Field Helicity

Helicity is motion along a helical, or corkscrew path. The equation for helicity is similar to the Lamb vector but involves the inner product of the vector field and its curl [18]. As a result, the helicity of an order ℓ tensor field is an order 2ℓ-2 tensor field. For a second-order tensor field, the helicity tensor is given by(3.17)Hml=Dmiεijk∂jDkl.Sample helicity plots are shown in the bottom row of Figure2. We use the trace of the helicity tensor to quantify the magnitude of helicity. Extrema of tr(H) indicate the location of centers of helices. In the appendix, we present an analytical justification of this choice.


## 4. The Helmholtz Decomposition

The Helmholtz decomposition [19] of a vector field, v, is given by(4.1)v=∇ϕ+∇×ψ+h,
where ∇ϕ is the gradient of a scalar potential field ϕ, ∇×ψ is the curl of a vector stream field ψ, and h is a harmonic vector field. Note that ∇ϕ is irrotational, so it is useful for isolating features such as local maxima and minima of divergence (foci of sources and sinks) in v without interference from curl-based features. Likewise, ∇×ψ is solenoidal and is useful for isolating centers of vortices in v. The harmonic vector field, h, is both solenoidal and irrotational and typically is of small magnitude.The Helmholtz decomposition is known to extend from vector fields to differential forms (where it is often referred to as the Hodge decomposition [20]), and therefore also to antisymmetric tensors.


### 4.1. Helmholtz Decomposition of Tensor Fields

We utilize the generalized Helmholtz decomposition as described by McGraw et al. [2]. Using the previously defined operators, the Helmholtz decomposition of second- and fourth-order tensor fields is given as(4.2)Dij=∂iϕj+ɛimn(∂mψnj)+Hij,Dijk=∂iϕjk+ɛimn(∂mψnjk)+Hijk,Dijkl=∂iϕjkl+ɛimn(∂mψnjkl)+Hijkl.
Just as in the vector field case, we have div(curlψ)=0 and curl(gradϕ)=0. The formulation can be made for tensors of arbitrary order, but we present the order 2, 3, and 4 decompositions since those are the basis for the experiments in Section 5. The second- and fourth-order tensors describe a wide range of physical phenomena. Tensors can represent diffusivity [21, 22], the fiber orientation distribution function (ODF) in DT-MRI [23], mechanical stress and strain [24], and electromagnetic quantities [25]. More abstractly, tensors may represent covariance, skew, kurtosis, and higher-order moments of multivariate probability distributions as well as homogeneous polynomials. Third-order tensors have been used to describe the apparent bidirectional reflectance distribution function (BRDF) in face relighting applications [26].For numerical computational purposes, we will be reshaping tensor fields into column vectors. For each tensor component, the elements of the field are vectorized in lexical order of the spatial coordinates(x,y,z). The components are then ordered within the vector according to lexical order of indices. An input second-order tensor field, D, with spatial dimensions n×m×p is then vectorized as [DxxDxyDxzDyxDyyDyzDzxDzyDzz]T which has 9mnp elements.We will represent the discretized operators as block matrices where the blocks correspond to finite difference operators applied to a single tensor component. For 3D fields, the multidimensional difference matrices are given by(4.3)Δx=Ip×p⊗Im×m⊗Δn×n,Δy=Ip×p⊗Δm×m⊗In×n,Δz=Δp×p⊗Im×m⊗In×n,
where In×n is an n×n identity matrix, ⊗ is the Kronecker product and, Δn×n is an n×n finite difference matrix. We use central differences for approximating derivatives, in which case Δ is given by(4.4)Δ=12[0-10⋯0+10-1⋱⋮0+10⋱0⋮⋱⋱⋱-10⋯0+10].
This definition of this matrix may be modified as needed to impose boundary conditions on the tensor field.We can approximate the curl of the second-order tensor fieldψij as(4.5)[CxxCxyCxzCyxCyyCyzCzxCzyCzz]=C[ψxxψxyψxzψyxψyyψyzψzxψzyψzz],
where C is the sparse block matrix(4.6)C=[000-Δz00Δy000000-Δz00Δy000000-Δz00ΔyΔz00000-Δx000Δz00000-Δx000Δz00000-Δx-Δy00Δx000000-Δy00Δx000000-Δy00Δx000].
Similarly, the gradient of the first-order tensor field ϕi is given by(4.7)[GxxGxyGxzGyxGyyGyzGzxGzyGzz]=G[ϕxϕyϕz],
where G is the sparse block matrix(4.8)G=[Δx000Δx000ΔxΔy000Δy000ΔyΔz000Δz000Δz].
Although we do not need the divergence operator to compute the Helmholtz decomposition, we will use it to identify critical points. The divergence of the second-order tensor field, D, is given by(4.9)[dxdydz]=D[DxxDxyDxzDyxDyyDyzDzxDzyDzz],
where D is the sparse block matrix(4.10)D=[Δx00Δy00Δz000Δx00Δy00Δz000Δx00Δy00Δz].The discretized operators for fourth-order tensors are not given here since they contain 81 rows each, but they are easily generated from the equations in the previous sections.To perform the generalized Helmholtz decomposition, we solve the least squares problem(4.11)min⁡ψ,ϕ‖D-Cψ-Gϕ‖F2,
where ∥·∥F denotes the Frobenius norm of the tensor ∥Xik∥F=trace(XijXjk).Using the fact thatCTG=GTC=0, we implement this numerically by alternately solving the normal equations(4.12)CTCϕ=CTD,(4.13)GTGψ=GTD,
using a stabilized biconjugate gradients (BiCG-stab) method until convergence is reached. Although the matrices on the left-hand sides of (4.12) are symmetric, they are not positive definite, so the standard conjugate gradients method cannot be used. No preconditioner was used for the results presented in this paper. The derivatives of all tensor components are constrained to be zero across each boundary. The biconjugate gradients method can be implemented in terms of basic matrix arithmetic and requires very little temporary storage. If the size of the dataset does not permit the matrices C and G to be stored in memory, the matrix multiplication can be implemented procedurally. This is the approach we have used for all of the fourth-order tensor data processed in this paper. We do not explicitly solve for H, the harmonic part of the field, but instead let H=D-gradϕ-curlψ. As a result there is no fitting error, but the quality of the decomposition can be assessed by checking the criteria which define the decomposition that gradϕ is irrotational, curlψ is solenoidal, and that H is harmonic (both irrotational and solenoidal).Synthetic second-tensor fields were generated (as described in the appendix) from 5 fields each containing one type of critical point (helix, saddle, source, spiral, and vortex). The tensor field plot shown in Figure5 was generated by plotting the radial surfaces r(x)=Dijxixj for unit vectors x.Helmholtz decomposition results for synthetic second-order tensor fields.
(a)
Helix(b)
Saddle(c)
Source(d)
Spiral(e)
Vortex(f)
Helix solenoidal part(g)
Saddle solenoidal part(h)
Source solenoidal part(i)
Spiral solenoidal part(j)
Vortex solenoidal part(k)
Helix irrotational part(l)
Saddle irrotational part(m)
Source irrotational part(n)
Spiral irrotational part(o)
Vortex irrotational part(p)
Helix harmonic part(q)
Saddle harmonic part(r)
Source harmonic part(s)
Spiral harmonic part(t)
Vortex harmonic partThe results of the generalized Helmholtz decomposition of the tensor fields are shown in the second and third rows of Figure5. The harmonic field, which is typically of small magnitude for vector field decompositions, can be substantial in terms of the tensor trace, but it is extremely smooth nearly constant in all of our synthetic field experiments as shown in the fourth row of Figure 5. In the decomposed fields, there seems to be a correspondence between sources of positive-definite tensors and vortices of negative-definite tensors and vice verse. For example, the source is decomposed as a sum of a negative vortex and a positive source (ignoring the harmonic part). The vortex is a sum of a positive vortex and a negative source. When decomposing the spiral, the components differ not only in sign, but in handedness. The right-handed spiral is decomposed as a right-handed positive spiral and a left-handed negative spiral. We quantify the positiveness/negativeness of tensors by using the tensor trace: negative-definite tensors have tr(D)<0 and positive-definite tensors have tr(D)<0. In the results section, we will use this face to distinguish sources, vortices, and spirals.The decompositions of the helix and saddle are also shown for reference, but the location and classification of these critical points do not depend on the Helmholtz decomposition.When analyzing the results of the Helmholtz decomposition, it is important to keep in mind that the fields produced are not unique. You may add any constant tensor field tocurlψ or gradϕ and subtract that field from H and obtain new fields which satisfy the conditions of the decomposition. This fact also allows us to generate a decomposition which does maintain positive (or negative) definiteness of the resulting fields. For example, we can make gradϕ and curlψ positive by adding a constant isotropic tensor field to each and subtracting the isotropic tensor fields from H. However, our subsequent visualizations do not depend on either technique.Differential operators allow us to classify critical points, but the Helmholtz decomposition gives us a method of classifying arbitrary points in the field by using the difference of errors,d, given by(4.14)eψ=‖D-Cψ‖F,eϕ=‖D-Gϕ‖F,d=eψ-eϕ.
This function can be used to segment the field into two regions. When d<0, the field is better approximated by curlψ, and when d>0, it is better approximated by gradϕ. Li et al. [6] used a similar decision criterion based on a ratio of errors for segmentation of vector fields. The contour line d=0 separates the field into relatively low-divergence and low-curl regions.


## 4.1. Helmholtz Decomposition of Tensor Fields

We utilize the generalized Helmholtz decomposition as described by McGraw et al. [2]. Using the previously defined operators, the Helmholtz decomposition of second- and fourth-order tensor fields is given as(4.2)Dij=∂iϕj+ɛimn(∂mψnj)+Hij,Dijk=∂iϕjk+ɛimn(∂mψnjk)+Hijk,Dijkl=∂iϕjkl+ɛimn(∂mψnjkl)+Hijkl.
Just as in the vector field case, we have div(curlψ)=0 and curl(gradϕ)=0. The formulation can be made for tensors of arbitrary order, but we present the order 2, 3, and 4 decompositions since those are the basis for the experiments in Section 5. The second- and fourth-order tensors describe a wide range of physical phenomena. Tensors can represent diffusivity [21, 22], the fiber orientation distribution function (ODF) in DT-MRI [23], mechanical stress and strain [24], and electromagnetic quantities [25]. More abstractly, tensors may represent covariance, skew, kurtosis, and higher-order moments of multivariate probability distributions as well as homogeneous polynomials. Third-order tensors have been used to describe the apparent bidirectional reflectance distribution function (BRDF) in face relighting applications [26].For numerical computational purposes, we will be reshaping tensor fields into column vectors. For each tensor component, the elements of the field are vectorized in lexical order of the spatial coordinates(x,y,z). The components are then ordered within the vector according to lexical order of indices. An input second-order tensor field, D, with spatial dimensions n×m×p is then vectorized as [DxxDxyDxzDyxDyyDyzDzxDzyDzz]T which has 9mnp elements.We will represent the discretized operators as block matrices where the blocks correspond to finite difference operators applied to a single tensor component. For 3D fields, the multidimensional difference matrices are given by(4.3)Δx=Ip×p⊗Im×m⊗Δn×n,Δy=Ip×p⊗Δm×m⊗In×n,Δz=Δp×p⊗Im×m⊗In×n,
where In×n is an n×n identity matrix, ⊗ is the Kronecker product and, Δn×n is an n×n finite difference matrix. We use central differences for approximating derivatives, in which case Δ is given by(4.4)Δ=12[0-10⋯0+10-1⋱⋮0+10⋱0⋮⋱⋱⋱-10⋯0+10].
This definition of this matrix may be modified as needed to impose boundary conditions on the tensor field.We can approximate the curl of the second-order tensor fieldψij as(4.5)[CxxCxyCxzCyxCyyCyzCzxCzyCzz]=C[ψxxψxyψxzψyxψyyψyzψzxψzyψzz],
where C is the sparse block matrix(4.6)C=[000-Δz00Δy000000-Δz00Δy000000-Δz00ΔyΔz00000-Δx000Δz00000-Δx000Δz00000-Δx-Δy00Δx000000-Δy00Δx000000-Δy00Δx000].
Similarly, the gradient of the first-order tensor field ϕi is given by(4.7)[GxxGxyGxzGyxGyyGyzGzxGzyGzz]=G[ϕxϕyϕz],
where G is the sparse block matrix(4.8)G=[Δx000Δx000ΔxΔy000Δy000ΔyΔz000Δz000Δz].
Although we do not need the divergence operator to compute the Helmholtz decomposition, we will use it to identify critical points. The divergence of the second-order tensor field, D, is given by(4.9)[dxdydz]=D[DxxDxyDxzDyxDyyDyzDzxDzyDzz],
where D is the sparse block matrix(4.10)D=[Δx00Δy00Δz000Δx00Δy00Δz000Δx00Δy00Δz].The discretized operators for fourth-order tensors are not given here since they contain 81 rows each, but they are easily generated from the equations in the previous sections.To perform the generalized Helmholtz decomposition, we solve the least squares problem(4.11)min⁡ψ,ϕ‖D-Cψ-Gϕ‖F2,
where ∥·∥F denotes the Frobenius norm of the tensor ∥Xik∥F=trace(XijXjk).Using the fact thatCTG=GTC=0, we implement this numerically by alternately solving the normal equations(4.12)CTCϕ=CTD,(4.13)GTGψ=GTD,
using a stabilized biconjugate gradients (BiCG-stab) method until convergence is reached. Although the matrices on the left-hand sides of (4.12) are symmetric, they are not positive definite, so the standard conjugate gradients method cannot be used. No preconditioner was used for the results presented in this paper. The derivatives of all tensor components are constrained to be zero across each boundary. The biconjugate gradients method can be implemented in terms of basic matrix arithmetic and requires very little temporary storage. If the size of the dataset does not permit the matrices C and G to be stored in memory, the matrix multiplication can be implemented procedurally. This is the approach we have used for all of the fourth-order tensor data processed in this paper. We do not explicitly solve for H, the harmonic part of the field, but instead let H=D-gradϕ-curlψ. As a result there is no fitting error, but the quality of the decomposition can be assessed by checking the criteria which define the decomposition that gradϕ is irrotational, curlψ is solenoidal, and that H is harmonic (both irrotational and solenoidal).Synthetic second-tensor fields were generated (as described in the appendix) from 5 fields each containing one type of critical point (helix, saddle, source, spiral, and vortex). The tensor field plot shown in Figure5 was generated by plotting the radial surfaces r(x)=Dijxixj for unit vectors x.Helmholtz decomposition results for synthetic second-order tensor fields.
(a)
Helix(b)
Saddle(c)
Source(d)
Spiral(e)
Vortex(f)
Helix solenoidal part(g)
Saddle solenoidal part(h)
Source solenoidal part(i)
Spiral solenoidal part(j)
Vortex solenoidal part(k)
Helix irrotational part(l)
Saddle irrotational part(m)
Source irrotational part(n)
Spiral irrotational part(o)
Vortex irrotational part(p)
Helix harmonic part(q)
Saddle harmonic part(r)
Source harmonic part(s)
Spiral harmonic part(t)
Vortex harmonic partThe results of the generalized Helmholtz decomposition of the tensor fields are shown in the second and third rows of Figure5. The harmonic field, which is typically of small magnitude for vector field decompositions, can be substantial in terms of the tensor trace, but it is extremely smooth nearly constant in all of our synthetic field experiments as shown in the fourth row of Figure 5. In the decomposed fields, there seems to be a correspondence between sources of positive-definite tensors and vortices of negative-definite tensors and vice verse. For example, the source is decomposed as a sum of a negative vortex and a positive source (ignoring the harmonic part). The vortex is a sum of a positive vortex and a negative source. When decomposing the spiral, the components differ not only in sign, but in handedness. The right-handed spiral is decomposed as a right-handed positive spiral and a left-handed negative spiral. We quantify the positiveness/negativeness of tensors by using the tensor trace: negative-definite tensors have tr(D)<0 and positive-definite tensors have tr(D)<0. In the results section, we will use this face to distinguish sources, vortices, and spirals.The decompositions of the helix and saddle are also shown for reference, but the location and classification of these critical points do not depend on the Helmholtz decomposition.When analyzing the results of the Helmholtz decomposition, it is important to keep in mind that the fields produced are not unique. You may add any constant tensor field tocurlψ or gradϕ and subtract that field from H and obtain new fields which satisfy the conditions of the decomposition. This fact also allows us to generate a decomposition which does maintain positive (or negative) definiteness of the resulting fields. For example, we can make gradϕ and curlψ positive by adding a constant isotropic tensor field to each and subtracting the isotropic tensor fields from H. However, our subsequent visualizations do not depend on either technique.Differential operators allow us to classify critical points, but the Helmholtz decomposition gives us a method of classifying arbitrary points in the field by using the difference of errors,d, given by(4.14)eψ=‖D-Cψ‖F,eϕ=‖D-Gϕ‖F,d=eψ-eϕ.
This function can be used to segment the field into two regions. When d<0, the field is better approximated by curlψ, and when d>0, it is better approximated by gradϕ. Li et al. [6] used a similar decision criterion based on a ratio of errors for segmentation of vector fields. The contour line d=0 separates the field into relatively low-divergence and low-curl regions.


## 5. Results

The generalized Helmholtz decomposition was implemented in Matlab and run on a system with Intel Quad Core QX6700 2.66 GHz CPU and 4 GB RAM. The algorithm was applied to the synthetic datasets as described below. The fitting quality of the decomposition can be analyzed in terms of the magnitude of the harmonic term∥D-gradϕ-curlψ∥, and ∥divcurlψ∥ and ∥curlgradϕ∥. All three of these terms should be very near to zero. Fitting results are presented in Tables 1 and 2, and timing results are given in Figures 8.Table 1
Random second-order tensor fitting error.
∥divcurlψ∥∥divD∥∥curl∇ϕ∥∥divD∥∥divH∥∥divD∥∥curlH∥∥divD∥mean5.73E-0173.09E-0176.07E-0165.11E-016variance7.12E-0385.46E-0384.23E-0353.99E-035max5.77E-0173.11E-0176.17E-0165.21E-016Table 2
Random fourth-order tensor fitting error.
∥div  curlψ∥∥divD∥∥curl∇ϕ∥∥divD∥∥divH∥∥divD∥∥curlH∥∥divD∥mean2.66E-0081.35E-0081.69E-0071.48E-007variance2.80E-0215.11E-0222.13E-0198.02E-020max2.67E-0081.35E-0081.70E-0071.48E-007We assessed the fitting quality by generating 100 random fields, with uniformly distributed tensor components in the range [−1, 1]. The tables show the mean, variance, and maximum of the each fitting parameters. Examples of the random fields and their error difference plots are shown in Figures6 and 7.A sample random second-order field and the error difference.
(a)
Order 2(b)
Order 4A sample random fourth-order field and the error difference.
(a)
Order 2(b)
Order 4For all datasets, the decomposition obeys the expected Helmholtz properties. The divergence of the solenoidal part is many orders of magnitude smaller than the divergence of the input field. Likewise, the curl of the irrotational part of the field is many orders of magnitude smaller than the input field. The resulting harmonic part has very small divergence and curl. The same observations hold when considering many randomly generated field. The variance of the fitting parameters is very small, suggesting that the decomposition has consistent behavior over a large number of fields.It is clear from these tables that the fitting quality of the generalized Helmholtz decomposition for the fourth-order tensor fields is worse than the second-order tensor fields but still shows a great reduction in divergence and curl compared to the input data.The runtime of our visualization approach is dominated by the time it takes to compute the generalized Helmholtz decomposition. A summary of timing results is presented in Figure8. The red line shows the time to compute the decomposition of a second-order tensor field by forming the sparse matrices C and G in memory and solving the normal equations using BiCG. The blue line shows the computation time for the same data when implementing matrix multiplications by C and G procedurally. For the fourth-order tensor field, we do not attempt to construct the matrices since they are so large.Figure 8
Run time of Helmholtz decomposition for tensor fields of various orders and sizes.We analyzed the source, vortex, and spiral features by generating many synthetic tensor fields. We start with a vector field containing a source (as described in the appendix) and generate different spiral by rotating each vector in the field. When the rotation angle is90∘ the source is transformed into a vortex. We also varied the strength of the spiral feature by adding an isotropic component to each voxel in the field. Several of the generated tensor fields are shown on the left side of Figure 9. This method of generating the data reflects the fact that these features (source, spiral, and vortex) lie on a continuum. We can quantify the position along this continuum using the trace of the irrotational part as given by the Helmholtz decomposition. As shown in the center column of Figure 9, the trace is positive when the spiral is more vortexlike, and negative when the spiral is more sourcelike. The location and strength of the feature may be quantified by the divergence of the irrotational part of the field, as shown in the right column of Figure 9. The solenoidal part of the tensor field contains similar information.Spiral tensor field features with varying spiral angle (left), the trace of the irrotational part ofD (middle), and the divergence of the irrotational part of D (right).
(a)
Spiral 1(b)
Spiral 2(c)
Spiral 3(d)
Spiral 4(e)
Spiral 5(f)
Spiral 6(g)
Spiral 7Our next experiment involves fields containing multiple critical points. Synthetic second- and fourth-order tensor fields were generated (as described in the appendix) from the sources and vortices shown in Figure10 and then summing them, D=D1+D2+D3+D4. The tensor fields in all following figures are visualized by plotting the radial surfaces r(x)=Dijxixj for unit vectors x in the second-order case, and plotting r(x)=Dijkxixjxk, and r(x)=Dijklxixjxkxl in the third- and fourth-order cases. The surface is colored blue when r is positive and red when r is negative.Vortices and sources used to construct the synthetic field.
(a)
D1(b)
D2(c)
D3(d)
D4Several interesting observations can be made from the results of the Helmholtz decomposition shown in Figure11. The critical points in the original field Figure 10(a) are not clearly visible. Summing the 4 fields together made it very difficult to visually discern the critical points in the original field, but in the decomposed fields, they are quite evident.Helmholtz decomposition results for fourth-order synthetic tensor fieldD=curlψ+gradϕ+H.
(a)
Synthetic fieldD=(D1+D2+D3+D4)2(b)
curlψ(c)
gradϕ(d)
harmonic termHThe error difference function,d, for the second- and fourth-order tensor decomposition, is shown in Figure 12 as a filled contour plot. Note that the “hot”, colored regions represent the solenoidal part of the field, and the “cool” colors represent the irrotational part of the field. The smaller contour curves encircle the critical point in the field. The contour at d=0 segments the field into two regions and even separates the nearby critical points in the center of the field. This function can be seen as a simple classifier for the critical points, separating nodes of sources/sinks from centers of vortices. The d=0 isocontour should not be interpreted as a hyperstreamline, but instead as the boundary between solenoidal and irrotational regions in the field. The generality of the decomposition with respect to tensor order is reflected in the similarity of these 2 functions. In both cases, the error difference function segments the field into disjoint regions of irrotational and solenoidal behavior.Error difference,d, for second- and fourth-order synthetic fields.
(a)
Order 2(b)
Order 4The previous experiment involved data generated by summing even-order tensor fields. The same experiments applied to third-order tensor fields did not have comparable results since directions whereDijkxixjxk are negative result in destructive interference when the fields are added together. This effect causes the locations of critical points in the resulting field to be significantly different from the original fields. To overcome this problem, we designed a new experiment in which a set of synthetic fields were generated from the vector field v=[sin2x+y,cos(x+y2),0] by repeated vector contraction, as described in the appendix. The fields are displayed as glyphs in Figure 13, and the resulting error differences are shown in Figure 14. The error differences all show similar results although the exact location of the d=0 isocontour (shown as a thick black line) varies with tensor order.Synthetic tensor fields.
(a)
Order 1(b)
Order 2(c)
Order 3(d)
Order 4Error difference of Helmholtz decomposition.
(a)
Order 1(b)
Order 2(c)
Order 3(d)
Order 4Five tensor fields of order 2 and 4 were generated by transforming and summing the fields shown in Figure5. The ground-truth locations of critical points in the original fields are shown in Figure 15. The second-order tensor field is shown in Figure 16, and the error difference is shown in Figure 17 with the d=0 isocontour being plotted as a thicker black line. Figures 18 and 19 show the results of field visualization using filled contour plots and crosses marking the locations of critical points. The plots contain the error difference contours, but the color map has been changed to differentiate these images from previous results. These plots also contain extrema and filled contours of the trace of the helicity tensor and the divergence of the Lamb tensor. Helicity contours are filled with purple, and local maxima are denoted by cyan-colored crosses. Saddle contours are filled with green, and local maxima are denoted by yellow-colored crosses. Spiral points are drawn in green but may occur inside of the red or blue contours of sources or vortices.Ground-truth location of critical points.
(a)(b)(c)(d)(e)Resulting tensor field.
(a)(b)(c)(d)(e)Error difference.
(a)(b)(c)(d)(e)Field visualization.
(a)(b)(c)(d)(e)Order 4 field visualization.
(a)(b)(c)(d)(e)Some spurious contours do occur, but some of these are due to boundary effects, as in Figure18(b), and some are very small, such as in Figure 18(a). Some unexpected contours occur due to the process of generating the tensor fields. For example, regions surrounded by multiple vortices resemble saddle points, and this is reflected in Figure 18(c). In all cases, we were able to identify and correctly classify all of the critical points which were present in the original fields.Also note the similarity between the order 2 and 4 results in Figures18 and 19. This suggests that in future work these same differential operators may be applied to tensor fields of order 5 and higher.


## 6. Conclusions

Differential operators can provide intuitive and useful information about the structure of tensor fields. Specifically, local peaks in magnitude of divergence and curl correspond to critical points in the tensor field. We have introduced generalizations of helicity and the Lamb vector for identifying helices and saddle points in tensor fields. The usefulness of these measures was validated analytically and demonstrated on synthetic tensor fields of order 2, 3, and 4. We have also presented a method of segmenting tensor fields in a meaningful way by finding a separating curve which partitions the field into approximately irrotational and solenoidal regions.The formulations we have presented are general with respect to tensor order and do not require eigenvalues to be computed. However, processing time for fourth-order tensor fields is high, and future work will involve exploiting constraints, such as symmetry, to reduce the computational complexity of the problem.Using these operators, we have developed skeleton-like visualizations which provide concise and intuitive representations of tensor fields. The separating curve and local extrema of differential operators form a sparse skeleton representation of the field which is still much information about the global structure of the field. These methods for analyzing and visualizing tensor fields hold promise for simplifying these rich and complex datasets. In contrast, direct visualization of tensor fields may obscure critical points due to the resulting visual clutter.In future work, we plan to improve the performance of the Helmholtz decomposition by exploiting CPU or GPU parallelism, and by using local (rather than global) optimization methods, such as moving least squares.



---

*Source: 142923-2011-07-02.xml*