Methods of modelling and characterising heart fiber geometry

Abstract

The identification and determination of aspects of the construction of a patient's heart is important for cardiologists and cardiac surgeons in the diagnosis, analysis, treatment, and management of cardiac patients. For example minimally invasive heart surgery demands knowledge of heart geometry, heart fiber orientation, etc. While medical imaging has advanced significantly the accurate three dimensional (3D) rendering from a series of imaging slices remains a critical step in the planning and execution of patient treatment. Embodiments of the invention construct using diffuse MRI data 3D renderings from iterating connections forms derived from arbitrary smooth frame fields to not only corroborate biological measurements of heart fiber orientation but also provide novel biological views in respect of heart fiber orientation etc.

Description

FIELD OF THE INVENTION

This invention relates to heart fiber geometry and more particularly to improved methods of optimizing connection forms in local neighbourhoods and employing these to statistically analyse heart fiber geometry using diffuse magnetic resonance imaging data.

BACKGROUND OF THE INVENTION

Cardiac myofibers are densely packed in the heart wall and are locally aligned to helical curves. Helices act as geodesics between points in the myocardium and mathematical analyses support the view that this alignment is mechanically optimal. As a result, geometric descriptions of cardiac fibers using the helix angle, taken to be the projected angle between the fiber direction and the short-axis plane (see first image 100A in FIG. 1), are popular within the prior art. Several accounts from both small-scale histology and voxel-scale studies based on Diffusion MRI (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°. The range of the transverse angle, which is the angle formed by a fiber moving away from a plane perpendicular to the transmural direction, is much smaller, about ±10°, and is therefore often ignored in the literature,

The analysis of myofibers from histological slices is cumbersome and their invasiveness does not easily admit an association with the original intact three-dimensional geometry. Thus, many modern analysis methods work with cardiac fiber orientation data derived from dMRI measurements. However, the scale at which current dMRI measurements are made is at least one order of magnitude larger than the length of individual cardiomyocytes The measured signal therefore reflects the composite behaviour of large groups of cardiac muscle cells within the collagen matrix (see third image 100C in FIG. 1). Within the prior art a promising characterization of the collective geometrical variation of cardiac fibers, using a method derived from texture flow analysis, has been reported. This work concluded that that the cardiac fiber directions across three mammalian species, the rat, the canine and the human, are locally described by a particular minimal surface, the generalized helicoid model (GHM).

However, a limitation of the GHM is that its streamlines lie on a planar manifold in spite that the heart wall is curved (see first image 100A in FIG. 1). The GHM thus captures the variability of cardiac fibers in a plane tangent to the local cardiac wall but not orthogonally to it. Moreover, experimental results have shown that the GHM is only accurate in the immediate neighborhood of a voxel, with fitting errors growing rapidly as the neighborhood in which the fits are applied is increased.

Accordingly it would be beneficial to provide an alternate model for modeling cardiac fibers within heart walls that removes the limitations of GHM.

SUMMARY OF THE INVENTION

It is an object of the present invention to heart fiber geometry and more particularly to improved methods of optimizing connection forms in local neighbourhoods and employing these to statistically analyse heart fiber geometry using diffuse magnetic resonance imaging data.

The inventors have attached a local frame field to the fiber and transmural directions and study the full differential geometry of this moving frame through the Maurer-Cartan connection one-forms. Further, the inventors demonstrate that form-based models on which the streamlines lie on ellipsoidal shells or homeoids give lower fitting errors than the generalized helicoid model fits.

More precisely within the method of moving frames, the Maurer-Cartan form is an operator that measures the differential structure of a manifold which is typically applied in a forward manner to study the geometrical characteristics of the manifold under consideration. However, the inventors have applied the theory of moving frames in the reverse direction such that the Maurer-Cartan connection forms are manipulated to generate manifolds or embeddings based on certain assumptions of their differential structure. Beneficially, this approach allows characterization of a smooth frame field in three dimensions as a parameterization on the space of Maurer-Cartan connection forms. By introducing various connection form embeddings for studying frame fields, the inventors show that GHM can fully describe locally using a combination of connection forms. Further, the inventors introduce a fitting energy for estimating connection forms from arbitrary smooth frame fields. Based upon applying these frame fields to dMRI data and studying its differential geometry, not only can the inventive methods described corroborate biological measurements of heart fiber orientation but they also provide novel biological views in respect of heart fiber orientation.

In accordance with an embodiment of the invention there is provided a method

acquiring data relating to a heart;

establishing a model relating to the heart;

establishing a fitting energy associated with the model;

iterating the model in dependence upon a mathematical process and the acquired data;

calculating a new fitting energy associated with the iterated model; and

determining whether the new fitting energy meets a predetermined threshold, wherein

when the predetermined threshold is met the iterative process is terminated; and

when the predetermined threshold is not met the iterative process continues.

In accordance with an embodiment of the invention there are provided computer executable instructions stored in a non-volatile, non-transient computer memory, the computer executable instructions for execution by a microprocessor and relating to a process comprising the steps of:

acquiring data relating to a heart;

establishing a model relating to the heart;

establishing a fitting energy associated with the model;

iterating the model in dependence upon a mathematical process and the acquired data;

calculating a new fitting energy associated with the iterated model; and

determining whether the new fitting energy meets a predetermined threshold, wherein

when the predetermined threshold is met the iterative process is terminated; and

when the predetermined threshold is not met the iterative process continues.

Other aspects and features of the present invention will become apparent to those ordinarily skilled in the art upon review of the following description of specific embodiments of the invention in conjunction with the accompanying figures.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the present invention will now be described, by way of example only, with reference to the attached Figures, wherein:

FIG. 1 depicts the definition of the helix angle, the transmural change in the helix angle in directions orthogonal to the heart wall tangent plane and a histological slice of cardiac tissue showing individual elongated cardiomyocytes and their nuclei (dark) in the intercellular collagen network;

FIG. 2 depicts histogram of c_ijk≡c_ijF_k for rat, human, and with measures normalized such that the horizontal axis is in millimeters;

FIG. 3 depicts the one-form extrapolation error for neighborhoods of size |N_i=3³, 5³, 7³;

FIG. 4 depicts the one-form approximation to the GHM for the neighborhood N₃ together with the error of each model a function of N_i;

FIG. 5 depicts error of fit for the different models analyzed for the human heart where the results are shown for isotropic voxel neighbourhoods of size 3³, 5³, 7³;

FIG. 6 depicts log-normal fits of the human heart extrapolation error;

FIG. 7 depicts the rotational angle φ of the different models analyzed for the human heart in isotropic voxel neighbourhoods of size 3³, 5³, 7³;

FIG. 8 depicts joint histograms of the angular errors e and φ for the human heart for different models for isotropic 3³ neighborhoods;

FIG. 9 depicts the geometry characterized by connection forms;

FIG. 10 depicts parameters after 200 Nelder-Mead iterations including convergence plots for c₁₂₃ and ε₁ and computational timings for different estimation techniques;

FIG. 11 depicts solutions for 3 neighboring spheres ρ=1, 1.15, 1.3 where K₁ models the turning towards the tangent vector f₁ and each column represents a distinct tuple of (K₂, K₃) which expresses the turning of the tangent vector towards f₂ when moving along the f₂ (spreading) and f₃ (turning) directions;

FIG. 12 depicts solutions to Equations (41) and (42) on the unit sphere, with initial frames are shown with column 1 depicts a single flow line whilst columns 2-4 show neighbouring flow lines from an initial starting point for various (K₁, K₂) tuples;

FIG. 13 depicts examples of ellipsoidal solutions integrated from a flat patch (first row), a spherical patch (second row), and a cylindrical patch (third row) for various values of c₁₃₁, c₁₃₂, c₂₃₁, c₂₃₂;

FIG. 14 depicts the helix angle definition employed within the specification;

FIG. 15 depicts the principal fiber direction and cardiac frame fields for a rat heart obtained from dMRI data;

FIG. 16 depicts endocardium (left) and epicardium (right) Euclidean distance, ranging from zero to 20 voxels and distance gradient vectors;

FIG. 17 depicts flow lines in the myocardium, found by tracking the first eigenvector of diffusion wherein a frame indicates the axial plane of the heart, the transmural direction, and the direction from the base to the apex;

FIG. 18 depicts volume histograms of c_ijk connections (radians/voxel) and fitting errors (radians) for Ω₃ by fusing all rats in the dataset using optimized parameter computations together with the minimum and maximum dataset envelope for each value;

FIG. 19 depicts selected dataset volumes of rat connection forms c_ijk (radians/voxel) and fitting errors ε_i (radians) in χ₃ for the optimized and direct computations;

FIG. 20 depicts transmural histograms of combined volumes of c_ijk≡c_ijf_k and fitting errors in Ω₃ for the rat data sets using optimized parameter computations;

FIG. 21 depicts normalized distributions of the extrapolation error ε_i (in radians) for different neighborhoods Ω_i;

FIG. 22 depicts a method comparison and effect of f₁ filtering on mean c_ijk and ε_i in the  ₁ for different standard deviations of the vector field f₁;

FIG. 23 depicts selected fiber directions f₁ for the rat A, before (left) and after (right) applying iterative filtering;

FIG. 24 depicts the effect of f₁ filtering on selected connection forms (rads/voxel) and extrapolation error ε_i (rads) wherein selected short-axis slices are shown in columns 1 to 3 for the rat subject A, and a full volume histogram is shown for the species in column 4;

FIG. 25 depicts the individually optimized connection form and the mean errors of fit using the constant Maurer-Cartan embedding as a baseline wherein the boxes illustrate the mean and twice the standard deviation; and

FIG. 26 depicts the mean extrapolation error as a function of neighborhood size for various embeddings.

DETAILED DESCRIPTION

The present invention is directed to heart fiber geometry and more particularly to improved methods of optimizing connection forms in local neighbourhoods and employing these to statistically analyse heart fiber geometry using diffuse magnetic resonance imaging data.

The ensuing description provides exemplary embodiment(s) only, and is not intended to limit the scope, applicability or configuration of the disclosure. Rather, the ensuing description of the exemplary embodiment(s) will provide those skilled in the art with an enabling description for implementing an exemplary embodiment. It being understood that various changes may be made in the function and arrangement of elements without departing from the spirit and scope as set forth in the appended claims.

A1. The Maurer-Cartan Form

The inventors characterize the differential geometry of fibers in the heart wall by measuring the manner in which they turn locally. Accordingly, a frame field is constructed F₁, F₂, F₃ε³, F_i·F_j=δ_ij, where δ_y is the Kronecker delta, in such a manner that the turning of the frame field characterizes the turning of the heart wall fibers. The frame field is then expressed as a rotation of the Cartesian frame [e₁e₂e₃] in Equation (1) where the attitude matrix AεSO(3) is a smoothly varying orthonormal matrix, and where the basis vectors e_j are treated as symbols such that F_i=a_ije_j. As a result the differential geometry of the fibers is new directly characterized by the attitude transformation. Its differential structure is given by Equation (2) where d is the differential operator, A⁻¹=A^T, C=(dA)A⁻¹ is the Maurer-Cartan form, and where for simplicity the notation dF_i=Σ_jc_ijF_j is used.

$\begin{array}{cc}{\left[F_1{F}_2{F}_3\right]}^T={A\left[e_1,e_2,e_3\right]}^T& \left(1\right)\\ d\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\end{array}\right]=\left[\begin{array}{c}d{F}_1\\ d{F}_2\\ d{F}_3\end{array}\right]=\left(dA\right)\left[\begin{array}{c}{e}_x\\ e_y\\ e_z\end{array}\right]=\left(dA\right)A^{-1}\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\end{array}\right]=C\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\end{array}\right]& \left(2\right)\end{array}$

The Maurer-Cartan matrix is skew symmetric, i.e. C=−C^T. Hence it has at most 3 independent, non-zero elements: c₁₂, c₁₃, and c₂₃. Each c_ij is a one-form in ³ that can be contracted on a vector v=[v₁, v₂, v₃]^Tε³ to yield the initial rate of turn of F_i towards F_j when moving in the direction of v. The inventors denote this contraction c_ij(v), which is found to be c_ij(v)=∇_vF_i·F_j|_x, where xε³ is a point in the fiber field and ∇_vF_i is the covariant derivative of F_i in the direction v. Accordingly, c_ij(v) is defined by Equation (3) where the components of the frame vectors are enumerated as F_i=[F_i1, F_i2, F_i3]^T where δ_x=δ/δx is used to denote partial derivatives. Since we are interested in studying the change of the frame field in the direction of its basis vectors we study the contractions c_ijk≡c_ijF_k. It should be noted that the frame field F₁, F₂, F₃ has 3 degrees of freedom. Since this field roams a 3-dimensional space, a linear model of the spatial change of the frame field must have 9 degrees of freedom, which are embodied in e_ijk.

$\begin{array}{cc}{c}_{\mathrm{ij}}\langle v\rangle =\left[F_{j1}{F}_{j2}{F}_{j3}\right]\left[\begin{array}{c}{\delta }_x{F}_{i1}{\delta }_y{F}_{i1}{\delta }_z{F}_{i1}\\ {\delta }_x{F}_{i2}{\delta }_y{F}_{i2}{\delta }_z{F}_{i2}\\ {\delta }_x{F}_{i3}{\delta }_y{F}_{i3}{\delta }_z{F}_{i3}\end{array}\right]\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]& \left(3\right)\end{array}$

The abstraction and the comprehensiveness in the one-form description of the geometrical behavior of a frame field can be harnessed to develop models that are descriptive of the variability of cardiac fiber orientations across multiple species. Accordingly, in section 3 the inventors introduces a class of fiber models based on one-forms and re-introduce the GHM as a planar approximation to the complete one—form parameterization.

A2. Measures on a Discrete Fiber Frame Field

The inventors analyze hearts represented as diffusion MRI volumes embedded in 3D rectangular lattices with coordinates x=xe₁+ye₂+ze₃=[x, y, z]^TεZ³. A tangent vector T is identified as the principal eigenvector of the diffusion tensor field. Consistency in T amongst voxel neighbours is enforced by adopting an adaptive cylindrical coordinate system. The centroid c_z of the chamber within each short-axis slice, s, is first determined. T(m) is then made to turn clockwise with respect to that centroid through Equation (4) where l(s_z) is the local approximation of the heart's long-axis.T(x)→sign((T×(x−c_z))·l(s_z))T(x)  (4)

For all the hearts that we consider, l(s_z) approximately coincides with the world's z axis. In accordance with the spirit of the GHM method the heart transmural direction B is estimated as the gradient vector of a distance transform produced as follows

• • the binary image (mask) of the heart is closed using mathematical morphological operations;
  • the closest distance to the heart wall is evaluated at every point;
  • the gradient of the distance transform is computed, and
  • the skeletal points of colliding fronts are removed and interpolated by thresholding the magnitude of the gradient vectors.

The normals {circumflex over (B)} are then aligned to point from outer to inner wall. With T and {circumflex over (B)} we specify a local frame as given by Equations (5A) to (5C) where B is the part of {circumflex over (B)} orthogonal to T. From here on, the inventors use the symbols T, N, and B interchangeably with the corresponding symbols F_j Further, the inventors also refer to the local plane spanned by T and N as the tangent plane.

$\begin{array}{cc}{F}_1=\frac{T}{T}& \left(5A\right)\\ F_2=N=F_3\times F_1& \left(5B\right)\\ F_3=B=\frac{\left(\hat{B}-\left(\hat{B}·T\right)T\right)}{\hat{B}-\left(\hat{B}·T\right)T}& \left(5C\right)\end{array}$

A2.1 One-Form Intuition

One-form contractions c_ijk can be interpreted as the amount of turning of F_i towards F_j, in the direction F_k. For example, c_TNB describes a transmural rotation of T towards N, as shown in second image 100B in FIG. 1. c_ijk were computed at each voxel from the discrete fiber frame field combined with Equation (3). The resulting histograms for three species (rat, canine, human) are shown in FIG. 2 and illustrate statistics on the local turning of the frame field. The rotations of T towards N, c_TNBF_k, are intuitively linked to the curvature parameters of the GHM method in the prior art. These rotations intuitively describe the manner in which fibers turn in the tangent plane of the heart, c_TNT in first image 200A in FIG. 2, describes their tangential curvature, c_TNN in first image 200B in FIG. 2, their fanning in the tangent plane, and c_TNN in second image 100B in FIG. 1, their transmural turning or equivalently the rate of change of the helix angle. This is arguably the most salient variation of the cardiac fibers. The rotations of T towards B, c_TBF_k, expresses the turning of the fibers towards the inner wall, c_TBT effectively measures the first the local curvature of the heart, c_TBN describes a twisting of the tangent plane and the rate of change of the transverse angle, and c_TBB measures a fanning or thickening of the local fiber population away from the tangent plane, towards the inner wall. As shown in seventh to ninth images 200G to 200I in FIG. 2 the remaining rotations of N towards B, c_NBF_k are an order of magnitude smaller in all three directions which indicates that the frame field axis N is constrained within the local tangent plane. c_NBT measures a twisting of the tangent plane, c_NBN the second curvature of the heart wall and c_NBB a transmural fanning or thickening.

A2.2 The One-Form Model

The Maurer-Cartan form extrapolates the local shape to first order as given by Equations (6) and (7) where h is an offset from the point at which the frame is expressed and {tilde over (T)}_h represents the predicted direction of this neighbor by the one-form extrapolation to first-order approximation. As within the GHM method the inventors construct an error measure by computing the average angular difference between the measured and predicted directions in an isotropic neighborhood N as given by Equation (8) where |N_i|=i³ for odd iεZ and T_h is the true neighbours measured direction.

$\begin{array}{c}{\tilde{T}}_h=T+c_{\mathrm{TN}}\langle h\rangle N+c_{\mathrm{TB}}\langle h\rangle B\left(6\right)\\ =T+\left(c_{\mathrm{TNT}}h·T+c_{\mathrm{TNN}}h·N+c_{\mathrm{TNB}}h·B\right)N+\\ \left(c_{\mathrm{TBT}}h·T+c_{\mathrm{TBN}}h·N+c_{\mathrm{TBB}}h·B\right)B\left(7\right)\end{array}$ $e\left(N_i\right)=\frac{1}{N_i}\sum _{h\in N_i}\arccos \left(T_h·\frac{{\tilde{T}}_h}{{\tilde{T}}_h}\right)\left(8\right)$

The associated errors of fit for different species are shown in FIG. 3. Diffusion MRI noise and resolution, heart size, and underlying fiber geometry are factors that account for the error disparity across different species. We delay further analysis of these errors until Section 3, where they will be compared against those of the other models the inventors introduce next.

A2.3 The Generalized Helicoid as a Subset of the One-Form Model

The generalized helicoid model of the GHM within the prior art expresses the local fiber direction in a plane tangent to the heart wall. Within the local coordinate frame, the fiber direction at a point x=x_TT+x_NN+x_BBε³ is given as the angle defined in Equation (9) where K_*εR are the GHM curvature parameters.

$\begin{array}{cc}\theta \left(x,K_T,K_N,K_B\right)=\arctan \left(\frac{K_T{x}_T+K_N{x}_N}{1+K_N{x}_T-K_T{x}_N}\right)+K_B{x}_B& \left(9\right)\end{array}$

Direct calculations show that a frame field spanned by T(θ), N(θ), B(θ) has instantaneous turning given by c_TNT(θ)=K_T, c_TNN(θ)=K_N, and c_TNB(θ)=K_B with the remaining one-forms all being zero. The GHM parameters can thus be estimated directly using Equation (3) and the GHM model may be evaluated directly using Equation (7) and central differences as an alternative to the generative model. To compare these two representations, the parameter vector K=(K_T, K_N, K_B) of the GHM was estimated at each voxel using a standard Nelder-Mead optimization scheme. The problem was formulated as the selection of the parameters K which minimize an extension of Equation (8) where {tilde over (T)}_h→{tilde over (T)}_h(K)=(cos θ, sin θ, 0) and θ=θ(K) as given by Equation (9.

The results for the fitting error and a comparison with the one-forms are shown in FIG. 4. These results indicate that the one-form model is able to capture the GHM's parameterization accurately and that it consistently yields lower errors. Note that in addition to an improved fitting method which is continuous rather than discrete the inventive method estimates heart wall normals slightly differently than is done within the GHM prior art. Consequently, the inventors parameter estimates are more precise and still support the overall shape distribution reported within the prior art. Accordingly, within the subsequent description the inventors when considering the GHM methodology will be referring to and using its one-form approximation, which they refer to as the ghm-form.

Accordingly, the inventors then proceeded to introduce a differential model, the homeoid, that can also be expressed using a subset of the one-forms, and has the advantage that it is intuitively connected to the large-scale structure of the heart by enforcing the ellipsoidal topology of the local tangent plane.

A2.4 The Generalized Helicoid on an Ellipsoid is a Homeoid

The calculations within Section 2.3 can be applied to model fibers with smoothly varying fiber orientations, such that the differential operations are well defined. Motivated by evidence that fibers wind around the heart wall while remaining approximately parallel to the tangent plane to the wall at each location the inventors consider a specialization to the case where the fibers lie locally on thin homeoids, which are shells composed of two concentric and similar ellipsoids.

As introduced in Section 1, the Maurer-Cartan form has only 3 independent one-forms: c_TN•, c_TB• and c_NB• each with 3 associated spatial degrees of freedom, for a total of 9 possible combinations. Working with the intuition given in Section 2.2 of each c_ijk, this is a convenient space to develop models of fiber geometry. For example, in Section 2.3 the inventors showed that for the GHM only c_TNT, c_TBT and c_NBT are non-zero. Based on a general description of the cardiac fiber architecture as collections of fibers that (i) vary smoothly and (ii) are locally constrained to the tangent space of smooth and orthogonal surfaces to the heart wall, the following one-form contractions given in Equation (10) must occur.c_TNT=α;c_TNN=β; c_TNB=γ; c_TBB≈0; c_NBB≈0  (10)

Locally, these fibers lie in the tangent plane of a thin homeoid. The parameter fields α, β and γ are introduced as the curvature parameters of the fibers. c_NBB must be zero otherwise the fibers could move in and out of the local tangent plane and the hypothesis (ii) would not be satisfied. The remaining contractions specify the shape of the homeoid and accordingly are given by Equation (11) where ρ_i and ρ₂ are the radii fields of the oscillating ellipsoid. Using Equation (7), the model can be employed to extrapolate the orientation of fibers in the neighborhood of a point x. Constraints given by Equations (10) and (11) are satisfied by enforcing the nullity of c_TBN and c_TBB such that we obtain Equation (12).

$\begin{array}{cc}{c}_{\mathrm{TN}}\langle T\rangle =\frac{1}{{\rho }_1};c_{\mathrm{TN}}\langle N\rangle =\frac{1}{{\rho }_2};c_{\mathrm{TB}}\langle N\rangle =0;c_{\mathrm{NB}}\langle T\rangle =0& \left(11\right)\\ {\tilde{T}}_h=T+c_{\mathrm{TN}}\langle h\rangle N+\left(c_{\mathrm{TBT}}h·T\right)B& \left(12\right)\end{array}$

A3. Model Space Comparison

The analytical models of fiber geometry described so far vary in their parametric complexity. The one-form, homeoid and generalized helicoid models respectively have 9, 5, and 3 parameters. We introduce the constant model which will serve as a base-line to which the remaining models can be compared. This parameter-free model simply assumes {tilde over (T)}_h=T in Equation (7). To compare the different models in terms of their fitting accuracy, we have evaluated each c_ijk on the human data set using first-order central differences on 3³ neighbors. We then used these one-forms to extrapolate each model using Equations (7) and (8) in isotropic neighbourhoods where N_i, where |N_i|=i³ for i=3, 5, 7, 9. FIG. 5 shows a distribution of the error of fit in the human dataset for the different models across increasingly large voxel neighborhoods. Error generally increases with neighborhood size but the relative performance of each model is difficult to assess. Accordingly the inventors therefore fitted a log-normal distribution to each error plot and show the resulting log-normal mode e^μ−σ2 and mean

$e^{\mu +\frac{1}{2}{\sigma }^2}$ as a function of neighborhood size in FIG. 6. As expected, the constant model provides an upper bound on the error of fit and is a measure of the smoothness in the data. The one-form model has the lowest error of fit when the neighborhood size reflects the scale at which central differences were computed but behaves poorly for larger neighborhood sizes, which we attribute to local over-fitting. On the other hand, the ghm-form and the homeoid models are well-behaved for all neighborhood sizes and differ only slightly from one another.

In Section 2.4 the inventors showed that their extrapolated {tilde over (T)}_h axis only differs by the c_TBT one-form which is a measure of the curvature of the heart wall. For the human, this value is small and therefore the two models should be very similar. The rat hearts is smaller in size and therefore has larger per voxel curvature. In this case Table 1 shows that the homeoid is a better fit.

TABLE 1
Extrapolation Errors in Radians for Each Species,
Differential Model and Neighbourhood Ni
$\mathrm{fr}\mathrm{om}\mathrm{the}\mathrm{mode}{e}^{\mu -{\sigma }^2}\mathrm{and}\mathrm{the}\mathrm{mean}{e}^{\mu +\frac{1}{2}{\sigma }^2}\mathrm{of}\log \text{-}\mathrm{normal}\mathrm{fits}.$
		Rat	Dog	Human
|Ni|		(mode, mean)	(mode, mean)	(mode, mean)
33	one-form	0.093, 0.113	0.076, 0.092	0.039, 0.054
	homeoid	0.094, 0.115	0.085, 0.102	0.046, 0.063
	ghm-form	0.099, 0.118	0.089, 0.106	0.050, 0.066
	Constant	0.163, 0.181	0.110, 0.128	0.070, 0.090
53	one-form	0.146, 0.175	0.157, 0.181	0.085, 0.112
	homeoid	0.145, 0.173	0.152, 0.176	0.085, 0.112
	ghm-form	0.150, 0.175	0.147, 0.171	0.085, 0.111
	Constant	0.271, 0.286	0.156, 0.178	0.109, 0.135
73	one-form	0.202, 0.239	0.220, 0.251	0.124, 0.161
	homeoid	0.200, 0.234	0.202, 0.235	0.117, 0.152
	ghm-form	0.206, 0.237	0.189, 0.221	0.114, 0.147
	Constant	0.389, 0.402	0.189, 0.215	0.141, 0.172

A moving frame in ³ has 3 degrees of freedom of which 2 are captured by the error vector T−{tilde over (T)}_h: the angular difference e(N) between dMR1 orientations and extrapolations specified by Equation (8), and a rotation φ about T. The third DOF is the rotation ψ of N about T. ψ strongly depends on the calculations of B and much less on the direct measurements. In contrast to the GHM, the homeoid model considers this angle. However, since we focus on the direct measurement of the fiber geometry given by T, then further investigation of ψ is left outside of this description. φ can be obtained by projecting T−{tilde over (T)}_h, onto the local NB plane and measuring its angle with respect to the frame axis N:

$\varphi =\arctan \frac{\left(T-{\tilde{T}}_h\right)·B}{\left(T-{\tilde{T}}_h\right)·N},i.e.$ φ=(0, π) and

$\left(\frac{\pi }{2},\frac{3\pi }{2}\right)$ respectively indicate alignment with the frame vectors N and B respectively. FIG. 5 showed the marginal distributions of e(x,N) for various neighborhoods sizes and FIG. 7 shows the marginal distributions of φ. In FIG. 8 the inventors show the joint histogram of these marginal distributions for the human heart and for the 4 different models in a 3³ extrapolation neighborhood. It is evident that the error along φ is negligible when e is small. The spread of φ measures the rotation of the local NB plane, which for the one-form model results in over-fitting for larger neighborhood sizes. Accordingly, the inventors have established a one-form model yielding the lowest fitting error for small neighborhoods and that the homeoid model is as accurate as or better than the GHM depending on the per-voxel curvature of the heart wall.

B1. Structure of the Frame Fields

Let a point x=Σ_ix_ie_iε³ be expressed in terms of the natural orthonormal coordinate system e₁, e₂, e₃. A differential orthonormal frame field embedded in ³ is denoted as F=[f₁, f₂, f₃]^T: ³→³ and defined by f₁·f₂=δ_ij, where δ_ij is the Kronecker delta and where · is the inner product, and f₁×f₂=f₃, where x is the 3-dimensional cross product. Since the frame E=[e₁, e₂, e₃]^T forms an orthonormal basis for ³, F can be expressed as the rotation f_i=Σ_ja_ije_j, where the elements of the attitude matrix A A={a_i,j}ε^3×3 are differentiable, and A⁻¹=A^T. Treating f_i and e_j as symbols then this rotation can be written in matrix form as given by Equation (13).

$\begin{array}{cc}{\left[f_1,f_2,f_3\right]}^T={A\left[e_1,e_2,e_3\right]}^T\iff F=\mathrm{AE}& \left(13\right)\\ \mathrm{dF}=\left(\mathrm{dA}\right)E+A\left(\overrightarrow{d}{\overrightarrow{E}}^0\right)=\left(\mathrm{dA}\right)A^{-1}F=\mathrm{CF}& \left(14\right)\\ {\mathrm{df}}_i=\sum _j{c}_{\mathrm{ij}}{f}_j& \left(15\right)\end{array}$

Since each e_i is constant, the differential geometry of the frame field is completely characterized by the attitude matrix A. Now taking the exterior derivative on both sides of Equation (13) we arrive at Equation (14) wherein d denotes the exterior derivative operator, and C=(dA)A⁻¹={c_ij}ε^3×3 is the Maurer-Cartan matrix of one-forms c_ij. Now writing f_i as symbols then Equation (14) can be understood as Equation (15).

The Maurer-Cartan matrix is skew symmetric and accordingly Equation (16) applies such that there are at most 3 independent, non-zero elements: c₁₂, c₁₃ and c₂₃. The differential of the elements a_ij(x): ³→ of A is expressed in terms of dx^k, the natural basis for one-forms in R³, and is given by Equation (17).

$\begin{array}{cc}C=\left[\begin{array}{ccc}0& c_{12}& c_{13}\\ -c_{12}& 0& c_{23}\\ -c_{13}& -c_{23}& 0\end{array}\right]& \left(16\right)\\ {\mathrm{da}}_{\mathrm{ij}}=\sum _k\frac{\delta a_{\mathrm{ij}}}{\delta x_k}d{x}^k& \left(17\right)\end{array}$

One-forms as noted supra are differential operators that may be applied to tangent vectors through a process denoted contraction, written as dw(v)ε for a general one-form dw on ³ and tangent vector vε³. The contraction of dw=Σ_kw_kdx^k on v is written as a bilinear operation on the canonical projection of v: dw(v)=Σ_iw_ide_i(Σ_jv_je_j)=Σ_iw_iv_i. Hence, when contracted, the elements of C in Equation (14) are given by Equation (18) and they express the amount of turning of f_j(x) when x moves in the direction v. This analysis naturally applies to frame fields embedded in ² by setting f₃ to be constant. In that case, C is a 2×2 matrix with a single one-form c₁₂ which describes the in-plane turning of the frame; hence, if f₁ is the tangent of a planar curve and f₂ its normal, then c₁₂(f₁) is the curves curvature and c₁₂(f₂) its fanning.

$\begin{array}{cc}{c}_{\mathrm{ij}}\langle v\rangle =\sum _k{a}_{\mathrm{jk}}\sum _l\frac{\delta a_{\mathrm{ik}}}{\delta x_l}{v}_l& \left(18\right)\\ {{f_i\langle v\rangle }_{x_0}=f_1}_{x_0}+d{f}_i\langle v\rangle +O\left({v}^2\right)& \left(19\right)\\ =f_1{|}_{x_0}+\sum _{j\ne 1}{f}_j\langle v\rangle \sum _k{v}_k{c}_{\mathrm{ijk}}{|}_{x_0}+O\left({v}^2\right)& \left(20\right)\end{array}$

The space of linear models for smooth frame fields is fully characterized by the elements of c_ij of C, and can be generated from the first order terms of a Taylor series of f_i centered at a point x₀. The first order approximation for the motion of a frame vector f_i in a direction v=Σ_kv_kf_k can be expressed as Equations (19) and (20) where we use the short-hand c_ijk≡c_ijf_k for the natural connections of the local frame. Since only 3 unique non-zero combinations of i, j are possible in ³ due to the skew symmetry of C, there are only 9 unique non-zero combinations of c_ijk possible for k=1, 2, 3. FIG. 9 illustrates the behavior of a frame field where each connection form predominates.

Considering the model at a point x₀ for all frame vectors and one-forms, an approximation of the frame from F(x₀) to F(x₀+v) can be obtained by the normalized model given in Equation (21) which in matrix notation yields Equation (22) where CⁱεR^3×3 denotes the matrix of connection forms with c_klⁱ=ε_jkc_kil and ε_ik−sgn(k−i). In general {tilde over (f)}_i(v), {tilde over (f)}₂(v), and {tilde over (f)}₃(v) will not form an orthogonal basis. To establish we can combine Equation (15) with Equation (19) to yield Equation (23) for k≠i, k≠j.

$\begin{array}{cc}{\tilde{f}}_i\left(x_o+v\right)=\frac{f_i+{\mathrm{df}}_i\langle v\rangle }{f_i+{\mathrm{df}}_i\langle v\rangle }& \left(21\right)\\ {\tilde{f}}_i\left(x_o+v\right)=\frac{f_i+{\mathrm{FC}}^i{F}^{T}v}{f_i+{\mathrm{FC}}^i{F}^{T}v}& \left(22\right)\\ {f_i\langle v\rangle }_{x_0}·f_i\langle v\rangle {|}_{x_0}=\{\begin{array}{cc}1+c_{\mathrm{ij}}^2+c_{\mathrm{ik}}^2& \text{:}i=j\\ c_{\mathrm{ik}}{c}_{\mathrm{jk}}& \text{:}i\ne j\end{array}& \left(23\right)\end{array}$

B2. Estimation of Connection Forms

The full space of linear models for smooth frame fields was introduced in Equation (20). It was found that this space has at most 9 independent parameters c_ijk which fully characterize the local geometry of a frame field. Accordingly, we need to compute these parameters for a smooth frame field.

B2.1 Direct Computation of Connection Forms

The connection one-forms c_ij can be directly obtained from the vector fields f_i and their differentials df_i using Equation (14) to yield Equation (24) and therein Equations (25) and (26). The differential df_i can be computed by applying the exterior derivative of a function g: Rⁿ→R as given in Equation (27) which yields Equations (28) to (30) where

$=\left[\frac{\delta x^i}{\delta x_j}\right]\in R^{3\times 3}J$ is the Jacobian matrix of partial derivatives and ds=(dx¹, dx², dx³)^T. Setting v=f_k then we finally obtain Equation (31). The Jacobian matrix J can be approximated to first order using for example finite differences such that

$\frac{\delta f_{\mathrm{ij}}}{\delta x_k}\left(x\right)\approx \frac{1}{2}\left(f_{\mathrm{ij}}\left(x+e_k\right)-f_{\mathrm{ij}}\left(x-e_k\right)\right).$

$\begin{array}{cc}{\mathrm{df}}_i·f_k=\left(\underset{j}{\sum ^3}{c}_{\mathrm{ij}}{f}_j\right)·f_k& \left(24\right)\\ =\stackrel{3}{\sum _j}{c}_{\mathrm{ij}}{\delta }_{\mathrm{jk}}\mathrm{since}{f}_i·f_k={\delta }_{\mathrm{ik}}& \left(25\right)\\ =c_{\mathrm{ik}}& \left(26\right)\\ dg=\sum _i^n\frac{\delta g}{\delta x_i}{\mathrm{dx}}^i& \left(27\right)\\ {\mathrm{df}}_i·f_j\langle v\rangle =\left(\sum _l^3\left(\sum _k^3\frac{\delta f_{\mathrm{ik}}}{\delta x_k}{\mathrm{dx}}^k\right)f_{\mathrm{jl}}\right)\langle v\rangle & \left(28\right)\\ =f_j^{T}J\left(f_i\right)\mathrm{dx}\langle v\rangle & \left(29\right)\\ =f_j^{T}J\left(f_i\right)v& \left(30\right)\\ c_{\mathrm{ijk}}={\mathrm{df}}_i·f_j\langle f_k\rangle =f_j^{T}J\left(f_i\right)f_k& \left(31\right)\end{array}$

B2.2 Computation Via Energy Minimization

The connection forms c_ijk at a point x₀ can also be found as the minimizer of the energy contained within a neighbourhood Ω as established in Equation (32) where ε_i is a function of the error for each frame axis, as given by Equation (33) where f_i(x₀+v) is the frame data term and {tilde over (f)}_i(x₀+v) is its approximation using Equation (22).

$\begin{array}{cc}{c}_{\mathrm{ijk}}^{*}\left(x_0\right)=\arg \underset{c_{\mathrm{ijk}}}{\min }\frac{1}{\Omega }\sum _{v\epsilon \Omega }\stackrel{3}{\sum _i}{\varepsilon }_i\left(x_0+v\right)& \left(32\right)\\ {\varepsilon }_i\left(x_0+v\right)=\arccos \left(f_i\left(x_0+v\right)·{\tilde{f}}_i\left(x_0+v\right)\right)& \left(33\right)\end{array}$

Here, Ω can take any shape. We will denote a cubic isotropic neighborhood in ³ of radius i as Ω_2i+1, and 6 nearest neighbors will be denoted as Ω₃⁺. The optimization energy in Equation (33) can be solved for using standard algorithms such as Nelder-Mead and BOBYQA.

B2.3 Method Comparison

Now referring to FIG. 10 there are compared 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 Ω₃⁺ and then evaluated in Ω₃. Finite difference computations of Section B2.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 3 seconds of computation time. In comparison, BOBYQA optimization takes about 1 minute 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 c_ijk 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 axis differential df_i=d(f_i1e₁+f_i2e₂+f_i3e₃) in f_i is bounded due to Equation (34).

$\begin{array}{cc}{\frac{\delta f_i}{\delta x_k}}_x\approx \frac{f_{\mathrm{ij}}\left(x-e_k\right)-f_{\mathrm{ij}}(x+e_{k)}}{2}\le 1& \left(34\right)\end{array}$

Using ∥·∥₁ to denote the Euclidean one-norm, then the bounds defined by Equations (35A) to (35C) are obtained.|c_i,jv|=|f_jTf_i(x₁,x₂, x₃)v|  (35A)≦f_jT·[∥v∥₁∥v∥₁∥v∥₁]^T  (35B)≦∥∥₁≦α for vεΩ_2α+1  (35C)

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 may be preferable. For the remainder of this specification however, a Nelder-Mead optimization scheme running for 300 iterations will be assumed. The upper bound of Equations (35A) to (35C) will be enforced by discarding volume elements that fall outside of permitted values.

The Pointcaré-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. Whereas in theory characterizations of open sections on manifolds can be made free from singularities, 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 discussed below in Section B4, singularities generate a sharp turning in the frame vectors. This behavior cannot be captured by the first order model of Equation (20) which will therefore only yield a coarse approximation. Using the direct computations of Equations (31), Equation (34) 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 Equation (32) or a regularization term can be added. These strategies should be seen as heuristics, and will in general not yield a good fit to the data close to singularities.

B3. Maurer-Cartan Embeddings

The geometrical intuition conveyed by the c_ijk parameters in Equation (8) and illustrated in FIG. 9 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. Given a seed frame f_i(x₀) with known parameters c_ijk, the solution to the manifold f_i(x) for all x can then be found by integration. In some cases, such as with the GHM as discussed below, there exist closed form solutions. In general, models may be solved for numerically by using Equation (22) in a standard integration scheme.

Now the inventors demonstrate 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.

B3.1. Generalized Helicoids as Connection Forms

The model for representing the geometry of 3D smooth streamline flows within the prior art known as the Generalized Helicoid Model (GHM) is based upon modeling texture flows and has been used to measure axonal geometry in white matter fiber tracts and to characterize muscle directions in the heart by modeling the local variation of a smooth frame field. The inventors show in this section that the GHM model is a subset of Equation (22), where all but 3 parameters are zero.

$\begin{array}{cc}{f}_1\left(x_0+v\right)=\cos \theta f_1\left(x_0\right)+\sin \theta f_2\left(x_0\right)& \left(36A\right)\\ \theta \left(x_0,v\right)=\arctan \left(\frac{K_1\left(x_0\right)v_1+K_2\left(x_0\right)v_2}{1+K_2\left(x_0\right)v_1-K_1\left(x_0\right)v_1}\right)+K_3\left(x_0\right)v_3& \left(36B\right)\\ A\left(\theta \right)=\left[\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -\sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]& \left(37\right)\end{array}$

f₁ is given by the GHM vector field in Equations (36A) and (36B) where θ(x₀, v) parameterizes a minimal surface over parameter fields K_i and where v=Σ_iv_if_i(x₀). The 9 parameters c_ijk of Equation (22) can then be evaluated for f₁. The trivial basis e=f₁(x₀) is chosen in order to construct the attitude matrix A, Equation (37), where θ=θ(x₀,v). Since f₁ stays parallel to the tangent f₁−f₂ plane at x₀ and f₃ remains unchanged, then {tilde over (f)}₂=f₃×f₁ will also stay parallel to that tangent plane at x₀, Hence, we must have c₁₃=c₂₃=0 which is confirmed by straightforward computations. As a result the one-form is found to be defined by Equation (38) since the derivatives are to be evaluated at the origin, and where K_i and df_i are evaluated at x₀. Accordingly we arrive at Equation (39) and conclude that the GHM is a subset of Equation (22), where all but c₁₂₁, c₁₂₂, and c₁₂₃ are zero. As a result, in frame fields where the 6 connection forms c_13k, c_23k|kε{1,2,3} are non-zero, the GHM will not be able to comprehensively characterize the underlying geometry of the frame field. Further other issues can arise with the GHM since the rotation angle θ(x₀, v) yields a frame field singularity when Equation (30) is satisfied.

$\begin{array}{cc}{c}_{12}=K_1{\mathrm{df}}_1+K_2{\mathrm{df}}_2+K_3{\mathrm{df}}_3& \left(38\right)\\ c_{12k}=c_{12}\langle f_k\rangle =K_k,k\in \left\{1,2,3\right\}& \left(39\right)\\ v={\left(1+{\left(\frac{K_1}{K_2}\right)}^2\right)}^{-1}\left(-K_2,K_1,0\right)& \left(40\right)\end{array}$

B3.2 A Model on Spheres

Connection forms can be constrained to generate thin spherical shell embeddings. Given a smooth surface where f₁ and f₂ span the local tangent plane, then we can generate a family of subspaces that are non-intersecting along f₃, by setting c₁₂f_kk, kε{1,2,3} as given in Equations (41A) and (41B) where K_i: ³→ are any smooth functions. The remaining parameters control the motion of frames within parallel spherical surfaces. For a point x on a spherical shell with center at

$x_0-\frac{1}{\rho }{f}_3$ we then obtain Equations (42A) and (42B) respectively. This completely specifies the local linear model for spherical shells using Equation (20). By numerical integration we can generate flow lines of f₁ in a local neighborhood, as illustrated in FIG. 11.c₁₂f₁=K₁c₁₂f₂=K₂c₁₂f₃=K₃  (41A)c₁₃f₃=0c₂₃f₃=0  (41B)

$\begin{array}{cc}{c}_{13}\langle f_1\rangle =\frac{1}{\rho }{c}_{13}\langle f_2\rangle =0& \left(42A\right)\\ c_{23}\langle f_1\rangle =0c_{23}\langle f_2\rangle =\frac{1}{\rho }& \left(42B\right)\\ \left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}\rho \cos \left(\varphi \right)\sin \left(\psi \right)\\ \rho \sin \left(\varphi \right)\sin \left(\psi \right)\\ \rho \cos \left(\psi \right)\end{array}\right]& \left(43\right)\end{array}$

Curves with f₁ as their tangent lie on a single shell and each forms a circle, as shown in FIG. 12. To prove that flow-lines of f_i form circles on a sphere, consider a sphere in polar representation as given by Equation (43) centered at 0 and with radius ρ, and consider the intersection with a plane orthogonal to the z-axis of height −ρ≦z₀ρ. Since z₀=ρ cos(ψ), using Equation (43), the curve of intersection is found to be

$\rho ={\left[\sqrt{{\rho }^2-z_0^2}\cos \varphi ,\sqrt{{\rho }^2-z_0^2}\sin \varphi ,z_0\right]}^T,$ which is a circle of radius R=√{square root over (ρ²−z₀²)} with the z-axis as the axis of rotation. Evaluating the Maurer-Cartan form of a frame field, where

$f_1=\frac{{\vartheta }_{\varphi }p}{{\vartheta }_{\varphi }p},f_3=\frac{{\left(x,y,z\right)}^T}{{\left(x,y,z\right)}^T},$ and f₁=f₃×f₁, we find that c₁₂f₁=z₀/√{square root over (ρ²−z₀²)}, c₁₃f₁=−1/ρ and c₂₃f₁=0 which exactly corresponds to the constraints given by Equations (41A)-(41B) and (42A)-(42B) when K₁=z₀/√{square root over (ρ²−z₀²)}. To align the two coordinate systems, we first set the tangent of the circle to be f_i, and set the change of f_i in the direction f_i to be the normal to the circle. The normal is then derived directly from Equations (22), (41A)-(41B) and (42A)-(42B) as the unit vector parallel to df₁f₁=K₁f₂−(1/ρ)f₃. Finally, the radius of the circle is found by solving for z₀ in terms K₁, z₀=±K₁ρ²/√{square root over (1+K₁²ρ²)} and inserting this into R which yields R=ρ/√{square root over (1+K₁²ρ²)}=1/|df₁f₁.

B3.3 A Model on Ellipsoids

The spherical embedding of Section B3.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₁ and f₂ span the tangent plane which is organized in thin shells, then Equations (41)-(41B) and (42A)-(42B) generalize to Equations (44A) and (44B).

$\begin{array}{cc}{c}_{13}\langle f_1\rangle =\frac{1}{{\rho }_1}{c}_{13}\langle f_3\rangle \ne 0& \left(42A\right)\\ c_{23}\langle f_3\rangle \ne 0c_{23}\langle f_2\rangle =\frac{1}{{\rho }_2}& \left(42B\right)\end{array}$

Here, c₁₃f₃ and c₂₃f₃ are only zero when f₁ and f₂ are aligned with the principal directions of the underlying surface. This completely specifies the local linear model for elliptical shells using Equation (22). As for the spherical shell model, flow lines can be generated in a local neighbourhood by numerical integration. However, these yield a more complex geometry, as illustrated in FIG. 13, for frame integrations on various surfaces: a flat, a spherical, and an ellipsoidal surface.

B4. 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

Accordingly, in this section the inventors 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. Various models have been proposed to explain this organization, based on both physical and mathematical considerations. Also, more recently it has become possible to study this organization within the intact myocardium using diffusion magnetic resonance imaging (dMRI). 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. 14), is ubiquitous. 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° in mammals. 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. 15. 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 and thus the diffusion signal reflects properties of a fibrous composite. Modeling myofiber geometry within the prior art using the GHM as described supra to parametrize myofiber orientation in the heart wall has been performed in three mammalian species, the rat, the canine, and the human. A limitation of using the GHM is that its flow lines are constrained to lie on planar manifolds, in spite of the curvature of the heart wall evident in FIG. 14. This limitation allows only three available geometrical parameters, resulting in a reduced ability for characterizing increasingly larger and complex myofiber neighborhoods. Accordingly, the inventors have instead established a methodology of fitting a frame field to cardiac dMR1 data and to then examine its connection forms.

The inventors describe below this analysis to a diffusion MRI database of healthy rat hearts, containing 8 rat subjects, which will be labeled as subjects A to H, with (0.25 mm)³ voxel resolution and a dimension of 64×64×128 voxels.

B4.1 Selecting a Cardiac Frame Field

Cardiac diffusion MRI volumes are sampled on 3D rectangular lattices with coordinates x=Σx_ie_i, x_iεZ³. In order to apply connection form analysis, a local cardiac frame needs to be defined. The first frame vector of this frame is chosen to be parallel to the orientation of the principal direction of diffusion u₁, i.e. f₁\\u₁. Because the direction of u₁ is locally ambiguous, our first task is to enforce consistency in the direction of f₁ among voxel neighbors. This is done using an adaptive cylindrical coordinate system. The centroid c_z of the chamber within each short-axis slice mask ξ_z is first determined using Equation (45). f₁ is then locally made to turn clockwise with respect to that centroid for each slice ξ_z, and for all xεξ_z as given by Equations (46) to (48) where l_z is the local approximation to the heart's long-axis. For most hearts under consideration, l_z approximately coincides with the volume's z axis.

$\begin{array}{cc}{c}_z=\frac{\sum _{x\in {\xi }_z}x}{{\xi }_z}& \left(45\right)\\ l_z=\frac{c_{z+1}-c_{z-1}}{c_{z+1}-c_{z-1}}& \left(46\right)\\ f_1\left(x\right)=\mathrm{sign}\left(\left(u_1\left(x\right)\times \left(x-c_z\right)·l_z\right)\right)u_1\left(x\right)& \left(47\right)\\ f_1\left(x\right)=v_x\left(u_1\right)u_1& \left(48\right)\end{array}$

We now need one additional orthogonal frame axis to f₁ in order to define a frame field. In fact, since we require f_iεR³ to have unit length, only one additional degree of freedom is required to fully describe a frame field: a rotation about f₁. In principle, any model for this rotation could be adopted, based on geometrical intuition, biological relevance, or computational considerations. However, in common with conventional cardiac spatial analysis, the inventors 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 {circumflex over (f)}₃ using the gradient vector of a distance transform produced as follows:

• • the binary image (mask) of the heart is closed using mathematical morphological operations;
  • the Euclidean distance to the heart wall is evaluated at every point to both the epicardium and endocardium (see FIG. 16), and
  • the endocardium gradient is negated and the average of the two gradients is then computed.

The normals {circumflex over (f)}₃ are then aligned to point from outer to inner wall. With f₁ and {circumflex over (f)}₃, a local frame is specified at x in Equations (49A)-(49C) where f₃ is the part of {circumflex over (f)}₃ orthogonal to f₁. The local plane spanned by f₁ and f₂ will be referred to as the tangent plane. FIG. 15 shows a sample of this frame field within two cross-sectional 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 and is due to the inherent smoothing of the diffusion signal. Various heuristics could be used to identify and mask out these locations, for example using Equation (35).

B4.2 Cardiac Maurer-Cartan Connections

The inventors now employ the Maurer-Cartan connection forms to analyse the cardiac frame field. To get a preview of the biological relevance of this methodology FIG. 17 illustrates the correspondence between three large-scale geometrical features exhibited by cardiac fibers and connection forms, namely, the variation of the helix angle is related to c₁₂₃, the short-axis curvature of the heart wall to c₁₃₁, and the long-axis curvature to C₂₃₂.

B4.2.1 Cardiac Moving Frame Intuition

Considering the frame at x₀, one-form contractions c_ijk 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₁ towards f₂, c₁₂f_k, are intimately linked to the curvature parameters of the GHM. These rotations intuitively describe the manner in which myofibers turn in the tangent plane of the heart, c₁₂₁ in first image 900A in FIG. 9, describes their tangential curvature, c₁₂₂, shown in first image 900B in FIG. 9, their fanning in the tangent plane, and c₁₂₃, shown in third image 900C in FIG. 9, their transmural turning or equivalently the rate of change of the helix angle. The rotations of f₁ towards f₃, c₁₃f_k, express the turning of the fibers towards the inner wall:

• • c₁₃₁ is related to a sectional curvature of the heart wall in the tangential direction, fourth image 900D in FIG. 9;
  • c₁₃₂ describes an upwards twisting of the tangent plane and the rate of change of the transverse angle, fifth image 900E in FIG. 9; and
  • c₁₃₃ describes the degree of transmural fanning of the local fiber population away from or towards the tangent plane, sixth image 900F in FIG. 9.

For the remaining rotations:

• • c₂₃₁ measures twisting of the tangent plane, seventh image 900G in FIG. 9;
  • c₂₃₂ is related to a second sectional curvature of the heart wall, eighth image 900H in FIG. 9; and
  • c₂₃₃ measures myofiber fanning (in or out) in the transmural direction, ninth image 900I in FIG. 9.

B4.2.2 Full Volume Histograms

Full volume histograms for optimized c_ijk, computations in Ω₃ are shown in FIG. 18 for each subject in the dataset. The histograms show that connection forms, in particular the variation of the helix angle c₁₂₃, 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. For example, in hypertrophic hearts the myocardium becomes more spherical globally; using our framework, this should result in sharper histograms for the sectional descriptors c₁₃₁ and c₂₃₂. 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 Ω₃ is shown in FIG. 19. As expected, optimized computations converge near values obtained by the direct computations while lowering the fitting error. Further spatial analysis in Section B4.3 (Table 3) supports these observations.

B4.2.3 Transmural Histograms

Now referring to FIG. 20 there are depicted the results of sampling the 9 connection forms c_ijk 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. 16. 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 c₁₂₃ 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 is not the subject of this specification. 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,

B4.3 Scale and Spatial Dependence

The direct computations of connection forms described in Section B2.1 depend only on the nature of the differential kernel used. On the other hand, the optimized parameter computations of Section B2.2 depend on the size and shape of the neighborhood Ω in which the energy (Equation (32)) is computed. Ω can be adjusted, and a filtering of the principal direction of 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 c_ijk connection forms using different isotropic neighborhoods Ω_i. In doing so, we can measure the smoothness of the frame field, and obtain a preliminary spatial localization of its geometrical variation.

B4.3.1 Neighborhood Shape

Table 2 shows the mean and standard deviation for all voxels of the dataset, for each connection form and various neighborhoods Ω_i. Fitting errors ε_i of Equation (33) 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 θ(∥v∥²) term in Equation (19). This indicates that a first-order 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.

TABLE 2
Effect of Neighbourhood Size on Optimized
Connection Forms and Fitting Errors
Ω	3	5	7
c121	0.008 ± 0.083	0.010 ± 0.113	0.013 ± 0.202
c122	0.015 ± 0.091	0.022 ± 0.113	0.035 ± 0.244
c123	−0.284 ± 0.237	−0.304 ± 0.236	−0.310 ± 0.313
c131	0.039 ± 0.043	0.040 ± 0.030	0.040 ± 0.029
c132	0.016 ± 0.047	0.016 ± 0.032	0.012 ± 0.034
c133	−0.017 ± 0.063	−0.017 ± 0.047	−0.014 ± 0.038
c231	0.000 ± 0.032	0.001 ± 0.022	−0.000 ± 0.018
c232	0.031 ± 0.046	0.028 ± 0.031	0.025 ± 0.024
c233	0.013 ± 0.074	0.010 ± 0.053	0.005 ± 0.042
ε1	0.014 ± 0.030	0.029 ± 0.042	0.047 ± 0.053
ε2	0.013 ± 0.028	0.027 ± 0.040	0.046 ± 0.051
ε3	0.005 ± 0.014	0.009 ± 0.014	0.013 ± 0.015

FIG. 21 explores the relationship between neighborhood size and the various embeddings of Section B3. Extrapolation errors increase with increasing Ω_i, 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 c₁₂₃ 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 Equation (33), a relationship that will be investigated in Section B4.4.

B4.3.2 Filtering Diffusion Directions

Filtering of diffusion directions after adopting cylindrical consistency using Equation (36) can be applied to compensate for the effect of noise in the diffusion volumes and to investigate the scale of fine muscle cardiac structures. 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 and their relationship to biology is not fully established. Thus, in the following description the inventors 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 G_σ+, with standard deviation σ⁺ on the vector field f₁ with a renormalization step, before computing the connection forms. The following update equation is f₁⁰≡f₁, f₁ⁿ⁺¹=G_σ+*f₁ⁿ/∥G_σ+*f₁ⁿ∥. The remaining fields f₂ and f₃ are then obtained using Equation (45). The inventors make use of the notation σ=σ⁺√{square root over (n)} to denote n iterations of standard deviation σ⁺.

As an example, FIG. 23 shows the effect of filtering on f₁ for σ=0.4√{square root over (10)}. FIG. 24 shows a volumetric helical sampling of selected connection forms and fitting errors, computed with respect to the unfiltered data, for various filtering magnitudes σ_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 σ_i is increased values in these regions are propagated to their neighbors and are largely responsible for the overall increase in fitting error. In particular, for c₁₂₃ at σ₄=0.4√{square root over (10)} there is a false suggestion that the voxels turn in a clockwise direction in these mislabeled regions. FIG. 22 and Table 3 summarize these results as means over all volumes for each connection form and fitting error.

TABLE 3
Effect of Iterative Gaussian Smoothing Applied to f1 on Extrapolation
Error for Direct an Optimised Parameter Computations in Ω3
	σ = 0	σ = 0.2{square root over (10)}	σ = 0.4{square root over (10)}
	Direct	Optimised	Direct	Optimised	Direct	Optimised
c121	.006 ± .078	.008 ± .083	.002 ± .076	.006 ± .097	−.001 ± .047	−.001 ± .053
c122	.013 ± .086	.015 ± .091	.008 ± .086	.011 ± .103	.044 ± .046	.003 ± .047
c123	−.235 ± .183	.284 ± .237	−.178 ± .232	−.205 ± .271	−.133 ± .163	−.135 ± .169
c131	.037 ± .054	.039 ± .043	.039 ± .040	.040 ± .036	.044 ± .026	.044 ± .028
c132	.015 ± .055	.016 ± .047	.012 ± .044	.013 ± .040	.010 ± .029	.009 ± .031
c133	−.017 ± .072	.017 ± .063	−.013 ± .067	−.013 ± .056	−.008 ± .045	−.008 ± .044
c231	.000 ± .034	.000 ± .032	.001 ± .032	.001 ± .029	.002 ± .029	.002 ± .031
c232	.028 ± .048	.031 ± .046	.027 ± .047	.029 ± .043	.024 ± .043	.025 ± .041
c233	.013 ± .068	.013 ± .071	.013 ± .062	.015 ± .056	.012 ± .045	.012 ± .047
ε1	.024 ± .053	.014 ± .030	.038 ± .082	.035 ± .087	.041 ± .066	.041 ± .066
ε2	.021 ± .052	.013 ± .028	.035 ± .081	.033 ± .086	.038 ± .064	.038 ± .065
ε3	.007 ± .021	.005 ± .014	.008 ± .023	.006 ± .020	.008 ± .020	.008 ± .021

B4.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 c_ijk only affects the respective ε_i and ε_j errors. FIG. 25 demonstrates that the most important connection form in lowering both ε₁ and ε₂ is c₁₂₃, since it describes a rotation of the helix angle, which has a large magnitude throughout the volume. c₁₂₁ and c₁₂₂ also significantly lower the error, which supports the use of low dimensional embeddings such as the GHM described in Section B3.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. ε₃ 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 predicted directly from FIG. 25.

Embeddings that offer greater complexity in capturing the variation of a particular frame axis will lower the frame error associated with it. FIG. 26 shows the three mean errors ε_i for various embeddings as a function of neighborhood size Ω_i, using optimized computations. In relation to the frame error ε₁, all embeddings with an unconstrained c₁₂₃ will perform similarly, as shown in FIG. 26, since it is the connection form that dominates most of the frame field variation. On the other hand, the third frame vector f₃ captures the extent to which the principal direction of diffusion remains normal to the heart wall normal. Since fiber directions u₁ generally run tangentially to the heart wall, the resulting error is low for all embeddings.

Accordingly, the inventors have shown that measurements of the geometry of cardiac myofibers can be performed 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) Helix Angle: c₁₂₃ 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 7 voxels. Integrating c₁₂₃ throughout this distance produces a total change of −116.3°, which is in close correspondence with values of 120° reported in the literature.

2) Wall Curvatures: c₁₃₁ and c₁₃₂ reflect the sectional curvatures of the heart wall as projected onto the local osculating ellipsoid to f₁. 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) Myocyte Fanning: c₁₂₂ and c₁₃₃ 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.

4) Myocyte Twisting: c₁₂₃, c₁₃₂ and c₂₃₁ are measures of twisting in the collagen matrix that contains cardiac myocytes. Whereas c₁₂₃ corresponds to the variational component of the helix angle, c₁₃₂ describes a turning that is directed in an upward fashion from base to apex, and c₂₃₁ one that rotates the tangent fiber plane.

Accordingly, the inventors have shown that the method of moving frames can be applied to dMRI data and employed to explore and compute the variation of a smooth frame field. The method allows the development of a selection of geometrical embeddings which imposed certain constraints—via thresholding (spherical and ellipsoidal manifolds) or by assuming a functional dependency (e.g. GHM) on connection forms. Although the embodiments of the invention have been presented with respect to simple cases, it would be evident that 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 to characterize the geometry of cardiac myocytes. 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.

C. Tissue Engineering

Heart wall myofibers are arranged in minimal surfaces to optimize organ function. Based upon the analysis performed by the inventors they have established that the orientation of myofibers within the heart are locally arranged in a very special manner which can described by a class of minimal surfaces called generalized helicoids. By describing these surfaces locally with a small number of parameters the inventors have been able to generate mathematical fits to myofiber orientation data measured using diffusion magnetic resonance imaging of rat, dog, and human hearts. The computer model shows how fibers should be oriented locally to give the heart wall its mechanical function, and thus could be used in tissue engineering applications which require the regeneration of heart wall tissue where it has been damaged, as in the case of infarctions. The model can also be used to provide atlases of normal fiber geometry to be used in clinical applications. Previous studies have described the shape of individual fibers as pieces of helical curves, but not their collective volumetric structure in the heart wall.

It would be evident therefore that the computer models of the human heart, for example, can thereby serve as a scaffold for artificial muscle construction. For example, the scaffold may be formed from a polymer using the mathematical model according to embodiments of the invention. Accordingly, implementing a material exploiting polymer fibers aligned and orientated as with a human heart muscle (wall) allows for the construction of a new kind of composite material with flexible, expandable and contractible properties, which could have any number of uses. Further, it would be evident that the models defined according to embodiments of the invention have benefits within heart tissue engineering and the diagnosis of heart muscle diseases.

Whilst the embodiments of the invention described supra have been described with respect to a heart it would be evident that other biological organs and/or elements may be similarly modeled and analysed exploiting data such as diffuse magnetic resonance imaging, for example. Accordingly, the mathematical model may provide the surfaces and orientation of fibers forming the biological structure.

The foregoing disclosure of the exemplary embodiments of the present invention has been presented for purposes of illustration and description. It is not intended to be exhaustive or to limit the invention to the precise forms disclosed. Many variations and modifications of the embodiments described herein will be apparent to one of ordinary skill in the art in light of the above disclosure. The scope of the invention is to be defined only by the claims appended hereto, and by their equivalents.

Further, in describing representative embodiments of the present invention, the specification may have presented the method and/or process of the present invention as a particular sequence of steps. However, to the extent that the method or process does not rely on the particular order of steps set forth herein, the method or process should not be limited to the particular sequence of steps described. As one of ordinary skill in the art would appreciate, other sequences of steps may be possible. Therefore, the particular order of the steps set forth in the specification should not be construed as limitations on the claims. In addition, the claims directed to the method and/or process of the present invention should not be limited to the performance of their steps in the order written, and one skilled in the art can readily appreciate that the sequences may be varied and still remain within the spirit and scope of the present invention.

Claims

1. A method comprising: acquiring diffuse magnetic resonance imaging data relating to a heart;establishing a model relating to the heart based on a generalized helicoid model that expresses a frame field of local fiber directions in a plane tangent to a heart wall;fitting each data point of the data as a local orthogonal frame expressed as (i, j, k) to represent the local fiber directions from the model; anddetermining connection forms cijk at each data point to represent a rotation of the local orthogonal frame in a spatial neighborhood of the data point in accordance with the data.

2. The method of claim 1, wherein determining connection forms cijk comprises computing the connection forms cijk as a minimizer of an energy contained within the spatial neighborhood.

3. The method of claim 2, wherein the energy is solved for using Nelder-Mead iterations.

4. The method of claim 2, wherein determining the connection forms cijk comprises enforcing bounds on variables used for minimizing the energy.

5. The method of claim 1, wherein at least one of the local fiber directions for the model corresponds to a local normal to the heart wall which is an approximate direction in which an endocardium moves when the heart beats.

6. The method of claim 1, wherein the connection forms cijk are Maurer-Cartan connection forms.

7. The method of claim 6, wherein (i, j, k) vary from 1 to 3, and wherein c123 is related to a rate of change of a helix angle of cardiac fibers within the heart; c131 and c132 are related to sectional curvatures of a heart wall of the heart; c122 and c133 are related to a rate of fanning out for cardiac myocytes within the heart; and c123, c132 and c231 are related to measures of twisting in a collagen matrix that contains cardiac myocytes within the heart.

8. The method of claim 1, wherein establishing a model relating to the heart comprises imposing constraints on the local fiber directions to control a shape and complexity of the frame field.

9. The method of claim 1, wherein determining connection forms cijk at each data point comprises limiting a size of the spatial neighborhood to minimize a fitting error.

10. A non-transitory computer readable medium having stored thereon program instructions executable by a processor for: acquiring diffuse magnetic resonance imaging data relating to a heart;establishing a model relating to the heart based on a generalized helicoid model that expresses a frame field of local fiber directions in a plane tangent to a heart wall;fitting each data point of the data as a local orthogonal frame expressed as (i, j, k) to represent the local fiber directions from the model; anddetermining connection forms cijk at each data point to represent a rotation of the local orthogonal frame in a spatial neighborhood of the data point in accordance with the data.

11. The computer readable medium of claim 10, wherein determining connection forms cijk comprises computing the connection forms cijk as a minimizer of an energy contained within the spatial neighborhood.

12. The computer readable medium of claim 11, wherein the energy is solved for using Nelder-Mead iterations.

13. The computer readable medium of claim 11, wherein determining the connection forms cijk comprises enforcing bounds on variables used for minimizing the energy.

14. The computer readable medium of claim 10, wherein at least one of the local fiber directions for the model corresponds to a local normal to the heart wall which is an approximate direction in which an endocardium moves when the heart beats.

15. The computer readable medium of claim 10, wherein the connection forms cijk are Maurer-Cartan connection forms.

16. The computer readable medium of claim 15, wherein (i, j, k) vary from 1 to 3, and wherein c123 is related to a rate of change of a helix angle of cardiac fibers within the heart; c131 and c132 are related to sectional curvatures of a heart wall of the heart; c122 and c133 are related to a rate of fanning out for cardiac myocytes within the heart; and c123, c132 and c231 are related to measures of twisting in a collagen matrix that contains cardiac myocytes within the heart.

17. The computer readable medium of claim 10, wherein establishing a model relating to the heart comprises imposing constraints on the local fiber directions to control a shape and complexity of the frame field.

18. The computer readable medium of claim 10, wherein determining connection forms cijk at each data point comprises limiting a size of the spatial neighborhood to minimize a fitting error.