Topological analysis for 3D real, symmetric second-order tensor fields using Deviatoric Eigenvalue Wheel

Computers & Graphics 54 (2016) 28–37

Contents lists available at ScienceDirect

Computers & Graphics journal homepage: www.elsevier.com/locate/cag

Special Issue on CAD/Graphics 2015

Topological analysis for 3D real, symmetric second-order tensor fields using Deviatoric Eigenvalue Wheel Xianyong Liu, Lizhuang Ma n School of Computer Science and Technology, Shanghai Jiao Tong University, Shanghai, China

art ic l e i nf o

a b s t r a c t

Article history: Received 22 March 2015 Received in revised form 11 July 2015 Accepted 12 July 2015 Available online 21 July 2015

Tensor fields are of particular interest not solely because they represent a wide variety of physical phenomena but also because crucial importance can be inferred from vector and scalar fields in terms of the gradient and Hessian, respectively. This work presents our key insights into the topological structure of 3D real, symmetric second-order tensor fields, of which the central goal is to show a novel computation model we call Deviatoric Eigenvalue Wheel. Based on the computation model, we show how the eigenanalysis for tensor fields can be simplified and put forward a new categorizer for topological lines. Both findings outperform existing approaches as our computation depends on the tensor entries only. We finally make a strict mathematical proof and draw a conclusion of R ¼ K cos ð3αÞ, wherein R is the tensor invariant of determinant, K is a tensor constant, and α is called Nickalls angle. This conclusion allows users to better understand how the feature lines are formed and how the space is divided by topological structures. We test the effectiveness of our findings with real and analytic tensor fields as well as simulation vector fields. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Tensor fields Visualization Topology Interaction technique

1. Introduction Visualizing 3D real, symmetric second-order tensor fields has found uses in various engineering contexts. The stress tensors for bone implantation in orthopaedics [1], the strain tensors for exhibiting isotropic responses in materials [2], the real, symmetric traceless tensors for simulating molecules in Nematic Liquid Crystals [3], and the diffusion tensors for monitoring the white matter structures in Magnetic Resonance Imaging [4,5] are a few examples. Differing from conventional means such as Hyperstreamlines [6], Tensorlines [7], HyperLIC [8], HOTlines [9], and Superquadric glyphs [10], our work centers around the topological analysis of 3D real, symmetric second-order tensor fields (referred to here simply as “tensor fields”, unless otherwise stated). Despite recent progress towards topology-based visualization of tensor fields [11–14,5,15–17], the study aiming at bettering our understanding about topological structures and features of underlying tensor fields remains a challenge while being meaningful. The reasons are threefold. First, in scientific data visualization the topology-based techniques, pursing the idea of showing less can be more, are proven to yield simplified, yet compact, depictions. Second, trained users can even reconstruct the data fields by

n

Corresponding author. E-mail address: [email protected] (L. Ma).

http://dx.doi.org/10.1016/j.cag.2015.07.009 0097-8493/& 2015 Elsevier Ltd. All rights reserved.

looking at the topological structures in most cases. Third, the comprehension to scalar/vector fields can be improved through the topological analysis of tensor fields. For example, by applying a data transformation, vector and scalar fields are mapped into tensor fields with respect to the gradient and Hessian, respectively. As a result, many appealing representations for the datum, which are hardly speculated from the preimages directly, can be inferred from the images of tensor fields. As shown in Fig. 1 (right), where the simulated 3D dynamical Lorenz attractor system is visualized using our approach in the transformed tensor field, the degenerate curves are symmetric and behave like the pivots of two vortices, which the flow trajectories are around. Early work on topology-based methods for visualizing tensor fields has laid an important background of this research project. Nevertheless, none of those techniques are flawless. It was showed that in 3D tensor fields double degeneracies are stable and form degenerate curves that could be detected by the “Extracting-andTracing” framework, however a fact that deviatoric (tensor) fields are topology-preserving was not sufficiently appreciated in [19,14,13]. Ref. [5] borrowed the result from [2] to regard degenerate lines of a tensor field as crease lines of the tensor invariant mode, unfortunately a strict mathematical proof to the result was not given by both. Moreover, the topological analysis about the essence of deviatoric fields was not explored either. Recent work of [17] went deeper and gave many interesting results, such as the traceless surface and the neutral surface, however its proposal of the eigenvalue manifold was not as

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

29

Fig. 1. Left: the Lorenz attractor (a.k.a. the Lorenz butterfly) of a simulated 3D dynamical system that exhibits chaotic flow [18] (rho ¼ 28; sigma ¼ 10; beta ¼ 8:0=3:0; t i ¼ ½0; 25). Right: The topological lines of the tensor field transformed from the flow field (ti ¼ 0, and the glyphs describe the vector field direction).

straightforward as the eigenvalue wheel in [20,5]. Additionally, the approach for detecting degenerate curves in [17] not only was timeconsuming but also encountered undesirable zig-zags. The central goal of this paper is to show a novel computation model we call Deviatoric Eigenvalue Wheel (DEW, demonstrated in Fig. 3), with which this paper makes the following contributions. Eigenanalysis: Conventional eigendecomposition has been widely used in visualization, yet it is expensive and sensitive to numerical computations. The presented eigenanalysis is faster and more robust, as its computation relies on the tensor entries only (see Section 8.1). Categorization: The categorization of feature lines accounts for reflecting how the fields are divided by topological structures. The concept of eigen difference was employed by [19,14,13], however interpolating colors failed to inform the real type of feature lines. Regarding that triple degeneracies in tensor fields are rare and unstable, [5] and [17] suggested to categorize topological lines by applying the extremal values of the mode and the LP Factor, respectively. In the paper, we put forward a new categorizer that depends on the sign of tensor determinant. Tensor invariant: Tensor invariants are closely related to analyzing the topological features of tensor fields. Ref. [5] pointed out that the degenerate curves are actually the crease (ridge and valley) lines of the tensor invariant mode, where the computations took advantage of the result in [2] while lacking a strict mathematical proof. Such a pitfall is overcome in this paper.

2. Previous work This section mentions the previous work highly relevant to this paper. For general tensor fields, we refer the interested readers to [21,6,22,15,16]. The salient topological features in 2D tensors fields are degenerate points [23], i.e., trisectors and wedges, whose tensor indices are half-integers due to the indeterminacy of eigendirections that is not seen in regular vector fields. The understanding of topology of 3D tensor fields remained vague until [24] provided some useful findings, e.g., the “double degeneracy” and “triple degeneracy”. Unfortunately, [24] did not recognize that in 3D cases the topological formalisms are continuous lines instead of isolated points. Rewardingly, using the theory of dimensional analysis, Zheng and Pang [19] proved that the triple degeneracy in 3D tensor fields is actually unstable and typically absent from practical datasets, conversely the double degeneracy is stable and constitutes one dimensional degenerate lines (a.k.a. degenerate curves, feature lines,

or topological lines). The local behaviors nearby degenerate points were unclear until [14] attempted this issue and experimentally showed that transition points are exactly the points where the degenerate tensors switch types between trisectors and wedges. While figuring out analytic equations for the transition points is still an open problem in tensor field analysis. Another local behavior is the tangent along degenerate curves which was fully discussed in [13] by using the Taylor's expansion of the well-known tensor discriminant, showing that the tangent around a degenerate point is just the direction from which the discriminant's Hessian vanishes, so do itself and its gradient. It was complained that the algorithms in [19,13] produced results that were sensitive to noise and essentially meaningless in the context of Diffusion Tensor Imaging (DTI) [4]. To overcome this difficulty, [5] changed our mindset by regarding the degenerate curves as the crease (ridge and valley) lines of the invariant mode defined as pffiffiffi P P pffiffiffiffiffi 2μ3 = μ2 3 , wherein μ1 ¼ i λi =3; μ2 ¼ i ðλi  μ1 Þ2 =3, and P μ3 ¼ i ðλi  μ1 Þ3 =3 are the central tensor moments, and λi are the tensor eigenvalues.

3. Fundamentals of tensor topology A 3D real, symmetric second-order tensor field is usually defined by a smooth tensor-valued function, associating each point in the field with a value T ¼ T ij ei  ej ð1 r i; jr 3Þ, wherein T ij are the tensor entries in a space spanned by the basis of dyads fei  ej g. In tensor field visualization, it was acknowledged that a tensor filed is fully equivalent to three eigenvector fields ν1  3 satisfying AðλÞνi ¼ 0 with three eigenvalues λ1  3 representing magnitudes. As is known, AðλÞ, i.e., λI T, is the characteristic matrix. The roots of the polynomial pðλÞ ¼ detðAðλÞÞ correspond to tensor eigenvalues that are ordered in a sequence of λ1 Z λ2 Z λ3 , in which the values λi ð1 r i r 3Þ are individually ranked as major, medium, and minor. So are the corresponding tensor eigenvectors νi ð1 r ir 3Þ. The coefficients of pðλÞ are defined as P ¼ trðTÞ, J ¼ ðtrðTÞ2  trðT2 ÞÞ=2, and R ¼ detðTÞ. They are all tensor invariants [26], among which P and R quantify the trace and determinant, respectively. Quantities of minor Q ¼ λ1 λ2 þ λ2 λ3 þλ1 λ3 , norm J T J ¼ ðλ21 þ λ22 þ λ23 Þ1=2 , and tensor discriminant D3 ¼ ðλ1  λ2 Þ2 ðλ2  λ3 Þ2 ðλ3  λ1 Þ2 are other important tensor invariants, which can be calculated via tensor entries [19]. The topological structure of tensor fields can be specified by tensor degeneracies and separatrices. Tensor degeneracies are formed by degenerate points that are the only places where hyperstreamlines can cross each other, and separatrices are the boundaries of hyperbolic sectors where trajectories sweep past degenerate points. As

30

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

Fig. 2. In 2D tensor fields degenerate points (a), including trisector and wedge marked by spheres in yellow and cyan, are physically isolated and theoretically convertible in space, doing a role similar to critical points (c) in vector fields. In 3D tensor fields degenerate points form continuous lines (b) to the contrary. (a) and (c) used the contrast correction technique in [25]. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

such, the degenerate point, of which two or more tensor eigenvalues are equal (i.e., D3 ¼ 0), is usually considered as a launching point to understanding the topology of tensor fields. The degenerate points of 2D tensor fields (see Fig. 2(a)), which are physically isolated and theoretically convertible although have different types of winding numbers, play a role similar to the critical points of 2D/3D vector fields (see Fig. 2(c)). On the contrary, in 3D tensor fields degenerate points form continuous lines. An instance is demonstrated in Fig. 2(b), where the planar (Type-P) degenerate points are red colored in the case of λ1 ¼ λ2 4 λ3 , and the linear (Type-L) degenerate points are blue colored in the case of λ1 4 λ2 ¼ λ3 . It is worth mentioning that tensor invariants are closely related to tensor topology, since they are independent of tensor “orientation” [27]. Therefore, they are often defined via tensor eigenvalues in visualization.

4. Proposal of the Deviatoric Eigenvalue Wheel Tensor decomposition [21,28,5,15,17] has been proven to be effective in tensor field visualization. For instance, a tensor value T can be uniquely decomposed into a sum of an isotropic part A and a deviatoric part D as 0 T 00 þ T 11 þ T 22 B T ¼ AþD ¼ B @

0

0

0

T 00 þ T 11 þ T 22 3

0

0

0

T 00 þ T 11 þ T 22 3

0 2T 00  T 11  T 22 B þB @

C C A

3

T 01

T 02

T 01

2T 11  T 00  T 22 3

T 12

T 02

T 12

2T 22  T 00  T 11 3

trðTÞ I 3

0

1

λ1

C A;

λ2 λ3

then it follows that Θ  1 TΘ ¼ Λ. Thus, we have Θ  1 DΘ ¼ Λ trðTÞ 3 I, where I is the identity tensor.□

1 C C: A

Based on the acknowledgement that the pattern of tensor “shape” [27] is determined by the hyperstreamlines following a specific eigenvector, it has been well recognized that deviatoric (tensor) fields defined via the mapping τ1 : T⟼D; s:t ., D ¼ T 

and B Λ¼@

1

3

Fig. 3. The Deviatoric Eigenvalue Wheel (DEW) computation model is distinguishable from the eigenvalue wheel model used in [5] by three aspects: (1) the origin O of the wheel is translated from μ1 to 0; (2) the radius r of the wheel is scaled from pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi 2μ2 to 2=3; (3) the medium eigenvalue λ2, whose range is highlighted by the purple bar, functions the same as the mode m in identifying tensor anisotropy. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

ð1Þ

are essentially topology-preserving [28,5]. In other words, the eigenspace of a tensor field T is the same as its deviatoric field D. To be complete, a proof is provided as follows. Proof. Suppose that Tνi ¼ λi νi where νi and λi are the eigenvectors and the eigenvalues, respectively. We name matrices Θ ¼ ðν1 ν2 ν3 Þ

The computation model of Deviatoric Eigenvalue Wheel (DEW, demonstrated in Fig. 3) is designed to be a two-phase process that is similar to [2,17]. In the first phase, given an input tensor field T the tensor transformation formulated by Eq. [1] is done. Thus, in the deviatoric field D we have P ¼ trðDÞ ¼ 0:

ð2Þ

Unlike the previous eigenvalue wheel model in [5] where the origin O of the eigenvalue wheel was μ1, in DEW the origin is translated to 0. In the second phase, we normalize the field D via the mapping as τ2 : D⟼Φ; s:t ., Φ ¼

D : JDJ

ð3Þ

Therefore, the pffiffiffiffiffiffiffiffiin DEWp ffiffiffiffiffiffiffiffiradius r of the eigenvalue wheel is rescaled from 2μ2 [5] to 2=3. As a result, the point here is to calculate the norm of tensor values in the deviatoric field. According to the proof

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

aforementioned, it follows that

0  λ1  P  3 B 1 1 B J Θ DΘJ ¼ J Θ J  J D J  J ΘJ ¼  @  

λ2  P3

through first multiplying their counterparts of Φ by J D J and then adding the results to trðTÞ=3. Tensor eigenvectors: By applying the analysis above to computing tensor eigenvalues, eventually tensor eigenvectors can be solved without having eigendecompositions.

1   C C : A  λ3  P3 

Further, considering J Θ  1 J ¼ J ΘJ  1 , we draw a conclusion of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J DJ ¼ ðλ1  P3Þ2 þ ðλ2  P3Þ2 þ ðλ3  P3Þ2 . Rearranging the equation of the conclusion, it holds that qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J DJ ¼ 23 P 2  2Q ;

31

ð4Þ

wherein P and Q are defined in Section 3 and computed based on the given tensor field T. As depicted at the top of Fig. 3, as the arm in black rotates from π=3 to 0 within the sector in the grey shadow, the wheel slides horizontally all the way from left to right along the green line. Empirically, throughout this paper the ∠α is called Nickalls angle [20] for the sake of simplicity. The cubic curve, moving along with the wheel, represents the polynomial pðλÞ defined in Section 3. The major, medium, and minor eigenvalues amount to the projection of the arm in black, blue, and deep red, respectively. The angle between any two neighboring arms, however, is a constant of 2π=3. Thus, with respect to DEW computation 2,ffiffiffiand λp 3 ffiffiffihave a range of pffiffiffithe p ffiffiffi pffiffiffi pffiffiffi model λ1, λp ½ 6=6; 6=3; ½  6=6; 6=6, and ½  6=6;  6=3, respectively, and they are subject to λ1 Zλ2 Z λ3 . In the middle of Fig. 3, the five smaller wheels are close-up shots where the black arm rotations with a step size of π=12. By observing the movements of the wheel, we find out that there exist three oneto-one and onto mappings defined as ζ i : m⟶λi ( 1 r m r  1 and 1 r i r 3), among which ζ2 is the most elegant since the range of λ2 is as symmetric as the mode m, but λ1 and λ3 are not. An important insight into the DEW computation model is that a natural relationship between the mode m and medium eigenvalue λ2 can be drawn. That is, λ2 behaves the same as m in identifying tensor anisotropy. As shown at the bottom pffiffiffi of Fig. 3, when m reaches 1 (i.e., the tensor is type-L) λ2 reaches  6=6; p analogously, when m reaches ffiffiffi 1 (i.e., the tensor is Type-P) λ2 reaches 6=6. What's more, m and λ2 reach 0 simultaneously just as the Nickalls angle reaches π=6, leading the tensor to be neutral. Such a result is consistent with what the eigenvalue manifold [17] delineated, but it looks more intuitive. Further, as we will be explaining in the next section, the eigenanalysis of tensor fields becomes much easier via the proposal of the DEW computation model.

Theorem 5.1. Let T represent a 3D real, symmetric second-order tensor field. λ is a known eigenvalue of T, and T  λI is rearranged as ðω1 ω2 ω3 ÞT . Then, the eigenvector ν corresponding to λ can be solved as ν ¼ ω1  ω2 þ ω2  ω3 þ ω3  ω1 . Proof. It follows that two of ωi must be linearly independent since T  λI is of rank two as it holds that ðω1 ω2 ω3 ÞT  ν ¼ 0.□ 5.2. Categorization of degenerate lines We put forward a new categorizer via the following theorem without offering a proof. Like the eigenanalysis, it is advantageous that the computation of our categorizer depends on the tensor entries only. Theorem 5.2. Given a 3D real, symmetric tensor field, Φ denotes the field generated by the Deviatoric Eigenvalue Wheel computation model. At a double degenerate point p the type is decided by the tensor determinant sign of its value, i.e., sgnðRðΦðpÞÞÞ, as follows: (i) R is positive if and only if p is linear; (ii) R is negative if and only if p is planar.

5.3. Tensor invariant of determinant Recalling the DEW computation model, in the output tensor field Φ we have 8 when i ¼ 1 > < λ1 ¼ r cos ðαÞ λi ¼ λ2 ¼ r cos ðα 2π=3Þ when i ¼ 2 > : λ ¼ r cos ðα þ2π=3Þ when i ¼ 3; 3

pffiffiffiffiffiffiffiffi where λi are the tensor eigenvalues, r ¼ 2=3, and α is the Nickalls angle that varies within ½0; π=3. In addition, it is known that the tensor invariant of determinant RðΦÞ is calculated as RðΦÞ ¼ λ1  λ2  λ3 ¼  ðλ2 þλ3 Þ  λ2  λ3 : Substituting the λi into the equation leads to RðαÞ ¼ r 3 cos ðαÞ  ð0:75  cos 2 ðαÞÞ:

5. Topological analysis for tensor fields 5.1. Eigenanalysis Tensor eigenvalues: Differing from [5] that suggested to compute tensor eigenvalues based on the input tensor field T directly, we propose to simplify the computation indirectly based on the output tensor field Φ ¼ τ2 ðτ1 ðTÞÞ, where τ1 and τ2 are the mappings introduced by the proposal of DEW computation model. Intuitively, it can be proven that the tensor eigenvalues of Φ are λ1 ¼ r cos ðαÞ λ2 ¼ r cos ðα  σÞ λ3 ¼ r cos ðα þ σÞ; ð5Þ pffiffiffi pffiffiffiffiffiffiffiffi wherein r ¼ 2=3, α ¼ arccosð3 6RÞ=3, and σ ¼ 2π=3. Notice that the relationship between tensor determinant R and Nickalls angle α will be discussed in depth in Section 5.3. Compared to those analytic equations in [5], ours are more straightforward. The eigenvalues of the input tensor field T, thus, can be obtained

Next, we prove that RðαÞ is a monotonically decreasing function over its support by verifying its first derivative with respect to the Nickalls angle as ∂α1 ðRðαÞÞ ¼ r 3 ½  sin 3 ðαÞ þ0:25 sin ðαÞ þ 2 sin ðαÞ  ð1  sin 2 ðαÞÞ ¼  r 3 ð2:25 sin ðαÞ  3 sin 3 ðαÞÞ ¼  r 3 sin ðαÞ  ð2:25  3 sin 2 ðαÞÞ r0: As a result, it follows that the extrema of RðαÞ have opposite signs, just on the boundaries of its support. Recalling the intermediate value theorem (a.k.a. Bolzano's theorem), it holds that RðαÞ ¼ 0 if and only if α ¼ π=6 (i.e., the tensor is neutral). We continue taking the second derivative and reach ∂α2 ðRðαÞÞ ¼ r 3 cos ðαÞ  ð0:75  cos 2 ðαÞÞ  ð  9Þ: Interestingly, we come to a conclusion of ∂α2 ðRðαÞÞ ¼  9  RðαÞ (in geometry such affine transformations including a scaling and a reflection in Euclidean space are illustrated by Fig. 4). This conclusion is a classic Ordinary Differential Equation ODE), and it can be solved analytically. The solution is RðαÞ ¼ C 1 cos ð3αÞ þ C 2 sin ð3αÞ, where C1

32

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

and C2 are constants unknown. Taking the two known conditions, i.e., RðαÞ is maximum when α ¼ 0, and the tensor pffiffiffi is neutral when α ¼ π=6, into account, it can be verified that C 1 ¼ 6=18, and C2 ¼0. Therefore, in the tensor field Φ we achieve pffiffiffi RððΦðαÞÞ ¼ 6 cos ð3αÞ=18 ¼ K cos ð3αÞ; ð6Þ where α is the Nickalls angle, and K is a tensor constant in tensor fields.

The findings of topological analysis of the fourth module have been presented in Section 5, yet they are about to be verified in Section 7. When the data sources of the first module are vector and scalar fields, the corresponding tensor fields can be generated by computing their gradient and Hessian, respectively. Taking a vector field v as an example (see Fig. 6(b)), a tensor field T is obtained as T ¼ 12 ð∇v þ ð∇vÞT Þ:

6. Implementation details Fig. 5 is an illustration of the implementation workflow that is comprised of four modules. This section focuses on the first and the third modules, and the second one has been elaborated in Section 4.

Fig. 4. The zero, first and second order derivatives of the tensor invariant of determinant with respect to the Nickalls angle elaborated in Fig. 3.

ð7Þ

Discrete tensor fields, in practice, are generated by sampling smooth tensor fields with particular data structures, such as a 3D grid. Then the values of the fields are interpolated by either linear interpolations (bilinear or trilinear), which guarantee C1 continuity only, or cubic interpolations (bicubic or tricubic), which guarantee C2 continuity but are more involved. Fig. 6 demonstrates the format of various data sources by glyphs. From now on, we start to narrate about the third module in the workflow. The “Extracting-and-Tracing” framework proposed by Zheng and others in [19,14,13] has laid the foundation for extracting the topological lines in tensor fields, upon which we are able to verify our findings about topological analysis. We, here, mostly focus on describing the modifications made to the existing framework in our implementation. Meanwhile, we refer the interested reader to the reference for more details. In the degenerate point extraction stage, our DEW computation model guarantees a criterion of T 22 ¼  T 00 T 11 , which is adopted for rewriting the discriminant constraint functions in [19]. Hence, as suggested by [14], the number of variables of the functions can be ! reduced to 5. Then, an array CF DCFs is compounded by the new constraint functions. The degenerate points to be extracted, actually, are a set of points satisfying ! D3 ðpÞ ¼ 0; s:t ., CF DCFs ðpÞ ¼ 0:

ð8Þ

Therefore, the modified Newton–Raphson method is used to assist this Least-Squares problem. In other words, given a point with an initial guess of p0 we update its coordinates by numerical iterations as !  ; ð9Þ pt þ 1 ¼ pt  ðχ T  χÞ  1 ðχ T  CF DCFs Þ X ¼ pt

Fig. 5. An illustration of the implementation workflow. Given an input data source (a), which can be a scalar, vector, or tensor field, it is first formatted (b) and then processed by the Deviatoric Eigenvalue Wheel (DEW) computation model, such that the trained users can analyze the topological features (d) of the underlying tensor field (c).

! wherein χ denotes ∂CF DCFs =∂X. We noticed that the previous extraction process in [19] was somewhat problematic because of the risk of “False Negative” (see Section 8.2). In short, the Newton–Raphson method suffers from occasionally missing degenerate points, leading discontinuities to the final feature curves in the next stage. We, accordingly, fix this issue by setting a sanity check, which serves to catch two exceptions: (1) the ðχ T  χÞ is not invertible; (2) an updated point

Fig. 6. The format of data sources used in the implementation. (a) The smooth vector/tensor fields are sampled by a 3D grid; (b) Cone glyphs: the orientation and color illustrate the vector direction and magnitude, respectively; (c) Traceless superquadric glyphs [17]: the axis of the cylinder indicates the non-degenerate eigenvector direction (green and red represent linear and planar, respectively), and the colored band intersects the glyph faces corresponding to eigenvalues with the same sign; (d) Disc glyphs: the glyphs are color-coded by medium eigenvalues, and if the values are positive discs are scaled by the major and medium eigenvalues and perpendicular to the minor eigenvector; otherwise, discs are scaled by the medium and minor eigenvalues and perpendicular to the major eigenvector. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

pt ðt Z 1Þ is out of the range of the face, on which the numerical iterations are acting. In the feature line tracing stage, the degenerate points extracted in the previous stage are assembled into a bunch of degenerate lines. A key modification we make in this stage is interpreted by Fig. 7(b).

Fig. 7. An illustration of the Neighboring-Connection algorithm [19] (a) and the Prediction-Correction algorithm [13] (b). In [13], the planes where the “PredictionCorrection” actions took place were perpendicular to the tangents during the tracing stage. Whereas in our implementation we always align the planes with an axis in space by dynamically checking collisions between the rays following the tangents and the grid faces. This change not only enables us to reuse the computation results in the previous extraction stage but also simplifies the Rayplane collision test.

33

7. Experimental results All the experimental results illustrated in this work are collected from a notebook PC with 64-bit OS, 2.4 GHz CPU, and 8 GB RAM. The first test case we carry out is related to the classical Lorenz attractor, which was first studied by Ed N. Lorenz, a meteorologist, around 1963 and derived from a simplified model of convection in the earth's atmosphere [18]. Also, it arises naturally in models of lasers and dynamos. The system is most commonly expressed by non-linear differential equations as 8 _ > < x ¼ sigma  ðy  xÞ vðx; y; zÞ ¼ y_ ¼ xðrho  zÞ  y > : z_ ¼ xy beta  z where rho ¼28, sigma¼ 10, and beta ¼ 8:0=3:0. Sometimes, rho is known as the Prandtl number, and sigma is the Rayleigh number. Given an initial particle positioned at the sphere in Fig. 1 (left), as the time ti goes from 0 to 25 the trajectory of the particle is patterned in the form of butterfly (a.k.a. the Lorenz butterfly). However, given a single time step visualizing the data without self-occlusions gets difficult as the system forms a 3D velocity vector field. As shown by Fig. 1 (right), this difficulty is overcome by visualizing the topological feature lines in the formatted tensor field. Note that the lines (both Type-P) are categorized by the theorem presented in Section 5.2. The absence of Type-L feature lines in the final results will be instructive for post-processing in

Fig. 8. An illustration of topological analysis to the Arnold–Beltrami–Childress Flows in Fluid Dynamics. Converting the normal “A : B : C ¼ 1 : 1 : 1” ABC Flows field (a) into a tensor field, the topological lines (b) are detected by [19] and categorized by Theorem 5.2. The existing algorithm [13] will benefit from our eigenanalysis in computing the tangents (c) along feature lines.

34

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

the future. For example, the flow visualization through Tensorlines [7] is suggested to choose the minor eigenvector. Next, we verify the effectiveness of the eigenanalysis presented in Section 5.1 by applying our findings to the algorithm in [13], where the tangent at a degenerate point is computed by the tensor gradient ∇T and a pair of orthogonal eigenvectors sitting in the degenerate plane. Therefore, computing the non-degenerate eigenvector perpendicular to the degenerate plane is critical. In our implementation, if a degenerate point is Type-L we compute the major eigenvector using Theorem 5.1 and treat it as the non-degenerate eigenvector, otherwise the minor eigenvector is computed. A test case (see Fig. 8(c)) is carried out upon a famous kinematic system known as Arnold–Beltrami–Childress (ABC) Flows. It is useful for studying the linear and nonlinear dynamo properties of a class of periodically forced flows [29]. Commonly, it is expressed as 8 _ > < x ¼ C  sin ðkzÞ þ B  cos ðkyÞ vðx; y; zÞ ¼ y_ ¼ A  sin ðkxÞ þ C  cos ðkzÞ > : z_ ¼ B  sin ðkyÞ þ A  cos ðkxÞ

where A, B, and C are constant coefficients, k is the wavenumber of the flow, and the domain is ½0; 2π3 . Particulary, for the normal ABC Flows (see Fig. 8(a)) whose coefficients are subject to “A : B : C ¼ 1 : 1 : 1”, the topological lines are detected in the tensor field and categorized by Theorem 5.2 (see Fig. 8(b)). We showed the results to domain experts in Fluid Dynamics. They were satisfied with the way in which the data was visualized and gave us some feedback: the topological structures of the tensor field provided pretty compact depictions. As highlighted by the small yellow circles in Fig. 8(b), the eight stagnation points, which are critical locations where the fluid velocity vanishes (see Fig. 1 in [30]), were clearly and accurately located. According to that in ABC Flows vorticity vectors are parallel to velocity vectors everywhere, those points indicate the places where the flow vorticity disappears as well. Furthermore, the categorized degenerate curves themselves are very much meaningful in physics that the planar and linear curves are markers of regions where the stretching rate of the flow is the minimum and maximum, respectively. We make another test with the double point load stress tensors, which are real-world data and commonly used as a benchmark for evaluating the effectiveness of relevant methods. Since the dataset

Fig. 9. The glyph representation of a double point load stress tensor field (left), where the topological lines (right) obtained by our approach are remarkably similar to those reported in [19,5]. Note that the two black arrows point to regions of triple degeneracies.

Fig. 10. (a) An illustration of the tensor field generated from the Sullivan vortex model [31]. As the Nickalls angle α ¼ 8π=27 the topological structures are displayed in (d), which is compounded by the feature lines (b) and the isocontour (c) of tensor determinant R.

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

tested by previous work [8,5] is not public, we turn to the open source VTK (http://www.vtk.org) that informs how to simulate single point load stress tensors such that double point load data can be obtained by applying a small displacement to the point where the first external force acts. The images in Fig. 9 display the glyph representation and topological lines of the field, respectively. The degenerate lines are categorized by the approach presented in Section 5.2. Our binary classifier makes discerning the type of feature lines as easy as [5,17] but saves computation cost. The topological structure of the stress tensor field, moreover, is physically significant in that the degenera-

35

cies are locations of pure pressure, while the two viscous stresses vanish. In addition, the linear and planar lines express tension dominant and compression dominant regions, respectively. Through our experiments, it was verified that the vertical loop, which lies directly under the planar feature line and connects the two triple degenerate points, is a stable topological feature in the sense that it persists even as the magnitudes of the two point loads are varied [19]. Finally, we start to verify our topological analysis about tensor invariant presented in Section 5.3. In the formatted tensor field Φ it has been proved that the tensor determinant R is a continuous

Fig. 11. Clockwise from top left: an illustration of the tensor field formulated in Eq. (10), the topological lines, and the tangents at degenerate points.

Fig. 12. A gallery of the deformations of isocontours of tensor determinant with respect to the Nickalls angle. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

36

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

and monotonically decreasing function with respect to the Nickalls angle α. To verify the result, we compute the isocontours of R, i.e., a series of level surfaces extracted by Marching Cubes [32]. According to the DEW computation model, as the Nickalls angle α goes to the boundaries (i.e., 0 and π=3) of its support, consequently the isocontours approach the topological lines of the tensor field. A test is taken in the flow field that simulates the Sullivan vortex model given by Sullivan (1959) [31] (see Fig. 10(a)). When α ¼ 8π=27 the level surface (see Fig. 10(c)) is composed of two patches that closely wrap the red-colored planar feature lines of the tensor field (see Fig. 10(d)). We proceed testing with an analytic tensor field (see Fig. 11), which is formulated as 0 1 10:0 0:0 x B z C T ¼ @ 0:0  20:0 ð10Þ A: x z 10:0 A gallery of the isocontours is shown in Fig. 12. When the Nickalls angle α ¼ π=27 (which is nearby the boundary of 0), the level surface consists of two “tubes” wrapping the blue-colored linear degenerate lines. As the angle increases, the surface starts to deform, intertwine, and merge gradually. Once the angle gets closer to π=6, the surface separates again, dividing the linear curves from the planar curves. As the angle is nearby the other boundary of π=3, the surface begins to shrink and eventually approaches the red-colored planar degenerate lines. Visualizing this process enables a better understanding of how the space is formed and divided by topological structures.

8. Algorithm complexity and comparison 8.1. Complexity analysis Let ðni ; nj ; nk Þ be the resolution of a 3D grid (see Fig. 6), with which tensor volume data is defined and numerical computation events are Table 1 A performance comparison between our algorithm and the state-of-the-art work [19,17] in terms of the extraction of topological features (the time is in seconds, and the symbols IL and IC denote interpolations being linear and natural cubic, respectively). Data set

Lorenz attractor ABC flows Stress tensor Sullivan Analytic field Tornado

Ours

[19]

[17]

IL

IC

IL

IC

IL

0.15 0.33 2.87 2.66 1.65 1.57

0.71 1.81 27.91 22.14 12.32 13.62

0.22 0.42 3.54 3.45 2.06 2.16

0.92 2.23 35.12 28.05 15.39 17.16

– 13.17 – 38.43 – 20.45

taken place. It is known that the number of cell faces where degenerate points to be examined equals 3ni nj nk þ ni nj þ nj nk þni nk . Hence, the complexity of nested loops in our algorithm is Θðn3 Þ with linear interpolations (Θðn7 Þ with natural cubic interpolations). Though the algorithm in [5] was of the same time complexity, it doubled the computation overhead because it extracted the ridge lines of mode associated with value þ1 as well as the valley lines of mode associated with value  1 separately. Inside a voxel, both [5] and [17] conformed to a similar rule in extracting topological lines, i.e., spatial curves can be solved by looking into the intersections between two algebraic surfaces. The efficiency of [17], however, was estimated to be poorer than [5], given that extra faces had been added through the Delaunay tetrahedralization. Table 1 shows a comparison between the presented algorithm and state-of-the-art work. It is indicated that our method achieves much higher performance than [17]. The statistics in the table, again, have evidenced that our algorithm is about 20% faster than [19] in terms of topological feature extraction. Substantial gains within the innermost loop of our implementation are attributable to the topological findings in Section 5 and the improvements and strategies in Section 6. 8.2. Experimental comparison As stated in Section 6, the sanity check is devised to tackle the challenge of “False Negative” that potentially leads discontinuities to feature lines during extracting degenerate points. Conducting natural cubic tensor interpolations is helpful in attenuating the problem to some extent (see the green boxes in Fig. 13), but it does not fix that issue satisfactorily. The sanity check, nevertheless, is rather effective in addressing the problem (see the yellow box in Fig. 13). Besides, through our investigation we notice that the algorithm in [17] not only is more time-consuming than ours but also suffers from an issue that probably is attributed to a risk called “oversegmentation” [33], which, likewise the known side-effect in image processing, tends to divide an object into different pieces that actually belong to one. In [17], the isolation processing forced the degenerate lines to go across certain faces produced by iteratively tetrahedralizing voxels. As a result, when the subdivisions were overdone the curves could be shifted away from their real positions, causing undesirable artifacts of zig-zags to the final visualization. Such a drawback, nevertheless, is not seen in our experimental results (see Fig. 14).

9. Conclusions and future work The prior work lays an important background for this project, while none of them is unflawed. The central goal of this paper is to show a novel computation model, which is called Deviatoric

Fig. 13. An illustration of comparison with the ABC Flows field, where the topological lines are extracted via linear interpolation without the sanity check (left), natural cubic interpolation without the sanity check (middle), and linear interpolation with the sanity check (right). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

X. Liu, L. Ma / Computers & Graphics 54 (2016) 28–37

Fig. 14. An illustration of comparison between our method (top) and the state-ofthe-art work [17] (bottom) with the Sullivan field. Note that in the bottom image the blue surface is the neutral surface. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Eigenvalue Wheel (DEW) and presented in Section 4. It evolves from the model of eigenvalue wheel [20,5], aiming to overcome the pitfall of the eigenvalue manifold model in [17]. Using the DEW computation model, in turn the existing “Extracting-and-Tracing” framework [19,14,13] can be improved. The eigenanalysis in Section 5.1 speeds up the computation of tangents along topological lines, and the categorizer in Section 5.2 is much easier to use. In Section 5.3, a strict mathematical proof to the result in [2,5] is given, and a relationship between the tensor determinant and the Nickalls angle is drawn to reveal how the space is divided by topological structures. Whilst the findings we present are impressive, they are not all we can see. We are investigating how the computation model we propose would facilitate our future work, such as the topological analysis of 3D asymmetric (or general) tensor fields that are more universal than the fields where we focus at the moment. Acknowledgements This work is supported by the NSFC Fund (No. 61133009) and National 973 Program of China (No. 2011CB302203). The authors would like to thank Weijie Song for his help on earlier versions. The comments by the anonymous reviewers helped to improve this paper. References [1] Dick C, Georgii J, Burgkart R, Westermann R. Stress tensor field visualization for implant planning in orthopedics. IEEE Trans Vis Comput Graph 2009;15 (6):1399–406. [2] Criscione JC, Humphrey JD, Douglas AS, Hunter WC. An invariant basis for natural strain which yields orthogonal stress response terms in isotropic hyperelasticity. J Mech Phys Solids 2000;48(12):2445–65. [3] Jankun-Kelly T, Mehta K. Superellipsoid-based, real symmetric traceless tensor glyphs motivated by nematic liquid crystal alignment visualization. IEEE Trans Vis Comput Graph 2006;12(5):1197–204.

37

[4] Schultz T, Theisel H, Seidel HP. Topological visualization of brain diffusion MRI data. IEEE Trans Vis Comput Graph 2007;13(6):1496–503. [5] Tricoche X, Kindlmann G, Westin C. Invariant crease lines for topological and structural analysis of tensor fields. IEEE Trans Vis Comput Graph 2008;14 (6):1627–34. [6] Delmarcelle T, Hesselink L. Visualizing second-order tensor fields with hyperstreamlines. Comput Graph Appl 1993;426(4):25–33. [7] Weinstein D, Kindlmann G, Lundberg E. Tensorlines: advection-diffusion based propagation through diffusion tensor fields. In: Proceedings of the conference on visualization: celebrating ten years, VIS '99; 1999. p. 249–53. [8] Zheng X, Pang A. HyperLIC. In: Proceedings of the 14th IEEE visualization 2003, VIS '03; 2003. p. 249–56. [9] Hlawitschka M, Scheuermann G. HOT-lines: tracking lines in higher order tensor fields. In: Proceedings of the conference on Visualization, 2005, VIS 05. IEEE; 2005. p. 27–34. [10] Schultz T, Kindlmann G. Superquadric glyphs for symmetric second-order tensors. IEEE Trans Vis Comput Graph 2010;16(6):1595–604. [11] Tricoche X, Garth C, Kindlmann G, Deines E, Scheuermann G, Ruetten M, Hansen C. Visualization of intricate flow structures for vortex breakdown analysis. In: Proceedings of the conference on Visualization, 2004. IEEE; 2004. p. 187–94. [12] Laramee RS, Hauser H, Zhao L, Post FH. Topology-based flow visualization, the state of the art. In: Topology-based methods in visualization; 2007. p. 1–19. [13] Zheng X, Parlett B, Pang A. Topological lines in 3D tensor fields and discriminant hessian factorization. IEEE Trans Vis Comput Graph 2005;11 (4):395–407. [14] Zheng X, Parlett B, Pang A. Topological structures of 3D tensor fields. In: Visualization, 2005, VIS 05; 2005. p. 551–8. [15] Zhang E, Yeh H, Lin Z, Laramee R. Asymmetric tensor analysis for flow visualization. IEEE Trans Vis Comput Graph 2009;15(1):106–22. [16] Chen G, Palke D, Lin Z, Yeh H, Vincent P, Laramee R, et al. Asymmetric tensor field visualization for surfaces. IEEE Trans Vis Comput Graph 2011;17 (12):1979–88. [17] Palacios J, Yeh HH-j, Wang W, Laramee RS, Vincent P, Zhang E. Robust extraction and analysis towards complete 3d tensor field topology. Technical Report. School of Electrical Engineering and Computer Science, Oregon State University; 2013. [18] Strogatz SH. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. [19] Zheng X, Pang A. Topological lines in 3D tensor fields. In: Proceedings of the conference on Visualization. IEEE; 2004. p. 313–20. [20] Nickalls RWD. A new approach to solving the cubic: Cardan's solution revealed. Math Gaz 1993;77(480):354–9. [21] Delmarcelle T, Hesselink L. Visualization of second order tensor fields and matrix data. In: IEEE Visualization; 1992. p. 316–23. [22] Zheng X, Pang A. 2D asymmetric tensor analysis. In: Proceedings of the conference on Visualization, 2005, VIS'05. IEEE; 2005. p. 3–10. [23] Delmarcelle T, Hesselink L. The topology of symmetric, second-order tensor fields, in: Pro. of the Conference on Visualization, VIS '94; 1994. p. 140–7. [24] Hesselink L, Levy Y, Lavin Y. The topology of symmetric, second-order 3D tensor fields. IEEE Trans Vis Comput Graph 1997;3(1):1–11. [25] Palacios J, Zhang E. Interactive visualization of rotational symmetry fields on surfaces. IEEE Trans Vis Comput Graph 2011;17(7):947–55. [26] Basser PJ, Pierpaoli C. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. J Magn Reson 1996;111 (3):209–19. [27] Ennis DB, Kindlmann G. Orthogonal tensor invariants and the analysis of diffusion tensor magnetic resonance images. Magn Reson Med 2006;55: 136–46. [28] Zheng X, Tricoche X, Pang A. Degenerate 3D tensors. In: Visualization and processing of tensor fields; 2006. p. 241–56. [29] Bouya I, Dormy E. Revisiting the ABC flow dynamo. Phys Fluids 2013; 25(3). [30] Dorch SBF. On the structure of the magnetic field in a kinematic ABC flow dynamo. Phys Scr 2000; 61(6). [31] Gibbon J, Fokas A, Doering C. Dynamically stretched vortices as solutions of the 3d Navier–Stokes equations. Physica D: Nonlinear Phenom 1999;132 (4):497–510. [32] Lorensen WE, Cline HE. Marching cubes: a high resolution 3D surface construction algorithm. In: SIGGRAPH '87; 1987. [33] Freixenet J, Muñoz X, Raba D, Martí J, Cufí X. Yet another survey on image segmentation: region and boundary information integration. In: Computer vision, ECCV, lecture notes in computer science, vol. 2352; 2002. p. 408–22.