Scientific Computing · Scientific Visualization

Tensor Fields and the Ambiguity of Eigenvector Glyphs

How a symmetric rank-two field is decomposed into three directions, why those directions have no sign, and what happens to the picture at isotropic and planar degeneracies.

A symmetric second-order tensor field assigns a $3 \times 3$ matrix to every point, and the matrix has six independent components. Diffusion imaging, stress analysis, and strain computations all produce fields of this kind, and none of them is visualizable by the means that work for scalars: there is no single number per point to color, and there is no single direction per point to draw an arrow along. The matrix describes how a quantity varies with direction, and every visualization of it is a projection of that six-dimensional object onto two or three perceivable dimensions.

The projection with the most structure behind it is the eigendecomposition. A real symmetric matrix has real eigenvalues and an orthonormal set of eigenvectors, so the field can be described at each point by three orthogonal directions with associated magnitudes. That description is a coordinate frame that follows the data, and it is what makes eigenvalue-derived glyphs and line integral convolution of eigenvector fields meaningful.

It is also where the trouble starts, because an eigenvector is defined only up to sign, and at degenerate points the frame is not defined at all.

What the decomposition gives and what it refuses to give

For a symmetric matrix $\mathbf{T}$, the eigendecomposition $\mathbf{T} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\mathsf{T}}$ exists with $\mathbf{Q}$ orthogonal and $\boldsymbol{\Lambda}$ diagonal with real entries. Sorting the eigenvalues, $\lambda_1 \ge \lambda_2 \ge \lambda_3$, gives three labelled directions: the major eigenvector for $\lambda_1$, the medium, and the minor.

Each eigenvector is a solution of $(\mathbf{T} - \lambda \mathbf{I})\mathbf{v} = 0$, and if $\mathbf{v}$ is a solution then $-\mathbf{v}$ is too. The magnitude of the eigenvalue is meaningful; the sign of the eigenvector is not, because the tensor is a quadratic form and $f(\mathbf{v}) = \mathbf{v}^{\mathsf{T}}\mathbf{T}\mathbf{v}$ depends only on the line through the origin. There is no physical measurement that distinguishes $\mathbf{v}$ from $-\mathbf{v}$, so a visualization that draws signed directions is drawing a choice made by the decomposition routine, not a property of the data.

The consequence for glyph-based rendering is that a glyph built from the three eigenvectors has four equally valid orientations, from flipping any subset of the axes, and a glyph built from them and colored by eigenvector sign will show a discontinuity across a surface where the routine’s arbitrary choice flipped. Any sign-dependent coloring therefore has to be post-processed into a consistent frame: the standard technique is to walk the field in a fixed order and flip each eigenvector whenever it makes an obtuse angle with the previously processed neighbour. The result is continuous except along surfaces where the flipped and unflipped choices are equally consistent, and those surfaces are real features of the field rather than artefacts. They cannot be removed by more careful bookkeeping.

Degeneracies are the second refusal. Where $\lambda_1 = \lambda_2 \neq \lambda_3$, the major eigenvector is any direction in the plane perpendicular to the minor eigenvector, so the eigenvector field has no unique value on a whole plane of directions. Where all three eigenvalues coincide, as they do in an isotropic region, the tensor is a multiple of the identity and every direction is an eigenvector: the frame is entirely undefined. A glyph that draws three axes in such a region is drawing an arbitrary rotation.

The practical response is to describe the tensor by measures that remain continuous where the eigenvectors do not. Westin’s shape measures interpret the eigenvalues as weights on three idealized geometries: a linear case ($\lambda_1 \gg \lambda_2 \approx \lambda_3$), a planar case ($\lambda_1 \approx \lambda_2 \gg \lambda_3$), and a spherical case ($\lambda_1 \approx \lambda_2 \approx \lambda_3$). Presented as barycentric coordinates, they form a continuous map into an equilateral triangle, and the triangle is the natural legend for them. The anisotropy fraction, usually written as fractional anisotropy, is a scalar summary,

$$\mathrm{FA} = \sqrt{\tfrac{3}{2}},\frac{\sqrt{(\lambda_1 - \bar{\lambda})^2 + (\lambda_2 - \bar{\lambda})^2 + (\lambda_3 - \bar{\lambda})^2}}{\sqrt{\lambda_1^2 + \lambda_2^2 + \lambda_3^2}},$$

which is zero for isotropy and approaches one as the tensor becomes rank one. It is invariant to the ordering of the eigenvalues and to their sign, so it is defined everywhere, including at the points where the frame is not.

Glyph design as a projection problem

The classical glyph is an ellipsoid whose semi-axis lengths are derived from the eigenvalues. For a diffusion tensor the quantity being visualized is the diffusivity along a direction, which is the inverse of the quadratic form, and the convention is therefore to draw a surface whose radius along the eigenvector direction is inversely proportional to the square root of the eigenvalue. The convention is a choice, not a mathematical identity, and a figure that changes it without saying so conveys the opposite of the intended anisotropy. A tensor that diffuses quickly in one direction draws a prolate glyph in every convention, but the aspect ratio is not preserved between them.

Ellipsoids have the advantage of encoding the full eigenvalue triple with no additional channels and the disadvantage of an ambiguous silhouette. Viewed from a direction that makes the projection of a prolate ellipsoid nearly circular, the shape loses the information that motivated it. Superquadrics, constructed as generalized spheres with two exponent parameters, sharpen the shape so that linear, planar, and spherical cases are distinguishable by silhouette at the cost of a shape whose axes are no longer literally proportional to the eigenvalues.

Glyph density is a separate design axis. A tensor field sampled on a grid has as many glyphs as voxels, and their projections overlap. Rendering a subset, chosen on a regular stride or by placing glyphs until a minimum screen-space separation is met, trades coverage for legibility, and the same trade appears in every dense glyph visualization. Where the field is noisy, glyph orientation inherits the noise directly: the eigenvector of a nearly isotropic tensor is determined by the small differences between nearly equal eigenvalues, so the glyph’s orientation is essentially random in those regions while its shape is nearly a sphere. This pairing of a meaningless orientation with a stable shape is the reason the shape measures exist as the primary encoding.

import numpy as np

def decompose(T, previous=None):
    """T: (..., 3, 3) symmetric tensor field. Returns eigenvalues (sorted desc)
    and eigenvectors with the sign chosen consistently where possible."""
    w, v = np.linalg.eigh(T)                 # ascending eigenvalues; v's columns are vectors
    w = w[..., ::-1]                         # descending
    v = v[..., :, ::-1]

    if previous is not None:
        # choose the sign that keeps the field continuous with the previous sample
        dots = np.sign(np.sum(v * previous, axis=-2))
        v = v * np.where(dots == 0, 1.0, dots)[..., None, :]

    mean = w.sum(axis=-1, keepdims=True) / 3.0
    num = np.sqrt(((w - mean) ** 2).sum(axis=-1))
    den = np.sqrt((w ** 2).sum(axis=-1))
    fa = np.sqrt(1.5) * np.where(den > 0, num / np.where(den > 0, den, 1.0), 0.0)
    return w, v, fa

np.linalg.eigh is the right routine rather than a general eigensolver: it is restricted to symmetric matrices, and for those it returns real eigenvalues and an orthonormal basis by construction, so no complex arithmetic and no separate deflation step is needed.

The sign-matching step is applied only against one neighbour because the problem is not solvable globally: on any closed loop where the frame rotates by a half turn, no consistent assignment exists, exactly as with a Möbius strip. The loop’s existence is a topological property of the field, and the residual sign discontinuity marks it rather than being eliminated by it.

Eigenvector fields as direction fields

Once a consistent sign convention is fixed locally, the major eigenvector defines a direction field, and everything that works on direction fields becomes available. Line integral convolution produces a dense texture whose streaks follow the field, which shows structure everywhere rather than at glyph positions. The visualization is then a scalar image rather than a set of objects, and its properties are those of a texture: the kernel length sets the streak length, a very long kernel washes out variation across the field, and a very short one gives noise.

Hyperstreamlines take a different approach, tracing the path of the major eigenvector while sweeping a cross-section whose shape is determined by the other two eigenvalues. The cross-section is what makes the method tensor-specific: a field that is locally planar produces a flat ribbon, so the twist and the change of aspect ratio along the line carry the medium and minor information. This is the analogue of a stream surface, and it shares the same failure modes, self-intersection in the projection and arbitrary termination when the line leaves the domain.

Degenerate points are where these methods need explicit handling. At an isotropic point, the major eigenvector direction changes discontinuously in every neighbourhood, so a line traced through it changes direction arbitrarily. The standard way to present this is to treat degenerate points as the points where lines must not be seeded, and to render their neighbourhoods through the shape measures instead. The classification of degenerate points in a direction field, by the number of lines entering and leaving, mirrors the classification of critical points in a vector field, and the same index arithmetic applies.

The depth of the projection

Every choice in this article is lossy. Glyphs at sampled positions reduce a continuous field to a set of markers; the shape measures reduce three eigenvalues to two independent numbers; fractional anisotropy reduces them to one; and the sign convention reduces an eigenvector field to something that can be drawn. None of these is a simplification of convenience, and the information that each discards is recoverable by the reader only if the figure states which projection was used.

The color channel in these images is itself subject to constraints that are easy to violate when the color is doing classification rather than encoding magnitude. Ternary shape measures mapped to hue run the risk of a rainbow ramp with no perceptual ordering, as described in Perceptually Uniform Colormaps and the Failure of Rainbow Ramps. When the tensor field is used to steer an integration rather than to be displayed directly, as in fibre tracking, the interpolated tensor at non-grid positions is a reconstruction with the same error structure as any other field, discussed in Scattered Data Interpolation and the Conditioning of Radial Basis Functions, and the tracker’s paths are integral curves of that reconstruction rather than of the measurement: the case analysed in Streamline Integration in Vector Fields and Where the Error Accumulates.

References

  1. Superquadric Tensor Glyphs — Gordon Kindlmann, Eurographics/IEEE Symposium on Visualization, 2004
  2. Visualizing Second-Order Tensor Fields with Hyperstreamlines — Delmarcelle & Hesselink, IEEE Computer Graphics and Applications, 1993
  3. Processing and Visualization for Diffusion Tensor MRI — Westin, Maier, Mamata, Nabavi, Jolesz & Kikinis, Medical Image Analysis, 2002
  4. MR Diffusion Tensor Spectroscopy and Imaging — Basser, Mattiello & LeBihan, Biophysical Journal, 1994