Abstract
The challenge of tensor field visualization is to provide simple and comprehensible representations of data which vary both directionally and 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.
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 Section 2, 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.
(a) Helix
(b) Saddle
(c) Source
(d) Spiral
(e) Vortex
3.1. Tensor Notation
The order 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 -dimensional space, each index can take one of different values. In 3 dimensions, an order- tensor then has 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,
Partial derivatives will be denoted by the symbol, where
This notation permits tensor equations to be expressed in a very compact manner. For example, the familiar vector field divergence operator in 3D which is given as in vector calculus notation, would be given as in index notation.
3.2. Tensor Field Divergence
In general, the divergence of an order tensor field is an order tensor field. For , the divergence is given in Einstein notation as This notation indicates that for all possible values of index , 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 Figure 3 and the velocity field . The mass flux into the left face is , and the flux out of the right face is , where is the material density. If mass is to be conserved, it is required that , or . For more general flows with nonzero and components, if we take the limit as the volume, , of the box approaches 0, this conservation property can be written in terms of the integral which is the definition of the divergence of .
Now consider the case of a flow field resulting from the action of a second-order tensor field on . The flux into the left side of the volume will be given by , and the flux out of the right side is given by . Conservation of mass for the one-dimensional case leads to the condition . By applying the product rule for differentiation, we get terms which depend on the rate of change of , and terms which depend on the rate of change of . The terms which depend on the rate of change of , are characterized by the divergence of .
In index notation, applying the product rule for differentiation to , The first term on the right hand side, , depends on the divergence of , and the second term depends on the gradient of . This can also be rewritten in vector calculus terms as .
In the context of diffusion, the divergence appears in Fick's second law which is a statement of conservation of mass for the diffusion process governed by concentration gradient and diffusion tensor field . 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 where is the permutation tensor The permutation tensor is often used to define the vector cross product . This leads to the vector calculus notation . The tensor curl can be analyzed just as the divergence was in (3.6). In this case, the result is that characterizes the part of which is due to the spatially varying tensor field ,
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 tensor field. For , the gradient is given by For order-0 tensors, the vector calculus notation is , 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 where 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 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 . In order to reduce the storage requirements for , we perform order reduction to obtain a lower-order tensor field which approximates . Analysis of is done in terms of its divergence, which again results in reduction of storage requirements. We have found that reduction of to its order-one (vector) approximation is useful for identifying saddle points in tensor fields.
(a)
(b)
(c)
(d)
(e)
(f) Helix
(g) Saddle
(h) Source
(i) Spiral
(j) Vortex
(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, 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 where denotes the th 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 tensor field. For a second-order tensor field, the helicity tensor is given by
Sample helicity plots are shown in the bottom row of Figure 2. We use the trace of the helicity tensor to quantify the magnitude of helicity. Extrema of 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, , is given by where is the gradient of a scalar potential field , is the curl of a vector stream field , and 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 without interference from curl-based features. Likewise, is solenoidal and is useful for isolating centers of vortices in . The harmonic vector field, , 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 Just as in the vector field case, we have and . 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 . The components are then ordered within the vector according to lexical order of indices. An input second-order tensor field, , with spatial dimensions is then vectorized as which has 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 where is an identity matrix, is the Kronecker product and, is an finite difference matrix. We use central differences for approximating derivatives, in which case is given by 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 as where is the sparse block matrix Similarly, the gradient of the first-order tensor field is given by where is the sparse block matrix 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, , is given by where is the sparse block matrix
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 where denotes the Frobenius norm of the tensor .
Using the fact that , we implement this numerically by alternately solving the normal equations 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 and 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 , the harmonic part of the field, but instead let . 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 is irrotational, is solenoidal, and that 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 Figure 5 was generated by plotting the radial surfaces for unit vectors .
(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 part
The results of the generalized Helmholtz decomposition of the tensor fields are shown in the second and third rows of Figure 5. 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 and positive-definite tensors have . 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 to or and subtract that field from 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 and positive by adding a constant isotropic tensor field to each and subtracting the isotropic tensor fields from . 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, , given by This function can be used to segment the field into two regions. When , the field is better approximated by , and when , it is better approximated by . Li et al. [6] used a similar decision criterion based on a ratio of errors for segmentation of vector fields. The contour line 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 , and and . 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.
We 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 Figures 6 and 7.
(a) Order 2
(b) Order 4
(a) Order 2
(b) Order 4
For 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 Figure 8. The red line shows the time to compute the decomposition of a second-order tensor field by forming the sparse matrices and 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 and procedurally. For the fourth-order tensor field, we do not attempt to construct the matrices since they are so large.
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 is 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.
(a) Spiral 1
(b) Spiral 2
(c) Spiral 3
(d) Spiral 4
(e) Spiral 5
(f) Spiral 6
(g) Spiral 7
Our 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 Figure 10 and then summing them, . The tensor fields in all following figures are visualized by plotting the radial surfaces for unit vectors in the second-order case, and plotting , and in the third- and fourth-order cases. The surface is colored blue when is positive and red when is negative.
| (a) |
| (b) |
| (c) |
| (d) |
Several interesting observations can be made from the results of the Helmholtz decomposition shown in Figure 11. 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.
| (a) Synthetic field |
| (b) |
| (c) |
| (d) harmonic term |
The error difference function, , 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 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 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.
(a) Order 2
(b) Order 4
The 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 where 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 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 isocontour (shown as a thick black line) varies with tensor order.
(a) Order 1
(b) Order 2
(c) Order 3
(d) Order 4
(a) Order 1
(b) Order 2
(c) Order 3
(d) Order 4
Five tensor fields of order 2 and 4 were generated by transforming and summing the fields shown in Figure 5. 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 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.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
Some spurious contours do occur, but some of these are due to boundary effects, as in Figure 18(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 Figures 18 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.
Appendices
A.
Consider the vector fields , , , , and . These represent a source, vortex, saddle, helix, and spiral, respectively. Tensor fields can be constructed by repeated contraction of these vector fields, for example, . We applied our experiments to combinations of these 5 basic fields and applied discrete numerical approximations of the differential operators to the data. In this appendix, we present a closed-form analytical analysis of the Lamb tensor (for saddle detection) and helicity tensor (for helix detection) when applied to the source, vortex, saddle, and helix. Since the differential operators used are linear operators, and are invariant to rotation, translation and scale, similar analyses can be performed on fields comprising combinations of these five basic fields.
A.1. Detecting Saddle Points
In Section 3, we describe saddle detection using the divergence of a vector field computed by reducing the order of the Lamb tensor of the original tensor field. The expression for each component of the Lamb vector for second-order tensor fields is given in Table 3. Expressions for the divergence of the reduced Lamb tensor for second-, third- and fourth-order fields are given in Tables 4, 5, and 6, respectively. Note that in these tables, the constants through are positive. For the second-order case, it is trivial to differentiate the saddle from the other features since the divergence is zero everywhere for the saddle, and a quadratic polynomial in all other cases.
For the order 3 and 4 fields, the divergence is no longer zero, but the saddle can be differentiated from other features by concavity of the divergence, as demonstrated in Figure 20. Therefore, we locate and classify saddle points as local maxima of the Lamb tensor divergence.
(a) Order 3 Source
(b) Order 3 Vortex
(c) Order 3 Saddle
(d) Order 4 Source
(e) Order 4 Vortex
(f) Order 4 Saddle
A.2. Detecting Helices
Expressions for the trace of helicity of the second-, third- and fourth-order fields are given in Tables 4, 5, and 6, respectively. In all cases, it is trivial to differentiate the helix from other features since is only nonzero when is helical.