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.

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; and determining connection forms c.sub.ijk 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 c.sub.ijk comprises computing the connection forms c.sub.ijk 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 c.sub.ijk 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 c.sub.ijk are Maurer-Cartan connection forms.

7. The method of claim 6, wherein (i, j, k) vary from 1 to 3, and wherein c.sub.123 is related to a rate of change of a helix angle of cardiac fibers within the heart; c.sub.131 and c.sub.132 are related to sectional curvatures of a heart wall of the heart; c.sub.122 and c.sub.133 are related to a rate of fanning out for cardiac myocytes within the heart; and c.sub.123, c.sub.132 and c.sub.231 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 c.sub.ijk 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; and determining connection forms c.sub.ijk 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 c.sub.ijk comprises computing the connection forms c.sub.ijk 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 c.sub.ijk 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 c.sub.ijk are Maurer-Cartan connection forms.

16. The computer readable medium of claim 15, wherein (i, j, k) vary from 1 to 3, and wherein c.sub.123 is related to a rate of change of a helix angle of cardiac fibers within the heart; c.sub.131 and c.sub.132 are related to sectional curvatures of a heart wall of the heart; c.sub.122 and c.sub.133 are related to a rate of fanning out for cardiac myocytes within the heart; and c.sub.123, c.sub.132 and c.sub.231 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 c.sub.ijk at each data point comprises limiting a size of the spatial neighborhood to minimize a fitting error.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

(2) 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;

(3) FIG. 2 depicts histogram of c.sub.ijkc.sub.ijF.sub.k for rat, human, and with measures normalized such that the horizontal axis is in millimeters;

(4) FIG. 3 depicts the one-form extrapolation error for neighborhoods of size |N.sub.i=3.sup.3, 5.sup.3, 7.sup.3;

(5) FIG. 4 depicts the one-form approximation to the GHM for the neighborhood N.sub.3 together with the error of each model a function of N.sub.i;

(6) 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.sup.3, 5.sup.3, 7.sup.3;

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

(8) FIG. 7 depicts the rotational angle of the different models analyzed for the human heart in isotropic voxel neighbourhoods of size 3.sup.3, 5.sup.3, 7.sup.3;

(9) FIG. 8 depicts joint histograms of the angular errors e and for the human heart for different models for isotropic 3.sup.3 neighborhoods;

(10) FIG. 9 depicts the geometry characterized by connection forms;

(11) FIG. 10 depicts parameters after 200 Nelder-Mead iterations including convergence plots for c.sub.123 and .sub.1 and computational timings for different estimation techniques;

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

(13) 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.sub.1, K.sub.2) tuples;

(14) 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.sub.131, c.sub.132, c.sub.231, c.sub.232;

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

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

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

(18) 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;

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

(20) FIG. 19 depicts selected dataset volumes of rat connection forms c.sub.ijk (radians/voxel) and fitting errors .sub.i (radians) in .sub.3 for the optimized and direct computations;

(21) FIG. 20 depicts transmural histograms of combined volumes of c.sub.ijkc.sub.ijf.sub.k and fitting errors in .sub.3 for the rat data sets using optimized parameter computations;

(22) FIG. 21 depicts normalized distributions of the extrapolation error .sub.i (in radians) for different neighborhoods .sub.i;

(23) FIG. 22 depicts a method comparison and effect of f.sub.1 filtering on mean c.sub.ijk and .sub.i in the .sub.1 for different standard deviations of the vector field f.sub.1;

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

(25) FIG. 24 depicts the effect of f.sub.1 filtering on selected connection forms (rads/voxel) and extrapolation error .sub.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;

(26) 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

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

DETAILED DESCRIPTION

(28) 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.

(29) 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.

(30) A1. The Maurer-Cartan Form

(31) 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.sub.1, F.sub.2, F.sub.3.sup.3, F.sub.i.Math.F.sub.j=.sub.ij, where .sub.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.sub.1e.sub.2e.sub.3] in Equation (1) where the attitude matrix ASO(3) is a smoothly varying orthonormal matrix, and where the basis vectors e.sub.j are treated as symbols such that F.sub.i=a.sub.ije.sub.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.sup.1=A.sup.T, C=(dA)A.sup.1 is the Maurer-Cartan form, and where for simplicity the notation dF.sub.i=.sub.jc.sub.ijF.sub.j is used.

(32) $\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)\\ \left[\begin{array}{c}{F}_1\\ F_2\\ F_3\end{array}\right]=\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\end{array}\right]=\left(A\right)\left[\begin{array}{c}{e}_x\\ e_y\\ e_z\end{array}\right]=\left(A\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}$

(33) The Maurer-Cartan matrix is skew symmetric, i.e. C=C.sup.T. Hence it has at most 3 independent, non-zero elements: c.sub.12, c.sub.13, and c.sub.23. Each c.sub.ij is a one-form in .sup.3 that can be contracted on a vector v=[v.sub.1, v.sub.2, v.sub.3].sup.T.sup.3 to yield the initial rate of turn of F.sub.i towards F.sub.j when moving in the direction of v. The inventors denote this contraction c.sub.ij(v), which is found to be c.sub.ij(v)=.sub.vF.sub.i.Math.F.sub.j|.sub.x, where x.sup.3 is a point in the fiber field and .sub.vF.sub.i is the covariant derivative of F.sub.i in the direction v. Accordingly, c.sub.ij(v) is defined by Equation (3) where the components of the frame vectors are enumerated as F.sub.i=[F.sub.i1, F.sub.i2, F.sub.i3].sup.T where .sub.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.sub.ijkc.sub.ijF.sub.k. It should be noted that the frame field F.sub.1, F.sub.2, F.sub.3 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.sub.ijk.

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

(35) 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 oneform parameterization.

(36) A2. Measures on a Discrete Fiber Frame Field

(37) The inventors analyze hearts represented as diffusion MRI volumes embedded in 3D rectangular lattices with coordinates x=xe.sub.1+ye.sub.2+ze.sub.3=[x, y, z].sup.TZ.sup.3. 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.sub.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.sub.z) is the local approximation of the heart's long-axis.
T(x).fwdarw.sign((T(xc.sub.z)).Math.l(s.sub.z))T(x)(4)

(38) For all the hearts that we consider, l(s.sub.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.

(39) 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.sub.j Further, the inventors also refer to the local plane spanned by T and N as the tangent plane.

(40) $\begin{array}{cc}{F}_1=\frac{T}{.Math.T.Math.}& \left(5A\right)\\ F_2=N=F_3{F}_1& \left(5B\right)\\ F_3=B=\frac{\left(\hat{B}-\left(\hat{B}.Math.T\right)T\right)}{.Math.\hat{B}-\left(\hat{B}.Math.T\right)T.Math.}& \left(5C\right)\end{array}$

(41) A2.1 One-Form Intuition

(42) One-form contractions c.sub.ijk can be interpreted as the amount of turning of F.sub.i towards F.sub.j, in the direction F.sub.k. For example, c.sub.TNB describes a transmural rotation of T towards N, as shown in second image 100B in FIG. 1. c.sub.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.sub.TNBF.sub.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.sub.TNT in first image 200A in FIG. 2, describes their tangential curvature, c.sub.TNN in first image 200B in FIG. 2, their fanning in the tangent plane, and c.sub.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.sub.TBF.sub.k, expresses the turning of the fibers towards the inner wall, c.sub.TBT effectively measures the first the local curvature of the heart, c.sub.TBN describes a twisting of the tangent plane and the rate of change of the transverse angle, and c.sub.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.sub.NBF.sub.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.sub.NBT measures a twisting of the tangent plane, c.sub.NBN the second curvature of the heart wall and c.sub.NBB a transmural fanning or thickening.

(43) A2.2 The One-Form Model

(44) 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)}.sub.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.sub.i|=i.sup.3 for odd iZ and T.sub.h is the true neighbours measured direction.

(45) $\begin{array}{c}{\tilde{T}}_h=T+c_{\mathrm{TN}}.Math.h.Math.N+c_{\mathrm{TB}}.Math.h.Math.B\left(6\right)\\ =T+\left(c_{\mathrm{TNT}}h.Math.T+c_{\mathrm{TNN}}h.Math.N+c_{\mathrm{TNB}}h.Math.B\right)N+\\ \left(c_{\mathrm{TBT}}h.Math.T+c_{\mathrm{TBN}}h.Math.N+c_{\mathrm{TBB}}h.Math.B\right)B\left(7\right)\end{array}$$e\left(N_i\right)=\frac{1}{.Math.N_i.Math.}\underset{h{N}_i}{.Math.}\arccos \left(T_h.Math.\frac{{\tilde{T}}_h}{.Math.{\tilde{T}}_h.Math.}\right)\left(8\right)$

(46) 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.

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

(48) 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.sub.TT+x.sub.NN+x.sub.BB.sup.3 is given as the angle defined in Equation (9) where K.sub.*R are the GHM curvature parameters.

(49) $\begin{array}{cc}\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}$

(50) Direct calculations show that a frame field spanned by T(), N(), B() has instantaneous turning given by c.sub.TNT()=K.sub.T, c.sub.TNN()=K.sub.N, and c.sub.TNB()=K.sub.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.sub.T, K.sub.N, K.sub.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)}.sub.h.fwdarw.{tilde over (T)}.sub.h(K)=(cos , sin , 0) and =(K) as given by Equation (9.

(51) 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.

(52) 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.

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

(54) 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.

(55) As introduced in Section 1, the Maurer-Cartan form has only 3 independent one-forms: c.sub.TN, c.sub.TB and c.sub.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.sub.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.sub.TNT, c.sub.TBT and c.sub.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.sub.TNT=;c.sub.TNN=; c.sub.TNB=; c.sub.TBB0; c.sub.NBB0(10)

(56) 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.sub.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 .sub.i and .sub.2 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.sub.TBN and c.sub.TBB such that we obtain Equation (12).

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

(58) A3. Model Space Comparison

(59) 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)}.sub.h=T in Equation (7). To compare the different models in terms of their fitting accuracy, we have evaluated each c.sub.ijk on the human data set using first-order central differences on 3.sup.3 neighbors. We then used these one-forms to extrapolate each model using Equations (7) and (8) in isotropic neighbourhoods where N.sub.i, where |N.sub.i|=i.sup.3 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.sup..sup.2 and mean

(60) ${}^{+\frac{1}{2}{}^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.

(61) In Section 2.4 the inventors showed that their extrapolated {tilde over (T)}.sub.h axis only differs by the c.sub.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.

(62) TABLE-US-00001 TABLE 1 Extrapolation Errors in Radians for Each Species, Differential Model and Neighbourhood N.sub.i $\mathrm{fr}\mathrm{om}\mathrm{the}\mathrm{mode}{e}^{-{}^2}\mathrm{and}\mathrm{the}\mathrm{mean}{e}^{+\frac{1}{2}{}^2}\mathrm{of}\log \text{-}\mathrm{normal}\mathrm{fits}.$ Rat Dog Human |N.sub.i| (mode, mean) (mode, mean) (mode, mean) 3.sup.3 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 5.sup.3 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 7.sup.3 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

(63) A moving frame in .sup.3 has 3 degrees of freedom of which 2 are captured by the error vector T{tilde over (T)}.sub.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)}.sub.h, onto the local NB plane and measuring its angle with respect to the frame axis N:

(64) $=\arctan \frac{\left(T-{\tilde{T}}_h\right).Math.B}{\left(T-{\tilde{T}}_h\right).Math.N},i.e.$
=(0, ) and

(65) 0$\left(\frac{}{2},\frac{3}{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.sup.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.

(66) B1. Structure of the Frame Fields

(67) Let a point x=.sub.ix.sub.ie.sub.i.sup.3 be expressed in terms of the natural orthonormal coordinate system e.sub.1, e.sub.2, e.sub.3. A differential orthonormal frame field embedded in .sup.3 is denoted as F=[f.sub.1, f.sub.2, f.sub.3].sup.T: .sup.3.fwdarw..sup.3 and defined by f.sub.1.Math.f.sub.2=.sub.ij, where .sub.ij is the Kronecker delta and where .Math. is the inner product, and f.sub.1f.sub.2=f.sub.3, where x is the 3-dimensional cross product. Since the frame E=[e.sub.1, e.sub.2, e.sub.3].sup.T forms an orthonormal basis for .sup.3, F can be expressed as the rotation f.sub.i=.sub.ja.sub.ije.sub.j, where the elements of the attitude matrix A A={a.sub.i,j}.sup.33 are differentiable, and A.sup.1=A.sup.T. Treating f.sub.i and e.sub.j as symbols then this rotation can be written in matrix form as given by Equation (13).

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

(69) Since each e.sub.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.sup.1={c.sub.ij}.sup.33 is the Maurer-Cartan matrix of one-forms c.sub.ij. Now writing f.sub.i as symbols then Equation (14) can be understood as Equation (15).

(70) The Maurer-Cartan matrix is skew symmetric and accordingly Equation (16) applies such that there are at most 3 independent, non-zero elements: c.sub.12, c.sub.13 and c.sub.23. The differential of the elements a.sub.ij(x): .sup.3.fwdarw. of A is expressed in terms of dx.sup.k, the natural basis for one-forms in R.sup.3, and is given by Equation (17).

(71) $\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}}=\underset{k}{.Math.}\frac{a_{\mathrm{ij}}}{x_k}{x}^k& \left(17\right)\end{array}$

(72) 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 .sup.3 and tangent vector v.sup.3. The contraction of dw=.sub.kw.sub.kdx.sup.k on v is written as a bilinear operation on the canonical projection of v: dw(v)=.sub.iw.sub.ide.sub.i(.sub.jv.sub.je.sub.j)=.sub.iw.sub.iv.sub.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.sub.j(x) when x moves in the direction v. This analysis naturally applies to frame fields embedded in .sup.2 by setting f.sub.3 to be constant. In that case, C is a 22 matrix with a single one-form c.sub.12 which describes the in-plane turning of the frame; hence, if f.sub.1 is the tangent of a planar curve and f.sub.2 its normal, then c.sub.12(f.sub.1) is the curves curvature and c.sub.12(f.sub.2) its fanning.

(73) $\begin{array}{cc}{c}_{\mathrm{ij}}.Math.v.Math.=\underset{k}{.Math.}{a}_{\mathrm{jk}}\underset{l}{.Math.}\frac{a_{\mathrm{ik}}}{x_l}{v}_l& \left(18\right)\\ {{f_i.Math.v.Math..Math.}_{x_0}=f_1.Math.}_{x_0}+d{f}_i.Math.v.Math.+O\left({.Math.v.Math.}^2\right)& \left(19\right)\\ =f_1{|}_{x_0}+\underset{j1}{.Math.}{f}_j.Math.v.Math.\underset{k}{.Math.}{v}_k{c}_{\mathrm{ijk}}{|}_{x_0}+O\left({.Math.v.Math.}^2\right)& \left(20\right)\end{array}$

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

(75) Considering the model at a point x.sub.0 for all frame vectors and one-forms, an approximation of the frame from F(x.sub.0) to F(x.sub.0+v) can be obtained by the normalized model given in Equation (21) which in matrix notation yields Equation (22) where C.sup.iR.sup.33 denotes the matrix of connection forms with c.sub.kl.sup.i=.sub.jkc.sub.kil and .sub.iksgn(ki). In general {tilde over (f)}.sub.i(v), {tilde over (f)}.sub.2(v), and {tilde over (f)}.sub.3(v) will not form an orthogonal basis. To establish we can combine Equation (15) with Equation (19) to yield Equation (23) for ki, kj.

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

(77) B2. Estimation of Connection Forms

(78) 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.sub.ijk which fully characterize the local geometry of a frame field. Accordingly, we need to compute these parameters for a smooth frame field.

(79) B2.1 Direct Computation of Connection Forms

(80) The connection one-forms c.sub.ij can be directly obtained from the vector fields f.sub.i and their differentials df.sub.i using Equation (14) to yield Equation (24) and therein Equations (25) and (26). The differential df.sub.i can be computed by applying the exterior derivative of a function g: R.sup.n.fwdarw.R as given in Equation (27) which yields Equations (28) to (30) where

(81) $=\left[\frac{x^i}{x_j}\right]R^{33}J$
is the Jacobian matrix of partial derivatives and ds=(dx.sup.1, dx.sup.2, dx.sup.3).sup.T. Setting v=f.sub.k then we finally obtain Equation (31). The Jacobian matrix J can be approximated to first order using for example finite differences such that

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

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

(84) B2.2 Computation Via Energy Minimization

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

(86) $\begin{array}{cc}{c}_{\mathrm{ijk}}^{*}\left(x_0\right)=\arg \underset{c_{\mathrm{ijk}}}{\min }\frac{1}{.Math..Math.}\underset{v}{.Math.}\stackrel{3}{\underset{i}{.Math.}}{\mathrm{.Math.}}_i\left(x_0+v\right)& \left(32\right)\\ {\mathrm{.Math.}}_i\left(x_0+v\right)=\arccos \left(f_i\left(x_0+v\right).Math.{\tilde{f}}_i\left(x_0+v\right)\right)& \left(33\right)\end{array}$

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

(88) B2.3 Method Comparison

(89) 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 .sub.3.sup.+ and then evaluated in .sub.3. 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.sub.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.sub.i=d(f.sub.i1e.sub.1+f.sub.i2e.sub.2+f.sub.i3e.sub.3) in f.sub.i is bounded due to Equation (34).

(90) $\begin{array}{cc}{.Math.\frac{f_i}{x_k}.Math.}_x.Math.\frac{f_{\mathrm{ij}}\left(x-e_k\right)-f_{\mathrm{ij}}(x+e_{k)}}{2}.Math.1& \left(34\right)\end{array}$

(91) Using .Math..sub.1 to denote the Euclidean one-norm, then the bounds defined by Equations (35A) to (35C) are obtained.
|c.sub.i,jv|=|f.sub.jTf.sub.i(x.sub.1,x.sub.2, x.sub.3)v|(35A)
f.sub.jT.Math.[v.sub.1v.sub.1v.sub.1].sup.T(35B)
.sub.1 for v.sub.2+1(35C)

(92) 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.

(93) 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.

(94) B3. Maurer-Cartan Embeddings

(95) The geometrical intuition conveyed by the c.sub.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.sub.i(x.sub.0) with known parameters c.sub.ijk, the solution to the manifold f.sub.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.

(96) 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.

(97) B3.1. Generalized Helicoids as Connection Forms

(98) 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.

(99) 0$\begin{array}{cc}{f}_1\left(x_0+v\right)=\cos f_1\left(x_0\right)+\sin f_2\left(x_0\right)& \left(36A\right)\\ \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(\right)=\left[\begin{array}{ccc}\cos & \sin & 0\\ -\sin & \cos & 0\\ 0& 0& 1\end{array}\right]& \left(37\right)\end{array}$

(100) f.sub.1 is given by the GHM vector field in Equations (36A) and (36B) where (x.sub.0, v) parameterizes a minimal surface over parameter fields K.sub.i and where v=.sub.iv.sub.if.sub.i(x.sub.0). The 9 parameters c.sub.ijk of Equation (22) can then be evaluated for f.sub.1. The trivial basis e=f.sub.1(x.sub.0) is chosen in order to construct the attitude matrix A, Equation (37), where =(x.sub.0,v). Since f.sub.1 stays parallel to the tangent f.sub.1f.sub.2 plane at x.sub.0 and f.sub.3 remains unchanged, then {tilde over (f)}.sub.2=f.sub.3f.sub.1 will also stay parallel to that tangent plane at x.sub.0, Hence, we must have c.sub.13=c.sub.23=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.sub.i and df.sub.i are evaluated at x.sub.0. Accordingly we arrive at Equation (39) and conclude that the GHM is a subset of Equation (22), where all but c.sub.121, c.sub.122, and c.sub.123 are zero. As a result, in frame fields where the 6 connection forms c.sub.13k, c.sub.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.sub.0, v) yields a frame field singularity when Equation (30) is satisfied.

(101) $\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}.Math.f_k.Math.=K_k,k\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}$

(102) B3.2 A Model on Spheres

(103) Connection forms can be constrained to generate thin spherical shell embeddings. Given a smooth surface where f.sub.1 and f.sub.2 span the local tangent plane, then we can generate a family of subspaces that are non-intersecting along f.sub.3, by setting c.sub.12f.sub.k.sub.k, k{1,2,3} as given in Equations (41A) and (41B) where K.sub.i: .sup.3.fwdarw. 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

(104) $x_0-\frac{1}{}{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.sub.1 in a local neighborhood, as illustrated in FIG. 11.
c.sub.12f.sub.1=K.sub.1c.sub.12f.sub.2=K.sub.2c.sub.12f.sub.3=K.sub.3(41A)
c.sub.13f.sub.3=0c.sub.23f.sub.3=0(41B)

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

(106) Curves with f.sub.1 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.sub.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.sub.0. Since z.sub.0= cos(), using Equation (43), the curve of intersection is found to be

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

(108) $f_1=\frac{{}_{}p}{.Math.{}_{}p.Math.},f_3=\frac{{\left(x,y,z\right)}^T}{.Math.{\left(x,y,z\right)}^T.Math.},$
and f.sub.1=f.sub.3f.sub.1, we find that c.sub.12f.sub.1=z.sub.0/{square root over (.sup.2z.sub.0.sup.2)}, c.sub.13f.sub.1=1/ and c.sub.23f.sub.1=0 which exactly corresponds to the constraints given by Equations (41A)-(41B) and (42A)-(42B) when K.sub.1=z.sub.0/{square root over (.sup.2z.sub.0.sup.2)}. To align the two coordinate systems, we first set the tangent of the circle to be f.sub.i, and set the change of f.sub.i in the direction f.sub.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.sub.1f.sub.1=K.sub.1f.sub.2(1/)f.sub.3. Finally, the radius of the circle is found by solving for z.sub.0 in terms K.sub.1, z.sub.0=K.sub.1.sup.2/{square root over (1+K.sub.1.sup.2.sup.2)} and inserting this into R which yields R=/{square root over (1+K.sub.1.sup.2.sup.2)}=1/|df.sub.1f.sub.1.

(109) B3.3 A Model on Ellipsoids

(110) 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.sub.1 and f.sub.2 span the tangent plane which is organized in thin shells, then Equations (41)-(41B) and (42A)-(42B) generalize to Equations (44A) and (44B).

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

(112) Here, c.sub.13f.sub.3 and c.sub.23f.sub.3 are only zero when f.sub.1 and f.sub.2 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.

(113) B4. Application to Myofiber Geometry

(114) 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

(115) 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.

(116) 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.

(117) 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.

(118) 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).sup.3 voxel resolution and a dimension of 6464128 voxels.

(119) B4.1 Selecting a Cardiac Frame Field

(120) Cardiac diffusion MRI volumes are sampled on 3D rectangular lattices with coordinates x=x.sub.ie.sub.i, x.sub.iZ.sup.3. 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.sub.1, i.e. f.sub.1\\u.sub.1. Because the direction of u.sub.1 is locally ambiguous, our first task is to enforce consistency in the direction of f.sub.1 among voxel neighbors. This is done using an adaptive cylindrical coordinate system. The centroid c.sub.z of the chamber within each short-axis slice mask .sub.z is first determined using Equation (45). f.sub.1 is then locally made to turn clockwise with respect to that centroid for each slice .sub.z, and for all x.sub.z as given by Equations (46) to (48) where l.sub.z is the local approximation to the heart's long-axis. For most hearts under consideration, l.sub.z approximately coincides with the volume's z axis.

(121) $\begin{array}{cc}{c}_z=\frac{\underset{x{}_z}{.Math.}x}{.Math.{}_z.Math.}& \left(45\right)\\ l_z=\frac{c_{z+1}-c_{z-1}}{.Math.c_{z+1}-c_{z-1}.Math.}& \left(46\right)\\ f_1\left(x\right)=\mathrm{sign}\left(\left(u_1\left(x\right)\left(x-c_z\right).Math.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}$

(122) We now need one additional orthogonal frame axis to f.sub.1 in order to define a frame field. In fact, since we require f.sub.iR.sup.3 to have unit length, only one additional degree of freedom is required to fully describe a frame field: a rotation about f.sub.1. 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.

(123) Specifically, we first estimate a transmural direction {circumflex over (f)}.sub.3 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.

(124) The normals {circumflex over (f)}.sub.3 are then aligned to point from outer to inner wall. With f.sub.1 and {circumflex over (f)}.sub.3, a local frame is specified at x in Equations (49A)-(49C) where f.sub.3 is the part of {circumflex over (f)}.sub.3 orthogonal to f.sub.1. The local plane spanned by f.sub.1 and f.sub.2 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).

(125) B4.2 Cardiac Maurer-Cartan Connections

(126) 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.sub.123, the short-axis curvature of the heart wall to c.sub.131, and the long-axis curvature to C.sub.232.

(127) B4.2.1 Cardiac Moving Frame Intuition

(128) Considering the frame at x.sub.0, one-form contractions c.sub.ijk can be interpreted as the amount of turning of f.sub.i towards f.sub.j, when considering neighboring frames in the direction f.sub.k. The rotations of f.sub.1 towards f.sub.2, c.sub.12f.sub.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.sub.121 in first image 900A in FIG. 9, describes their tangential curvature, c.sub.122, shown in first image 900B in FIG. 9, their fanning in the tangent plane, and c.sub.123, 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.sub.1 towards f.sub.3, c.sub.13f.sub.k, express the turning of the fibers towards the inner wall: c.sub.131 is related to a sectional curvature of the heart wall in the tangential direction, fourth image 900D in FIG. 9; c.sub.132 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.sub.133 describes the degree of transmural fanning of the local fiber population away from or towards the tangent plane, sixth image 900F in FIG. 9.

(129) For the remaining rotations: c.sub.231 measures twisting of the tangent plane, seventh image 900G in FIG. 9; c.sub.232 is related to a second sectional curvature of the heart wall, eighth image 900H in FIG. 9; and c.sub.233 measures myofiber fanning (in or out) in the transmural direction, ninth image 900I in FIG. 9.

(130) B4.2.2 Full Volume Histograms

(131) Full volume histograms for optimized c.sub.ijk, computations in .sub.3 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.sub.123, 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.sub.131 and c.sub.232. Other cardiomyopathies that are more localized could be investigated from the local geometrical signature of the cardiac frame field.

(132) A comparison of selected volume histograms for the direct and optimized parameter estimation methods in .sub.3 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.

(133) B4.2.3 Transmural Histograms

(134) Now referring to FIG. 20 there are depicted the results of sampling the 9 connection forms c.sub.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.sub.123 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,

(135) B4.3 Scale and Spatial Dependence

(136) 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.sub.ijk connection forms using different isotropic neighborhoods .sub.i. In doing so, we can measure the smoothness of the frame field, and obtain a preliminary spatial localization of its geometrical variation.

(137) B4.3.1 Neighborhood Shape

(138) Table 2 shows the mean and standard deviation for all voxels of the dataset, for each connection form and various neighborhoods .sub.i. Fitting errors .sub.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.sup.2) 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.

(139) TABLE-US-00002 TABLE 2 Effect of Neighbourhood Size on Optimized Connection Forms and Fitting Errors 3 5 7 c.sub.121 0.008 0.083 0.010 0.113 0.013 0.202 c.sub.122 0.015 0.091 0.022 0.113 0.035 0.244 c.sub.123 0.284 0.237 0.304 0.236 0.310 0.313 c.sub.131 0.039 0.043 0.040 0.030 0.040 0.029 c.sub.132 0.016 0.047 0.016 0.032 0.012 0.034 c.sub.133 0.017 0.063 0.017 0.047 0.014 0.038 c.sub.231 0.000 0.032 0.001 0.022 0.000 0.018 c.sub.232 0.031 0.046 0.028 0.031 0.025 0.024 c.sub.233 0.013 0.074 0.010 0.053 0.005 0.042 .sub.1 0.014 0.030 0.029 0.042 0.047 0.053 .sub.2 0.013 0.028 0.027 0.040 0.046 0.051 .sub.3 0.005 0.014 0.009 0.014 0.013 0.015

(140) FIG. 21 explores the relationship between neighborhood size and the various embeddings of Section B3. Extrapolation errors increase with increasing .sub.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.sub.123 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.

(141) B4.3.2 Filtering Diffusion Directions

(142) 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.sub..sub.+, with standard deviation .sup.+ on the vector field f.sub.1 with a renormalization step, before computing the connection forms. The following update equation is f.sub.1.sup.0f.sub.1, f.sub.1.sup.n+1=G.sub..sub.+*f.sub.1.sup.n/G.sub..sub.+*f.sub.1.sup.n. The remaining fields f.sub.2 and f.sub.3 are then obtained using Equation (45). The inventors make use of the notation =.sup.+{square root over (n)} to denote n iterations of standard deviation .sup.+.

(143) As an example, FIG. 23 shows the effect of filtering on f.sub.1 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 .sub.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 .sub.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.sub.123 at .sub.4=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.

(144) TABLE-US-00003 TABLE 3 Effect of Iterative Gaussian Smoothing Applied to f.sub.1 on Extrapolation Error for Direct an Optimised Parameter Computations in .sub.3 = 0 = 0.2{square root over (10)} = 0.4{square root over (10)} Direct Optimised Direct Optimised Direct Optimised c.sub.121 .006 .078 .008 .083 .002 .076 .006 .097 .001 .047 .001 .053 c.sub.122 .013 .086 .015 .091 .008 .086 .011 .103 .044 .046 .003 .047 c.sub.123 .235 .183 .284 .237 .178 .232 .205 .271 .133 .163 .135 .169 c.sub.131 .037 .054 .039 .043 .039 .040 .040 .036 .044 .026 .044 .028 c.sub.132 .015 .055 .016 .047 .012 .044 .013 .040 .010 .029 .009 .031 c.sub.133 .017 .072 .017 .063 .013 .067 .013 .056 .008 .045 .008 .044 c.sub.231 .000 .034 .000 .032 .001 .032 .001 .029 .002 .029 .002 .031 c.sub.232 .028 .048 .031 .046 .027 .047 .029 .043 .024 .043 .025 .041 c.sub.233 .013 .068 .013 .071 .013 .062 .015 .056 .012 .045 .012 .047 .sub.1 .024 .053 .014 .030 .038 .082 .035 .087 .041 .066 .041 .066 .sub.2 .021 .052 .013 .028 .035 .081 .033 .086 .038 .064 .038 .065 .sub.3 .007 .021 .005 .014 .008 .023 .006 .020 .008 .020 .008 .021

(145) B4.4 Model-Space Exploration

(146) 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.sub.ijk only affects the respective .sub.i and .sub.j errors. FIG. 25 demonstrates that the most important connection form in lowering both .sub.1 and .sub.2 is c.sub.123, since it describes a rotation of the helix angle, which has a large magnitude throughout the volume. c.sub.121 and c.sub.122 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. .sub.3 reduction is more variable although considerably small because of the stability of the heart wall normal.

(147) 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.

(148) 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 .sub.i for various embeddings as a function of neighborhood size .sub.i, using optimized computations. In relation to the frame error .sub.1, all embeddings with an unconstrained c.sub.123 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.sub.3 captures the extent to which the principal direction of diffusion remains normal to the heart wall normal. Since fiber directions u.sub.1 generally run tangentially to the heart wall, the resulting error is low for all embeddings.

(149) 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.

(150) 1) Helix Angle: c.sub.123 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.sub.123 throughout this distance produces a total change of 116.3, which is in close correspondence with values of 120 reported in the literature.

(151) 2) Wall Curvatures: c.sub.131 and c.sub.132 reflect the sectional curvatures of the heart wall as projected onto the local osculating ellipsoid to f.sub.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.

(152) 3) Myocyte Fanning: c.sub.122 and c.sub.133 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.

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

(154) 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 constraintsvia 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.

(155) 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.

(156) C. Tissue Engineering

(157) 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.

(158) 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.

(159) 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.

(160) 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.

(161) 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.