Methods

Hue-balls

The first part of our method uses color to visually distinguish regions exhibiting different diffusion tensors. It is based on a simple way to visualize unit-length vectors: a colormap on the sphere that gives a continuous mapping from direction to color. A closely analogous approach in vector visualization used a two dimensional hue-saturation colormap on the sphere to visualize perturbation velocity in an application of direct volume rendering to computational fluid dynamics [19]. What turns this into a tensor visualization technique is the application of the tensor matrix as a linear transform. At all locations in the tensor field, a single user-specified input vector is multiplied by the diffusion tensor matrix to create an output vector. The tensor is assigned color by using the direction of the output vector as the lookup into the spherical colormap. We use the term hue-ball to describe a spherical colormap used in this way. One motivation for the hue-ball is that throughout a region of the tensor field with high spatial coherence, multiplying by the tensor will tend to give the same result, and the assigned color will be nearly uniform. Discerning coherent structures in the tensor field becomes a task of visually detecting color coherence in the rendered image.

The properties of the tensor's matrix representation are important for understanding how the hue-ball functions. It is useful to consider the input vector in the basis formed by the eigenvectors. Given a tensor matrix $\mathrm{M}$ with unit-length eigenvectors $\mathbf{e}_1$, $\mathbf{e}_2$, and $\mathbf{e}_3$, an input vector $\mathbf{v}$ can be expressed as

$$\mathbf{v} = (\mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_1\hspace{-1p... ...pace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_3\hspace{-1pt}) \hspace{1pt}\mathbf{e}_3$$

$$= ( \mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_1\hspace{-1... ...pace{-1pt}, \mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_3\hspace{-1pt})$$ (1)

Then the output vector $\mathrm{M} \mathbf{v}$ can be expressed as

$$\mathrm{M} \mathbf{v} = \mathrm{M} (\mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_1... ...pace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_3\hspace{-1pt}) \hspace{1pt}\mathbf{e}_3$$

$$= \lambda_1 (\mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_1\... ...pace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_3\hspace{-1pt}) \hspace{1pt}\mathbf{e}_3$$

$$= ( \lambda_1 \mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_1... ..., \lambda_3 \mathbf{v}\hspace{-1pt}\cdot\hspace{-1pt}\mathbf{e}_3\hspace{-1pt})$$ (2)

where $\lambda_i$ is the eigenvalue corresponding to eigenvector $\mathbf{e}_i$. The coordinates of the output vector in the eigenvector basis are the input vector's coordinates, scaled by the corresponding eigenvalues. We term the change in direction between the input and output vectors the deflection caused by the tensor.

Equation 2 indicates that the vector is always deflected towards the principal eigenvector, since the coordinate of the input vector in the principal eigenvector direction will by definition grow proportionally larger than the components along the other eigenvectors. There is also a relationship between the amount of deflection and the tensor's anisotropy. Because the anisotropy of a tensor is in general related to the disparity among its three eigenvalues, multiplying a vector by a tensor with high anisotropy will cause a greater relative change among its coordinates, and hence a greater deflection. But there is a limit to how much the vector can be deflected. Since the diffusion tensor matrix has non-negative eigenvalues, multiplying by it can not change the sign of any of the vector's coordinates. Both the input and output vectors will be in the same octant of the eigenvector basis, so the angle between input and output vectors cannot exceed 90 degrees.

Figure 2: The diffusion tensor deflects the input vector by an amount related to its orientation and anisotropy

Figure 2 illustrates this deflection with two sequences of tensors, represented with their ellipsoidal form. Within each row, the anisotropy of the tensor is fixed, but the orientation is changing gradually. The upper row exhibits more anisotropy than the lower, but the orientations are matched between the two rows. The input vector is always pointing straight downward, above the ellipsoid, and the output vector is below the ellipsoid, pointing away from it. The deflection is always towards the principal eigenvector, but there is less deflection in the less anisotropic case. This illustrates the hue-ball's ability to de-emphasize isotropic regions (having low anisotropy). If the input vector points to a darker or more neutral color than the rest of the hue-ball, then isotropic regions will be inconspicuous, since in these regions there is little deflection of the input vector. Where the tensor field is more anisotropic, multiplying by the tensor has the potential to deflect the vector more, and the hue-ball can make these regions more noticeable.

Figure 3: Comparison within one dataset slice of the amount of deflection (a) from multiplication by the tensor matrix, with anisotropy (b) calculated from eigenvalues.

Figure 3 shows how well the amount of deflection can correspond to anisotropy in real data. Figure 3(a) was generated with a gray-scale hue-ball which was black at one pole, with brightness increasing linearly with angular deviation from the pole. The input vector points into the page, towards the black pole. Figure 3(b) shows the result of evaluating a specific anisotropy measure (called the ``anisotropy index'') everywhere in the same slice of data. While the details of the anisotropy calculation are given in Section 3.2, the relevant point is that the anisotropy is calculated from the eigenvalues of the tensor at each point in the field. This requires finding the roots of the tensor's cubic characteristic polynomial, demanding much more computation than is needed to do the simple matrix multiplication used for the hue-ball. Yet, the results are comparable. Both have a large bright ``X'' shape in the center, indicating the corpus collosum (connecting the two sides of the brain), and other finer structures also appear in both images.

One important difference, however, is the dark seam which extends down the middle of the corpus collosum in Figure 3(a). This occurs when the input vector is perpendicular to the tensor's principal eigenvector, becoming aligned with either of the tensor's minor eigenvectors. In this configuration, there is no deflection towards the principal eigenvector. Using a different input vector can move the problem to another region of the dataset, but does not eliminate it. Also, if the tensor's major eigenvector is aligned with the input vector, here again there will be little deflection. These are basic limits on how well the deflection of the hue-ball's input vector can approximate the anisotropy.

The free parameters in the hue-ball method of tensor visualization are the color assignment on the sphere, and the input vector to use for multiplication with the diffusion matrix. For the sake of simplicity we have used only the hue-ball mapping shown in Figure 8(a) (colorplate). The sphere has a band of saturated colors around its equator, with saturation decreasing to the top and bottom poles, which are a medium gray. All the colors have the same ``lightness'' in the HSL color space [8], since for the sake of coloring a tensor field, it is visually less confusing if the hue-ball varies only color, letting the shading model control intensity². All the hues appear twice on the hue-ball so as to create 180 degree rotational symmetry. As seen in Figure 2, the input vector can be deflected towards either end of the ellipsoid, depending on its orientation. The rotational symmetry insures that in either case the same hue is assigned.

In the same way that lights illuminating a scene can either be fixed in world or view coordinates, the hue-ball input vector can either be fixed relative to the dataset or relative to the viewpoint. If it is fixed relative to the dataset, changing the viewpoint allows one to inspect a single hue-ball mapping from different angles. When the hue-ball input vector is fixed in view space, the element of interactivity or motion can provide additional cues about the directionally dependent nature of the diffusion tensor. Image interpretation can be simplified by consistently orienting the hue-ball so that the input vector points to its neutral pole, and so that some hue is always aligned with an ``up'' direction. Regions in the volume which didn't significantly deflect the input vector will stay a neutral color, and a consistent mapping between color and deflection direction is maintained.

To illustrate how the hue-ball colors measured tensor data, the same dataset slice which was shown in Figures 1 and 3 has been mapped by the HSL hue-ball described above and is shown in Figure 8(b) (colorplate). Some previous techniques for colormapping diffusion tensor data assign color based on the direction of the principle eigenvector, and then modulate the color by the anisotropy, so that isotropic regions do not stand out. Using an appropriately chosen hue-ball mapping, with the input vector pointing to a neutral color, this happens automatically, so no anisotropy calculation is needed.

Barycentric Opacity Mapping

Another basic ingredient in direct volume rendering is an opacity assignment that makes uninteresting regions transparent so they do not occlude important structures. We term such a mapping from a diffusion tensor to opacity an opacity function. In our application, the objects of interest (white matter fiber tracts) are by definition very anisotropic, so high opacity should be assigned to them. Conversely, low or no opacity should be assigned to the gray matter on the exterior of the brain. The literature provides various metrics for anisotropy based on the tensor matrix's three sorted eigenvalues $\lambda_1 \geq \lambda_2 \geq \lambda_3$ [13,20]. We have chosen to use the ones by Westin et al. [20] due to the simple geometric motivation behind them. Metrics for three different kinds of anisotropy are given:

$$c_l = \frac{\lambda_1 - \lambda_2}{\lambda_1 + \lambda_2 + \lambda_3} \hspace{16pt}\raisebox{-15pt}{\epsfig{file=eps/c_l.eps, width=40pt}}$$ (3)

$$c_p = \frac{2(\lambda_2 - \lambda_3)}{\lambda_1 + \lambda_2 + \lambda_3... ...e{-30pt} \hspace{15pt}\raisebox{-10.5pt}{\epsfig{file=eps/c_p.eps, width=40pt}}$$ (4)

$$c_s = \frac{3\lambda_3}{\lambda_1 + \lambda_2 + \lambda_3} \vspace{-30pt} \hspace{15pt}\raisebox{-15pt}{\epsfig{file=eps/c_s.eps, width=40pt}}$$ (5)

It can be shown that all the metrics fall in the range $[0,1]$, and that they sum to unity: $c_l + c_p + c_s = 1$. The ellipsoids drawn next to the anisotropy metrics indicate the shape of diffusion tensor for which that metric will be high; it will be near zero for the other two shapes. Where only $c_l$ is high, the tensor field is said to be linearly anisotropic, where only $c_p$ is high, the tensor field is planarly anisotropic. The last metric, $c_s$ is actually for isotropy; $c_s = 1$ only when all the eigenvalues are equal. Therefore, a single anisotropy metric called the ``anisotropy index'' is defined as:

$$c_a = 1 - c_s = c_l + c_p = \frac{\lambda_1 + \lambda_2 - 2\lambda_3} {\lambda_1 + \lambda_2 + \lambda_3}$$ (6)

Figure 4: Different anisotropies measured on slice.

To see how the anisotropy can vary in measured data, Figure 4 shows the metrics $c_l$, $c_p$, and $c_a$ evaluated over the same dataset slice seen in previous figures, with brighter areas indicating higher anisotropy.

These anisotropy measures provide a simple way of assigning opacity to diffusion tensors for direct volume rendering. A function whose domain and range are both $[0,1]$ can map the values from an anisotropy measure (such as $c_l$) to opacity, assuming that the chosen anisotropy measure is high exactly in those regions that are deemed important for visualization.

A more flexible approach uses multiple anisotropy measures to facilitate more complex opacity functions. In light of the normalization built into $c_l$, $c_p$, and $c_s$, we propose the use of barycentric coordinates to depict the space of possible anisotropies, as shown in

Figure 5: Barycentric space of anisotropies

Figure 5. For every point in the triangle, there is a corresponding ellipsoid for which the anisotropy measures ($c_l$, $c_p$, and $c_s$) evaluate to the point's barycentric coordinates. In the figure, the three ellipsoids accompanying the corners of the triangle are representative of the ellipsoids which correspond to those corners. At each vertex of the triangle, one of the anisotropy measures is one, while the two others are both zero. Along the sides of the triangle, one of the anisotropy measures is zero, and the other two measures sum to one.

Figure 6: Examples of barycentric opacity maps and resulting volumes

The same two-dimensional barycentric space of anisotropies can serve as the domain of an opacity function for direct volume rendering a diffusion tensor field. This class of opacity functions can accentuate the range of anisotropy types in a very flexible way, while including the previously described opacity functions (which depend on a single anisotropy metric) as special cases. Though space does not permit an exploration of the idea, we also note here that the barycentric space can also serve as the domain of a color mapping (``transfer function'') applied to the tensor field, allowed different anisotropies to be discerned by varying their color.

Figure 6 demonstrates some barycentric opacity maps. Each opacity map is depicted by gray-scale representation: brighter regions in the triangle correspond to higher opacity assignment. For the purposes of this figure, the effect of the opacity map is demonstrated by applying the map to the the tensor dataset, resulting in a scalar volume of opacity values. This new scalar volume is visualized with a linear opacity function, and shaded according to the gradient of opacity values. A single white light is coming from the upper right. One can see that analogous to Figure 4, appropriately chosen opacity functions allow one to see the form of structures in the dataset which have one predominant type of anisotropy. It is these structures which will be colored by the action of the hue-ball, and illuminated according to lit-tensors.

Lit-Tensors

Streamlines used in vector visualization are sometimes hard to interpret because they lack the shading cues which we are accustomed to seeing on surfaces. However, illuminated streamlines (``lit-lines'') have shading and highlights that give information about their direction and curvature, creating the appearance of shiny filaments [22]. In the case of diffusion tensor visualization, we have made tensors opaque based on their anisotropy, but no other information about the tensor is indicated. We aim to show that a means of illuminating diffusion tensors, called lit-tensors, can help disambiguate the type and orientation of anisotropy. The constraints we followed for designing such a scheme are:

1.

In regions of complete linear anisotropy, the lighting model should be identical to that of illuminated streamlines. Complete linear anisotropy means movement by diffusion is constrained to one dimension, so it is sensible for the lighting model to degenerate to one already developed for vector fields.

2.

In regions of complete planar anisotropy, the lighting model should be the same as with traditional surface rendering. The obvious choice for the ``surface normal'' for a planar anisotropic tensor is the third eigenvector, perpendicular to the plane formed by the span of the first two eigenvectors (associated with the largest two eigenvalues).

3.

There has to be a smooth interpolation between these two extremes. Since tensor data can exhibit a wide variety of anisotropies, allowing small variations in anisotropy to lead to large changes in shading will probably create a more confusing image.

This can be seen as a problem of how to interpolate illumination between different codimensions. The codimension of the diffusion tensor's representative ellipsoid is two in the linear anisotropy case, and one with planar anisotropy. Previous work [1] has rigorously developed illumination methods for general manifold dimension and codimension, but did not cover cases part-way between different codimensions. Unlike that work, no claim to physical accuracy or plausibility is made for the model presented here; it is just one simple way of satisfying the constraints above.

We take as our starting point the Blinn-Phong lighting model [3]:

$$I = I_{\mathrm{ambient}} + I_{\mathrm{diffuse}} + I_{\mathrm{specular}}$$

$$= k_\mathrm{a} A_\lambda O_\lambda + I_\lambda (k_\mathrm{d} O_\lambda \mathbf{L}\cdot\mathbf{N} + k_\mathrm{s} (\mathbf{H}\cdot\mathbf{N})^n)$$ (7)

$k_\mathrm{a}$, $k_\mathrm{d}$, and $k_\mathrm{s}$ control the contributions of ambient, diffuse, and specular reflection to the final image. Following Foley et al. [8], we add the subscript $\lambda$ to those variables which vary according to color. For example, there are separate values $I_\mathrm{r}$, $I_\mathrm{g}$, $I_\mathrm{b}$, for the red, green, and blue components of the directional light source. The ambient light color is $A_\lambda$. Instead of representing the intrinsic object color with different $k_\mathrm{a}$ and $k_\mathrm{d}$ for red, green, and blue, we use $O_\lambda$ for object color and keep $k_\mathrm{a}$ and $k_\mathrm{d}$ as separate controls. In our case, the intrinsic object color is determined by the output of the hue-ball. $\mathbf{L}$ is the vector pointing towards the directional light source, $\mathbf{V}$ points towards the eye, and $\mathbf{N}$ is the surface normal. Note that instead of using $(\mathbf{V}\cdot\mathbf{R})^n$ for the specular component, where $\mathbf{R}$ is the reflection of $\mathbf{L}$ across $\mathbf{N}$, we are using the ``half-way'' vector $\mathbf{H}$ in $(\mathbf{H}\cdot\mathbf{N})^n$. $\mathbf{H}$ is the normalized average of $\mathbf{L}$ and $\mathbf{V}$, and $n$ is the shininess exponent.

Because a streamline is one-dimensional, at any given point along it there is an infinite set of normals, all perpendicular to the tangent direction $\mathbf{T}$, radiating outwards in a circle. If naively using Equation 7 to illuminate a streamline, one must find the normal which is in the plane spanned by $\mathbf{L}$ and $\mathbf{T}$ to evaluate $\mathbf{L}\cdot\mathbf{N}$. Similarly, another specific normal must be found to evaluate $\mathbf{H}\cdot\mathbf{N}$. The insight which makes lit-lines simple is that one does not need to actually find a specific normal in order to evaluate a dot product with it. With Pythagoras's theorem, the dot product can be expressed in terms of the tangent $\mathbf{T}$:

$$\mathbf{U}\cdot\mathbf{N} = \sqrt{1 - (\mathbf{U}\cdot\mathbf{T})^2}$$ (8)

where $\mathbf{U}$ is either $\mathbf{L}$ or $\mathbf{H}$, for the diffuse and specular terms, respectively.

The relevant property of Equation 8 is that the lighting calculation depends on a vector (the tangent) which gives the object's direction, instead of a depending on a normal. The direction and orientation of a diffusion tensor is determined by not one, but two vectors: the first and second eigenvectors³. Both of these could be interpreted as tangents, but their relative importance is determined by the magnitudes of the corresponding eigenvalues. To characterize the relative importance of the first two eigenvectors in determining the tensor's orientation, we introduce a parameter $c_\theta$. Assuming that the eigenvalues are ordered $\lambda_1 \geq \lambda_2 \geq \lambda_3$, we define

$$c_\theta = \frac{\pi}{2}\frac{c_p}{c_a}= \frac{\pi(\lambda_2 - \lambda_3)} {\lambda_1 + \lambda_2 - 2\lambda_3}$$ (9)

As anisotropy varies from completely linear ($c_l = 1$; $c_p = 0$) to completely planar ($c_l = 0$; $c_p = 1$), $c_\theta$ varies from 0 to $\frac{\pi}{2}$. The role of $c_\theta$ is to control how much the second eigenvector contributes to the lighting of the diffusion tensor. In the linear case, only the first eigenvector determines the tensor orientation, and in the planar case, both the first and second eigenvectors matter equally.

The expression to be used in lieu of dot products with $\mathbf{N}$ is:

$$\sqrt{1 - (\mathbf{U}\cdot\mathbf{e}_1)^2 - (\mathbf{U}\cdot\mathbf{e}_2\,\sin(c_\theta))^2}$$ (10)

Note that in the case of linear anisotropy, $\sin(c_\theta) = \sin(0) = 0$, so the contribution from $\mathbf{e}_2$ vanishes, and the expression reduces to the formula for lit-lines (Equation 8), with the principal eigenvector $\mathbf{e}_1$ taking the role of the tangent $\mathbf{T}$. This is appropriate, since in linear anisotropy, the principal eigenvector points in the direction of movement, as does a streamline's tangent vector.

In planar anisotropy, $\sin(c_\theta) = \sin(\frac{\pi}{2}) = 1$, and the contributions of the two dot products are equal. This means that for any other vector $\mathbf{W}$ such that

$$(\mathbf{W}\cdot\mathbf{e}_1)^2 + (\mathbf{W}\cdot\mathbf{e}_2)^2 = (\mathbf{U}\cdot\mathbf{e}_1)^2 + (\mathbf{U}\cdot\mathbf{e}_2)^2$$ (11)

Equation 10 will have the same value. Therefore, in planar anisotropy the lighting model is rotationally symmetric around $\mathbf{e}_3$. Rotational symmetry in this case is actually an important feature of the lighting model. In planar anisotropy, the diffusion tensor ellipsoid degenerates to a disc, and any vector in the plane spanned by the disc is an eigenvector. Because of this numerical instability, the calculated directions of the first and second eigenvectors will be essentially random. The illumination should not be sensitive to this arbitrary orientation, and should only be a function of the third eigenvector. In fact, one can use Pythagoras's theorem to show that if $c_\theta = \frac{\pi}{2}$, Equation 10 gives an exact formula for $\mathbf{U}\cdot\mathbf{e}_3$. Interpreting both $\mathbf{e}_1$ and $\mathbf{e}_2$ as surface tangents, then the surface normal $\mathbf{N}$ is aligned along $\mathbf{e}_3$. Therefore the model contains standard surface shading as a special case.

Figure 7: Sequence of volumes of differing anisotropy, rendered with lit-tensors. Anisotropy varies gradually between the nine volumes, going in scanline order.

To demonstrate lit-tensors, Figure 7 shows nine different synthetic diffusion tensor datasets which were direct volume rendered with a fixed viewpoint and light. The anisotropy index $c_a$ of the sphere is also constant in every case, but $c_\theta$ is changing. The dataset in the upper left has complete linear anisotropy in a concentric circular pattern (along lines of latitude). The dataset in the middle has complete planar anisotropy (and hence looks just like a standard surface rendering). The dataset in the lower left has complete linear anisotropy along lines of longitude, going from pole to pole. The images provide a convincing sense of surface anisotropy, which is not a typical trait in direct volume renderings. Figure 8(c) on the colorplate shows the same sequence of volumes colored with a hue-ball.

Footnotes

... intensity²

If one is seeking a truly constant luminance colormap, HSL colorspace is too simplistic.

... eigenvectors³

Because the eigenvectors always form an orthogonal basis, and because we are adopting two-sided lighting, the third eigenvector does not contribute any additional information.