A global approach to cardiac tractography

A GLOBAL APPROACH TO CARDIAC TRACTOGRAPHY Carole Frindel, Jo¨el Schaerer, Pierre Gueth, Patrick Clarysse, Yue-Min Zhu and Marc Robini CREATIS-LRMN, CNRS UMR 5220, INSERM U630, INSA of Lyon and University of Lyon1 69621 Villeurbanne Cedex, France ABSTRACT Cardiac myofibrilles bundles in the human heart can be located by iteratively tracing the local water diffusion direction inferred from diffusion weighted MRI images. This well known process, called streamlining, needs to be initialized by a seed. In this paper, a global, graph modeling approach for cardiac myofibrille tractography is presented. Seed points are no longer needed : it predicts the fibers in one shot for the whole DT-MRI volume without initialization artifacts. The main merits of the presented methodology are its ability to give a unique estimation of the heart architecture (independent of seed initialization), the fact that it has no hyperparameters and that it provides an optimal balance between the density of fibers and the amount of available data. Index Terms— Diffusion tensor MRI, human heart fiber architecture, tractography, graph theory 1. INTRODUCTION In cardiac tractography, myofibrille paths are estimated by tracing iteratively the direction of maximal water diffusion estimated from diffusion weighted MRI images. Despite the potentially multi-directional environment within a voxel, water diffusion is highly anisotropic in most regions of the cardiac muscle and thus the direction of maximal water diffusion aligns with the myofibrille bundle direction. Such estimated fiber paths can be used to investigate and better understand the fiber architecture of the human heart and its muscular organization. Classical tractography methods use simple line propagation techniques, in which a single trajectory is propagated bidirectionally from a manually defined seed point. These local approaches estimate 3D curves by integration of the major eigenvector field (eg. [1]). Euler or higher-order Runge-Kutta schemes are typically used with interpolation of the diffusion tensor field to reconstruct smooth and high density curves. More recent work can be divided into approaches based on Bayesian models [2, 4] and approaches based on a geometrical method, such as level set or fast marching methods [6]. This work was supported by Siemens, the French ANRT and the French Research Unit (GDR) STIC-Sant´e from CNRS-INSERM.

978-1-4244-2003-2/08/$25.00 ©2008 IEEE

883

However, all these methods are local approaches and are subject to the specification of a starting point for the fiber tracts. In this paper, we depend on same hypotheses, but propose a completely new point of view on the tractography problem : our approach is based on graph modeling, where nodes are voxels of the DT-MRI tensor volume and edges represent the connection of two neighbooring voxels by a fiber bundle. Our graph model is thus directly based on a discrete framework and hence provides optimal balance between fiber density and the amount of available information. This choice prevents interpolation artifacts in regions where diffusion is badly defined and where several fiber populations cohabit. Once the graph is defined, our goal is to select the edges that best correspond to the available data and remove all others. The main merit of this global approach to tractography is that the initialization of the tractography process with a seed is not necessary anymore, thus preventing initialization artifacts. It results in an estimation of the principal fiber bundles which capture the global organization of the muscle.

2. MATERIALS AND METHODS 2.1. Data acquisition We studied a set of sixteen ex vivo human hearts from healthy to severely diseased ones (ischemic cardiomyopathy). Each heart was placed in a plastic container and filled with hydrophilic gel to maintain a diastolic shape. The data were acquired with a Siemens Avanto 1.5T MR Scanner. Diffusion images were obtained using an echo planar imaging (EPI) pulse sequence with the following parameters : image reconstructed to 128 × 128 matrix, resolution of 2 × 2 × 2 mm3 , 52 axial contiguous slices. The DT-MRI acquisition protocol is generally defined by the number Nd of diffusion sensitizing directions and the number Ne of excitations used for signal averaging: (Nd ,Ne ). In this work, the acquisition protocol (Nd ,Ne ) = (12, 4) was selected among other combinations in agreement with the acquisition protocol quality comparison conducted in [5]. The acquisition time for a 3-D dataset was 7 37 .

ISBI 2008

Fig. 1. Edge based data fidelity term. Left : configuration where the edge matches the local fiber direction. Right : opposite configuration.

Fig. 2. Tangent-based data fidelity term Tangent based data fidelity term

2.2. Graph model We represent a diffusion tensor volume as a boolean-weighted undirected graph. Nodes are voxels of the diffusion tensor volume and edges represent the connection of two neighbooring voxels by a fiber bundle. Our method selects the edges that best fit the available data by minimizing an energy describing the tractography problem globally. Edges with weight equal to 1 corresponds to an existing fiber connection and vice versa. In order to select the best edges in the graph, we minimize a cost functional made of three terms: a data fidelity term, a topological criterion, and a geometrical criterion. Compatibility of edge and tensor : data fidelity term The role of the data fidelity term is to favor edges that connect two voxels which are likely to be connected by fiber bundles. There are two ways to define this term : edgewise or pointwise. We investigate both options and discuss the pros and cons of the two approaches. Edge based data fidelity term In order to determine whether two neighboring voxels are connected, one must compare their tensors. We introduce a tensor metric, Medge (e), which indicates whether the current edge e matches with the local fiber direction or not. Medge is defined as : Medge (e) = 1 − |1 .e| × |2 .e|

(1)

where 1 and 2 are the principal diffusion directions of the two tensors defining the edge e and . denotes the standard Euclidean inner product. With this expression, the metric Medge (e) will be close to 0 if an edge corresponds to a fiber direction (i.e 1 , e and 2 are collinear) and close to 1 when 1 .e or 2 .e equals 0 (both cases are illustrated in Fig.1). Using this metric, we define the data fidelity term as follows: Fedge (e) = (1 − Medge (e))we + (1 − we )Medge (e)

(2)

where e is the current edge, we is a boolean weight which equals 1 if the edge e is selected and 0 otherwise.

884

One problem with the previous data fidelity term is that it restricts fiber directions to edge orientations. This is highly unsatisfactory since fibers can be of arbitrary orientation, and edge orientations are very limited. To overcome this problem, we propose a tangent-based data fidelity term. We consider the two neighboring edges of a point to jointly define the tangent of the fiber bundle at this point (illustrated in Fig.2). This approach allows many more possible fiber directions. Furthemore, it would be possible to consider higher-order neighborhoods in order to allow even more fiber orientations. However, since the quality of human heart typical DT-MRI data is limited, we do not consider this to be necessary. The new data fidelity term, Ftgt (v), is defined as:  2 )−1 we1 we2 Mtgt (Tv , e1 , e2 ) (3) (CE(v) {e1 ,e2 }⊂E(v) e1 =e2

where E(v) is the set of edges adjacent to node v, Ck2 = k! 2(k−2)! and Tv is the tensor at node v. The orientation metric associated with two adjacent edges e1 and e2 is defined as :  Tv e¯2   (4) Mtgt (Tv , e1 , e2 ) = 1 − ¯ e2 Tv 2 where e¯ = e1 − e2 represents the mean orientation defined by adjacent edges e1 , e2 to node v (the   edges are oriented to point in the opposite direction), and .2 is the spectral norm induced by the Euclidean norm   A = sup Ax2 , 2 x=0 x2

(5)

which reduced to the largest eigenvalue in modulus when A is symmetric. The metric Mtgt favors the case where the mean orientation of a pair of adjacent edges is aligned with the node tensor. Topological criterion Fibers are tubular objects. There are however many configurations of the graph that do not correspond to fibers, such as multiple branches at one node. Histological studies have shown that crossing doesn’t occur in the human left ventricule [7]. It is thus necessary to add a criterion to specify that the graph is composed of tubular

structures. Efficient way to enforce this, at a local level, is to force each node to be connected to at most two neighbors. Let E(v) be the set of edges adjacent to node v and V be the set of nodes of the graph. We consider the topological cost functional to be:  Jtopo (v) (6) Jtopo = v∈V

where,  Jtopo (v) =

(d(v) − 2)2 0

if d(v) > 2 otherwise.

(7)

where d(v) is the degree of the node v. Geometrical criterion In the human heart, it is assumed that the centers of connected neighboring voxels are joined by fiber tracts forming circular arcs. Simple geometrical considerations can be used to discriminate (high energy) arcs whose radius of curvature are not compatible with the human heart architecture (fiber bend) ; for instance : G(e1 , e2 ) = 1 + e1 . e2

(8)

where e1 and e2 are two edges adjacent to the same node (the edges are oriented to point in the opposite direction). Note that this criterion is only needed with edge based data fidelity term. In order to enforce geometrical constraints with the tangent-based term, we simply ignore combinations of edges with high curvature (with an angle smaller or equal to π/2). Cost functional All the previous criteria can be gathered in a global cost functional characterizing our problem. We have two possibilities depending on chosen data fidelity term :  Fedge (e) + αJtopo + β Jedge = e∈E   2 (CE(v) )−1 G(e1 , e2 )we1 we2 v∈V

Jtgt =

{e1 ,e2 }⊂E(v) e1 =e2



Fig. 3. Fiber tracts estimation from linear and circular artificial data. Left: Principal direction fields associated with artificial datasets, Middle: Fiber tracts estimation with the edge-based data fidelity term, Right: Fiber tracts estimation with the tangent-based data fidelity term of all its neighbors, in the purpose of minimizing the global cost functional. Since the energy of an edge can be estimated from its neighborhood, we can make a decision for the flip by evaluating only this local energy. This algorithm is applied sequentially and repeated until no more flip occurs. Fiber trajectory extraction We extract the different connected components of the graph. Each of these components correspond to the topology and approximate geometry of a fiber bundle: we know voxels through which each fiber runs, but we don’t know exactly where in each voxel. Since cardiac fibers should be reasonably smooth, we chose to constrain the geometry with smoothing B-Splines. 3. RESULTS 3.1. Artificial data

(9)

Ftgt (v) + αJtopo

v∈V

where V is the nodes set, E(v) is the set of edges adjacent to node v and E is the edges set of the graph. Note that α is not an empirically fixed hyper-parameter : it is gradually increased during the minimization process so that the connexity of each voxel is equal to at most two in the final configuration. Minimization of the cost functional We use iterated conditional mode (ICM) [3]. This minimization algorithm iteratively estimates the state (0 or 1) of each edge given the states

885

Fig. 3 shows results obtained on artificial data. We note that, in certain cases, the edge-based data fidelity term gives innacurate results. When fibers are straight (Fig.3 up), the fiber population in the right region has a direction that does not follow edge directions anymore : edge-based criterion is unable to follow fibers properly. The tangent data fidelity term shows better performance. In the circular data example (Fig.3 down), the two terms give very similar results : the directional limitation observed with the edge-based term is only critical when fiber curvature is low. 3.2. Real human heart data Fig. 4 shows results on human data, as described in section 2.1. Results show good agreement between predicted

fiber and histological knowledge about the structure of the human myocardium. The computations took approximately 2 hours for the tractography of a whole heart (128 × 128 × 52 nodes), with an unoptimized implementation written in the Python language, on a standard PC. We estimate that with a fast implementation written in C or Fortran, only a few minutes would be necessary to capture the fiber architecture of a whole heart.

The graph model we propose, together with the tangentbased data fidelity term is thus a new and promising approach to the cardiac tractography problem. The approach does not require the specification of seed points, nor does it rely on the choice of hyper-parameters. The method is able to accurately capture the hearts principal fiber tracts in one shot, and provides optimal balance between the density of fibers and the amount of available data. This is particularly interesting for a better understanding of the architecture and organization of the cardiac muscle continuum. Future work will include further validation against other imaging techniques and/or histological cuts, and the improvement of the data fidelity term. 5. ACKNOWLEDGEMENTS We would like to thank J. Dardenne for his help for the development of the visualization tool, and L. Guigues for constructive discussions on the graph model. We would also like to thank P. Croisille, E. Stephant and S. Rapacchi who provided us with the DT-MRI data. 6. REFERENCES [1] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi. In vivo fiber tractography using DT-MRI data. Magn. Reson. Med., 44:625–632, 2000. [2] T. Behrens, M. Woolrich, M. Jenkinson, H. JohansenBerg, R. Nunes, S. Clare, P. Matthews, J. Brady, and S. Smith. Characterization and propagation of uncertainty in diffusion-weighted MR imaging. Magn. Reson. Med., 50(5):1077–1088, 2003. [3] J. Besag. On the statistical analysis of dirty pictures. Journal of the Royal Statistical Society. Series B (Methodological), 48(3):259–302, 1986. [4] O. Friman, G. Farneb¨ack, and C.F. Westin. A bayesian approach for stochastic white matter tractography. IEEE Trans. Med. Imaging, 25(8):965–978, 2006.

Fig. 4. Global fiber prediction using tangent-based data fidelity term, (A) seen from the top, (B) detailled through the left ventricular wall and (C) cut along long axis.

4. CONCLUSIONS AND PERSPECTIVES We have presented and studied two different global functionals for cardiac tractography. It clearly apprears from the results that the tangent-based data fidelity term is superior: it gives better results and requires no hyper-parameters.

886

[5] C. Frindel, M. Robini, S. Rapacchi, E. Stephant, Y. Zhu, and P. Croisille. Towards in vivo diffusion tensor MRI on human heart using edge-preserving regularization. Proc. 29th Int. Conf. IEEE EMBS, pages 6007–6010, 2007. [6] G. J. Parker, C. A. Wheeler-Kingshott, and G. J. Barker. Estimating distributed anatomical brain connectivity using fast marching methods and diffusion tensor imaging. IEEE Trans. Med. Imaging, 21(5):505–512, 2002. [7] D. F. Scollan, Alex Holmes, Raimond Winslow, and John Forder. Histological validation of myocardial microstructure obtained from diffusion tensor magnetic resonance imaging. Am. J. Physiol. Heart Circ. Physiol., 275:H2308–H2318, 1998.