2492

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 37,

NO. 12,

DECEMBER 2015

Maurer-Cartan Forms for Fields on Surfaces: Application to Heart Fiber Geometry Emmanuel Piuze, Jon Sporring, and Kaleem Siddiqi, Senior Member, IEEE Abstract—We study the space of first order models of smooth frame fields using the method of moving frames. By exploiting the Maurer-Cartan matrix of connection forms we develop geometrical embeddings for frame fields which lie on spherical, ellipsoidal and generalized helicoid surfaces. We design methods for optimizing connection forms in local neighborhoods and apply these to a statistical analysis of heart fiber geometry, using diffusion magnetic resonance imaging. This application of moving frames corroborates and extends recent characterizations of muscle fiber orientation in the heart wall, but also provides for a rich geometrical interpretation. In particular, we can now obtain direct local measurements of the variation of the helix and transverse angles, of fiber fanning and twisting, and of the curvatures of the heart wall in which these fibers lie. Index Terms—Differential geometry, moving frames, Maurer-Cartan form, diffusion MRI, heart wall myofibers

Ç 1

INTRODUCTION

T

HE

method of moving frames, as refined by Cartan[1], describes the geometrical structure and the invariants of Lie groups or equivalently, of smooth manifolds. The method generalizes the familiar concepts of curvature and torsion for space curves and of the shape operator for surfaces embedded in euclidean space R3 to higher dimensions. Its intuitive conceptualization produces a variety of useful geometrical descriptors for studying three-dimensional entities. In fact, any smooth manifold in R3 can be represented via a set of orthonormal frames which can then be analyzed using the method of moving frames. The study of smooth frame fields in three-dimensions has applications in many areas of pattern analysis, including fluid flow analysis in computational mechanics, solving equivalence and symmetry problems in geometry in physics, and the development of invariant geometrical features in computer vision [2]. Examples of the latter include Faugeras’ study of the geometry and evolution of curves in the euclidean, affine and projective planes [3], Flash and Handzel’s analysis of the differential geometry of motion paths [4] and Boutin and Bazin’s consideration of the classical structure from motion problem from a differential geometric viewpoint [5]. Related work in modeling streamlines include Ben-Shahar and Zucker’s analysis of texture flows arising in natural scenes [6], Savadjiev et al.’s use of tensor connections in diffusion MRI [7], Calabi’s development of object recognition signatures [8], and Olver’s joint invariant signatures [9].

The first efforts towards the methods presented herein were developed in [10], [11], [12], [13], where experiments showed that the method of moving frames and a minimal surface model, the generalized helicoid model (GHM), could be leveraged to characterize cardiac fiber geometry. We extend these findings in several ways. First, we have a more rigorous exposition of the method of moving frames. This leads to a detailed discussion on the generation of various geometrical embeddings, of which the heart is a specialization. More precisely, in the method of moving frames, the Maurer-Cartan form is an operator that measures the differential structure of a manifold. The operator is typically applied in a forward manner to study the geometrical characteristics of the manifold under consideration. In the present article we also apply the theory of moving frames in the reverse direction: the Maurer-Cartan connection forms are manipulated to generate manifolds or embeddings based on certain assumptions of their differential structure. In fact, this idea can be used to characterize a smooth frame field in three dimensions as a parameterization on the space of Maurer-Cartan connection forms. We also extend prior work by developing and comparing a number of fitting methods for estimating various geometrical embeddings. Lastly, we include extensive experimental work demonstrating the potential of this analysis as well as provide a deeper interpretation of these novel geometrical probes and their relation to past literature on heart wall fibre orientation. Our key contributions are listed below.





E. Piuze and K. Siddiqi are with the School of Computer Science and the Centre for Intelligent Machines, McGill University, Montreal, QC H3A 2A7, Canada. E-mail: {epiuze, siddiqi}@cim.mcgill.ca.
J. Sporring is with the Department of Computer Science, University of Copenhagen, Universitetsparkein 1, Copenhagen, DK-2100, Denmark. E-mail: [email protected]. Manuscript received 25 Feb. 2014; revised 15 Jan. 2015; accepted 22 Jan. 2015. Date of publication 2 Mar. 2015; date of current version 6 Nov. 2015. Recommended for acceptance by R. Vidal. For information on obtaining reprints of this article, please send e-mail to: [email protected], and reference the Digital Object Identifier below. Digital Object Identifier no. 10.1109/TPAMI.2015.2408352





We introduce various connection form embeddings for studying frame fields, and we show that the GHM used in [12] can be fully described locally using a combination of connection forms. We introduce a fitting energy for estimating connection forms from arbitrary smooth frame fields. In application to modeling cardiac fiber geometry, we attach a frame field to diffusion MRI data and study its differential geometry, thereby corroborating and also providing novel biological measurements of heart fiber orientation.

0162-8828 ß 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

PIUZE ET AL.: MAURER-CARTAN FORMS FOR FIELDS ON SURFACES: APPLICATION TO HEART FIBER GEOMETRY

2493

The article is organized as follows. The method of moving frames is introduced in Section 2. Techniques are then developed to compute the Maurer-Cartan parameterization for describing the geometrical behaviour of a frame field to first order approximation in Section 3. Section 4 describes various manifold embeddings as constraints on connection forms. Finally, in Section 5, we employ these methods to characterize the geometry of cardiac muscle fibers.

2

THE STRUCTURE OF FRAME FIELDS

P Let a point x ¼ i xi e i 2 R3 be expressed in terms of the natural orthonormal coordinate system e 1 ; e 2 ; e3 . A differential orthonormal frame field embedded in R3 will be denoted as F ¼ ½ff 1 ; f 2 ; f 3 T : R3 ! R3 , and defined by f i  f j ¼ dij , where dij is the Kronecker delta and where  is the inner product, and by f 1  f 2 ¼ f 3 , where  is the three-dimensional cross product. Since the frame E ¼ ½ee1 ; e2 ; e3 T forms an orthonormal basis for R3 , F can be P expressed as the rotation f i ¼ j aij ej , where the elements of the attitude matrix A ¼ faij g 2 R33 are differentiable, and A 1 ¼ AT . Treating f i and ej as symbols, we can write this rotation in matrix form: T

T

½f 1 f 2 f 3  ¼ A½e 1 e 2 e3  ()F F ¼ AE :

(1)

Since each e i is constant, the differential geometry of the frame field is completely characterized by the attitude matrix A. Taking the exterior derivative on both sides of (1), we have [14]

% E ¼ ðdA AÞA A1 F ¼ C F ; dF F ¼ ðdA AÞE E þ A dE 0

(2)

where d denotes the exterior derivative operator, and C ¼ ðdA AÞA A1 ¼ fcij g 2 R33 is called the Maurer-Cartan matrix of 1-forms cij . Writing f i as symbols, (2) is to be understood as X dff i ¼ cij f j : (3) j

The Maurer-Cartan matrix is skew symmetric [14], hence we have 2 3 0 c12 c13 (4) 0 c23 5; C ¼ 4 c12 c13 c23 0 such that there are at most three independent, non-zero elements: c12 , c13 , and c23 . xÞ : R3 ! R of A is The differential of the elements aij ðx k expressed in terms of dx , the natural basis for 1-forms in R3 : daij ¼

X @aij k

@xk

dxk :

(5)

1-forms are differential operators that may be applied to tangent vectors through a process denoted contraction, written as dwhvvi 2 R for a general 1-form dw on R3 and tangent P vector v 2 R3 . The contraction of dw ¼ k wk dxk on v is

Fig. 1. Geometry characterized by connection forms.

written as aP bilinear operation onPthe canonical projection of P v : dwhvvi ¼ i wi deei h j vj e j i ¼ i wi vi . Hence, when contracted, the elements of C in (2) are given by cij hvvi ¼

X

ajk

k

X @aik l

@xl

vl ;

(6)

and they express the amount of turning of f i ðx xÞ towards xÞ when x moves in the direction v . This analysis natuf j ðx rally applies to frame fields embedded in R2 by setting f 3 to be constant. In that case, C is a 2  2 matrix with a single 1-form c12 which describes the in-plane turning of the frame; hence, if f 1 is the tangent of a planar curve and f 2 its normal, then c12 hff 1 i is the curves curvature and c12 hff 2 i its fanning [6]. The space of linear models for smooth frame fields is fully characterized by the elements cij of C , and can be generated from the first order terms of a Taylor series of f i centered at a point x 0 . The first order approximation P for the motion of a frame vector f i in a direction v ¼ k vk f k can be expressed as    f i hvvix ¼ f i x þ dff i hvvix þ Oðkvvk2 Þ 0

0

0

X  X   f j x vk cijk x þOðkvvk2 Þ; ¼ f i x þ 0

j6¼i

0

k

0

(7) (8)

where we use the short-hand cijk  cij hff k i for the natural connections of the local frame. Since only three unique nonzero combinations of i; j are possible in R3 due to the skew symmetry of C , there are only nine unique non-zero combinations of cijk possible for k ¼ 1; 2; 3. Fig. 1 illustrates the behavior of a frame field where each connection form predominates. Considering the model at a point x 0 for all frame vectors and 1-forms, an approximation of the frame from F ðx x0 Þ to F ðx x0 þ v Þ can be obtained by the normalized model

2494

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 37,

NO. 12,

DECEMBER 2015

Fig. 2. (From left to right) parameters after 200 Nelder-Mead iterations (in radians/voxel for cijk and radians for i ), convergence plot for c123 and 1 , and computational timings, for different estimation techniques. i

f þ dff i hvvi ; x0 þ v Þ ¼ i f~i ðx kff i þ dff i hvvik

(9)

which in matrix notation yields

cijk ¼ dff i  f j hff k i ¼ f Tj Jðff i Þff k :

f þ FC F v x0 þ v Þ ¼ i ; f~i ðx kff i þ F C i F T v k i

@x where J ¼ ½@x  2 R33 is the Jacobian matrix of partial derivj atives and dx x ¼ ðdx x1 ; dx x2 ; dx x3 ÞT . Setting v ¼ f k , we finally obtain

T

(10)

where C i 2 R33 denotes the matrix of connection forms with C ikl ¼ "ik ckil and "ik ¼ sgnðk  iÞ. In general f~i ðvvÞ; f~2 ðvvÞ; f~3 ðvvÞ will not form an orthogonal basis. To see this, combine (3) and (7) to obtain    1 þ c2ij þ c2ik : i ¼ j (11) f i hvvix  f j hvvix ¼ 0 0 cik cjk : i 6¼ j; for k 6¼ i; k 6¼ j.

The Jacobian matrix J can be approximated to first order  @ff using, e.g., finite differences: @xij ðx xÞ  12 f ij ðx x þ ek Þ  k x  ek ÞÞ. f ij ðx

3.2 Computation via Energy Minimization The connection forms cijk at a point x 0 can also be found as the minimizer of the energy contained within a neighbourhood V: x0 Þ ¼ arg min c ijk ðx cijk

3

ESTIMATION OF CONNECTION FORMS

The full space of linear models for smooth frame fields was introduced in (8). It was found that this space has at most nine independent parameters cijk , which fully characterize the local geometry of a frame field. This section will be devoted to methods for computing these parameters for a smooth frame field.

3.1 Direct Computation of Connection Forms The connection 1-forms cij can be directly obtained from the vector fields f i and their differentials dff i : ! 3 X cij f j  f k using (2) (12) dff i  f k ¼ j

¼

3 X

since f i  f k ¼ dik

cij djk

(13)

j

¼ cik :

(14)

The differential dff i can be computed by applying the exterior derivative of a function g : Rn ! R, dg ¼

n X @g dx xi ; @x i i

(15)

which yields dff i  f j hvvi ¼

3 3 X X @ff ik l

k

@xk

! dx x

k

!

f jl hvvi

(16)

¼ f Tj Jðff i Þdx xhvvi

(17)

¼ f Tj Jðff i Þvv;

(18)

(19)

3 1 XX "i ðx x0 þ v Þ; jVj v 2V i

(20)

where "i is a function of the error for each frame axis,   "i ðx x0 þ v Þ ¼ arccos f i ðx x0 þ v Þ  f~i ðx x0 þ v Þ ;

(21)

where f i ðx x0 þ v Þ is the frame data term and f~i ðx x0 þ v Þ is its approximation using (10). Here, V can take any shape. We will denote a cubic isotropic neighborhood in R3 of radius i as V2iþ1 , and 6 nearest neighbors will be denoted as Vþ 3. The optimization energy in (21) can be solved for using standard algorithms such as Nelder-Mead [15] and BOBYQA [16].

3.3 Method Comparison Fig. 2 compares mean connection form estimation results and timings for direct and optimized computations using a sample of a smooth frame field of biological (cardiac) nature containing 50,000 nodes. Here, all connection forms are estimated in a neighborhood Vþ 3 and then evaluated in V3 . Finite difference computations of Section 3.1 are referred to as Finite, and Nelder-Mead computations as NM. Seeded computations make use of the results obtained using finite differences. In general, unseeded computations systematically yield a larger error than seeded ones, but both eventually converge to the same error, at around 150 iterations, which requires about 3s of computation time. In comparison, BOBYQA optimization takes about 1min to fully converge, and yields a slightly lower error than Nelder-Mead and finite computations. One major advantage of the BOBYQA optimizer is the possibility of enforcing bounds on the variables to fit, since it is a constrained optimization scheme. In general, bounds can be obtained on the values of the cijk terms, based on both empirical and theoretical observations. In the discrete case for a smooth orthonormal frame field, and using a central differentiation scheme, the frame

PIUZE ET AL.: MAURER-CARTAN FORMS FOR FIELDS ON SURFACES: APPLICATION TO HEART FIBER GEOMETRY

axis differential dff i ¼ dðff i1 e 1 þ f i2 e2 þ f i3 e 3 Þ in J f i is bounded since      @ff i  f ðx x  e k Þ  f ij ðx x þ e k Þ     ij  1: @x    2 k x

(22)

Using kk1 to denote the euclidean 1-norm, the following bound is thus obtained: jcij hvvij ¼ jff Tj J f i ðx1 ; x2 ; x3 Þ v j

(23a)

f Tj  ½ kvvk1

kvvk1

kvvk1 T

(23b)

kvvk1 a;

for v 2 V2aþ1 :

(23c)

In certain applications where computation time is not an issue and where there is a prior on connection form bounds, the seeded BOBYQA optimization scheme can be preferable. For the remainder of this article, a Nelder-Mead optimization scheme running 300 iterations will be assumed. The upper bound of (23) will be enforced by discarding volume elements that fall outside of permitted values. The Pointcare-Hopf theorem states the existence of at least one singular point for frame fields embedded in surfaces with a non-zero Euler characteristic. An example is the singularity that naturally arises in the GHM (see (28)). Whereas in theory characterizations of open sections on manifolds can be made free from singularities [17], in dealing with acquired or fitted data these could still arise due to the discretization of the underlying frame field. In applications such as the one we focus on in Section 5, singularities would generate a sharp turning in the frame vectors. This behavior cannot be captured by the first order model of (8) which will therefore only yield a coarse approximation. Using the direct computations of (19), (23) shows that connection forms will be bounded numerically near singularities. However, optimized computations may yield large values, in which case a hard threshold can be placed on (20) or a regularization term can be added, as in [11]. These strategies should be seen as heuristics, and will in general not yield a good fit to the data close to singularities.

4

MAURER-CARTAN EMBEDDINGS

The geometrical intuition conveyed by the cijk parameters in (8) and illustrated in Fig. 1 can be used to develop local embeddings or connection form models. By imposing constraints on certain parameters, the shape and the complexity of the resulting frame field can be controlled explicitly. x0 Þ with known parameters cijk , the Given a seed frame f i ðx solution to the manifold f i ðx xÞ for all x can then be found by integration. In some cases—as with the GHM to be discussed below [12]—there exist closed form solutions. In general, models may be solved for numerically by using (10) in a standard integration scheme. This section demonstrates embeddings based on generalized helicoids, and those which lie on spherical and ellipsoidal shells. The general case in which all connection forms are used will be referred to as the full connection form model, and the one where all connection forms are set to zero as the constant connection form model.

2495

4.1 Generalized Helicoids as Connection Forms A model for representing the geometry of 3D smooth streamline flows was proposed by Savadjiev et al. in [13], which built on earlier work by Ben-Shahar and Zucker [6] on modeling texture flows. The Generalized Helicoid Model was used to measure axonal geometry in white matter fiber tracts [13] and to characterize muscle directions in the heart [12] by modelling the local variation of a smooth frame field. In this section it is shown that the GHM model is a subset of (10), where all but three parameters are zero. In [12], f 1 is given by the GHM vector field f 1 ðx x0 þ v Þ ¼ cos u f 1 ðx x0 Þ þ sin u f 2 ðx x0 Þ   K1 ðx x0 Þ v1 þ K2 ðx x 0 Þ v2 uðx x0 ; v Þ ¼ arctan 1 þ K2 ðvv0 Þ v1  K1 ðx x 0 Þ v1

(24a)

(24b)

x 0 Þ v3 ; þ K3 ðx where uðx x0 ; v Þ parameterizes aPminimal surface over paramx0 Þ. The nine parameeter fields Ki and where v ¼ i vi f i ðx ters cijk of (10) can then be evaluated for f~1 . The trivial basis ei ¼ f i ðx x0 Þ is chosen in order to construct the attitude matrix A, which yields 2 3 cos u sin u 0 AðuÞ ¼ 4  sin u cos u 0 5; (25) 0 0 1 where u ¼ uðx x0 ; v Þ. Since f 1 stays parallel to the tangent f 1  f 2 plane at x 0 and f 3 remains unchanged, thus f~2 ¼ f 3  f 1 will also stay parallel to that tangent plane at x0 . Hence, we must have c13 ¼ c23 ¼ 0, which is confirmed by straightforward computations. The remaining 1-form is found to be c12 ¼ K1 dff 1 þ K2 dff 2 þ K3 dff 3 ;

(26)

since the derivatives are to be evaluated at the origin, and where Ki and dff i are evaluated at x 0 . Thus, c12k ¼ c12 hff k i ¼ Kk ; k 2 f1; 2; 3g;

(27a)

and we conclude that the GHM is a subset of (10), where all but c121 , c122 , and c123 are zero. As a result, in frame fields where the six connection forms c13k ; c23k jk 2 f1; 2; 3g are non-zero, the GHM will not be able to comprehensively characterize the underlying geometry of the frame field. Other issues can arise with the GHM since the rotation angle uðx x0 ; v Þ of (24b) yields a frame field singularity when [18] K1 v ¼ ð1 þ ð Þ2 Þ1 ðK2 ; K1 ; 0Þ: K2

(28)

4.2 A Model on Spheres Connection forms can be constrained to generate thin spherical shell embeddings. Given a smooth surface where f 1 and f 2 span the local tangent plane, we can generate a family of subspaces that are non-intersecting along f 3 , by setting c12 hff 1 i ¼ K1 ;

c12 hff 2 i ¼ K2 ;

c13 hff 3 i ¼ 0;

c12 hff 3 i ¼ K3 ;

c23 hff 3 i ¼ 0;

(29a) (29b)

2496

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 37,

NO. 12,

DECEMBER 2015

Fig. 3. Solution to (29) and (30) for three neighboring spheres r ¼ 1; 1:15; 1:3 (yellow, blue, red). K1 models the turning towards the tangent vector f 1 . Each column represents a distinct tuple of (K2 ; K3 ), which expresses the turning of the tangent vector towards f 2 when moving along the f 2 (spreading) and f 3 (turning) directions.

where Ki : R3 ! R are any smooth functions. The remaining parameters control the motion of frames within parallel spherical surfaces. For a point x 0 on a spherical shell with center at x 0  r1 f 3 , we have 1 c13 hff 1 i ¼ ; r

c13 hff 2 i ¼ 0;

(30a)

c23 hff 1 i ¼ 0;

1 c23 hff 2 i ¼ : r

(30b)

This completely specifies the local linear model for spherical shells using (8). By numerical integration we can generate flow lines of f 1 in a local neighborhood, as illustrated in Fig. 3. Curves with f 1 as their tangent lie on a single shell and each forms a circle, as shown in Fig. 4. To prove that flowlines of f i form circles on a sphere, consider a sphere in polar representation, 2 3 2 3 x r cos ðfÞ sin ðcÞ 4 y 5 ¼ 4 r sin ðfÞ sin ðcÞ 5; z r cos ðcÞ

inserting this into R, which yields R ¼ r= 1=jdff 1 hff 1 ij.

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ K12 r2 ¼

4.3 A Model on Ellipsoids The spherical embedding of Section 4.2 is a specialization of ellipsoidal geometry. In the general case, manifolds are embedded within thin ellipsoidal shells, which we refer to as homeoids. Given any smooth surface with anisotropic principal curvatures, and where f 1 and f 2 span the tangent plane which is organized in thin shells, (29) and (30) generalize to 1 (32a) c13 hff 1 i ¼ ; c13 hff 3 i 6¼ 0; r1 1 c23 hff 3 i 6¼ 0; c23 hff 2 i ¼ : (32b) r2 Here, c13 hff 3 i and c23 hff 3 i are only zero when f 1 and f 2 are aligned with the principal directions of the underlying surface. This completely specifies the local linear model for elliptical shells using (10). As for the spherical shell model,

(31)

centered at 0 and with radius r, and consider the intersection with a plane orthogonal to the z-axis of height r z0 r. Since z0 ¼ r cos c, using (31) the curve of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi intersection is found to be p ¼ ½ r2  z20 cos f; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2  z20 sin f; z0 T , which is a circle of radius qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ r2  z20 with the z-axis as the axis of rotation. Evaluating the Maurer-Cartan form of a frame field, where @ p

T

f1 ¼ jj@ff pjj, f3 ¼ ðx;y;zÞT , and f2 ¼ f3  f1 , we find that jjðx;y;zÞ jj ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q c12 hf1 i ¼ z0 =ðr r2  z20 Þ, c13 hf1 i ¼ 1=r, and c23 hf1 i ¼ 0, which exactly corresponds to the constraints given by (29) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and (30) when K1 ¼ z0 =ðr r2  z20 Þ. To align the two coordinate systems, we first set the tangent of the circle to be f 1 , and set the change of f 1 in the direction f 1 to be the normal to the circle. The normal is then derived directly from (10), (29), and (30) as the unit vector parallel to dff 1 hff 1 i ¼ K1 f 2  ð1=rÞff 3 . Finally, the radius of the circle is found by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi solving for z0 in terms of K1 , z0 ¼ K1 r2 = 1 þ K12 r2 and

Fig. 4. Solutions to (29) and (30) on the unit sphere, with initial frames shown in red, green, and blue. Column 1 shows a single flow line (yellow) with the yellow arrow pointing to the origin. Columns 2-4 shows neighbouring flow lines from an initial starting point (red) for various ðK1 ; K2 Þ tuples.

PIUZE ET AL.: MAURER-CARTAN FORMS FOR FIELDS ON SURFACES: APPLICATION TO HEART FIBER GEOMETRY

2497

Fig. 6. The helix angle aH is defined as the angle between the short axis of the heart and the projected myofiber direction in a local normal plane.

Fig. 5. Examples of ellipsoidal solutions to (29) and (32) integrated from a flat patch (first row), a spherical patch (second row), and a cylindrical patch (third row) for various values of c131 ; c132 ; c231 ; c232 .

flow lines can be generated in a local neighbourhood by numerical integration. However, these yield a more complex geometry, as illustrated in Fig. 5, for frame integrations on various surfaces: a flat, a spherical, and an ellipsoidal surface.

5

APPLICATION TO MYOFIBER GEOMETRY

The walls of the ventricles in the mammalian heart are composed of elongated muscle cells called cardiomyocytes, which are densely packed within a collagen matrix. This matrix forms the bulk of the heart, which has a truncated ellipsoidal shape. Cardiomyocytes stack approximately end on end, forming smoothly varying structures known as myofibers. The arrangement of myofibers is critical for normal heart function because it is the alternate contraction and relaxation along their length that determines pumping efficiency. In the following we apply connection form based modeling to characterize myofiber geometry. This could have many practical uses including differentiation between normal and pathological arrangements, integration of myofiber geometry into patient-specific cardiac models and monitoring changes in heart wall structure in studies of development and aging. Early work based on histology has shown that cardiac myofibers are locally aligned to helical curves which wrap around the ventricles [19]. Various models have been proposed to explain this organization, based on both physical and mathematical considerations [20], [21], [22]. More recently it has become possible to study this organization within the intact myocardium using diffusion magnetic resonance imaging (dMRI) [23], [24], [25]. In analyses of myofiber geometry, the helix angle—taken to be the angle between the fiber direction projected onto a plane orthogonal to the penetration direction and the short-axis plane (see Fig. 6)—is ubiquitous [25], [26], [27]. Several accounts from both small-scale histology and voxel-scale studies using dMRI report that along a transmural penetration line from

the heart’s outer to inner wall, the helix angle varies smoothly and regularly, undergoing a total change in orientation of about 120 degree in mammals [19], [25], [26], [28], [29]. The variation in helix angle provides a coarse description of fiber geometry, which in turn depends on the choice of penetration direction, as illustrated in Fig. 7. The methods developed in the previous section have the potential to provide more complete local characterizations of heart wall myofiber geometry. The scale of current dMRI measurements is at least one order of magnitude larger than the length of individual cardiomyocytes [30] and thus the diffusion signal reflects properties of a fibrous composite. Modeling myofiber geometry at this scale Savadjiev et al. [12] have used the GHM of Section 4.1 to parametrize myofiber orientation in the heart wall in three mammalian species, the rat, the canine, and the human. A limitation of using the model in [12] is that its flow lines are constrained to lie on planar manifolds, in spite of the curvature of the heart wall (Fig. 6). As described in Section 4.1, this limitation allows only three available geometrical parameters, resulting in a reduced ability for characterizing increasingly larger and complex myofiber neighborhoods. Motivated by this model we propose instead to fit a frame field to cardiac dMRI data and to then examine its connection forms. We shall apply this analysis to a diffusion MRI database of healthy rat hearts [11], [12], containing 8 rat subjects, which will be labelled as subjects A to H, with ð0:25 mmÞ3 voxel resolution and a dimension of 64  64  128 voxels.

5.1 Selecting a Cardiac Frame Field Cardiac diffusion MRI volumes arePsampled on 3D rectangular lattices with coordinates x ¼ i xi ei ; xi 2 Z3 . In order to apply connection form analysis, a local cardiac frame

Fig. 7. Principal fiber direction and cardiac frame field f 1 (red), f 2 (green), and f 3 (blue) for a rat heart, obtained from dMRI. The base is located upwards and the apex downwards. Color indicates the helix angle of the fiber direction, from 90 degree (green) to þ90 degree (red).

2498

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

Fig. 8. Endocardium (left) and epicardium (right) euclidean distance, ranging from zero (black) to 20 voxels (red), and distance gradient vectors (blue).

needs to be defined. The first frame vector of this frame is chosen to be parallel to the orientation of the principal direcu1 . Because the direction of u 1 tion of diffusion u 1 , i.e., f 1 nnu is locally ambiguous, our first task is to enforce consistency in the direction of f 1 among voxel neighbors. This is done using an adaptive cylindrical coordinate system. The centroid cz of the chamber within each short-axis slice mask zz is first determined using P cz ¼

x 2zz

jzz j

x

:

c zþ1  c z1 lz ¼ kcczþ1  c z1 k   f 1 ðx xÞ ¼ sign ðu 1 ðx xÞ  ðx x  cz ÞÞ  l z Þ u 1 ðx xÞ u1 Þu u1 ; ¼ nx ðu

(34) (35) (36)

where l z is the local approximation to the heart’s long-axis. For most hearts under consideration, lz approximately coincides with the volume’s z axis. We now need one additional orthogonal frame axis to f 1 in order to define a frame field. In fact, since we require f i 2 R3 to have unit length, only one additional degree of freedom is required to fully describe a frame field: a rotation about f 1 . In principle, any model for this rotation could be adopted, based on geometrical intuition, biological relevance, or computational considerations. Following [12] and conventional cardiac spatial analysis, we use a model for the frame field which is based on an estimate of the local normal to the heart wall. The normal is approximately the direction in which the endocardium moves when the heart beats, which is also naturally orthogonal to local myofiber orientations, and thus is biologically meaningful. This also leads to a consistent choice of the remaining frame field directions throughout the heart wall, and a smooth rotation of the frame field from one voxel to its neighbor. Specifically, we first estimate a transmural direction f^3 using the gradient vector of a distance transform produced as follows: a) the binary image (mask) of the heart is closed using mathematical morphological operations, b) the euclidean distance to the heart wall is evaluated at every point to both the epicardium and endocardium (Fig. 8), and c) the endocardium gradient is negated and the average of

NO. 12,

DECEMBER 2015

Fig. 9. Flow lines in the myocardium, found by tracking the first eigenvector of diffusion originating in the red boxes. A frame indicates the axial plane of the heart (red axis), the transmural direction (green axis), and the direction from the base to the apex (blue axis). The fiber geometry expressed in this coordinate system corresponds to the connection forms c123 (a), c232 (b) and c131 (c), using Fig. 1 for comparison.

the two gradients is then computed. The normals f^3 are then aligned to point from outer to inner wall. With f 1 and f^3 , a local frame is specified at x , f1 ¼

nx ðu u1 Þu u1 ; ku u1 k

(37a)

f3 ¼

ðf^3  ðf^3  f 1 Þff 1 Þ ; kf^3  ðf^3  f 1 Þff 1 k

(37b)

(33)

f 1 is then locally made to turn clockwise with respect to that centroid for each slice zz , and for all x 2 zz :

VOL. 37,

f 2 ¼ f 3  f 1;

(37c)

where f 3 is the part of f^3 orthogonal to f 1 . The local plane spanned by f 1 and f 2 will be referred to as the tangent plane. Fig. 7 shows a sample of this frame field within two crosssectional cardiac slices. Near structural boundaries, and due to various anomalies, discontinuities may arise in the frame field. This is typically not observed due to biology [27], and due to the inherent smoothing of the diffusion signal [24]. Various heuristics could be used to identify and mask out these locations, for example using (23).

5.2 Cardiac Maurer-Cartan Connections We shall now employ the Maurer-Cartan connection forms to analyse the cardiac frame field. To get a preview of the biological relevance of our approach Fig. 9 illustrates the correspondence between three large-scale geometrical features exhibited by cardiac fibers and connection forms: the variation of the helix angle is related to c123 , the short-axis curvature of the heart wall to c131 , and the long-axis curvature to c232 . 5.2.1 Cardiac Moving Frame Intuition Considering the frame at x 0 , 1-form contractions cijk can be interpreted as the amount of turning of f i towards f j when considering neighboring frames in the direction f k . The rotations of f 1 towards f 2 , c12 hff k i, are intimately linked to the curvature parameters of the GHM (see Section 4.1). These rotations intuitively describe the manner in which myofibers turn in the tangent plane of the heart: c121 , shown in Fig. 1a, describes their tangential curvature, c122 , shown in Fig. 1b, their fanning in the tangent plane, and c123 , shown in Fig. 1c, their transmural turning or equivalently the rate of change of the helix angle. The rotations of

PIUZE ET AL.: MAURER-CARTAN FORMS FOR FIELDS ON SURFACES: APPLICATION TO HEART FIBER GEOMETRY

2499

Fig. 11. Selected dataset volumes of rat connection forms cijk (radians/ voxel) and fitting errors i (radians) in V3 for the optimized and direct computations.

Fig. 10. Volume histograms of cijk connections (radians/voxel) and fitting errors (radians) for V3 by fusing all rats in the dataset using optimized parameter computations (shown in red) together with the min (red) and max (blue) dataset envelope for each value.

f 1 towards f 3 , c13 hff k i, express the turning of the fibers towards the inner wall: c131 is related to a sectional curvature of the heart wall in the tangential direction, c132 describes an upwards twisting of the tangent plane and the rate of change of the transverse angle, and c133 described the degree of transmural fanning of the local fiber population away from or towards the tangent plane. c231 measures twisting of the tangent plane, c232 is related to a second sectional curvature of the heart wall, and c233 measures myofiber fanning (in or out) in the transmural direction.

5.2.2 Full Volume Histograms Full volume histograms for optimized cijk computations in V3 are shown in Fig. 10 for each subject in the dataset. The histograms show that connection forms, in particular the variation of the helix angle c123 , are globally consistent among different healthy subjects of the same population. These connection histograms could be used for studying various cardiomyopathies in place of focusing solely on the helix angle. The latter cannot be used for detecting or surgically restoring anomalies such as myocardial infarction and hypertrophy [27]. For example, in hypertrophic hearts the myocardium becomes more spherical globally; using our framework, this should result in sharper histograms for the sectional descriptors c131 and c232 . Other cardiomyopathies that are more localized could be investigated from the local geometrical signature of the cardiac frame field.

A comparison of selected volume histograms for the direct and optimized parameter estimation methods in V3 is shown in Fig. 11. As expected, optimized computations converge near values obtained by the direct computations while lowering the fitting error. Further spatial analysis in Section 5.3 (Table 2) supports these observations.

5.2.3 Transmural Histograms Fig. 12 shows the result of sampling the nine connection forms cijk throughout the cardiac wall based on distance to the epicardium. A normalized marginal histogram is obtained for each distance by sampling all voxels that correspond to a given value of the distance transform in Fig. 8. Because of the exploratory and novel nature of these distributions, only a preliminary interpretation can be made of their shape. Any anatomical explanation for these variations is subject to errors in the segmentation (e.g. near the septum and within the right ventricle) as well as boundary effects near the walls which tend to spread out the distributions. In particular, the c123 histogram shows that the course of the cardiac helix angle undergoes a major transition near the mid-wall. In general, connection forms are quite variable transmurally. An exact spatial localization and biological explanation of these variations remains to be fully explored. By looking at spatial and scale dependence of the connection forms in the following section, we can however gain insight into the large-scale architecture of the cardiac frame field and answer some of these questions. 5.3 Scale and Spatial Dependence The direct computations of connection forms described in Section 3.1 depend only on the nature of the differential kernel used. On the other hand, the optimized parameter computations of Section 3.2 depend on the size and shape of the neighborhood V in which the energy (20) is computed. V can be adjusted, and a filtering of the principal direction of

2500

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 37,

NO. 12,

DECEMBER 2015

TABLE 1 Effect of Neighborhood Size on Optimized Connection Forms (Radians/Voxel) and on Fitting Errors (Radians) for the Rat Dataset as Mean  Standard Deviation V

3

c121 c122 c123 c131 c132 c133 c231 c232 c233 1 2 3

:008  :083 :015  :091 :284  :237 :039  :043 :016  :047 :017  :063 :000  :032 :031  :046 :013  :071 :014  :030 :013  :028 :005  :014

5 :010  :113 :022  :113 :304  :236 :040  :030 :016  :032 :017  :047 :001  :022 :028  :031 :010  :053 :029  :042 :027  :040 :009  :014

7 :013  :202 :035  :244 :310  :313 :040  :029 :012  :034 :014  :038 :000  :018 :025  :024 :005  :042 :047  :053 :046  :051 :013  :015

i of (21) are shown for the full connection form model embedding

diffusion can be performed to better target the scale of features that are to be extracted. This section thus investigates some of the possibilities in shaping the differentiation kernel for filtering and in selecting the energy neighborhood for computing the 9 cijk connection forms using different isotropic neighborhoods Vi . In doing so, we can measure the smoothness of the frame field, and obtain a preliminary spatial localization of its geometrical variation.

5.3.1 Neighborhood Shape Table 1 shows the mean and standard deviation for all voxels of the dataset, for each connection form and various neighborhoods Vi . Fitting errors i of (21) are shown for the full connection form embedding. The results indicate globally stable mean values, although a local spatial analysis should be investigated to draw further conclusions. All errors increase following neighborhood size, meaning that increasingly more information is lost to the nonlinear Oðkvvk2 Þ term in (7). This indicates that a firstorder linear model of the natural cardiac frame field is not practical to describe variability beyond a few voxels, given the current resolution at which the data was acquired.

Fig. 12. Transmural histogram of combined volumes of cijk  cij hff k i and fitting errors in V3 for the rat data sets using optimized parameter computations. The horizontal axis represents voxel depth from the epicardium and the vertical axis represents volume values. The mode of a Gaussian distribution fit to the marginal histogram for each distance is traced in blue.

Fig. 13 explores the relationship between neighborhood size and the various embeddings of Section 4. Extrapolation errors increase with increasing Vi , similarly for all embeddings. This linear relationship among different embeddings arises from the fact that most of the frame field variability is captured by the c123 connection, which is considered in each embedding (with the exception of the constant model), and that added complexity does not proportionally contribute to a diminishing energy measure in (21), a relationship that will be investigated in Section 5.4.

5.3.2 Filtering Diffusion Directions Filtering of diffusion directions after adopting cylindrical consistency using (36) can be applied to compensate for the

TABLE 2 Effect of Iterative Gaussian Smoothing Applied to f1 on Extrapolation Errors and Connection Forms for the Direct and Optimized Parameter Computations in V3 s¼0 direct c121 c122 c123 c131 c132 c133 c231 c232 c233 1 2 3

:006  :078 :013  :086 :235  :183 :037  :054 :015  :055 :017  :072 :000  :034 :028  :048 :013  :068 :024  :053 :021  :052 :007  :021

optimized :008  :083 :015  :091 :284  :237 :039  :043 :016  :047 :017  :063 :000  :032 :031  :046 :013  :071 :014  :030 :013  :028 :005  :014

pffiffiffiffiffi s ¼ 0:2 10 direct optimized :002  :076 :008  :086 :178  :232 :039  :040 :012  :044 :013  :067 :001  :032 :027  :047 :013  :062 :038  :082 :035  :081 :008  :023

:006  :097 :011  :103 :205  :271 :040  :036 :013  :040 :013  :056 :001  :029 :029  :043 :015  :056 :035  :087 :033  :086 :006  :020

pffiffiffiffiffi s ¼ 0:4 10 direct :001  :047 :004  :046 :133  :163 :044  :026 :010  :029 :008  :045 :002  :029 :024  :043 :012  :045 :041  :066 :038  :064 :008  :020

optimized :001  :053 :003  :047 :135  :169 :044  :028 :009  :031 :008  :044 :002  :031 :025  :041 :012  :047 :041  :066 :038  :065 :008  :021

PIUZE ET AL.: MAURER-CARTAN FORMS FOR FIELDS ON SURFACES: APPLICATION TO HEART FIBER GEOMETRY

2501

Fig. 13. Normalized distributions of the extrapolation error 1 (in radians) for different neighborhoods Vi .

effect of noise in the diffusion volumes and to investigate the scale of fine muscle cardiac structures. We refer the reader to [31] for an overview of diffusion based smoothing schemes for vector-valued images. A specialized literature also exists for the case of smoothing diffusion weighted images or smoothing the diffusion tensors directly or indirectly, (see, e.g., [32]). In the context of cardiac tissue, whereas the principal direction of diffusion is widely accepted as corresponding to the orientation of cardiomyocytes, the second and third eigenvectors exhibit a high degree of spatial variability [26] and their relationship to biology is not fully established. Thus, in the following we opt for a much simpler smoothing strategy, one that focuses on the first principal direction of diffusion only. The method we employ is based on an element-wise iterative normalized convolution with a Gaussian kernel Gsþ with standard deviation s þ on the vector field f 1 with a renormalization step, before computing the connection forms. The following ¼ Gs þ update equation is used with f 01  f 1 , f nþ1 1 f n1 =kGsþ f n1 k: The remaining fields f 2 ; f 3 are then obtained pffiffiffi using (37). We will make use of the notation s ¼ s þ n to denote n iterations of standard deviation s þ . As an example, Fig. 15 shows the effect of filtering on pffiffiffiffiffi f 1 for s ¼ 0:4 10. Fig. 16 shows a volumetric helical sampling of selected connection forms and fitting errors, computed with respect to the unfiltered data, for various filtering magnitudes s i . As connection forms become smoother following an increase in filtering magnitude, geometrical measurements are obtained which pertain to large scale myofiber arrangements rather than local ones. Although the vast majority of errors are small (less than 0.05 radians), the last row shows that the fitting error is higher near the right ventricle, which is in part due to an incorrect labeling of background or interior voxels as lying within the myocardium. As a consequence, as s i is increased values in these regions are propagated to their neighbors and are largely responsible for the overall pffiffiffiffiffi increase in fitting error. In particular, for c123 at s 4 ¼ 0:4 10 there is a false suggestion that the voxels turn in a clockwise direction in these mislabeled regions. Fig. 14 and Table 2 summarize these results as means over all volumes for each connection form and fitting error.

5.4 Model-Space Exploration We now examine the contribution of each individual connection form in lowering the mean fitting error, by settings all others to zero. As expected, each cijk only affects the respective i and j errors. Fig. 17 demonstrates that the most important connection form in lowering both 1 and 2 is c123 , since it describes a rotation of the helix angle, which has a large magnitude throughout the volume. c121 and c122 also significantly

Fig. 14. Method comparison and effect of f 1 filtering on mean ffiffiffi pffiffiffi cijkpand inffiffiffi V3 p for i p ffiffiffi the dataset, with ðs 0 ; s 1 ; s 2 ; s 3 ; s 4 Þ ¼ ð0; 0:2 5; 0:2 10; 0:4 5; 0:4 10Þ.

lower the error, which supports the use of low dimensional embeddings such as the GHM described in Section 4.1. The smaller contribution of all c13k forms in lowering the fitting error is likely due to the local scale and the isotropic neighborhoods in which the errors were measured, and to the small magnitude of the geometric features that they measure. 3 reduction is more variable although considerably small because of the stability of the heart wall normal. The ability of the various Maurer-Cartan embeddings to lower the fitting error using the cardiac frame field can be

Fig. 15. Selected fiber directions f 1 (in gray) for the p ratffiffiffiffiffiA, before (left) and after (right) applying iterative filtering (s ¼ 0:4 10). Hue is the azimuthal angle of the orientation, which represents directions upward perpendicular to (green), downward perpendicular to (red), and aligned (blue) with the short axis plane.

2502

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 37,

NO. 12,

DECEMBER 2015

Fig. 16. Effect of f 1 filtering on selected connection forms (rads/voxel) and extrapolation error 1 (rads). Selected short-axis slices are shown pffiffiffiin columns A, and a full volume histogram is shown for the species in column 4, where s 0 ¼ 0; s 1 ¼ 0:2 5; s 2 ¼ pffiffiffi pffiffiffi 1 to 3 for pffiffiffi the rat subject 0:2 10; s 3 ¼ 0:4 5; s 4 ¼ 0:4 10.

PIUZE ET AL.: MAURER-CARTAN FORMS FOR FIELDS ON SURFACES: APPLICATION TO HEART FIBER GEOMETRY

2503

Fig. 17. The individually optimized connection form and the mean errors of fit i in V3 using the constant Maurer-Cartan embedding as a baseline. Boxes illustrate the mean and twice the standard deviation.

predicted directly from Fig. 17. Embeddings that offer greater complexity in capturing the variation of a particular frame axis will lower the frame error associated with it. Fig. 18 shows the three mean errors i of various embeddings as a function of neighborhood size Vi , using optimized computations. In relation to the frame error 1 , all embeddings with an unconstrained c123 will perform similarly, as shown in Fig. 18, since it is the connection form that dominates most of the frame field variation. On the other hand, the third frame vector f 3 captures the extent to which the principal direction of diffusion remains normal to the heart wall normal. Since fiber directions u1 generally run tangentially to the heart wall, the resulting error is low for all embeddings.

5.5 Summary Measurements of the geometry of cardiac myofibers were made using the method of moving frames. These results corroborate and extend existing cardiac literature, most of which have not been reported before. More precisely, the following was observed. (1)

(2)

Helix angle. c123 measures the rate of change of the helix angle, and is in the order of 0:290 rad/voxel. In a typical heart from the dataset, the average transmural depth from apex to sub-atria amounts to about seven voxels. Integrating c123 throughout this distance produces a total change of 116:3 degree, which is in close correspondence with values of 120 degree reported in the literature. Wall curvatures. c131 and c232 reflect the sectional curvatures of the heart wall as projected onto the local osculating ellipsoid to f 1 . Their respective mean values of 0.039 and 0.031 rads/voxel imply radii of curvature of 26 and 32 voxels. These values are in the range of the two principal radii of an aligned ellipsoid to a typical heart in the dataset, which at the mid-region has a circular shape and a half-width of about 25 voxels.

(3)

(4)

6

Myocyte fanning. c122 and c133 are measures of how much cardiac myocytes fan out, and their mean values are 0:015 and 0:017 rads/voxel. Based on histological studies, these values are expected to be small, since myocytes form a largely homogeneous, parallel medium. Myocyte twising. c123 , c132 and c231 are measures of twising in the collagen matrix that contains cardiac myocytes. Whereas c123 corresponds to the variational component of the helix angle, c132 describes a turning that is directed in an upward fashion from base to apex, and c231 one that rotates the tangent fiber plane.

CONCLUSION

We employed the method of moving frames in Sections 2 and 3 to explore and compute the variation of a smooth frame field. In Section 4 the method was used to develop a selection of geometrical embeddings which imposed certain constraints—via thresholding (spherical and ellipsoidal manifolds) or by assuming a functional dependency (GHM) on connection forms. Although we limited ourselves to these simple cases, there are many additional possibilities for developing such embeddings. By carefully tailoring connection form constraints to the application at hand, one can design a powerful geometrical probe. This analysis can be performed whenever smooth flow lines or trajectories need to be interpreted geometrically. Different embeddings were employed in Section 5 to characterize the geometry of cardiac myocytes. It was shown that the resulting connection forms relate to established cardiac measurements, many of which were until now only determined from indirect empirical data. These measurements include the rate of change of the helix and transverse angles, measures of cardiomyocyte fanning and twisting, and sectional heart wall curvatures. This research yields new possibilities for differentiation between healthy and pathological cardiac tissue, and for generating models of synthetic cardiac orientations. A more thorough spatial analysis of the localization, clustering and co-dependence of connection forms, related to the rotation of the cardiac frame field in both healthy and diseased hearts is a fruitful direction for future work.

ACKNOWLEDGMENTS Fig. 18. Mean extrapolation error i as a function of neighborhood size V for various embeddings.

The authors are grateful to FQRNT (Quebec) and NSERC (Canada) for research funding. Emmanuel Piuze is the corresponding author.

2504

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

REFERENCES [1] [2] [3] [4] [5] [6] [7] [8]

[9] [10] [11] [12]

[13] [14] [15] [16] [17] [18] [19]

[20]

[21] [22] [23]

[24]

E. Cartan, La Theorie des Groupes Finis et Continus et la Geometrie Differentielle Traitees par la Methode du Repere Mobile, vol. 18, Paris, 1937. P. J. Olver and J. Pohjanpelto, “Maurer–Cartan forms and the structure of lie pseudo-groups,” Selecta. Math., vol. 11, no. 1, pp. 99–126, 2005. O. Faugeras, Cartan’s Moving Frame Method and Its Application to the Geometry and Evolution of Curves in the Euclidean, Affine and Projective Planes. New York, NY, USA: Springer, 1994. T. Flash and A. A. Handzel, “Affine differential geometry analysis of human arm movements,” Biol. Cybern., vol. 96, no. 6, pp. 577– 601, 2007. M. Boutin and P.-L. Bazin, “Structure from motion: A new look from the point of view of invariant theory,” SIAM J. Appl. Math., vol. 64, no. 4, pp. 1156–1174, 2004. O. Ben-Shahar and S. W. Zucker, “The perceptual organization of texture flow: A contextual inference approach,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 25, no. 4, pp. 401–417, Apr. 2003. P. Savadjiev, G. L. Kindlmann, S. Bouix, M. E. Shenton, and C.-F. Westin, “Local white matter geometry from diffusion tensor gradients,” NeuroImage, vol. 49, no. 4, pp. 3175–3186, 2010. E. Calabi, P. J. Olver, C. Shakiban, A. Tannenbaum, and S. Haker, “Differential and numerically invariant signature curves applied to object recognition,” Int. J. Comput. Vision, vol. 26, no. 2, pp. 107– 135, 1998. P. J. Olver, “Joint invariant signatures,” Found. Comput. Math., vol. 1, no. 1, pp. 3–68, 2001. E. Piuze, J. Sporring, and K. Siddiqi, “Moving frames for heart fiber geometry,” Inform. Process. Med. Imaging, vol. 23, pp. 524– 535, 2013. E. Piuze, H. Lombaert, J. Sporring, G. J. Strijkers, A. J. Backermans, and K. Siddiqi, “Atlases of cardiac fiber differential geometry,” in Proc. 7th Int. Conf. Functional Imag. Modeling Heart, 2013, pp. 442–449. P. Savadjiev, G. J. Strijkers, A. J. Bakermans, E. Piuze, S. W. Zucker, and K. Siddiqi, “Heart wall myofibers are arranged in minimal surfaces to optimize organ function,” Proc. Nat. Acad. Sci. USA., vol. 109, no. 24, pp. 9248–9253, 2012. P. Savadjiev, S. W. Zucker, and K. Siddiqi, “On the differential geometry of 3d flow patterns: Generalized helicoids and diffusion mri analysis.” in Proc. Int. Conf. Comput. Vis., 2007, pp. 1–8. H. Flanders, Differential Forms with Applications to the Physical Sciences. Mineloa, New York, USA: Dover, 2012. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed. Cambridge, U.K.: Cambridge Univ. Press, 2007. M. J. Powell, “The BOBYQA algorithm for bound constrained optimization without derivatives,” Cambridge NA Rep. NA2009/ 06, Univ. Cambridge, Cambridge, U.K., 2009. T. Needham, Visual Complex Analysis. Oxford, U.K.: Clarendon, 1998. E. Piuze, P. G. Kry, and K. Siddiqi, “Generalized helicoids for modeling hair geometry,” Comput. Graph. Forum, vol. 30, no. 2, pp. 247–256, 2011. D. D. Streeter, “Gross morphology and fiber geometry of the heart,” in Handbook of Physiology, Section 2. The Heart., R. M. Berne and N. Sperelakis, Eds. New York, NY, USA: Williams and Wilkins, 1979, pp. 61–112. S. H. Gilbert, A. P. Benson, P. Li, and A. V. Holden, “Regional localisation of left ventricular sheet structure: Integration with current models of cardiac fibre, sheet and band structure,” Eur. J. Cardio-Thoracic Surg., vol. 32, no. 2, pp. 231–249, 2007. A. Horowitz, M. Perl, and S. Sideman, “Geodesics as a mechanically optimal fiber geometry for the left ventricle,” Basic. Res. Cardiol., vol. 88, no. suppl. 2, pp. 67–74, 1993. C. S. Peskin, “Mathematical aspects of heart physiology,” Courant Inst. Math. Sci., New York Uni., New York, NY, USA, Tech. Rep., 1975. D. F. Scollan, A. Holmes, R. Winslow, and J. Forder, “Histological validation of myocardial microstructure obtained from diffusion tensor magnetic resonance imaging,” Am. J. Physiol. Heart Circulatory Physiol., vol. 275, no. 6, pp. H2308–H2318, 1998. E. Hsu, A. Muzikant, S. Matulevicius, R. Penland, and C. Henriquez, “Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation,” Am. J. Physiol. Heart Circulatory Physiol., vol. 274, no. 5, pp. H1627–H1634, 1998.

VOL. 37,

NO. 12,

DECEMBER 2015

[25] L. Geerts, P. Bovendeerd, K. Nicolay, and T. Arts, “Characterization of the normal cardiac myofiber field in goat measured with mr-diffusion tensor imaging,” Am. J. Physiol., Heart Circ. Physiol., vol. 283, pp. H139–145, 2002. [26] H. Lombaert, J. Peyrat, P. Croisille, S. Rapacchi, L. Fanton, F. Cheriet, P. Clarysse, I. Magnin, H. Delingette, and N. Ayache, “Human atlas of the cardiac fiber architecture: Study on a healthy population,” IEEE Trans. Med. Imag., vol. 31, no. 7, pp. 1436–1447, Jul. 2012. [27] J. C. Walker, J. M. Guccione, Y. Jiang, P. Zhang, A. W. Wallace, E. W. Hsu, and M. B. Ratcliffe, “Helical myofiber orientation after myocardial infarction and left ventricular surgical restoration in sheep,” J. Thoracic Cardiovascular Surg., vol. 129, no. 2, pp. 382–390, 2005. [28] J. Chen, W. Liu, H. Zhang, L. Lacy, X. Yang, S.-K. Song, W. A. Wickline, and X. Yu, “Regional ventricular wall thickening reflects changes in cardiac fiber and sheet structure during contraction: Quantification with diffusion tensor MRI,” Am. J. Physiol., Heart Circ. Physiol., vol. 289, pp. H1898–1907, 2005. [29] D. Streeter and D. Bassett, “An engineering analysis of myocardial fiber orientation in pig’s left ventricle in systole,” Anatomical Rec., vol. 155, no. 4, pp. 503–511, 2005. [30] I. J. LeGrice, B. H. Smaill, L. Z. Chai, S. G. Edgar, J. B. Gavin, and P. J. Hunter, “Ventricular myocyte arrangement and connective tissue architecture in the dog,” Am. J. Physiol., Heart Circ. Physiol., vol. 269, pp. H571–582, 1995. [31] R. Deriche and D. Tschumperle, “Diffusion PDE’s on vectorvalued images: Local approach and geometric viewpoint,” IEEE Signal Process. Mag., vol. 19, no. 5, pp. 16–25, 2002. [32] M. Martin-Fernandez, E. M. noz Moreno, L. Cammoun, J.-P. Thiran, C.-F. Westin, and C. Alberola-L opez, “Sequential anisotropic multichannel wiener filtering with Rician bias correction applied to 3d regularization of DWI data,” Med. Image Anal., vol. 13, pp. 19–35, 2009. Emmanuel Piuze-Phaneuf received the BSc degree in joint honors physics and computer science in 2009, and the MSc degree in computer graphics in 2011, both from McGill University. He is currently working toward the PhD degree in computer science at McGill University, and is a member of McGill’s Centre for Intelligent Machines. His research interests lie in computer vision, computer graphics and medical imaging.

Jon Sporring received the PhD degree in 1998 in computer science from the University of Copenhagen, Denmark. He has held positions at IBM Research Center, Almaden; Computer Vision and Robotics Lab at FORTH, Greece; 3D-Lab, School of Dentistry, University of Copenhagen; Nordic Bioscience a/s, Denmark; School of Computer Science, McGill University, Canada. He is currently an associate professor at the Department of Computer Science, University of Copenhagen, and a co-founder and CSO of DigiCorpus Aps. Kaleem Siddiqi received the BS degree from Lafayette College in 1988 and the MS and PhD degrees from Brown University in 1990 and 1995, respectively, all in the field of electrical engineering. He is currently a professor, William Dawson scholar, and an associate director research of the School of Computer Science at McGill University. He is also a member of McGill’s Centre for Intelligent Machines. Before moving to McGill in 1998, he was a postdoctoral associate in the Department of Computer Science at Yale University (1996-1998) and held a visiting position in the Department of Electrical Engineering at McGill University (1995-1996). His research interests include the areas of computer vision, image analysis, and medical imaging. He is a member of the Phi Beta Kappa, Tau Beta Pi, and Eta Kappa Nu. He is a senior member of the IEEE. " For more information on this or any other computing topic, please visit our Digital Library at www.computer.org/publications/dlib.