Regularized Eigenvalue Derivatives

This section will assume a tensor fieldD(t)which is differentiable in a single scalart ∈R and denote its derivative with respect to t by D^′(t). From perturbation theory, it is known that for a differentiable, symmetric tensor field, there exist differentiable eigenvalue functions (cf. Chapter two in [109]). They are given by the diagonal elements of [D^′(t)]E, which is obtained by applying the rotation matrix E whose rows are the orthonormal eigenvectors of D(t)to the derivative D^′(t):

[D^′(t)]E =ED^′(t)E^T (6.1)

However, these differentiable eigenvalue functions do not produce sorted eigenvalues, and ordering them generally introduces non-differentiable cusps at pointstat which D(t)

78

6.2. Regularized Eigenvalue Derivatives is degenerate. This is illustrated in Figure 6.1: Even though eigenvalues can be described by differentiable functions (red and blue lines in Figures 6.1 (a) and (d)), sorted eigenval-ues are non-differentiable where the original functions cross (Figure 6.1 (b) and (e)). We propose to solve this problem by regularizing the derivatives to make them everywhere continuous. This is achieved by smoothly blending the derivatives near degeneracies (Fig-ure 6.1 (f)).

Within a small neighborhood of a degeneracy, the repeated eigenvalue typically splits into different eigenvalues, which constitute its λ-group [109]. Even if a repeated eigen-value itself is not differentiable, the mean eigen-value of its λ-group is (dashed black lines in Figure 6.1). From [D^′]E, this mean derivative can be extracted as the average of the diagonal entries that belong to the duplicated eigenvalue. Our regularized derivatives ˜λ^′_i are defined to preserve this mean derivative (λ˜^′₁+ ˜λ^′₂ = (λ1+λ2)^′).

In Figure 6.1 (c), the regularized derivativesλ˜^′_i have been integrated numerically to ob-tain regularized eigenvalue functions λ˜i, which behave mostly like the sorted eigenvalues in Figure 6.1 (b), but are modified locally around degeneracies to be everywhere differ-entiable. These functions are shown purely for illustration. Our implementation blends the derivatives directly and does not make use of the resulting regularized eigenvalue functions.

In order to decide at which point we should start to blend the derivatives of two eigen-values λi and λj, we introduce the measure ρij of their relative distance:

ρ_ij = |λi−λj| λi+λj

(6.2) Since we aim at diffusion tensors, this definition assumes that λi > 0. Blending starts when ρ drops below a threshold ǫ, which was set toǫ= 0.05in our experiments.

Let wij = (ρij/ǫ−1)². With this, the exact formula we used to define λ˜^′₁ in terms of the discontinuous derivatives λ^′_i is:

λ˜^′₁ =

The artifacts in Figure 6.2 (a) become more pronounced when considering edge maps of non-linear functions of sorted eigenvalues, like the Westin measures.

An alternative way to blend eigenvalue derivatives is to find a rotation E¯ such that those diagonal entries of[D^′]E¯ which correspond to a repeated eigenvalue equal itsλ-group mean derivative. This is more complex to implement, but has the advantage of cleanly separating changes in shape (diagonal elements) from changes in orientation (off-diagonal) while preserving the total derivative magnitude, as measured by the rotationally invariant Frobenius norm. The key to this method is to observe that in case of a degeneracy, we are

6. Eigenvalue Derivatives for Edge Detection in DT-MRI

(a)k∇λ¹k (without regularization) (b)k∇˜λ¹k (with regularization)

Figure 6.2.: In edge maps, ignoring degeneracies leads to small, but noticeable artifacts (a), which are fixed by the proposed regularization (b).

free to choose any set of mutually orthogonal eigenvectors e_i which span the eigenspace of the repeated eigenvalue.

Let b_i be the vectors of the assumed basis, and let D_j,k^′(1) be entry (j, k) of the matrix representation D^′(1)= [D^′]E. Then, our algorithm works as follows:

1. If ρ1,2 < ǫ, create D^′(2) by rotating D^′(1) aroundb₃ such that D^′(2)_1,1 =D^′(2)_2,2. 2. If ρ2,3 < ǫ, create D^′(3) by rotating D^′(2) aroundb₁ such that D^′(3)_2,2 =D^′(3)_3,3.

3. If ρ1,3 < ǫ, identify i such that D_i,i^′(3) is in between the remaining two diagonal entries,D^′(3)_j,j ≤D_i,i^′(3) ≤D^′(3)_k,k. If D_i,i^′(3) is larger (smaller) thanµ= 0.5·(D_j,j^′(3)+D^′(3)_k,k), rotate aroundb_k (bj) such that D^′(4)_i,i =µ. Afterwards, rotate around b_i such that D_j,j^′(5) =D^′(5)_k,k.

The final matrix D^′(5) equals [D^′]E¯. The correct angles φ for the rotations are found by writing the desired elements of D^′(n+1) as trigonometric functions of elements from D^′(n) and φ and solving the specified equalities for φ. To avoid visible boundaries that would result from a fixed threshold ǫ, we perform a gradual transition between the non-degenerate and the non-degenerate case (ρ= 0) by scaling rotation anglesφby(1−ρ/ǫ). Due to periodicity, there are infinitely many values of φ which solve the given trigonometric equalities. However, scaling φ for interpolation only produces the expected result when selecting the smallest possible value ofφ.

6.3. Experimental Results

We used component-wise convolution with a cubic B-spline kernel [154] to obtain a dif-ferentiable tensor field from the discrete sample values. This method preserves positive definiteness and implies slight smoothing. Figure 6.3 (a) presents aclmap of a coronal sec-tion of the brainstem. It reveals several tracts, which have been annotated by an expert:

The pontine crossing tract (pct), the superior cerebellar peduncle (scp), the decussation

80

6.3. Experimental Results

(a) Annotatedclmap of the brainstem

(b) Edges incl, from eigenvalue derivatives

(c) Approximate edges from scalar cl samples

(d) Edges in FA (e) Edges incs (f) Plots ofcl (—), FA (– –), and cs(· · ·)

Figure 6.3.: Edges in cl, based on eigenvalue derivatives, separate adjacent fiber tracts in DT-MRI data (b). Neither evaluatingclat grid points (c) nor mapping edges in other shape measures (d+e) produces results of comparable quality.

of the superior cerebellar peduncle (dscp), thecorticopontine/corticospinal tract (cpt/cst), and the middle cerebellar peduncle (mcp). Figure 6.3 (b) is produced by computing the total derivative of cl,

c^′_l = λ1(−2λ^′₂−λ^′₃) +λ2(2λ^′₁+λ^′₃) +λ3(λ^′₁−λ^′₂)

(λ1+λ2+λ3)² (6.4)

and using regularized eigenvalue derivatives to evaluate it. Minima in edge strength nicely separate adjacent fiber bundles, which was confirmed by overlaying the edges onto a color coded direction map. Figure 6.3 (c) illustrates that due to the non-linearity of cl, it is not sufficient to evaluate this measure at grid points and to construct an edge map by computing gradients from the resulting scalar samples.

Similar results have previously been obtained by Kindlmann et al. [115], who employ valley surfaces of fractional anisotropy (FA) to reconstruct interfaces between adjacent tracts of different orientation. Figure 6.3 (d) presents an FA edge map of the same region.

Unlike cl, FA can be formulated directly in terms of tensor components, so exact edge maps do not require eigenvalue derivatives. A comparison to Figure 6.3 (b) suggests that cl produces more pronounced fiber path boundaries. In particular, the dscp is hardly separated in the FA edge map.

The observation that a shape measure like FA can be used to find boundaries in orien-tation has been explained by the fact that partial voluming and component-wise interpo-lation lead to more planar shapes in between differently oriented tensors [115]. Sincecl is

6. Eigenvalue Derivatives for Edge Detection in DT-MRI

more sensitive to changes between linearity and planarity than FA is, it is better suited to identify such boundaries.

Further evidence for this reasoning is given in Figure 6.3 (e), which presents edges in cs, an isotropy measure that completely ignores the difference between linearity and planarity, and which consequently is less effective at separating tracts than FA. The visual impression is confirmed by Figure 6.3 (f), which plotscl (solid line), FA (dashed line) and cs (dotted line, uses the axis on the right) against vertical voxel position along a straight line that connects the centers of dscp and scp. It exhibits a sharp minimum incl, a shallow minimum in FA, and no extremum in cs.

6.4. Invariant Gradients in Terms of Eigenvalue