Continuous Histograms for Anisotropy of 2D Symmetric Piece-wise Linear Tensor Fields

The analysis of contours of scalar fields plays an important role in visualization. For example the contour tree and contour statistics can be used as a means for interaction and filtering or as signatures. In the context of tensor field analysis, such methods are also interesting for the analysis of derived scalar invariants. While there are standard algorithms to compute and analyze contours, they are not directly applicable to tensor invariants when using component-wise tensor interpolation. In this chapter we present an accurate derivation of the contour spectrum for invariants with quadratic behavior computed from two-dimensional piece-wise linear tensor fields. For this work, we are mostly motivated by a consistent treatment of the anisotropy field, which plays an important role as stability measure for tensor field topology. We show that it is possible to derive an analytical expression for the distribution of the invariant values in this setting, which is exemplary given for the anisotropy in all details. Our derivation is based on a topological sub-division of the mesh in triangles that exhibit a monotonic behavior. This triangulation can also directly be used to compute the accurate contour tree with standard algorithms. We compare the results to a na\"ive approach based on linear interpolation on the original mesh or the subdivision.


Introduction
Contours or isosurfaces of scalar fields play a central role in visualization. We will only use the term contour in this chapter. There is a large body of work centered around contour computation and analysis with many applications. A prominent example is the contour tree, which provides a structural overview of the data. Complementary information is provided by some kind of contour statistics, sometimes also referred to as continuous histogram, which is an extension of the discrete histogram to encode the distribution of the scalar values of continuous functions. It can serve as a valuable quantitative signature [3] of a data set. Both structures are frequently used for interaction and filtering, for example, to determine interesting iso-values or for transfer function design [19]. Histograms have also been applied for tensor field visualization and exploration to display the distribution of scalar invariants derived from tensor fields [13]. The join tree, a subset of the contour tree, of the anisotropy builds the foundation for a stability measure of topological features in the tensor field [18]. There are many standard algorithms to compute and analyze contours mostly based on piece-wise linear or multi-linear behavior. This assumption is, however, often violated for tensor invariants when using component-wise tensor interpolation. In this chapter, we present the derivation of an exact closed-form expression for continuous histograms for quadratic tensor invariants that is consistent with a piece-wise linear interpolation of a tensor filed given over a two-dimensional triangulated domain. Thereby we are especially interested in the anisotropy, an important characteristic in many applications, defined as the difference between major and minor eigenvalue.
Our main contributions are: Framework for the analysis of scalar invariants that is consistent with a piece-wise linear interpolation of the tensor components. This includes topologically correct contours, the correct contour tree, and the exact contour histograms. A generic closed-form formulation of histograms for continuous data based on a generalization of the cumulative histogram. This approach neither requires explicit contour length computations nor the computation of the gradient. Explicit solution for the histogram for piece-wise quadratic functions with positive Hessian over a two-dimensional domain. The anisotropy over a piece-wise linear component interpolation provides a relevant example. Comparison to naïve approaches using a linear interpolation of the anisotropy and a discussion of the results.
In the remainder of this section, we first summarize the development of the concept of histograms from a signature of discrete data to continuous fields. Then we summarize some aspects of tensor interpolation and invariants, that motivated our work.

Continuous histograms and contour statistics
Originally histograms have been developed for discrete data. They go back to Pearson and refer to a graphical representation displaying the frequency of the data items as bars over the data range [10]. As such a histogram provides a summary of the data and they are also used to compare and characterize data sets. Formally the histogram for a discrete set of data values F = {fi|i ∈ 1, · · · , m} is a one-dimensional bar chart displaying the frequency of each distinct element f ∈ F by a bar with height where δ is the delta function which is equal to one if f = fi and zero otherwise. A related concept is the cumulative histogram, applicable to ordered data, counts the cumulative number of data values smaller than a given value f ∈ F . This results in values Today, the term is often not only used for the graphical representation but also for the concept capturing the distribution of function values by binning the function range and counting the frequency of data samples within these bins. The resulting histograms, however, are strongly dependent on the binning size and the sampling strategy of the function in the domain. The continuous nature of the function between the sample points is not represented. For these reasons several concepts to extend histograms to continuous functions have been proposed resulting in a distribution of the function values. Bajaj et al. [3] have introduced the contour spectrum as a data-signature to find interesting iso-values for volume rendering. The contour spectrum assigns the geometric measures of the contour length, respectively surface area, to each scalar value. In addition, they consider areas and volumes of regions below or above a given iso-value. They also propose a method to compute these values exactly under the assumption of a piece-wise linear interpolation and a piece-wise constant gradient. Later Carr et al. [5] investigated the relation between histograms and the contour spectrum based on a nearest-neighbor interpolation. Scheidegger et al. [16] completed this work and introduced a natural generalization of histograms to the continuous setting as the distribution given by the contour spectrum weighted by the "local isosurface density" expressed by the inverse gradient magnitude. Given a scalar field f and scalar value h the distribution is given as For the derivation of the distributions, they refer to Coarea formula that provides a relationship between the sum of area integrals and a global volume integral as used by Mullen et al. [14]. For the computation of the distribution, they propose an approximation based on the count of active cells weighted by the inverse of the cell span, the difference between max and min per cell, to approximate the gradient. For piece-wise linear function, this yields a more or less good approximation of the real distribution.
Similar concepts have also be considered for multivariate fields, for example, generalizations of scatterplotts [1] and parallel coordinates [9] to continuous data. In both cases, the mathematical model is based on mass conversation of a mapping from the function range to the scatterplot space, respective parallel coordinates space. The continuous histogram is a special case considering a one-dimensional range. While the mathematical model is generic, the proposed computational algorithms are limited to specific interpolation models in the spatial domain. The algorithm by Bachthaler et al. [1] assumes a linear interpolation in tetrahedral cells. The work by Heinrich et al. [9] introduces a splatting approach to generate the distributions by computing footprints for Gaussian input kernels. These limitations have also been discussed by Bachthaler et al. [2]. They propose an adaptive refinement approach that does not require a linear interpolation. The algorithm, however, still assumes that the interpolation is convex which is not the case for most tensor invariants based on a component-wise tensor interpolation.
In our work, we focus on invariants with a quadratic behavior and propose an exact solution pursuing a slightly different approach than most of the previous methods. Instead of directly generalizing the concept of histograms we start with the generalization of the cumulative histogram. This approach requires neither surface area or contour length nor gradient approximations and is thus much more flexible in terms of requirements for the input functions and interpolation schemes. The continuous scatterplot or contour distribution follows than as a derivative of the cumulative histogram.

Notes on tensor field interpolation
In this section, we briefly summarize some notes related to tensor field interpolation without going much into detail. This has been a frequently discussed topic in many applications and it is agreed that this is a challenging topic.
While many theories and visualization models assume data given as fields on continuous domains, the data-reality are data sets sampled at discrete locations. Depending on the origin of the data this can be voxel-based data from imaging or meshes for simulation data. In any case, it is necessary to approximate and reconstruct the field from these samples. Thereby the most commonly used methods for tensor fields are based on a component-wise interpolation. Similar to other applications and other data types these are linear, bi-, or tri-linear interpolation depending on the grid type. Concerning such interpolations, a frequently discussed topic is the non-linear dependency of most tensor invariants on the tensor components [12]. There have been many methods proposed which explicitly try to preserve certain tensor characteristics [4,17]. Recently there appeared a survey over interpolation methods for positive definite tensors [8].
We will, however, stay with the most simple interpolation methods, linear component-wise interpolation within tetrahedra, a method that is also often the basis for finite element simulations. The non-linear dependence of the tensor invariants can lead to a non-convex behavior of the invariants inside the cells. An example, we are especially interested in, is the behavior of the anisotropy, defined as the difference between the maximum and minimum eigenvalue, which is typically decreased inside the cells. Sometimes this is referred to as "swelling-artifact", but we rather argue that this is a fundamental property of tensor behavior. It reflects the fact that the anisotropy is linked to the directional behavior of the tensor. This is expressed by the fact that the anisotropy is always zero at degenerate points and can directly be used to measure the stability of the degenerate points [11]. Therefore the analysis of the anisotropy field must be handled consistently with the chosen interpolation method for the tensor field. This concerns the computation of contours, the contour tree and also histograms.
For our consideration in this work, we assume that we have a continuous tensor field given on a two-dimensional triangulated domain. We further assume a piece-wise linear component-wise linear interpolation of the tensor field within these triangles. The anisotropy has then a quadratic behavior in the triangle. One can observe similar behavior for other invariants, for example, the determinant too. All the above-proposed methods for histogram computations fail when the interpolation of the data is not convex. Similarly, typical contour or contour tree computation is also based on convex interpolation. Figure 1 For the scalar field on the left, the join tree, split tree and the contour tree are shown.

Contour trees, a topological summary of scalar functions
The contour tree keeps track of the changes of sub-and superlevel sets of a function and thus provides a valuable summary of the function. It can be assembled form the join tree tracing the changes in contours when the function value is increased from −∞ to ∞ and the split tree one where the function value is decreased from ∞ to −∞. In more detail, the join tree tracks the creation and merging of sublevel sets, which are recorded in a tree. At the local minima of the function, new branches are created. As the function value increases, branches are extended and merged at saddle points where two sublevel sets merge. The global maximum of the function becomes the root of the tree. The split tree is constructed accordingly. The most commonly used algorithm for contour tree computation assumes a piecewise linear interpolation of the scalar field where all critical points are located at the cell vertices [6]. In the case of the anisotropy, we are especially interested in the join tree since many of the minima correspond to the degenerate points in the tensor field. The join tree can be used to quantify the stability of these points [18]. For correct results an accurate join tree computation based on the quadratic behavior or the anisotropy is essential.

Problem statement and Solution Overview
Given a tensor field sampled over a 2D triangular mesh, our goal is to compute the accurate contour-tree and histogram of the tensor anisotropy consistent with linear interpolation of tensor components. Anisotropy is a quadratic function under this interpolation. Figure 2 gives an overview of our solution approach. For each triangle in the input mesh, we first determine the coefficients of the quadratic function for the anisotropy based on tensor values at the triangle vertices. We show that anisotropy is a special quadratic function which either has a single minimum equal to zero and elliptical contours, or in a degenerate case, minima along a line and the contours are pairs of parallel lines, refer to Section 4. The first step is the subdivision of the triangle into monotonous sub-triangles, which is the basis for the topologically correct contour extraction, the correct contour tree computation and the derivation of the histogram. We obtain the histogram as the derivative of the cumulative histogram, which is computed first from the area of the sub-level set for the isovalues. To compute the explicit area of the sub-level sets we apply a linear transformation to the coordinate system such that elliptical contours become circular contours centered at the origin, Section 4.1.

Histogram
Cumulative Histogram Input triangle with tensors at the vertices Quadratic function for anisotropy After translation

After rotation
After scaling Devision into monotonous triangles Figure 2 Solution overview: First, the coefficients of the quadratic function are determined. Then a transformation is applied such that the minimum is at the origin and the contours are circular, followed by a sub-division into monotonous triangles. For each monotonous triangle, the cumulative anisotropy histogram is computed using sub-level set areas, which are added to get the cumulative histogram for the original triangle. Finally, the derivative of the cumulative histogram yields the anisotropy histogram.

Background and Notations
In this section, we provide the notations and a brief mathematical background required for the proposed method discussed in sections later. First, in section 3.1 we set up the notations for tensors and the invariant of interest to us that is tensor anisotropy. Then in section 3.2, we describe the barycentric interpolation within a triangle, since this approach is used for linear interpolation of tensors within a triangle. Lastly, in section 3.3, we discuss the general bi-variate quadratic function, its critical point, and the shape of its contours.

Second order symmetric tensors and anisotropy
For a second order symmetric tensor, T = e f f g , the two eigenvalues are: Depending on the application different measures for anisotropy are used. For positive definite tensors, relative measures like the fractional anisotropy are common. In mechanical engineering a typical measure is the difference between the maximum and the minimum eigenvalues α(T ) = λ − µ = (e + g) 2 − 4(eg − f 2 ). It has been shown that this value is also related to the stability of degenerate points in tensor field topology [11] and is the measure we are mostly interested in. In the following, we will however consider the squared value of α(T ), which has the same topological characteristics but simplifies the computations a lot.

Barycentric coordinates and piece-wise linear interpolation
For all our considerations in this paper we use a component wise linear interpolation of the tensor field. We use barycentric coordinates for interpolation in the triangle. Since we require an explicit representation of the the tensor field for the derivation of the exact histogram we briefly review the main definitions and formulas in this section.
For a triangle with vertices p1 = (x1, y1), p2 = (x2, y2) and p3 = (x3, y3), the barycentric coordinates (β1, β2, β3) of an arbitrary point p = (x, y) within a triangle are given by We assume that some scalar function s is sampled at the vertices of the triangle with scalar values s1, s2 and s3 at vertices p1, p2 and p3. Then, the scalar value s at an arbitrary point p = (x, y) within the triangle can be determined as This function s(x, y) is linear in x and y, and can be alternatively written as

Bivariate quadratic functions and their critical points
As the anisotropy, as defined in equation (1), is a quadratic function we summarize some general facts about quadratic functions and define the notation that we will use later in the paper. A general quadratic function in two variables can be written as s(x, y) = Ax 2 + Bxy + Cy 2 + Dx + Ey + F or as quadratic form as The critical point of the scalar function s is a point where the gradient ∇s(x, y) = 0 is zero. The partial derivatives of s(x, y) with respect to the two variables is: The location of the critical point pc = (xc, yc) of s can be obtained after solving the linear equations ∂s/∂x = 0 and ∂s/∂y = 0 using equations (3) and (4). The resulting coordinates are The function s and its critical point can be classified based on the sign of the determinant of the Hessian, H If H > 0, then the critical point is either a maximum or minimum and the contours are ellipses, compare Figure 4(a,b). The type of the critical point depends on the sign of A and C. If A, C > 0, then the critical point is a minimum, otherwise it is a maximum. In the other case when H < 0, then the critical point is a saddle and the contours are hyperbolas.
For the case when H = 0, there are no critical points and the contours are either parabolas, parallel or coincident lines. The type of the contour when H = 0 depends on the invariant I defined as If I = 0, the contours are parabolic, compare Figure 4

Anisotropy for 2D piece-wise linear tensor fields
In this section, we have a closer look at the anisotropy function for a linear tensor field. The main observation is in the non-degenerate case, that the anisotropy always has exactly one minimum with a value of zero. This corresponds to the fact that the tensor field always has exactly one critical point. All contour lines are ellipses.
Given a triangle with vertices p1 = (x1, y1), p2 = (x2, y2) and p3 = (x3, y3). Let the tensors at these vertices be T1, T2 and T3, respectively. Let the components of the tensors be The barycentric coordinates can be used for linear interpolation of the tensor field within the triangle. Tensor T at an arbitrary point p = (x, y) within the triangle can be found as: As described earlier in Section 3.2, the tensor components e, f and g linearly interpolated within the triangle as can be written as f (x, y) = fxx + fyy + fc, g(x, y) = gxx + gyy + gc.
The function, we are interested in is the anisotropy ν, the explicit expression for which can be obtained by substituting the linear expressions for tensor components in equation (1) yielding Comparing equation (13) with the general quadratic equation (2), we obtain the coefficients of the quadratic function in dependence of the tensor components: Substituting the above equations in equation (6), we derive the determinant of Hessian for the anisotropy to be (for derivation refer to Appendix A.1): From equation (20), we can conclude that the contours of anisotropy function are never hyperbolic. In most cases, when H is strictly greater than 0, the contours are elliptical. Moreover, in that scenario, from equations (14) and (16), we can deduce that the type of the critical point is always a minimum.
Appendix A.1 contains the detailed analysis of ν. We show that the bi-variate quadratic function ν is a special function that has elliptical contours or in some cases a pair of parallel lines, ruling out the possibility of hyperbolic or parabolic contours.

Field normalization using coordinate transformations
In the following, we derive a coordinate transformation such that the bivariate quadratic scalar function determined for a triangle has a standard format. This allows for the application of a unified strategy for further analysis.
Equation (13) gives the expression of anisotropy in general bivariate quadratic form. We have also established that contours of this function will are elliptical. Therefore we first transform the coordinates such that elliptical contours are centered at the origin and aligned to the axes. This can be achieved by applying a translation and rotation, both rigid-body transformations preserving areas. Then, a transformation is applied such that contours of the bivariate quadratic function become circles rather than ellipses. This can be achieved by applying a non-uniform scaling, a linear transformation which distorts the area by a constant factor given by the determinant of the transformation matrix.
After translation such that the minimum pc = (xc, yc) falls into the origin the anisotropy becomes with the translated coordinates xt = x − xc, yt = y − yc and Ft = F . Using equations (14)- (19), it can be shown that Ft = 0 for ν. This implies that minimum value of ν is zero at the critical point and this point is a degenerate point in the tensor field. Equation The area distortion because of this scaling is given by the factor √ λ1λ2.

Subdivision in monotonous sub-triangles
We consider a triangle with vertices p1 = (x1, y1), p2 = (x2, y2) and p3 = (x3, y3) and the tensors at these vertices are T1, T2 and T3. The tensors are linearly interpolated within the triangle and the anisotropy is a quadratic function. For the extraction of contours, the computation of the contour tree as well as for the correct histogram we require piece-wise monotonous behavior inside the triangles, which is given in this chapter. Similar, subdivisions have also been proposed by Dillard et al. [7] and by Nucha et al. [15]. There are five different cases depending on the location of the global minima and the local minima at the edges. Especially the cases when the minimum pc lies within the triangle or outside the triangle need a different treatment for the computation of the anisotropy histogram.

Case a
The minimum pc is within the triangle In this case, the input triangle can be partitioned into six sub-triangles such that ν behaves monotonously within the triangle. See Figure 6b.

Case b
The minimum pc is outside the triangle Here, we have four possibilities: 1. No triangle edge has an edge minimum. See Figure 6i.    Figure 6c. In all these cases, the triangle can be sub-divided into an appropriate number of monotonous triangles. Although these sub-divided triangles are monotonous and look similar to the first case, an important difference is that in general none of the vertices lies at the origin.

6
Computation of the histogram for ν In this section, we derive the exact continuous histogram of the anisotropy for linearly interpolated tensor fields. We use the term histogram here not strictly as it has been originally introduced but in the sense of the distribution of function values. It can also be considered as the weighted contour spectrum using the terminology of Bajaj et al. [3] as introduced by Scheidegger et al. [16]. In contrast to previous methods, we approach the problem by first computing the cumulative histogram CH and then derive the histogram from it by computing its derivative. This makes the exact computation of the histogram much more feasible.
Where Areai(ν ≤ ν0) is the area of the sublevel set in triangle i. Specifically, we consider here the anisotropy with its quadratic behavior, however the method also directly applies for linear fields. In our derivation of Areai(ν ≤ ν0) we consider the following setting. Given is a triangle ABC with vertices A = p1 = (x1, y1), B = p2 = (x2, y2) and C = p3 = (x3, y3) and respective tensors T1, T2 and T3, which are linearly interpolated within the triangle. Further we ordered the vertices such that the anisotropy values are ν1 < ν2 < ν3. To obtain the contribution of the triangle to the cumulative histogram at the value ν we compute the area of the sublevel set at ν within the triangle. In the following, we assume that anisotropy is monotonous within the triangles resulting from the subdivision introduced in Section 5. We further assume that the transformations as described in Section 4.1 have been applied. This means that the anisotropy has the form as given in equation (25) with a global minimum pc = (0, 0) at the origin and the level sets that are circles. The area within the transformed triangle is distorted by a constant factor, which can be appropriately multiplied to get the exact areas for the original triangle. We consider two cases depending on whether the vertex A of the triangle vertex lies at the origin or not.

Case 1: The global minimum p c lies at one triangle vertex
Let's assume that the global minimum pc lies in vertex A = p1 of the triangle ABC. Then the shape of the sublevel set for ν can have two different types depending on whether ν is smaller or larger than ν2 = r 2 AB . If ν ≤ ν2 the sublevel set is a sector of a circle with radius r = √ ν and opening angle θA, as shown in Figure 7. The area of the sublevel set for ν ∈ [ν1, ν2] is then given by Where θA is angle between the edges AB and AC.
The rate of change of the area is given by which is a constant not depending on the specific value of ν. If the isovalue ν is greater than ν2 and less than ν3, the contour is composed of a circular sector and an additional more complex shape DBF E as illustrated in Figure 7. This region is enclosed by two circular segments DB and EF and two line segments DE and BF and can be computed as With r 2 = ν and r 2 AB = ν2 this results in In the expression above, the given value of ν determines the values of θ.

Case 2: The global minimum p c does not lie at any triangle vertex
These triangles are the more generic case, we still assume that the anisotropy is monotonous inside the triangle and the vertices such that the anisotropy values are ν1 < ν2 < ν3. The triangles look similar to Case 1, but none of the vertices lie at the origin, compare Figure 8. To compute the area of a sublevel set ABF E within these triangles, we considering the sublevel sets in two triangles where we can use the algorithm of Case 1. In Figure 8 they are the triangles DBC and DHC. For the final area, we have in addition to consider the area of triangle AHB which has to be added or subtracted depending on the exact position of A.

A brief note on implementation
We compute the cumulative histogram of anisotropy by computing the sublevel set area. The detailed algorithm is provided in Algorithm 1. For a given tensor mesh, we pre-process all the triangles to determine the coefficients of the quadratic function for anisotropy within the triangle, compare Section 4. Each triangle is appropriately subdivided into monotonous triangles, compare Section 5. Then for each anisotropy value in the histogram bin, we compute the sublevel set area by adding up the sublevel set area in all monotonous triangles. We also handle triangles with degenerate minimum appropriately.

Parallelization
The algorithm described in Algorithm 1 has lots of options for parallelization. The loops at Lines 4 and 38 are embarrassingly parallel because each triangle in the mesh can be pre-processed in parallel. Similarly, each bin in the histogram can be computed in parallel. Further, with atomic add operations, the computation of the sublevel set area for anisotropy value ν corresponding to a particular histogram bin can be parallelized over the set of monotonous triangles. These are the loops at Lines 41 and 44 in Algorithm 1.

Numerical issues
Since there are a few floating-point operations and checks involved in the Algorithm 1, floating-point errors can be introduced which may propagate to the output resulting in a wrong histogram. We handle these errors by introducing reasonable checks, for example making sure that the cumulative histogram thus obtained is always a monotonically increasing function and sums up to the total area of the domain for the largest value of anisotropy. However, a deeper study of the errors and more robust computation is left for future work. We are confident that using a multi-precision library for computation will remove the floating-point errors.

Results
We apply the proposed approach of computing anisotropy histograms to three different case studies. We show results for synthetic, simulation and experimental data. For all the case studies we compute the continuous histograms for three different interpolations of the anisotropy.

Interpolation [a]:
Linear interpolation of the anisotropy within the original triangulation.

Interpolation [b]:
Linear interpolation of the anisotropy within the sub-divided monotonous triangles, compare Section 5.

Interpolation [c]:
Anisotropy based on the linear interpolation of the tensor components, resulting in a quadratic behavior or the anisotropy, compare Section 4. The resulting anisotropy fields are directly shown using a color map where dark blue refers to an anisotropy value of zero and red assigned to the maximum value. We also show a set of contours in these images as white lines. In addition, we computed the join tree for all cases, which are always the same for interpolation [b] and [c]. Finally, we computed the continuous histogram where the results for one data set are plotted in one image, interpolation [a] displayed as a red curve, interpolation [b] as a green curve and interpolation [c] as a black curve.

Synthetic data
The first example is a synthetic data set where tensors with user specified tensor components are placed at grid locations of a 5x5 grid. This grid is triangulated to provide a mesh with 32 triangles. This simple example is well suited to demonstrate the differences between the accurate histogram . The upper row shows the resulting field as a color-map superimposed with contour lines. The ellipses in the vertices represent the tensors defining the field. The second row shows the respective join trees for the three anisotropy fields. The red square corresponds to the zeros inside of the original triangles. All zero-leaves in the tree represent degenerate points of the tensor field topology.
on the subdivided mesh with monotonous triangles [b], as shown in Figure 9. It can immediately be seen that for interpolation [a] the topology of the contours is not correct. While the exact contours differ between interpolation [b] and interpolation [c], they have the same topology. This is also reflected in the corresponding join trees, which are the same for both fields Figure 9.
In the next step, we randomly perturb the directions of the tensors at the vertices without changing their eigenvalues to generate an ensemble of four tensor fields. Note that since we do not change the eigenvalues, the anisotropy at the vertices of the mesh remains unchanged after perturbation. Hence, for interpolation [a] all the fields, the contours, the join trees, and their histograms will be the same as evident from  Figure 11. The contour tree for the interpolation [a] is the same for all 4 fields and already given in Figure 10a, so we only show the trees for the sub-division in monotonous triangles. The degenerate points of the tensor field where the anisotropy is zero appear also as minima in the join tree and are highlighted in red. The join trees provide an overview of the possible cancellations of degenerate points and thus their stability [18]. Comparing the trees it can be seen that the four data sets vary significantly for their topological structure and stability of its degenerate points. The join tree for the linear anisotropy, Figure 10a, has no zero and is thus not consistent with the direction fields given by the tensors.

Simulation data
The difference between the three different interpolations for the anisotropy field gave significantly different results for the synthetic data. In the next step, we want to investigate the impact of these differences on real-world data. At first, we consider a simulation data set from mechanical engineering. It is a well-known data set, often referred to as two-point load, that represents a stress field in a solid block resulting from the application of two external forces, Figure 12. The data is given on a cubic mesh. The anisotropy measure used so far corresponds to the squared von Mises stress, which plays an important role in failure analysis of mechanical components. The direction field in one slice is shown in Figure 12b. In our analysis, we consider a section of this slice, which is shown in Figure 12c displaying one eigenvector field. Color represents the anisotropy using a bi-linear interpolation in the original mesh. In our analysis, we use a triangulated version of the data. Figure 13 shows the anisotropy fields using the three different interpolations in comparison. In Figure 13a one can observe the asymmetry introduced by the triangulation but is similar to 12c. (1) (2) (3) (4) Figure 11 Ensemble of four tensor datasets (1)-(4) generated by random variation of the eigenvector directions at the vertices. The join trees summarize the hierarchical structure of the minima in the anisotropy field and the degenerate points in the tensor filed topology. Degenerate points inside the triangles are marked as red nodes in the tree. The anisotropy fields in 13(b) and (c) are both based on the monotonous subdivision and are very similar to each other but differ strongly from 13a. The asymmetry due to the triangulation is largely reduced. The minima inside the original triangles in the field capture the locations of the degenerate points of the direction field. The corresponding histogram can be seen in Figure 13(d) and (e). As expected they differ strongly for very small values but are rather similar for large values of the anisotropy.

Measurement data
As the last example, we examine two sections of a slice of 3D Diffusion Tensor Imaging data. Specifically, we compute and compare the anisotropy histograms of noisy regions outside the brain and a region inside the brain. The 2D slice from the data along with selected regions is shown in Figure 14. Figure 15 shows the histograms computed for the noisy region. With the interpolation approach [c] for computing histograms, we can observe a high contribution of anisotropy values near zero, which hints at the existence of a lot of degenerate points in the noisy region. This is not captured by the histogram computed with interpolation approach [a]. It should be noted that interpolation [b] (green plot) although yields better results than interpolation [a] (red plot), it is still quite different than the correct histogram (black plot). We did the same experiment for a region that is within the brain and hence contains less noise. The results are shown in Figure 16. Here we observe, that histograms are very similar whether we use the accurate quadratic function for anisotropy of linearly interpolate it within the triangles.

Figure 14
Two segments selected from Diffusion Tensor Imaging data. The 20x20 2D grid shown in red outline is chosen from noisy region of the data while the selection shown in black outline is chosen from the region occupied by the brain. We plot the histograms for these selections in Figures  15 and 16 respectively.

Conclusions
In this paper we explore the behavior of the anisotropy, as an example for a nonlinear derived tensor invariant, when applying a linear component wise interpolation of the tensor field. We demonstrate that a linear interpolation of the invariant itself, the interpolation approach [a], leads to an incorrect topology of contours as well as a bias in the histogram. With this analysis we want to emphasize the importance of being consistent with the chosen interpolation in all analysis steps. For our analysis we have chosen a linear interpolation of tensor components, which is the most commonly used method in simulations and provides a valid continuous field. An independent interpolation of the direction field and the anisotropy violates the preservation of topological invariants and does not result in a valid tensor field. More specifically we have presented a derivation for the computation of correct contours and histograms in this setting. Component wise linear interpolation of tensor components results in a quadratic function for anisotropy. The method is based on a subdivision of the mesh into triangles with monotonous behavior. This subdivision with a linear interpolation of the anisotropy, interpolation approach [b], already results in topologically correct contours. However, the histograms are not accurate and show a bias towards larger anisotropy values. This is especially prominent in regions with many degenerate points. In areas of high anisotropy, interpolation [b] provides a good approximation. The method described in this chapter, interpolation approach [c], can be used to compute correct continuous histogram for anisotropy.
All derivations in this chapter are given for the anisotropy defined as the difference between the major and minor eigenvalue, of 2D tensor fields. Although not trivial, extension of this work to 3D tensor fields is feasible. An extension to the determinant, which is also quadratic but not elliptic, would be possible in a similar way. A general extension to other non-linear tensor invariants, however might not be possible in a closed form and will require a good approximation schema. Therefore, we plan to explore methods for efficient approximations to the correct distributions with clear error bounds. Also, computing continuous scatterplots to visualize the space of multiple invariants at once is an interesting topic and will be subject of future work.

A.1 Detailed analysis of anisotropy and its contours
Proof of H ≥ 0 Substituting the equations (14) - (19) in equation (6), we can analyse the determinant of Hessian for anisotropy ν: Since H ≥ 0, we have shown that the contours of νare never hyperbolic.

Proof of
From equation (38), we conclude that when H is 0, I is also 0. This means that the contours of ν are never parabolic.
To conclude, based on equations (30) and (38), we deduce that the contours of anisotropy ν and hence α are either ellipses or a set of parallel lines.