Muscular fibre estimation

Abstract

A method of estimating muscle fibre architecture for a muscular structure is described. The method comprises receiving a plurality of images that represent the motion of the muscular structure during a motion sequence, receiving an initial model of muscle fibre architecture, receiving a model of mechanical coupling between muscular structure motion and muscle fibres, and generating a joint function using muscular structure motion as represented in the plurality of images, the initial model of muscle fibre architecture, and the model of mechanical coupling. The method further comprises applying an optimisation procedure to optimise the joint function and, from the optimisation procedure, determining a model of muscle fibre architecture consistent with the muscular structure motion indicated in the plurality of images.

Description

MUSCULAR FIBRE ESTIMATION TECHNICAL FIELD

[0001] The present invention relates to estimating and modelling fibre architecture for a muscular structure.Embodiments described in detail relate to myocardial fibre architecture.BACKGROUND

[0002] Modelling myocardial fibre architecture is of interest in the field of medical research. A realistic heartmodel provides a platform for simulating the effects of drugs or implanted devices on the heart, and thus helpin developing treatments for various heart defects and diseases. However, the myocardial fibre architecture is still not completely understood, and current methods of estimating and modelling myocardial architecture have limited accuracy. Similar problems apply to modelling of other complex muscle structures in the body.

[0003] Two aspects are considered when modelling myocardial fibre architecture: the fibre architecture itself,which includes the structure of the myocardial fibres, and the movement of the fibres during a heart cycle. Regarding the fibre architecture, the myocardium, the middle layer of the heart muscle, is made up of individual fibres (myocytes) which arrange themselves into sheetlets. This structure is shown in Figure 1, which containsa schematic illustration of myocardial fibre and sheetlet architecture in the heart 10, alongside a zoom-in ofsheetlets 12 showing the individual myocytes making up the sheetlet. As shown in Figure 1, the myocytes locally align with each other along a first orientation within the sheetlet, and the sheetlets then orientate in a second, perpendicular direction. The orientation is important as cardiomyocyte orientation determines the preferential electrical wave propagation and tissue contraction in the heart, and so proper orientation of the myocytes is essential for valid and accurate heart models.

[0004] Current methods for determining the orientation of myocardial fibres and creating heart models includesimplified geometrical descriptions such as rule-based models (RBMs), and imaging methods such asDiffusion Tensor Imaging (DTI).

[0005] Considering firstly RBMs, RBMs assign myocyte orientations using mathematical descriptions basedon rules derived from histological and experimental observations. Using the mathematical descriptions, the myocyte orientations are incorporated into a cardiac computational model. However, the resultant model isextremely simplified since the model is based on general biomedical observations. Additionally, this methodcannot provide a personalised model of myocardial fibre architecture for a patient, or properly represent any variation in a population, as again, the mathematical rules describing the myocardial fibres are based on general observations. Thus, RBMs essentially provide a general model of the myocardial fibre architecture.

[0006] An alternative method for determining myocardial fibre architecture is DTI, which is a widely usedimaging modality for assessing tissue microstructure, primarily in the neuroimaging field. DTI captures theanisotropy of water diffusion in a tissue, and so enables determination of individual cell orientation and sheetlet orientation, and thus tissue structure. In more detail, a medical image that is sensitised to diffusion in a particular direction is generated through application of magnetic field gradients to a subject. This process isrepeated multiple times, applying the diffusion-encoding magnetic field gradients in different directions, togenerate a three dimensional diffusion tensor (DT) for each voxel of the image. The more directions used, themore accurate the resultant tensor. The three eigenvectors (E1, E2, E3) of each diffusion tensor (DT) are then estimated, and the eigenvectors correspond to cell orientation. E1 is the principal eigenvector, that defines the direction of fastest diffusion and thus ideally corresponds to the direction of the fibres in the tissue. When this method is applied to imaging the heart, (cardio (c)DTI), these vectors, E1, E2, E3, correspond to the myocyteorientation, ^^, the sheetlet orientation, ^^, and the sheetlet-normal direction, ^^, respectively, as illustrated inFigure 1.

[0007] Within cDTI, both ex-vivo and in-vivo methods exist. While ex-vivo cDTI is used for estimating themyocardial architecture, and is more accurate than in-vivo cDTI, ex-vivo cDTI has the obvious disadvantageof not being applicable to patients or volunteers, reducing it to some post-mortem donors, and thus reducingthe available acquisitions to a very small number of cases. If an estimate of the myocardial fibre architecture in a population is required, the small number of cases available would lead to an estimate and model that is not representative of the population. A heart model generated using this method would therefore not be an appropriate model on which to base medical research. In addition, keeping an extracted heart in its natural in-vivo configuration corresponding to any typical cardiac phase is extremely challenging. It is thereforeimpossible to guarantee that the imaged myocardial tissue state corresponds to any in-vivo cardiac phase, orto the blood-substitute filling-fluid volume and pressure, which introduces inaccuracies in the interpretation ofthe estimated fibre architecture.

[0008] In-vivo cDTI has also emerged in recent years as an alternative method for estimating andsubsequently modelling myocardial fibre architecture for an individual patient. For example, Ferreira et al.(Ferreira et al. (2014), In vivo cardiovascular magnetic resonance diffuse tensor imaging shows evidence ofabnormal myocardial laminar orientations and mobility in hypertrophic cardiomyopathy, Journal ofCardiovascular Magnetic Resonance, 16(1), 87) describe performing cDTI at two cardiac phases in 22patients, to measure cross-myocyte diffusion, and identified that patients with hypertrophic cardiomyopathy,demonstrated impaired reorientation of sheetlets in the diastolic phase of the heart cycle. However in-vivocDTI is still an experimental modality, and so is not yet acquired regularly in clinics. Inherent cardiac motionand the resultant deformation of the fibres means that more complex image sequences are required, which leads to a long acquisition time, a high noise level, and a high degree of inaccuracy when acquiring in-vivo cardio diffuse tensor images. Indeed, a disadvantage of DTI in general is the low signal to noise ratio, which can affect the accuracy of this technique.

[0009] In addition to the above-described disadvantages, neither of the above methods (RBMs and cDTI)take into consideration the movement of the myocardium during the heart cycle as a source of information.Myocardial movement is extremely important as movement leads to deformation of the myocardial fibres,which affects and is affected by fibre structure. Therefore, the response and influence of the fibres onmovement is desirable for a realistic, accurate model of the heart. An imaging modality that shows movementis cine-Magnetic Resonance Imaging (MRI). Cine-MRI images are regularly acquired in clinics, and are dynamic images that are able to represent movement of the heart during a complete cardiac cycle. To acquirea cine-MRI image a series of MRI images are acquired in one slice or a stack of parallel slices, over severalcardiac cycles, and then the multiple images are rearranged into a series showing a single normalised cardiaccycle. Several cine-MRI image series are typically acquired in different anatomical plane orientations. Ingeneral, in the short-axis several MRI slices are acquired at different heights of the heart, which enables themto be visualised as a temporal sequence of 3D images:. Figure 9a is an MRI slice shown at a central heartheight for three different time points 920, 940, 960, and Figure 9b is the same cine-MRI image shown at asingle time point for a stack of images at different heart heights. In the long-axis, one slice is typically acquiredfor each of three cardiac anatomical planes, as illustrated in Figure 9c, which also shows the intersection of the long-axis with one slice of the short-axis cine-MRI image.

[0010] Therefore, while there are some methods available for estimating heart fibre architecture, andgenerating models of the myocardium, current methods are limited in their use and accuracy.

[0011] It is therefore an aim of the present invention to provide methods for estimating myocardial fibrearchitecture using readily available, generic cardiac images of the human heart, where the method overcomessome of the limitations discussed above. SUMMARY OF THE INVENTION

[0012] In a first aspect, the invention provides a method of estimating muscle fibre architecture for a muscularstructure, the method comprising receiving a plurality of images that represent the motion of the muscular structure during a motion sequence, receiving an initial model of muscle fibre architecture, receiving a model of mechanical coupling between muscular structure motion and muscle fibres, and generating a joint function using muscular structure motion as represented in the plurality of images, the initial model of muscle fibre architecture, and the model of mechanical coupling. The method further comprises applying an optimisation procedure to optimise the joint function and, from the optimisation procedure, determining a model of muscle fibre architecture consistent with the muscular structure motion indicated in the plurality of images.

[0013] The optimisation procedure may comprise parameterising the muscular structure motion indicated inthe plurality of images to determine a set of motion parameters, parameterising the fibre architecture indicated in the initial model of muscle fibre architecture to determine a set of fibre architecture parameters, andoptimising the joint function by alternatively or jointly optimising the set of motion parameters and the set offibre architecture parameters until convergence is reached.

[0014] Optimising the set of motion parameters and the set of fibre architecture parameters may comprisecalculating the gradient of the joint function with respect to the set of motion parameters and the set of fibre architecture parameters, and applying a gradient-based optimisation algorithm to the joint function.

[0015] Determining the model of muscle fibre architecture may comprise using the converged set of motionparameters and set of fibre architecture parameters.

[0016] Muscular structure motion may be parameterised using a continuous time-varying velocity field.

[0017] Parametrising the muscle fibre architecture may comprise parameterising the local fibre orientationand parameterising the spatial variation in fibre orientation.

[0018] The function describing motion, the initial model of muscle fibre architecture and the model ofmechanical coupling may be probabilistic models.

[0019] The method may additionally comprise generating a conditional probability density function torepresent muscular structure motion as indicated in the plurality of images, for the initial model of muscle fibrearchitecture, and for the model of mechanical coupling. The joint function may be generated using theprobability density functions.

[0020] The plurality of images may be received from a cine-MRI.

[0021] The initial model of muscle fibre architecture may be population-based fibre information provided fromrule-based models.

[0022] The model of mechanical coupling may be derived from strain properties of muscle fibres.

[0023] The method may additionally comprise receiving muscle fibre information for an individual, generatinga function describing the muscle fibre information for the individual, and generating the joint function usingmuscular structure motion as represented in the plurality of images, the initial model of muscle fibre architecture, the model of mechanical coupling, and the function describing the muscle fibre information for the individual.

[0024] The muscle fibre information for the individual may comprise at least one diffuse tensor image, andthe function describing the myocardial fibre information for the individual may be a conditional probabilitydensity function.

[0025] The muscular structure may be a heart, the muscle fibre architecture may be myocardial fibrearchitecture, and the motion sequence may be a cardiac cycle. In such embodiments, the method of estimatingmyocardial fibre architecture may be used to model the function of a heart.

[0026] The muscular structure may be a uterus.

[0027] In a second aspect, the invention provides a method of image enhancement in imaging a muscularstructure. The method comprises estimating muscle fibre architecture according to the first aspect of the present invention, determining the muscular structure motion during the motion sequence, applying an inverse motion parameter to at least some of the plurality of images to provide images compensated for the muscularstructure motion, and combining the compensated images using image fusion or image super resolution toprovide at least one higher resolution image of the muscular structure.

[0028] In a third aspect, the invention provides a system for estimating muscle fibre architecture using aplurality of images indicating motion of a muscular structure. The system comprises an input engine forreceiving each of the plurality of images indicating muscular structure motion, and a function creator configuredto generate a joint function using muscular structure motion as represented in the plurality of images, an initialmodel of muscular fibre architecture, and a model of mechanical coupling between muscular structure motion and muscular fibres. The system further comprises an optimisation engine configured to optimise the joint function, and a model generator configured to determine, using the optimised joint function, a model of muscular fibre architecture consistent with the muscular structure motion indicated in the plurality of images. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Specific embodiments of the invention will now be described, by way of example, with reference to theaccompanying figures, of which:

[0030] Figure 1 is a histological image of myocardial sheetlets, alongside a schematic illustration of sheetletsshowing the individual myocytes making up the sheetlet;

[0031] Figure 2 is a schematic block diagram showing the overview of an exemplary system for implementingan improved method of modelling myocardial fibre architecture according to embodiments of the present invention;

[0032] Figure 3 is a flow diagram showing an overview of a method used to generate a model of myocardialfibre architecture in accordance with embodiments of the present invention;

[0033] Figure 4 is a flow diagram showing the input processing of Figure 3 in greater detail;

[0034] Figure 5 is a flow diagram showing the optimisation procedure of Figure 3 in greater detail;

[0035] Figure 6 is a flow diagram showing a method of gradient calculation, in accordance with embodimentsof the present disclosure;

[0036] Figure 7 is a flow diagram showing an example gradient descent optimisation algorithm that may beused for the optimisation procedure of Figure 3;

[0037] Figure 8 is a schematic diagram showing an example application for the heart model generated usingthe method of Figure 3;

[0038] Figure 9a is a cine-MRI image acquired in the short axis shown at three time points for a slice at acentral heart height;

[0039] Figure 9b is the same cine-MRI image shown at a single time point for the stack of slices at differentheart heights;

[0040] Figure 9c is a cine-MRI image at one time point, showing MRI images acquired in the long-axis foreach of the three typical anatomical planes, and intersecting with one slice of the short-axis cine-MRI image.. DETAILED DESCRIPTION

[0041] Prior to setting forth the invention, a number of definitions are provided that will assist in theunderstanding of the invention. All references cited herein are incorporated by reference in their entirety. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs.

[0042] The terms ‘cardiomyocytes’, ‘cardiac fibres’, ‘myocardial fibres’, ‘cardiac fibres’ and similar terms areused interchangeably herein and refer to the individual cells which make up the myocardium (the middle layerof heart muscle) of the heart. Cardiomyocytes align with each other and then aggregate to form sheetlets. Thecardiomyocytes are responsible for contraction of the heart.

[0043] The term ‘fibre orientation’ does not refer to individual cardiomyocytes, but to the average orientationof the aggregate of cardiomyocytes in the local neighbourhood, which are approximately aligned in the same 3-dimensional orientation. The term ‘fibre orientation’ also implicitly includes the sheetlet orientation.

[0044] The term ‘fibre architecture’ comprises the local fibre orientation and sheetlet orientation, and theirglobal pattern in the whole myocardium. The term fibre architecture is therefore the structure (organisation)and orientation of fibres within a muscle structure.

[0045] The term ‘initial model’ of fibre architecture does not refer to a model for an initial phase of a motionsequence (e.g. the initial phase of a heart cycle), but to an initial estimate of fibre architecture. Accordingly,the initial model is an approximate model of fibre architecture, and may be any model of fibre architecture that is realistic.

[0046] The term ‘generic estimate of fibre architecture’ refers to a first guess of the fibre architecture. Thegeneric estimate is not personalised for a particular patient: it is an approximate, broad estimation of fibre architecture. The same generic estimate can be used for different patients. The method of the present invention refines the generic estimate, to create a personalised, patient specific, model of muscle fibre architecture.Overview of Method

[0047] The method of generating a heart model according to an embodiment of the present inventioncomprises estimating cardiac fibre architecture (the orientation of fibres and sheetlets in the myocardium) using 2 principal inputs: a plurality of images that represent the myocardial motion during a heart cycle, and anapproximation of fibre architecture, provided by a prior probability distribution representing fibre architecturefor a population. The method also uses an expression that approximates the mechanical coupling betweencardiac motion and myocardial fibres, in order to take into account the effect of fibre orientation on cardiacmotion. Embodiments described in detail here relate to modelling of a heart – however, it should be noted thatthis approach can be applied to other muscular structures in the body, particularly those where motion takesplace without the control of the subject (such as the uterus, for example). Where reference is made below tothe heart, to myocardial fibres, and to a cardiac cycle, the skilled person will appreciate that the teaching provided may be applied to another muscular structure with its own fibre types and its own motion sequence.

[0048] The block diagram in Figure 2 shows the overview of an exemplary system 20 for implementing animproved method of modelling myocardial fibre architecture according to embodiments of the presentinvention. The system includes an Input Engine 22 operatively coupled to a Probability Function Creator 24.The output of the Probability Function Creator 24 is provided to an Optimisation Engine 26, which comprisesa Gradient Engine 32 and an Optimiser 34. The Probability Function Creator 24 and the Optimisation Engine26 are both in communication with a datastore 30 containing a model for physical constraints 36, a model ofan approximate fibre architecture 38, and algorithms 40, each of which will be discussed in greater detail later.The Optimisation Engine 26 is operatively coupled to a Model Generator 28, which outputs a realistic modelof the myocardial fibre architecture.

[0049] In use, the Input Engine 22 receives subject-specific inputs: an image sequence describing cardiacmotion and, optionally, cardiac fibre information for an individual. The image sequence may be, for example,a cine-MRI or a 3D echocardiography, but other image sequences are possible. Several stacks of images or single image slices acquired in different planes (for example, short-axis or long-axis images showing 2, 3, or4 heart chambers), or acquisitions using different image modalities and temporal resolutions, may becombined. The optional subject-specific fibre information may be obtained, for example, from a cDTI (or fromcombining several cDTI images acquired at different phases during the heart cycle). It should be noted thatsubject-specific fibre information, that is a cDTI image for a subject, is an optional input, and in some embodiments this input is not included.

[0050] The Input Engine 22 provides the received inputs to the Probability Function Creator 24, which alsoretrieves the model of physical constraints 36 and the model of approximate fibre architecture 38 from thedatastore 30. The model of approximate fibre architecture is a first guess (a generic estimate) of the fibrearchitecture and is used to initialise the model. The Probability Function Creator 24 processes each of thereceived inputs to provide a log-probability function to the Optimisation Engine 26. The Optimisation Engine26 optimises the received log-probability function to generate optimised parameters for the fibre architecture.In greater detail, the Gradient Engine 32 calculates the gradient of the log-probability function, and providesthe gradient to the Optimiser 34. The Optimiser 34 then retrieves a gradient-based (e.g. gradient-descent)algorithm 40 from the datastore, and applies this algorithm 40 using the calculated gradients, to determine theoptimal parameters for fibre architecture, where the optimal parameters are compatible with both fibrearchitecture and motion (as described in the input image sequence). The optimised parameters are output tothe Model Generator 28, which uses the parameters to create a model of myocardial fibre architecture. Thus,the method of the present invention comprises optimising the initial, approximate model of fibre architecture todetermine an accurate model of fibre architecture. The resultant model of fibre architecture is robust tovariations in the initial model of fibre architecture.

[0051] Alongside the one or two subject-specific inputs (the image sequence describing motion and fibreinformation), the Probability Function Creator 24 incorporates two further subject-independent expressionswhen generating the log-probability function: an approximate model of fibre architecture 38 and a model ofphysical constraints 36. Considering firstly the approximate model of cardiac architecture 38, the approximatemodel 38 is constant across all subjects and may be obtained using theoretical (such as RBM) or statisticalinformation on the population. Turning to the second subject-independent expression, the physical constraints36 are an expression representing the mechanical coupling between cardiac motion and myocardial fibres.The physical constraints 36 are described analytically using equations deduced from the physical properties of the tissues constituting the myocardium, where the equations relate the cardiac deformation and the fibre and sheetlet orientation. Including equations that couple fibres with motion is key for the estimation of thepersonalised fibre architecture, exploiting that the orientation of myocardial fibres affects cardiac motion. Theexpression for the physical constraints 36 is the same for all subjects.

[0052] The method 300 carried out by the system 20 and used to generate a model of myocardial fibrearchitecture according to embodiments of the present invention is shown in the flowchart in Figure 3. Firstly,one or several image sequences describing cardiac motion (which is subject-specific) is received at Step 310,and, optionally, subject specific fibre information is received at Step 320. The approximate model 38 of fibrearchitecture and the expression for the physical constraints 36 is then retrieved at Step 330, and all inputs areprocessed at Step 340. As described with reference to Figure 2, the input processing results in a log-probabilityfunction, and an optimisation function is then applied, at Step 350, to the log-probability function. The optimisation function outputs optimised parameters for the fibre architecture, that are compatible with the cardiac motion observed in the input cine-MRI image sequence. The optimised parameters are used to generate, at Step 360, a model of the cardiac fibre architecture.

[0053] The method of input processing 400, as carried out by the Probability Function Creator 24, is describedin greater detail by the flowchart in Figure 4. After receiving the image sequence describing motion and (ifavailable) the subject-specific fibre information, and retrieving the model of approximate fibre architecture 38and expression for the physical constraints 36, a conditional probability density function (PDF) for each inputis calculated at Step 410. The three or four generated PDFs are then combined, at Step 420, into a singlecombined posterior probability function (PPF) for the joint distribution of motion and fibre architecture. A log-probability function, which is defined as the negative logarithm of the posterior probability function, is thencreated at Step 430, and output to the Optimisation Engine 26.Calculating Probability Density Functions

[0054] As detailed above, the method of generating a model of the heart fibre architecture according toembodiments of the present invention uses two subject-independent expressions (a distribution of fibrearchitecture, which represents fibre information for a population (from which an approximate model of fibrearchitecture 38 can be determined), and an expression coupling fibre architecture with cardiac motion) andtwo subject-specific inputs (an image sequence describing motion and subject-specific fibre information).However, as described above, in some embodiments not all of the inputs are required, and the method maybe applied when the subject-specific fibre information is not provided. Equally, in some embodiments, multiplerepresentations of the same input may be provided, for example any number of independent image sequences representing motion or multiple sources of subject-specific fibre information.

[0055] As illustrated in Figure 4, a conditional probability density function is calculated for each input, and itshould be appreciated any suitable PDF may be used. The derivation of example PDFs that may be used insome embodiments are described below.

[0056] In the following PDF calculations several equations show the proportionality of a probability to anexpression. The proportionality factor depends only on given inputs, and is independent on the motion andfibres which are meant to be estimated. This makes the proportionality factors irrelevant for the estimation ofthe motion and fibre parameters, and thus the proportionality factors can be discarded from the resultant PDFs.Image Sequence Describing Motion (subject-specific)

[0057] A required input for the method of the present invention is a sequence of images representing themotion of the heart during a cardiac cycle. A sequence of images ^(^, ^) is given for a sequence of times^^, ^^, … … , ^^^^ during the cardiac cycle. Typically these times are equally spaced (^^ = ^^ / ^, where ^ is themotion period and the times can be considered cyclic (^^ = ^^). The sequence of images ^(^, ^) reflect themotion during a cardiac cycle. Therefore, the conditional probability density function is modelled as having aparticular image sequence given a motion: ^(^│Motion). Any number and combination of image sequencesmay be input. In some embodiments, sequences of cine-MRI are considered, however, in other embodiments,other modalities that show motion, such as 3D echocardiography, may be used. An example outlining thederivation of an appropriate PDF ^(^│Motion) for an image sequence ^(^, ^) obtained from a cine-MRIaccording to an embodiment of the present invention is described in detail below.

[0058] Firstly an ideal image sequence ^^ (^, ^) (an image sequence without noise) is conceptualised as thedeformation with motion of an ideal undeformed image ^^^(^): ^^(^, ^) = ^^^(^(^, ^, 0))Equation 1where ^ is a function representing cardiac motion, so that any pair of images corresponding to times ^^ and^^ can be related by the motion between times ^^ and ^^:^^(^, ^^) = ^^(^(^, ^^ , ^^), ^^)Equation 2

[0059] The real image sequence ^(^, ^) is modelled as a noisy version of the ideal sequence ^^ (^, ^), and sothe real image sequence ^(^, ^) satisfies Equation 2 only approximately. A reasonable approximation for cine-MRI images is to model the noise in the real image sequence ^(^, ^) as additive white Gaussian noise, (otherimaging modalities may require a more complex modelling of noise) and thus, the probability density functionfor the image sequence, given the ideal image sequence and the motion isEquation 3 where ^^^^represents image noise, and ^^ ^^ = 0 describes an image with no noise.

[0060] Equation 3 is dependent on the ideal image ^^^, and since the ideal image ^^^ is unknown, thisdependency is removed from the conditional PDF. Assuming a naïve homogenous prior distribution for ^^^, thejoint distribution ^^^, ^^^^Motion^, is proportional to the conditional distribution in Equation 3: ^^^, ^^^^Motion^ ∝^(^│^^^, Motion). To remove ^^^, ^^^ is marginalised by integrating for ^^^. To integrate for ^^^, firstly a change invariables is required: the summation over the voxel positions ^ is re-expressed as a summation over theundeformed positions ^^, which is achieved by approximating the summation by a continuous integral andassuming that the motion is incompressible. It should be noted that without the assumption of incompressiblemotion, contributions from the Jacobian of the transformation ^(^^, 0, ^) (from the reference at time 0 to theposition at time t of the image, given by the motion between both times) should be incorporated in the probability densities. Expressing the summation over the undeformed voxel positions ^^,^^^, ^^^^Motion^ becomes:where ^^(^^) = ^(^^, 0, ^).

[0061] The integral over ^^^ follows a standard procedure for the marginalisation of multivariate Gaussiandistributions, where the result is also Gaussian:

[0062] Assuming again incompressible motion, the summation over points ^^ can be approximated by theequivalent summation over points ^ in any other time,^(^|Motion)Equation 6

[0063] The exponent is proportional to the typical sum of square error dissimilarity metric between each pairof images in the input image sequence.

[0064] It should be noted that in Equation 6 all pairs of images are used to define the PDF. While anembodiment of the present invention calculates the conditional PDF for the image sequence using Equation6, and each image is compared with every other image in the sequence, in other embodiments a reduced setof comparisons may be used, for example considering only consecutive times:{(^^ , ^^ ), , ^^ ), … , (^^^^ , ^^ )}, or the first ^ neighbours:{(^^ , ^^ , … , ^^), (^^ , ^^ , … , ^^^^), … , (^^^^ , ^^ , … , ^^^^)}.Approximate Model of Fibre Architecture (Population-based fibre information) - Subject-independent

[0065] Population-based prior information 38 on myocardial fibre architecture is also required for the methodaccording to the present invention. Population-based prior fibre information 38 is subject-independent (and sothe expression is constant across all subjects), and based on available information of the population. Thepopulation-based fibre information 38 can be estimated from any general biological observations on themyocardial fibre architecture, such as the biological or histological data integrated in RBMs or from a statistical atlas built from the fibre architecture of an independent set of representative subjects. The prior information onthe fibres provides an initial generic estimate of the myocardial architecture, and thus a probability distributionthat reflects the uncertainty in the estimate, and is dependent on the source of the initial estimate Source^, ismodelled ^^Fibre│Source^^. The derivation of an example PDF that may be used to model the fibre architectureprior distribution is outlined below, but again, any suitable PDF may be used.

[0066] As described previously, the myocardial fibre architecture comprises local myocytes oriented in anapproximately parallel direction, which then aggregate into planar sheetlets. Thus the myocardial fibrearchitecture can be modelled as a field of frames (orthonormal triads) {^^, ^^, ^^ }, which are defined in theundeformed myocardium (at the reference configuration, t = 0). The first vector, ^^, represents the orientationdirection of the myocardial fibres, the second vector, ^^, represents the orientation direction of the sheetlet,perpendicular to the fibres, and the third vector, ^^, represents the normal to the sheetlet (as shown in Figure1).

[0067] A fibre frame at a particular point may be parametrised in different ways, including by angle and axisof rotation, by quaternions, or by Euler angles. Representing a frame at a particular point using Euler angles is analogous to the parameterisation used in RBMs, where the fibre frame is given by the three angles, helical(^^^), transverse (^^^), and sheet (^). These three angles represent rotations aligned with and applied to aheart-adapted fixed reference frame, {^^, ^^, ^^} (as described in Bayer, J. D., Blake, R. C., Plank, G., &Trayanova, N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation tocomputational heart models. Annals of biomedical engineering, 40, 2243-2254), denoting the circular,longitudinal, and transmural directions following known conventions.

[0068] The spatial variation of the field of frames can also be parametrised in different ways: as a dense fieldwith some regularisation or with a lower dimensional set of parameters based, for example, on splines. In theexample discussed below, Euler angles and 3D b-splines on the adapted coordinates of the myocardium(circular, longitudinal, and transmural), (as described in Bayer, J., Prassl, A. J., Pashaei, A., et. al. (2018).Universal ventricular coordinates: A generic framework for describing position within the heart and transferringdata. Medical image analysis, 45, 83-93) are used to parameterise the field of frames.

[0069] Using Euler angles and B-splines, the field of frames, ^^, ^^ } are expressed as a function of aof spline parameters: {^^(^; {^^})}. The probability distribution of the fibres is therefore also expressed in termsof the spline parameters:^(Fibre|Source^) = ^({^^}|Source^)Equation 7As shown in Equation 7, the probability distribution depends on the information available from the Source^, anda typical and versatile model is a multivariate Gaussian distribution on the parameters:Equation 8where the mean parameters’ value, ^ ^^^^ and precision matrix, ^^^, depend on Source^.Physical Constraints – Subject Independent

[0070] The method according to the present invention also uses an expression for the set of physicalconstraints representing the coupling between cardiac motion and myocardial fibre architecture. Theexpression for the physical constraints 36 is subject-independent and based on an approximate modelling ofthe mechanical relationship between local fibre orientation and tissue deformation given by the motion. Thiscoupling is not deterministic but modelled as a conditional probability density, ^(Motion│Fibres). The derivationof an example PDF is outlined in detail below, including first the derivation of a set of physical constraints 36.

[0071] The physical constraints 36 are related to fibre-motion coupling, and thus firstly, the transformation ofmyocardial fibres with motion is considered. This transformation can be formulated in different ways, and oneexample is described below for illustrative purposes.

[0072] The field of frames (orthonormal triads) {^^, ^^, ^^)} representing the fibre architecture of theundeformed myocardium (at t=0) do not propagate with motion as vectors. The field of frames propagate in acomplex way, however, for simplicity, their vectorial propagation can be considered as an auxiliary field: ^^(^^ , ^) = ^(^^, ^)^^(^^)Equation 9where ^^ = ^^(^^) = ^(^^, 0, ^) denotes the motion from the undeformed configuration to the deformedconfiguration at time t, and ^(^^, ^) is the Jacobian tensor of the transition. To obtain the motion value at agiven point ^^ in the deformed tissue at time ^, the inverse transformation is required: ^ ^^^ = ^^ (^^) =^(^^ , ^, 0) .

[0073] The auxiliary field defined in Equation 9 is not orthonormal, and does not properly represent thedeformation of the fibre architecture. However, given the hierarchical definition of the orientation directions ofmyocardial fibres, an appropriate transformation can be obtained by the Gram–Schmidt orthonormalization ofthe auxiliary field in the hierarchy order, which is that the first vector ^^ represents fibre axis, the second vector^^ represents the sheetlet plane, and the third vector ^^ represents the sheetlet normal. Equations 10 and 11below show the orthonormalization method used:Equation 10 represent intermediate vectors which are orthogonal but not unit vectors (normalised). The fibreframes transformed by the motion ^^ are obtained by normalisation of the intermediate vectors:Equation 11

[0074] The physical constraints 36 used in embodiments of the present invention are derived from theimportant physical characteristic that there is highly rigid connectivity between fibres in the same sheetletcompared with the easy sliding between sheetlets. This property is part of the known mechanical properties ofthe myocardium, and is reflected in the very small shear strain expected on the sheetlet plane in comparisonwith any transversal direction. These physical constraints 36 can be expressed using the Green-Lagrangianstrain tensor E, and its projections ^^^to the undeformed fibre frame: ^^^(^^, ^) = ^^(^^)˕^(^^, ^)^^(^^)Equation 12Thus, the ideal deformation with shear between sheetlets occurs if ^^^(^^, ^) = 0. This equation implies aconstraint between motion and fibres.

[0075] Expressing the Green-Lagrangian strain tensor E in terms of the Jacobian tensor ^, the model forphysical constraints 36 described in Equation 12 becomes^^^(^^, ^) =1 ^^(2^ ^^^ − ^)^1 ^=^^^^1 ^^ = 2^ ^2^^(^^, ^)^^(^^)^ ∙ ^^(^^, ^)^^(^^)^Equation 13where ^^ ∙ ^^ = 0. The physical constraints 36 as described in Equation 13 can be expressed with the auxiliaryfield defined by the vectorial propagation (Equation 9 above):^^^(^^,Equation 14

[0076] In some embodiments, quasi-incompressibility of the myocardium may also be considered: the volumechange that arises from deformation may be encoded by the determinant of the Jacobian tensor, det ^(^^, ^) =1, or by the divergence of the Eulerian velocity field of the motion ^ ∙ ^(^, ^) = 0. However, incorporating thisproperty is optional, and is not discussed in detail.

[0077] As the model of the physical constraints 36 described by Equation 14 is an approximation, aprobabilistic description of the constraints is used. A Gaussian-like conditional probability is considered for themotion given the fibre architecture:^(Motion│Fibres)Equation 15

[0078] In Equation 15 ^^^ is the constraint standard deviation, and so indicates how precise the physicalconstraints 36 are expected to be. For an exact constraint, ^^^ → 0, but as the expression for the physicalconstraints 36 is an approximation, ^^^ > 0. The inclusion of the motion period T and the myocardium volume^^^^^ ensures that ^^^ is dimensionless. The ^^^ is expected to be small.Fibre Information – Subject-specific (Optional)

[0079] In some embodiments subject-specific fibre information, for example a cDTI, (or any similar imagemodality that allows extraction of the local fibre orientation) is also included, and used to generate the modelof myocardial fibre architecture. This input is optional, however it is advantageous to incorporate this input ifavailable as a subject-specific cDTI provides specific information on local fibre orientation for a subject, andthis information improves the precision and accuracy of the estimated fibre architecture. Any number of cDTI images can be combined for this input, including acquisitions in different slices, plane orientations and cardiacphases (such as atrial diastole, ventricular diastole, among many others).

[0080] The available fibre images are described by a conditional probability, ^^DTI│Fibre, Motion^, wherethe images are dependent on the fibre architecture and myocardial motion. It should be appreciated that if the DTI are acquired at a single cardiac phase taken as reference, there is no dependence on the motion. However, for DTI acquired at several cardiac phases, the dependence on the motion cannot be neglected,since, as described above (Equations 9 to11), the fibre frame at time t, ^^(^^ , ^), depends on the fibres at thereference configuration, ^^(^^), and the motion, ^^(^^).

[0081] Multiple DTI can be used in some embodiments of the present invention. When a series of DTI atdifferent cardiac phases are present, the estimation of the fibres incorporates the information present in all of the images in a single fibre model.

[0082] The derivation of an example PDF that may be used when DTI are acquired at several cardiac phases,and thus there is dependence on motion, is outlined below.

[0083] At each voxel / pixel of the DTI, the ideal diffusion tensor ^ is expected to have as eigenvectors thefibre frame: ^ = ∑^ ^^^ ^^ ^^ ^ ^^ , with diffusivity eigenvalues in decreasing order, ^^ ^ ^^ ^ ^^ ^ 0. However, theacquired image will be affected by noise, resulting in a probability distribution of the measured diffusion tensoraround the ideal one. This tensorial probability distribution is complex but a practical approximation used inthe method according to the present invention is a Gaussian distribution with isotropic variance ^^^^^for all tensor components:16where the second expression results from 2 characteristics: the frame {^^} is orthonormal, and ∑^ ^^^ = ^^^^^^.In addition, ^^^^^^ is independent of the fibre frame {^^}. Thus, the ^^^^^^ term does not affect the fibresestimation and can be disregarded as a constant factor, so that the relevant proportional part is17

[0084] Equation 17 is a matrix Bingham distribution (a normal distribution for the space of orthonormalframes). This is an antipodally symmetric probability distribution, in agreement with the property that fibres andsheetlets have no defined sense or sign, so that ±^^ represent the same direction.

[0085] Equation 17 models the local fibre orientation at one point and time {^^(^^ , ^)} corresponding to thediffusion tensor at one pixel in one DTI image. Therefore, the complete probability density involves thesummation for all DTI pixels:

[0086] It should be appreciated that the points over which the summation is performed are the DTI pixelswhere information on the fibres can be extracted. Thus, the points do not coincide, in general, with the voxelsover which summation was performed on the image sequence ^(^, ^) (Equation 6). For instance, if the fibreinformation is extracted using only a few cDTI slices, the summation should run through the pixels segmented as myocardium on those DTI slices.

[0087] The four probability density functions described above are example PDFs that may be used for eachinput, but as outlined, any other appropriate PDF may be used.Joint Posterior Probability Density Function

[0088] Returning to Figure 4, as illustrated Step 2 of input processing comprises combining the four (or threeif subject-specific fibre images are not included) conditional probability density functions into a singularposterior joint probability of motion and fibres, which is conditional to the input image sequence, DTI, and the source of population-based fibre information: ^(Motion, Fibres|I, DTI, Source^) ∝ ^(Motion, Fibres, I, DTI| Source^)= ^(I, DTI|Motion, Fibres) × ^(Motion, Fibres| Source^)= ^(I|Motion) × ^(DTI|Fibres, Motion) × ^(Motion|Fibres) × ^(Fibres| Source^)Equation 19

[0089] Equation 19 is simply the product of the four probability density functions described previously. Inderiving Equation 19, three properties involving conditional independence have been used. First, the imagesequences and DTI are conditionally independent once fibres and motion are given:^(^, DTI|Fibres, Motion) = ^(I|Fibres, Motion) × ^(DTI│Fibres, Motion)Equation 20 Second, the image sequence is conditionally independent of the fibres once the motion is given: ^(I|Fibres, Motion) × ^(I│ Motion)Equation 21Third, the motion is independent of the population-based prior fibre information once the fibres are given:^(Motion|Fibres, Source^) = ^(Motion│Fibres)Equation 22

[0090] In embodiments where no subject-specific DTI are available, and the method only uses three inputs(an image sequence describing motion, an expression for the physical constraints of the fibres 36, and anapproximate model of fibre architecture 38) the method is still applicable and the posterior probability functionis simplified to: ^(Motion, Fibres|I, Source^) ∝ ^(I|Motion) × ^(Motion|Fibres) × ^(Fibres│Source^)Equation 23

[0091] The two Posterior Probability Functions described in Equations 19 and 23 are simple cases ofhierarchical Bayesian models. Typically, Bayesian estimation uses Gibbs posterior sampling, however the present invention involves motion tracking (that is, heart movement), which is better modelled as an optimisation problem. Therefore, the algorithm for optimising the fibre architecture parameters is developed asa maximum posterior estimation. It should be appreciated however that in some embodiments, Gibbs posteriorsampling may be used to optimise the Posterior Probability Functions. Log-probability Function

[0092] Determining the maximum of the joint posterior probability function is equivalent to determining theminimum of the negative logarithm of the posterior probability function, so named the log-probability function.Thus, as illustrated in Figure 4, the final step of input processing comprises creating a log-probability function.The derivation of two example log-probability functions, using the posterior probability functions in Equations 19 and 23, are outlined below: when subject-specific DTI are provided: −log ^(Motion, Fibres | I, DTI, Source^) = − log ^(^ | Motion) − log ^(DTI | Fibres, Motion)− log ^(Motion | Fibres) − log ^(Fibres | Source^) + constantEquation 24 Or, when no subject-specific fibre information is provided −log ^(Motion, Fibres | I, Source^) = − log ^(^ | Motion) − log ^(Motion | Fibres) − log ^(Fibres | Source^)+ constantEquation 25

[0093] The constant in Equations 24 and 25 is dependent on the inputs but independent of the fibre andmotion parameters that are to be estimated. Thus, the constant is irrelevant for this estimation (the optimal isthe same irrespective of the value of this constant) and is therefore ignored when computing the negativelogarithm of the posterior probability function.

[0094] The resulting expression for the log-probability function is a summation of positive defined metrics,most of which are quadratic. Substituting Equation 24 with the negative logarithm of the PDFs described inEquations 6, 18,15, and 8, the log-probability function becomes:− log ^(Motion, Fibres | I, DTI, Source^)Equation 26The first term in Equation 26 relates to the image sequence representing motion, the second term to subject-specific fibre information, the third term to the set of physical constraints 36, and the final term to the population-based fibre information 38. In embodiments where no subject-specific fibre information is provided, the secondterm in the equation (−log ^(DTI | Fibres, Motion), is not included.

[0095] The precision of the maximum-posterior estimates (corresponding to the minimum of the log-probability function obtained using Equation 26) can be quantified by the second derivatives (Hessian matrix)of Equation 26, but this will not be described in detail here. Optimisation Procedure

[0096] The optimisation method 500 used in embodiments of the present invention, where the minimum ofthe log-probability function is determined, is illustrated by the flowchart in Figure 5.

[0097] Once the log-probability function is generated, the gradient of the log-probability function is calculated,at Step 510, with respect to both the fibre and motion parameters. The resultant gradients are then used to apply, at Step 520, a gradient descent algorithm, which uses an alternation of partial gradient descents for each parameter type (fibre and motion) in order to find the minimum of the log-probability function. Gradient descent is an iterative optimisation algorithm for finding a local minimum of a differentiable function, and so this optimisation method is well-suited to the present application. The output of the gradient descent algorithm is a set of optimised parameters representing motion and fibre architecture, which can be used to construct a heart model. The resultant optimised parameters for the fibre architecture will be compatible with the motion observed in the input image sequence, and thus will provide a heart model with increased accuracy compared to models generated using prior art methods. Gradient Calculation

[0098] As described above, to determine the minimum of the log-probability function, the gradient of the log-probability function must be calculated. Since the log-probability function (Equation 26) is the sum of different terms, the gradient is the sum of the gradient of each term, and thus the gradient of each term is calculated separately.

[0099] The method of gradient calculation 600 is described in the flowchart in Figure 6. Firstly, an expressionfor the gradient with respect to the motion parameters for each term in the log-probability function is determinedat Step 610. This includes the gradient of image dissimilarity with respect to motion in the input imagesequence, and the gradient of the fibre–motion coupling (physical constraints term). If subject-specific DTI isavailable in different cardiac phases, the gradient of the corresponding term with respect to motion is also included. It should be noted that at this stage, the motion parameters are not concreted, and are genericallydenoted by ^ = {^^}. The motion is then parameterised, at Step 620, using B-splines to obtain an iterativeformula that parameterises the motion using parameters ^^,^. It should be noted that while this is shown asoccurring after expressions for the gradient with respect to motion parameters are obtained, parameterisingthe motion can occur in parallel with or before determining expressions for the gradient. The motion parametersobtained at Step 620 are then used in conjunction with the gradient expressions determined at Step 610 to calculate, at Step 630, the gradient of each term with respect to motion explicitly using the motion parameters ^^,^. It should also be noted that explicit calculation using motion parameters may not be required for computing the gradient of every term in a function. For example, in Equation 45 (the gradient of motion with respect to motion parameters) only the first term requires explicit motion parameters ^^,^, This is useful when the method is repeated many times (for example, receiving inputs from multiple different subjects), as it means the method requires less computational power.[000100] An expression for the gradient with respect to the myocardial fibres for each term in the log-probability function is then determined at Step 640. At this stage the myocardial fibres are described generally,and in some embodiments this method step may occur before or in parallel with determining the gradients withrespect to motion. The fibre architecture is parameterised, at Step 650, using Euler angles with spatialdependency parameterised in turn by B-splines on myocardium adapted coordinates. Thus, the gradientsdetermined at Step 640 are calculated explicitly using the fibre B-spline parameters ^^,^ at Step 660. Once allgradients are calculated, the gradient descent algorithm is applied at Step 670 in order to determine theminimum of the negative log-probability function, and thus determine the optimal parameters for generating amodel of the myocardial architecture.[000101] Details on how each of these gradients may be calculated, as well as an example optimisationalgorithm, is described in detail below.Gradient with respect to motion[000102] As illustrated in Figure 6, an expression for the gradient of each term in the log-probabilityfunction with respect to motion parameters is determined. At this stage, the motion is considered to beparametrised by a set of general parameters ^ = {^^}. The dependence on motion appears in three terms ofthe log-probability function: log ^(^ | Motion), log ^(Motion | Fibres), and log ^(DTI | Fibres, Motion), if included(relating to the image sequence describing motion, the expression for physical constraints 36, and the subject-specific fibre information respectively).[000103] The first of these terms is the image dissimilarity, relating to the input sequence of images thatrepresent cardiac motion. Computing the gradient with respect to motion of the image dissimilarity term in thenegative log-probability function (Equation 26) gives:Equation 27[000104] Two new factors appear compared to Equation 26: the gradient of the image with respect toposition, ∇^^(^, ^) (which is computed numerically using finite differences or an appropriate convolutional filter,^^^as is typical in the field), and the gradient of the motion with respect to the motion parameters,.[000105] The term representing the image dissimilarity, −log ^(^ | Motion), is the most computationallyexpensive term when computing the gradient. For this term, all pairs of images are considered: ^= {(^, ^)|^, ^ ∈ [0, ^ − 1], ^ ≠ ^)}Equation 28 However, in some embodiments the method may begin with a relaxed objective function, considering a subsetof pairs, which can be progressively increased in a sequence, ^^ ⊂ ^^ ⊂ ⋯ ⊂ ^, defining stages of theoptimisation algorithm, analogous to a usual multiresolution cascade. This method improves both the computational efficiency and the stability of the function. For the initial subset of image pairs, only consecutive times in a cycle are considered ^^ = {(0,1), (1,2), … , (^ − 2), (^ − 1,0)}Equation 29[000106] Initially considering consecutive times in some embodiments is advantageous as consecutiveimages are the most similar, and so the gradients obtained from each of the pairs are more accurate estimates of the direction of optimisation. After several iterations of the optimisation algorithm, all warped images are expected to be close enough to work similarly well for the next iterations. The inclusion of more distant pairsof images is expected to refine the motion extracted from the plurality of images, increasing the stability toimage noise, and reducing cumulative errors in the full cardiac cycle motion from small over / under-estimations of the motion between consecutive pairs.[000107] The second term in the log-probability function that is dependent on motion is the expressionrepresenting the fibre-motion coupling (the physical constraints 36 of the fibre architecture)log ^(Motion| Fibres). Computing the gradient for this term with respect to motion gives:∂log ^(Motion| Fibres)− ^^^Equation 30[000108] Two gradients of motion derived quantities appear: the gradient of the Jacobian tensor,and the gradient of the divergence of the velocity field,^^^.[000109] The third term dependent on motion is, if included, the subject-specific DTI fibre information.The gradient of this term with respect to motion only appears when DTI is available for more than one cardiacphase. When only a single cardiac phase is available, this term is not included.[000110] Computing the gradient of the DTI term in Equation 26 with respect to motion gives:^(DTI | Fibres, Motion)− ^^^Equation 31 One new gradient factor has appeared: the gradient of the transformed fibreGiventhe frame definition above (Equations 10 and 11), by orthonormalizing the vector transformation of the framein the undeformed reference time, the gradient of the transformed frame with respect to motion can beexpressed in terms of the angular velocity tensorEquation 32 Equation 33 Equation 34usingEquation 35 Equation 36 Equation 37where, from Equation 9,^^^(^^,^)^^=^^(^ ). Accordingly, the gradient of the Jacobian Tensor is the sameas shown above for the gradient of the fibre–motion coupling.[000111] Three gradient factors have thus appeared when computing the gradient with respect tomotion for the log-probability function:To explicitly calculate these gradients, motion must first be parametrised. Motion parameterisation[000112] In the method according to embodiments of the present invention, parameterising the motioncomprises parametrising a velocity field, and then modelling the motion using the parameterised velocity fields.In embodiments, the motion ^ is parameterised by a continuous time-varying (Eulerian) velocity field ^(^, ^),which ensures that there is diffeomorphic motion when integrated:^^^(^^, ^^ , ^^) = ^^ + ^ ^(^(^ ^^ , ^^ , ^ ), ^^)^^′^^Equation 38[000113] Diffeomorphic motion means that the motion is bijective, differentiable, and has a differentiableinverse, which is a desirable property for most deformable image registration problems, especially for motion-tracking of deformable tissues, since it enforces tissue continuity along time, avoiding degeneratesuperposition of different regions on the same position.[000114] Turning to describing the velocity field, in some embodiments the velocity field may bedescribed using a dense representation, however this requires regularisation, which is typically in the form ofGaussian smoothing. An alternative representation is to describe the velocity field using a spatio-temporal B-spline diffeomorphic framework. The spatio-temporal B-spline framework parameterises the velocity field usingspatio-temporal (4D) B-splines:Equation 39where the multi-index notation ^ =is used to allow a compact notation (each index locates acontrol point, in four dimensions). Each index ^^ runs through each of the B-spline control points in thecorresponding dimension. The motion is then parameterised by the set of coefficients ^^^, where the number of coefficients is dependent on the B-spline resolution, that is, the number of B-spline control points. Theresultant B-spline parameters ^^^ can thus be used as motion parameters ^^.[000115] The B-spline description allows the continuity and smoothness to be controlled both spatiallyand temporally. In addition, the strain rate tensor (describing the rate of change of deformation), relating to thegradient of the velocity field, can be computed analytically and expressed linearly in the B-spline parameters ^^^:Equation 40The partial derivative of the cubic B-splines is separable and quadratic in the involved coordinate ^^. Thegradient of ^^(^, ^) is therefore:^ ^ ^^(^, ^) = ^^ ^^^(^^)^ ^^^ (^) ⇒ ∇^^^(^, ^) = ^ ^ ^^^ (^ ) ^ ^ ^^^(^^)^ ^^^ (^)^^^ ^^^,^^^Equation 41These derivatives are required for the velocity divergence, ∇^^^(^, ^), appearing in the negative log-probabilityfunction (Equation 26).[000116] Returning to Equation 38 and calculating the motion from the velocity field, the integral definingthe motion from the velocity field (Equation 38) is numerically integrated. This may be carried out using any ofthe Runge–Kutta methods of different orders. In the present invention, a first order Euler method is used, as outlined below, however other implementations are possible. Following the Euler method:Equation 42[000117] The time interval ∆^ is split into ^ intervals, so that ∆^ =^ . This can be formulated as the recursion: ^(^^ , ^^ , ^^ + ∆^) ≅ ^(^^ , ^^ , ^^) + ∆^ ^(^(^^ , ^^ , ^^), ^′)Equation 43where ^^ = ^^ + ^∆^ for any ^ = 0, … , ^ − 1. It should be noted that the recursion is equally valid for integralsforward in time (∆^ > 0) or backwards in time (∆^ < 0).[000118] Since time resolution of the cine-MRI image is sufficiently high, one interval per image can beused, so that ^ = ^ − ^ (it should be appreciated, however, that more refined integrations may also be used).Thus Equation 43 is simplified, and the iterative formula for describing motion is:^(^^ , ^^ , ^^^^) ≅ ^(^^ , ^^ , ^^)) + ∆^^^(^(^^ , ^^ , ^^), ^^)Equation 44∆^^ = ^^^^ − ^^, which for an equally spaced time sequence (which is typical for cine-MRI images) is ∆^^ =^ / ^. It should be appreciated that for backwards motion, the recursion is changed so that ∆^^ = ^^^^ − ^^ =−^ / ^.Gradients with respect to motion – Explicit Calculation[000119] As shown above, when first deriving the gradients with respect to the motion parameters, thethree gradient factors were defined with respect to a set of motion parameters ^ = {^^}. Once the motion andvelocity divergence are represented using the spatio-temporal B-spline diffeomorphic framework (Equations38 and 39 respectively), the gradient expressions for the motion features are computed explicitly in terms ofthe B-spline motion parameters, ^^,^. This includes the gradient of the motion function, ^^(^, ^^, ^^)^^^,^, the gradient an tensor,( )^^∙^(^,^) of the Jacobi^^ ^^,^^^^,^, and the gradient of the divergence of the velocity field,^^^,^.[000120] Firstly, considering the gradient of the motion function, from the iterative formula above(Equation 44), a recursive formula for the gradient can be obtainedEquation 45with the initial conditions ^(^ , ^^ , ^^) = ^^ ,^^^,^= 0, and the spatial gradient of the velocity givenEquation 40.[000121] Turning to the gradient of the Jacobian tensor, in order to calculate the gradient, the JacobianTensor ^^^(^^, ^) must first be calculated, where the Jacobian Tensor is given by:Equation 46 Multiple methods may be used to calculate the Jacobian Tensor, for example numerically, by finite differences,by an appropriate convolutional filter, or using recursion. An example of calculating the Jacobian Tensor andit’s gradient using recursion is described below.[000122] Using the equation of motion defined in Equation 38 the following integral formula for theJacobian can be derived:Equation 47Equation 47 can be approximated with a recursion, simultaneously or after computing motion:Equation 48[000123] Accordingly, the gradient of the Jacobian with respect to the motion parameters is then alsoobtained recursively:Equation 49 The above recursions, simultaneously with the iterative formula describing the motion (Equation 44) specialised for a = 0, can be performed independently per point ^^, without needing to store the full field and allowing its trivial parallelisation.[000124] Finally, considering the gradient of the velocity divergence with respect to the motionparameters ^^,^(^^∙^(^,^)^^^,^), the divergence of the Eulerian velocity field is first obtained analytically using thespatial gradient of the velocity (Equation 40):^ ∙ ^(^, ^) = ^^^∇^^^(^, ^).Equation 50 This is linear in motion parameters, and so the gradient is simplyEquation 51In this way, the three required gradients with respect to the motion parametershave been computed(Equations 45, 49, and 51). Parameterising the fibre architecture[000125] Returning to Figure 6, once the gradients have been computed with respect to the motionparameters ^^,^, the fibre architecture is parameterised. While parameterising the fibre architecture is shownas occurring after computing the gradients with respect to motion, in some embodiments parameterising thefibre architecture may occur before or in parallel with computing motion gradients.[000126] As discussed previously, at each point ^^ of the undeformed configuration of fibres, the localfibre orientation is described by an orthonormal triad, also known as a fibre frame {^^}. A fibre frame can berepresented and parametrised in different ways, for instance by angle and axis of rotation, by quaternions, orby Euler angles. Any of these representations may be used for the present invention.[000127] In an example of the present invention, Euler angles are used to represent fibre frames, sincethis is analogous to the parametrization used in RBMs, where the fibre frame is given by the three angles,helical (^^^), transverse (^^^), and sheet (^) angles, representing the sequence of rotations ^ aligned with andapplied to the heart-adapted fixed reference frame, {^^, ^^, ^^}, denoting the circular, longitudinal, andtransmural directions following the convention described in Bayer et al 2012 (paper cited previously).Equation 52 ^^ = ^(^^^, ^^^, ^)^^, ^^ = ^(^^^, ^^^, ^)^^, ^^ = ^(^^^, ^^^, ^)^^Equation 53 Equation 54 Equation 55Using Euler angles to parameterise fibre frames allows an intuitive interpretation of the resulting parametersand to reuse part of the method in RBMs.[000128] Expanding the rotations ^^, ^^, ^^ into components, the resulting fibre frame is then written interms of the myocardium adapted reference frames {^^, ^^, ^^} as:^^ = cos ^^^ cos ^^^ ^^ + sin ^^^ cos ^^^ ^^ + sin ^^^ ^^Equation 56^^ = (cos ^^^ sin ^^^ sin ^ − sin ^^^ cos ^) ^^ + (sin ^^^ sin ^^^ sin ^ + cos ^^^ cos ^) ^^ − cos ^^^ sin ^ ^^Equation 57^^ = −(cos ^^^ sin ^^^ cos ^ + sin ^^^ sin ^) ^^ + (−sin ^^^Equationwhich provides the fibre frames at each point,{^^(^^)}, from the angles at this same point,^^^(^^), ^^^(^^), ^(^^).[000129] The spatial variation of the field of frames can also be parametrised in different ways: as adense field with some regularisation (smoothing) or using a lower dimensional set of parameters based, forinstance, on splines. In an embodiment of the present invention, B-splines on the adapted coordinates of themyocardium (circular, longitudinal, transmural) are used to parameterise the spatial variation in the field offrames (consistent with the examples provided above for calculating a PDF for population-based fibreinformation).[000130] Each control point in the B-spline representation corresponds to a finite-support, basefunction, ^^(^^), defined in any point ^^ of the myocardium. Analogous to the B-splines parametrising themotion velocity field, the number of control points in each of the adapted coordinates (circular, longitudinal,transmural) controls the corresponding degree of variability of the angles along that coordinate. Thus, the setof parameters for the fibre frame (denoted above generically as {^^}) will be the set of B-spline control points, for each of the three Euler angles (^^^,^^^, and ^), where L is an index that runs through all controlpoints, and m is an index that runs through each of the angles:Equation 59 Equation 60 Equation 61[000131] To compute the gradient of the fibre frame (which describes the fibre orientation at each point)with respect to the fibre orientation parameters the derivative chain rule is applied:Equation 62 Equation 63 Equation 64To obtain the expressions of the appearing derivatives with respect to each of the three Euler angles is straight-forward from Equations 56, 57, and 58.Gradients with respect to fibre orientation[000132] Once the fibre-orientation is parameterised, the gradients with respect to the fibre orientationfor the log-probability function can be calculated. The dependence on the fibre orientation parameters appearsin three terms in log-probability function: log ^(DTI | Fibres, Motion), log ^(Motion | Fibres), andlog ^(Fibres | Source^) (corresponding to subject-specific fibre information, the expression for physicalconstraints 36, and the population-based fibre information respectively 38).[000133] The first of these terms relates to the subject-specific fibre information, provided, for example,by DTI. The gradient of the DTI term (described in Equation 26) with respect to the fibre orientation parameters ^^,^is: ^log ^(DTI | Fibres, Motion)2^( ) ^( )^^^ (^^ , ^)− = ^^^,^^^ ^ ^ ^ ^^ ^^ , ^ ^^ ^^ , ^ ^^^(^^ , ^)^^^^,^^ ^ ^ ^^ Equation 65[000134] As discussed above, this term only appears if subject-specific DTIs are used as an input.However, in contrast to when computing the gradient with respect to motion, the gradient with respect to fibreorientation is included even if DTI is available for only a single cardiac phase.[000135] Equation 65 includes a gradient with respect to the fibre parameters,^^^(^^,^)^^^,^. Computing thisgradient shares the same structure as for the equivalent term with respect to the motion parameters (^^^^(^^,^)^^^as shown in Equations 32, 33, and 34 above), and uses the angular velocity tensor ^^::Equation 66 Equation 67 Equation 68again, usingEquation 69 Equation 70 Equation 71where from Equation= ^(^^, ^)^^^(^^)^^^,^. Sincehas been determined above (Equations 62, 63,and 64), the gradient of the DTI with respect to the fibre orientation parameters is obtained.[000136] The second term for which a gradient with respect to fibre orientation parameters is calculatedis the fibre–motion coupling term (−log ^(Motion | Fibres) shown in Equation 26. Calculating the gradient ofthe fibre–motion coupling term with respect to the fibre orientation parameters gives:^ log ^(Motion | Fibres)− ^^^,^Equation 72Equation 72 is again dependent on^^^(^^)^^^,^, and is computed as shown in Equations 62, 63, and 64.[000137] The third term in the log-probability function for which the gradient with respect to fibreparameters is calculated is the population-based fibre information 38 term (term 4 in Equation 26). The genericfibre parameters {^^} are replaced by the particular ones, {^^,^} :− log ^(FibresEquation 73where denotes the precision matrix as in Equation 8 but expressed in the particularCalculating the gradient of this term:^ log ^(Fibres | Source )−^Equation 74 Gradient Descent Algorithm[000138] The gradients outlined above describe the two gradients required for computing the gradientof the log-probability function: the gradient with respect to motion and the gradient with respect to the cardiac fibres: ^log ^(Fibres, Motion | I, DTI, Source ) ^ log ^(Fibres, Motion | I, DTI, Sour )−^^^,and − ce^. ^^^^^,^Equation 75 Equation 76[000139] Using these two gradients, any gradient-based optimisation algorithm can be used to optimisethe function. Examples of gradient-based optimisation algorithms include gradient descent with constant learning rate, conjugate gradient descent, AdaGrad, stochastic gradient descent, or any of the quasi-Newtonmethods such as the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. A flowchart illustrating anexample optimisation algorithm 700 is shown in Figure 7.[000140] Firstly, the motion parameters are initialised at Step 710, using a null value ^^^ = 0,representing a trivial motion, ^(^, ^^ , ^^) = ^. It should be appreciated that other initialisation values may beconsidered if relevant information is available that allows a better initial guess. For instance, if values obtainedfrom a statistical atlas of heart motions are available. The fibre orientation parameters ^^,^ are then alsoinitialised, at Step 720, with a best estimate which maximises ^(Fibres | Source^). In an embodiment where aGaussian prior fibre distribution for a population is used as the input, the best prior estimate is the mean of thedistribution, ^ ^^^ ^,^ .[000141] The optimisation algorithm 40 is then applied, alternating for the optimisation of the motionand the fibre parameters. The fibre parameters are frozen in their current value, and the gradient descent isapplied to optimise the motion parameters at Step 730. Then, the motion parameters are frozen at their currentvalue, and the gradient descent is applied to optimise the fibre parameters at Step 740. The algorithm 40determines, at Step 750, if a convergence threshold is met, and if so, the algorithm 40 outputs the optimalparameters for the fibre architecture. If the convergence threshold is not met, Steps 730 and 740 are repeated,alternately applying the algorithm 40 to the motion and fibre parameters until the criterion of convergence issatisfied. An alternative approach to achieving convergence could be developed involving joint optimisation ofthe motion and the fibre parameters.[000142] In some embodiments, a gradient descent with golden section line search may be used as theoptimisation algorithm 40. The line search aims at optimising the learning rate for each step, and is a veryefficient and thus commonly used line search algorithm.[000143] Once the gradient-descent algorithm 40 is deemed to have minimised the log-probabilityfunction (which may be determined, for example, by a pre-determined threshold), the resultant fibre and motion parameters are used to generate a model of the cardiac fibre architecture. The resultant parameters are the optimised parameters, that is, the parameters that best describe the motion and the fibre architecture, satisfying their coupling. The optimised parameters are compatible with the motion observed with the input cine-MRI images and, if available, the subject-specific DTI. The parameterised motion and fibre architecture can then be used to provide a visualisation of cardiac motion, including fibre and sheetlet reorientation with motion.[000144] The resultant model may be used for many applications in the field of research anddevelopment. For example, the generated heart model may be used to determine the influence of implanteddevices on the myocardial fibres, or for testing the effects of drugs on the heart. A specific example of anapplication for the heart model generated using the method according to the present invention is illustrated inFigure 8. As shown, the heart model is input to a simulation process 80, which also receives a test input. Anexample of a simulation process 80 may be, for example, a simulation for examining impacts and likelihood of success of proposed reconstructive procedures in a diseased heart. The corresponding test inputs may beproposed reconstructive procedures for the heart. In response to receiving a heart model for a patient and atest input, the simulation outputs a result, which, in the example described, would relate to whether theproposed reconstructive procedure is likely to be successful. Using a more realistic heart model for a patientprovides more accurate results, and allows clinicians to make better choices when determining thereconstructive strategies to use for a patient, thus reducing the risk to patients. For instance, a more accurate motion and fibre estimate allows obtaining better estimates of myocardial local strains, which can be of clinical significance, in particular when pathological or anomalous local contractions are present, such in the case of infarcted tissue. In a related application, the generated optimum parameters for fibre architecture may be usedto monitor remodelling of myocardial fibres in a subject after, for example, a myocardial infarction. As moreheart models are generated in this manner a heart model that reflects a subset of the population (for example,a population over a certain age or with a certain disease) rather than an individual can be developed using a statistical model of the population. The heart model representing the subset of the population can again be used in research to more accurately test the effect of implantation devices or drugs on the particular subset of the population.[000145] Another potential application for the heart model generated using the method according to thepresent invention is to combine the heart model with methods that model cardiac muscle function. For example,the heart model may be used for method described in WO2022 / 106580, which describes a heart simulationmethod that incorporates muscular electrophysiology and electromechanical properties of the heart. Using theheart model of the present invention, which includes information on cardiac motion and fibre structure,alongside the method described in WO2022 / 106580, results in a better, more realistic simulation on which totest the effect of, for example, drugs.[000146] This method may also be applied to muscular structures, in particular bodily organs whichexperience motion that cannot be consciously controlled by a subject. All muscles include myocytes that are oriented according to some fibre architecture, and this fibre architecture will substantially influence their motion.A particularly relevant bodily organ is the uterus – an area of particular interest is contraction of the uterusduring birth or during a menstrual period. The same approach of capturing a plurality of images during a muscular motion cycle can be used, and initial models of fibre architecture and of mechanical coupling optimised according to the method to determine a model of fibre architecture consistent with the motion indicated in the plurality of images. This can be used, for example, to study the relationship between uterus motion and inflammation or pain (dysmenorrhea).[000147] While the primary use of this method described above is in providing a better model of thefibre architecture, a further use is to improve imaging – in the case of DTI image capture, DTI imaging. Oneexemplary case of this in the context of DTI cardiac imaging is described below.[000148] In the conventional case, the presence of multiple cardiac phases can be problematic forestablishing high resolution fibre images, as there is not an effective way to combine information from imagesobtained in different cardiac phases. In order to obtain higher resolution, it would be normal to devise an acquisition sequence in which all image slices were obtained in the same cardiac phase. That wouldnecessarily lead to an extended period of time being required for DTI image capture. Moreover, in addition tomaking image capture more time consuming, such an approach may introduce additional complexity if thereare progressive changes in cardiac behaviour over multiple cardiac cycles. However, using a method as described here, cardiac motion can be described from images obtained from multiple cardiac phases, and this information can be used to enhance their resolution even though the images are from different cardiac phases. As a model of cardiac motion is available, by applying the inverse motion to all DTI images at different cardiac phases, these images can be explicitly fused into a single 3D DTI image of higher resolution. In effect, the ability to transform the images for the determined motion enable the type of fusion and super resolution strategies customarily applied to scalar images, such as cine-MRI, to tensorial DTI images.

Claims

CLAIMS 1. A method of estimating muscle fibre architecture for a muscular structure, the muscle fibrearchitecture comprising muscle fibre orientation, the method comprising:receiving a plurality of images that represent the motion of the muscular structure during a motion sequence; receiving an initial model of muscle fibre architecture, the initial model comprising a genericestimate of muscle fibre orientation; receiving a model of mechanical coupling between muscular structure motion and muscle fibres; generating a joint function using muscular structure motion as represented in the plurality ofimages, the initial model of muscle fibre architecture, and the model of mechanical coupling; applying an optimisation procedure to optimise the joint function; andfrom the optimisation procedure, determining a model of muscle fibre architecture consistentwith the muscular structure motion indicated in the plurality of images.

2. The method of Claim 1, wherein the optimisation procedure comprises:parameterising the muscular structure motion indicated in the plurality of images to determine a set of motion parameters; parameterising the fibre architecture indicated in the initial model of muscle fibre architectureto determine a set of fibre architecture parameters; andoptimising the joint function by alternatively optimising the set of motion parameters and the set of fibre architecture parameters until convergence is reached.

3. The method of Claim 2, wherein optimising the set of motion parameters and the set of fibrearchitecture parameters comprises: calculating the gradient of the joint function with respect to the set of motion parameters andthe set of fibre architecture parameters; andapplying a gradient – based optimisation algorithm to the joint function.

4. The method of Claim 2 or 3, wherein determining the model of muscle fibre architecturecomprises using the converged set of motion parameters and set of fibre architectureparameters.

5. The method of any of Claims 2 to 4, wherein the muscular structure motion is parameterisedusing a continuous time-varying velocity field.

6. The method of any of Claims 2 to 5, wherein parametrising the muscle fibre architecturecomprises parameterising the local fibre orientation and parameterising the spatial variation infibre orientation.

7. The method of any of Claims 1 to 6, wherein the function describing motion, the initial model ofmuscle fibre architecture and the model of mechanical coupling are probabilistic models.

8. The method of any of Claims 1 to 7, wherein the method additionally comprises generating aconditional probability density function to represent muscular structure motion as indicated in theplurality of images, for the initial model of muscle fibre architecture and for the model ofmechanical coupling.

9. The method of Claim 8, wherein the joint function is generated using the probability densityfunctions.

10. The method of any of Claims 1 to 9, wherein the plurality of images are received from a cine-MRI.

11. The method of any of Claims 1 to 10, wherein the initial model of muscle fibre architecture ispopulation-based fibre information or fibre information provided from rule-based models.

12. The method of any of Claims 1 to 11, wherein the model of mechanical coupling is derived fromstrain properties of muscle fibres.

13. The method of any of Claims 1 to 12, wherein the method additionally comprises:receiving muscle fibre information for an individual;generating a function describing the muscle fibre information for the individual; andgenerating the joint function using muscular structure motion as represented in the pluralityof images, the initial model of muscle fibre architecture, the model of mechanical coupling, andthe function describing the muscle fibre information for the individual.

14. The method of Claim 13, wherein the muscle fibre information for the individual comprises atleast one diffuse tensor image, and the function describing the myocardial fibre information for the individual is a conditional probability density function.

15. The method of any preceding claim, wherein the muscular structure is a heart, the muscle fibrearchitecture is a myocardial fibre architecture, and the motion sequence is a cardiac cycle.

16. A method of modelling function of a heart, comprising a method of estimating myocardial fibrearchitecture as claimed in claim 15.

17. The method of any of claims 1 to 14, wherein the muscular structure is a uterus.

18. A method of visualizing muscular motion, the method comprising estimating muscle fibrearchitecture for a muscular structure according to any of claims 1 to 17 and using the optimisedset of motion parameters and the set of fibre architecture parameters to provide the visualization.

19. A method of image enhancement in imaging a muscular structure, comprising:estimating muscle fibre architecture as claimed in any of claims 1 to 14; determining the muscular structure motion during the motion sequence; applying an inverse motion parameter to at least some of the plurality of images to provide images compensated for the muscular structure motion; and combining the compensated images using image fusion or image super resolution to provide at least one higher resolution image of the muscular structure.

20. A system for estimating muscle fibre architecture using a plurality of images indicating motion ofa muscular structure, the muscle fibre architecture comprising muscle fibre orientation, thesystem comprising: an input engine for receiving each of the plurality of images indicating muscular structuremotion; afunction creator configured to generate a joint function using muscular structure motion asrepresented in the plurality of images, an initial model of muscular fibre architecture, the initial model comprising a generic estimate of muscle fibre orientation, and a model of mechanical couplingbetween muscular structure motion and muscular fibres;an optimisation engine configured to optimise the joint function; and a model generator configured to determine, using the optimised joint function, a model of muscular fibre architecture consistent with the muscular structure motion indicated in the plurality of images.

Images