NL2026924B1 - Spectro-dynamic magnetic resonance imaging - Google Patents

Abstract

The present patent disclosure relates to a method and a device for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, the method comprising performing optimization With an objective function based on a measurement model and a dynamical model for the time dependent displacements of the moving object and the dynamical parameters in order to determine an optimized or final set of the time dependent displacements of the moving object and the dynamical parameters.

Description

SPECTRO-DYNAMIC MAGNETIC RESONANCE IMAGING The present patent disclosure relates to a method and a device for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object, and a computer program product for performing the method.

Magnetic resonance imaging (MRI) is an imaging modality used for many applications and with many sequence parameters that can be tuned and many imaging parameters that can be IO observed to extract e.g. different kinds of biological information.

The human body is a dynamical system. Its performance is not only determined by its static anatomy, but also by its functioning in a dynamic setting. One specific example of a dynamic organ is the heart. More than 26 million people suffer from heart failure worldwide, and its prevalence is increasing. Left ventricular ejection fraction, measured using echocardiography, is the most widely used parameter for the diagnosis of heart failare. Myocardial strain can give important additional information in the prognosis of several cardiac diseases.

Another example of a biomechanical system is the musculoskeletal system. Conditions related to muscles and joints can cause long-term pain, and a large portion of the population suffers from musculoskeletal conditions.

Magnetic resonance imaging (MRI) is a medical imaging modality that provides excellent soft tissue contrast, making it an effective tool for imaging organs such as the heart or other muscles. Furthermore, MRI can obtain 3D information and does not involve ionizing radiation. However, MRI in its current form is a relatively slow modality, complicating the analysis of dynamics of organs and other tissue.

To obtain MR images, the magnetisation of the subject is sampled in the spectral domain, called k-space. Typically, the k-space is sampled line by line, where the acquisition of one line is called a readout. Image reconstruction by Fourier transformation requires a fully sampled k-space, for which many readouts and thus long acquisition times are needed (depending on the wanted resolution). When k-space is fully sampled, images can be reconstructed without artefacts, such as Nyquist artefacts. The opposite is true for undersampled data. In view of the above, when adding a third spatial or a temporal dimension, inferring time-resolved 3D dynamic information using MRI is very challenging.

Several techniques exist to accelerate MRI acquisitions for imaging dynamic systems such as the heart. In other words, it is attempted to reconstruct images at a high temporal resolution, and to derive the dynamical information from those reconstructed images. An example is to use so- called real-time imaging, a technique that acquires all data in a single heartbeat. However, the resulting images suffer from reduced image quality, and acquiring 3D images is not possible with this technique. Another approach is to use so-called cine imaging, where the acquisition is gated such that all data is acquired during a specific phase of the cardiac motion. This requires that the motion is periodic, which is problematic for the frequent aperiodic dynamics occurring within patients, such as patients with arrhythmia. Furthermore, using cine imaging to obtain respiratory motion requires additional cumbersome motion correction strategies.

An alternative to obtaining the dynamics from reconstructed images obtained at high temporal resolution, is to estimate dynamics directly from spectral (spatial frequency, k-space) data, removing the requirement of a fully sampled k-space at every point in time. A recently IO developed method called MR-MOTUS has shown that it is possible to reconstruct motion fields directly from undersampled k-space data. While this method does provide motion fields at a high temporal resolution, MR-MOTUS is limited, among other limitations, by the need for a reference image from a separate MR-scan, e.g. during breath hold by the patient.

It is an object, among objects, of the present patent disclosure to provide improved methods, devices and systems for obtaining dynamical information related to an object, e.g. a patient or an organ of the patient, using MR.

According to a first aspect, there is provided a computer implemented method for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object, the method comprising: providing a measurement model defining a relation between the time dependent displacements of the moving object in the spatial frequency domain and the measured SFDMR signal in the spatial frequency domain; providing a dynamical model dependent on a set of dynamical parameters modelling dynamics of the moving object in the spatial frequency domain in relationship to the time dependent displacements of the moving object in the spatial frequency domain; and performing optimization with an objective function based on the measurement model and the dynamical model for the time dependent displacements of the moving object and the predetermined parameter in order to determine an optimized or final set of the time dependent displacements of the moving object and the predetermined parameter.

The above method, which is also referred to as “spectro-dynamic MRI” is thus a method to probe time-resolved dynamic information using MRI data. Two models are combined to identify the dynamics of a system. The first model is the measurement model, providing a relation between the displacements of a moving object and the measured magnetization in the spectral domain. This model has the effect of acting as a data consistency term and can be used individually to reconstruct the dynamics when fully sampled data is available.

The provided dynamical model, when used together with the measurement model, allows the use of undersampled data. The dynamical model may include prior knowledge about the dynamics, Together with the measurement model, the displacements or motion fields can be reconstructed at a high spatial and temporal resolution, since undersampled data can be used.

Additionally, the dynamical model allows one to obtain the dynamical parameters, which can give information about the material properties of the scanned object at various points or locations thereof. The material properties, in turn, provide a way to determine or deduce tissue parameters such as myocardial strain referred to above.

The method allows to infer time-resolved dynamic information of objects such as a heart which would provide valuable information for the diagnosis and understanding of diseases such as arrhythmia and heart failure. Besides a high temporal resolution, the resulting set of final or optimized data also has a high spatial resolution. This makes the method highly suitable for time- resolved dynamic imaging, e.g. for investigating the loads and motion of joints.

The term “displacements” of the moving object, in an embodiment, is position information of one or more points in the spatial domain, or in the spatial frequency domain, associated with the object in reference to an origin or to a previous position, or denoted as generalized coordinates. Alternatives for displacements are movement information, motion information, displacement information, position information, and time dependent positions.

The measurement model and the dynamical model preferably both comprise at least one derivative of the time dependent displacements of the moving object. This interrelatedness of the two models allows for a faster solving of the optimization problem to be solved. The derivative can be any order, preferably first or second order. The orders are not required to be the same in the two models.

It will be understood that the providing of the dynamical model can more generally be implemented as providing a second model defining a predetermined relationship between a predetermined parameter in the spatial frequency domain and the time dependent displacements of the moving object in the spatial frequency domain.

Preferably the magnetic resonance excitation of the moving object is according to an applied pulse sequence.

In an embodiment, the dynamical model is (or is implemented as or comprises) a first differential equation defining a relation between the dynamical parameters, the displacements of the moving object and one or more temporal derivatives of the displacements. Additionally or alternatively, the relation defined by the differential equation comprises spatial derivatives of the displacements. Because in the method the dynamical model and the displacements are all defined in the spatial frequency domain, the signal can comprise small time steps for each value of k. Due to these small time steps, numerical differentiation methods such as finite differences can be used to approximate the temporal derivatives in the dynamical model.

In an embodiment, the measurement model is (or is implemented as or comprises) a second differential equation defining a relation between the measured SFDMR signal, the displacements of the moving object and one or more temporal derivatives of the displacements, wherein preferably the relation defined by the second differential equation further comprises temporal derivatives of the SFDMR signal. The same advantage as for the use of temporal derivatives in the dynamical model, in particular of using finite differences, applies here. In particular the use of temporal derivatives of the displacements in both the dynamical model and the measurement IO model is advantageous, as this allows for an effective route for solving the optimization problem.

The temporal derivatives may comprise one or more of a first order derivative, a second order derivative and a third derivative.

In an embodiment, the temporal derivatives of the displacements of the moving object are approximated using numerical differentiation, wherein preferably the temporal derivatives are {5 approximated from a plurality of finite differences in the SFDMR signal measured for a plurality of times. This allows for an efficient way of performing the optimization.

In an embodiment, the SFDMR signal comprises a plurality of signal values measured at a same spatial frequency and respective times, wherein a time difference between consecutive ones of the signal values measured at the same spatial frequency is smaller than a threshold time, or predetermined time period.

Preferably the threshold time is between 1 and 100 repetition times, more preferably between 1 and 50 repetition times (TR). In modern MRI systems, the repetition time can be in the order of 2-5 ms. When imaging, for instance, movement of lungs, a 100 repetition times would still be short enough to obtain useful displacements/displacement information (dynamical information) of the lungs when performing the method, in particular when using finite differences to approximate temporal derivatives.

In an embodiment, the SFDMR signal is measured for a plurality of spatial frequencies.

In an embodiment, the SFDMR signal is an undersampled SFDMR signal. The SFDMR signal may have been obtained using a sampling scheme or pattern wherein a certain spatial frequency is sampled at least twice within a predetermined time period, e.g. the time period described above. The sampling pattern may be such that the same spatial frequency (readout line, or k-point) is measured a plurality of times before different spatial frequencies are sampled. Alternatively, or more specifically, the sampling pattern may be an interleaved sampling pattern.

The method may comprise measuring the SFDMR signal for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object, preferably according to a sampling pattern, such as the ones described above.

In an embodiment, the dynamical parameters comprise, or are associated with, one or more of a mass matrix, a stiffness matrix and a damping matrix associated with the moving object. Dynamical models based on one or more of these parameters are available for many moving objects or organs, and can thus be used with no or relatively minor modification.

5 In an embodiment, the displacements of the moving object are represented by a displacement field describing the displacements of the moving object in the spatial frequency domain for each time of the plurality of times and for each of the plurality of spatial frequencies; or a generalized coordinate vector comprising generalized coordinates associated with the moving object for each time of the plurality of times.

It is noted that the generalized coordinate vector can be obtained directly from the displacement field and vice versa, using suitable basis functions, e.g. a set of spatially varying basis functions. The use of the generalized coordinate vector reduces the dimensionality of the optimization problem.

In general, a displacement field may comprise an assignment of displacement vectors for all points in a region or body (e.g. of or associated with the object) that is displaced from one state to another. In general, a displacement vector could specify the position of a point or a particle (e.g. of or associated with the object) in reference to an origin or to a previous position. For example, a displacement field may be used to describe the effects of deformation on a (solid) body such as the object.

In an embodiment, the dynamical model is implemented or defined as a differential equation comprising one or more of: d2 Mp(8) zz a(t), Lp(9) 3; q(t). and Kp(8)q(t), wherein Mp, is the mass matrix, Lp is the damping matrix, Kp is the damping matrix, q is the generalized coordinate vector, and t is the time, Preferably, the dynamical model is implemented as: d? d Mp (8) —=a() + Cp (8) 7790 + Kn(©)a(t) +7 = 0 wherein Z is a force term. Z can be non-zero or zero, depending on whether non- conservative and/or external forces are present or not, respectively. Also, Mp, Cp, and Kp, can each be zero; as long as one of these matrices is non-zero a preferred dynamical model of the object is described.

In an embodiment, the measurement model is defined as:

0 | d | 0 GM, UK t)= M+ 2m) IM * = Uj + >} M x Ie; == Uj =0

J J wherein M is the SFDMR signal in the spatial frequency domain, k a vector comprising spatial frequencies for each datapoint of the measured SFDMR signal, t the time, j the spatial dimension, k; the spatial frequency along the j-th dimension, and U, a displacement field describing displacements of the moving object in the spatial frequency domain for spatial dimension j. It is noted that *“*” in the present patent disclosure indicates a convolution operation.

In an embodiment, the step of performing optimization comprises solving an optimization problem based on the measurement model and the dynamical model, wherein the dynamical model is used as a regularisation term or as a constraint on the measurement model. Alternatively, the measurement model can be used as a constraint on the dynamical model. In both cases, the constraint may be an equality constraint or an inequality constraint.

In an embodiment, the optimization problem is defined as the minimization, for the displacements of the moving object and the dynamical parameters, of a sum of the norm of the measurement model and the dynamical model multiplied by a trade-off parameter that changes an amount of regularization by the dynamical model on the measurement model.

In an embodiment, the optimization problem is defined as: ple (Mk, 6), UK 6), Kk, OI + AF(Uk ©), 0, k, t)[|¢ wherein: G is the measurement model, F is the dynamical model, M is a vector or matrix comprising the measured SFDMR signal in the spatial frequency domain, U is a displacement field describing displacements of the moving object in the spatial frequency domain, k is a vector comprising spatial frequencies for each datapoint of the measured SFDMR signal, t is the time, © is a vector comprising the dynamical parameters, A is the trade-off-parameter being a scalar parameter that changes the amount of regularization of the dynamical model on the measurement model, and a is the power to which results of the respective norms (||..||) are raised, wherein preferably ais 1 or 2. The norms may be the Euclidian norm or any other suitable norm. It is noted that @ is position (or spatial frequency) dependent and/or time dependent.

In an embodiment, the optimization problem is defined as: minflG(M(k, 6), Uk, £), Kk, IS such that [F(UK, 1), 0,K ||? =O Here, ||F|| = 0 is an equality constraint. Alternatively, [|F|| can be used as an inequality constraint such as ||F|| < € or ||F|| < e. Alternatively, the minimization can be done for F with equality or inequality constraints on G.

In an embodiment, the optimization problem is defined as: minfiG (Mk, £), q(t). KOI + AIF (q(t), 8, kK, OI wherein the q is the generalized coordinate vector. Also here the norms may be the Euclidian norm or any other suitable norm.

According to a second aspect, there is provided a device or (computer) system for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, the device comprising one or more processors configured to: provide a measurement model defining a relation between the time dependent displacements of the moving object in the spatial frequency domain and the measured SFDMR signal in the spatial frequency domain; provide a dynamical model dependent on a set of dynamical parameters modelling dynamics of the moving object in the spatial frequency domain in relationship to the time dependent displacements of the moving object in the spatial frequency domain; and perform optimization with an objective function based on the measurement model and the dynamical model for the time dependent displacements of the moving object and the dynamical parameters in order to determine an optimized or final set of time dependent displacements of the moving object and the dynamical parameters.

It will be apparent that any advantage relating to the herein described methods) is readily applicable to the present device/system. Also any embodiments, options, alternatives, etc., described above for method are applicable to the present device/system.

According to a third aspect, there is provided, a computer program product comprising computer-executable instructions for performing the method described above and below, when the progran is run on a computer.

The accompanying drawings are used to illustrate presently preferred non-limiting exemplary embodiments of (the) method(s) and (the) device(s)/system(s) of the present disclosure.

The above and other advantages of the features and objects of the disclosure will become more apparent and the aspects and embodiments will be better understood from the following detailed description when read in conjunction with the accompanying drawings, in which: Fig. la is a schematic and simplified plot of a sampling pattern used in conventional MRI; Fig. 1b is a schematic and simplified plot of an example sampling pattern used in the presently disclosed embodiments; Fig, lc is a schematic and simplified plot of another, interleaved, example sampling pattern used in the presently disclosed embodiments; Fig. 2 shows schematic drawings of four different example dynamical systems of which example movement information is obtained using embodiments of the present disclosure, of which Fig. 2 a) shows a simple pendulum, b) coupled pendula, c) a spherical pendulum and (d) coupled spherical pendula;

Fig. 3 shows photos of the used experimental setup from three different view points; Fig. 4 shows plots of L-curves of root mean square errors obtained respectively for the four example dynamical systems of Fig. 2; Fig. 5 shows four plots of 1D FFT results for all readout lines in measurements of the four example dynamical systems of Fig. 2; Fig. 6 shows four plots of estimates displacements versus time obtained from measurements of the four example dynamical systems of Fig. 2; Fig. 7 shows four plots of estimated natural frequencies obtained from measurements of the four example dynamical systems of Fig. 2; Fig. 8 shows two plots of estimated displacements for the undersampled datasets, wherein plot a) shows results for a simple pendulum of Fig. 2{a) and plot b) shows results for a spherical pendulum of Fig. 2(c); and Fig. 9 shows a device according to an embodiment of the present patent disclosure. Theory Signal Acquisition In an MRI experiment, the signal intensity m(r, t) at spatial location r € R and time ¢ is proportional to the local magnetization of the object. The number of spatial dimensions is indicated by d. By applying linear gradient magnetic fields to the static magnetic field inside the scanner, this signal can be sampled in the spectral domain at the spatial frequencies k € R2. The spatial and spectral domains are related through the Fourier transform: MK, t) = F(m(r,t)) = [aa m(r,t)e” "8" dr (1) mr, t) = FHM 1) = [pa MK, )e?™ 7 dk (2) Note that variables in the spatial domain are indicated with lower case symbols, while upper case symbols are used for the spectral domain. Vectors are indicated in bold font. In conventional MRI, the measured data in the spectral domain M(K, t) is reconstructed to an image m(r, t) by using the inverse discrete Fourier transform, thereby approximating (2). Since the signal contribution of one specific voxel affects the whole k-space, and each sample in k-space contains information about the complete field of view in the image domain, the k-space must be fully sampled to prevent undersampling artifacts in the reconstructed image.

Fig. 1a shows sampling pattern as used in conventional imaging. Each black square indicates one readout, separated in time by the repetition time (TR). The vertical axis represents the different spatial frequencies that need to be sampled to get a fully sampled k-space (the phase- encoding steps for a Cartesian acquisition). With every readout, a different line in k-space is sampled. Thus, after N readouts (N = 6 in this schematic illustration), an image can be reconstructed. Typically, N is in the order of several thousand repetitions for a 3D acquisition. Fig. Ib shows a sampling pattern as can be used in spectro-dynamic MRI, i.e. in accordance with the present patent disclosure. The same readout line is repeated Nrep times (3 in this case) before different k-space frequencies are sampled. This greatly reduces the time difference At between the acquisition of the same data point in k-space (from N - TR to 1 - TR), thereby allowing the computation of temporal derivatives through finite differences. Figure lc shows another sampling pattern as can be used used in spectro-dynamic MRI, i.e. in accordance with the present patent disclosure. The same readout line is repeated Nrep times (3 also in this case), but in an interleaved manner, The time difference At between the acquisition of the same data point in k-space is reduced from from N - TR to 2 - TR in this interleaving example. Also with the sampling pattern of Fig. lc, the computation of temporal derivatives through finite differences is possible.

A goal of spectro-dynamic MRI is to infer the dynamics of a system as well as identify its dynamical parameters. To do this, the dynamics are modelled with spatial and temporal derivatives. As will be explained below, these temporal derivatives can be estimated using finite differences.

If one would use the measured data in the spatial domain, this would require the image values m(r, t) and m(r, t + At) to be available at some spatial location r. Since these images must first be reconstructed, this also means that two fully sampled k-spaces at these time points need to be available. Therefore, At cannot be smaller than the time it takes to sample a complete k-space (Fig. 1(a)), which for 3D acquisitions is in the order of seconds. This low temporal resolution leads to motion artifacts in the reconstruction.

Instead, here the data in the spectral domain is used directly. To estimate temporal derivatives using finite differences, measurements M(K, t) and M(K, t + At) in k-space are used. Since the intermediate reconstruction step is no longer required, At can be as short as the repetition time TR, which for a typical MRI system can be in the order of milliseconds (Fig. 1{b)). Thus, the temporal resolution of the measurements is improved by using (highly undersampled) k-spaces at successive points in time.

The dynamics of a moving and possibly deforming object can be described by a displacement field u(r, t) € R%, with its Fourier transform U(K, t). An embodiment of the measurement model G, which will be derived below, will be used to describe the relation between this displacement field and the measurements M(K, £) in the spectral domain: GMk t), Uk 6), K) =0 (3)

When the data is fully sampled, this measurement model is sufficient to reconstruct the unknown displacement field.

However, fully-sampled k-space data cannot be acquired for a typical moving object, thus the measured data will be undersampled.

To compensate for this, additional information about the system under investigation is provided.

The measurement model G will be 19 used in combination with a dynamical model F which provides a priori knowledge about the dynamics of the system, This dynamical model also allows for the reconstruction of a set of dynamical parameters, yielding more information about the system’s dynamical properties.

The dynamics of a system can be modelled with a differential equation, such as a partial differential equation (PDE). In general, an embodiment of the differential equation would be a function of the displacement field, its spatial and temporal derivatives, and a set of dynamical parameters 8(T), which can be spatially dependent.

Like the measurement model, this dynamical model can also be converted to the spectral domain, resulting in: F(U(K. t), 8(K),k) = 0(4)

The specific form for (4) depends on the dynamical system that is being investigated.

See below for more details about the dynamical systems used to show the methods provided herein give accurate results.

Given the provided measurement model G and the dynamical model F, the dynamics as described by the displacement field and the dynamical parameters can be reconstructed.

One option is to minimize the L2 norm of the error of one model, while using the other model as an equality constraint.

However, it is important to note that both (3) and (4) will be subject to imperfections.

Measurement noise will introduce errors in the measurement model, while model imperfections will result in a difference between the observed dynamics and the modelled dynamics.

Therefore, a single optimization problem will be solved, where we use the dynamical model as a penalty term (or regularizer) on the measurement model: min IGM UKDI+A1FUOKDHIE 6)

il In this embodiment, 4 is introduced as a trade-off parameter that changes the amount of regularization of the dynamical model on the measurement model.

If 4 = 0, only the measurement model is used.

This will allow for the reconstruction of the displacement field in case of a fully sampled acquisition, but it will be more sensitive to noise, and no dynamical parameters can be estimated.

Increasing 4 will increase the robustness to measurement noise, but it will create a bias towards the dynamics included in the dynamical model.

Ultimately, when A goes to infinity, the reconstructed solution is equivalent to the one obtained when the dynamical model is used as an equality constraint.

In this case, any model imperfections in the dynamical model will no longer be corrected, and any dynamics not included in the dynamical model will no longer be reconstructed.

Therefore, a suitable value for 4 is found which balances the residuals of both models, as is shown in further detail in relation to Fig. 4 below.

Measurement Model The measurement model describes how the measured magnetization changes due to the dynamics of the measured system.

If the magnetization is in steady-state, and the receive and excitation fields are homogeneous, it can be assumed that the total magnetization is conserved.

Therefore, the derivation of the measurement model is similar to the one used for the conservation of mass when modelling fluid mechanics.

We start by looking at the volume integral of the magnetization over an arbitrary volame element V, € R3. This control volume is fixed in space, while the object is moving.

The rate of change of the magnetization within V is equal to the inflow of magnetization across the surface Sy of this volume element: = [If mr, 6) dV-Hs, mr, V(r, t) :n dS = 0 (6)

Here v(T, t) is the velocity field, and n is the unit normal vector of the surface element.

Because our volume element is not moving, we can move the temporal derivative inside the integral: Jl, mae ty dv = fff, sme) dv (7) Furthermore, we can use the divergence theorem: s[m(r, Ov(r,t)] nds = fff, V-[m(r,)v(r,)] dV (8)

Combining (6), (7), and (8), we can write the measurement model in integral form:

J Ik, zm, 6 av + Jf, V- [m(r, )v(r t]dV = 0 (9) Since (9) must hold for any Vo, we can extract the integrand and rewrite our measurement model in the differential form: Zm(r, 6) + Vm(r,t) - v(r, £) + OV v(x, 9] = 0 (10) To couple this measurement model with the dynamical model, the displacements u(r, t) are introduced into our measurement model. To do this, it is assumed that the velocity field is the temporal derivative of the displacement field: a ; vir, t) = 5 u(r, t) (in) This assumption is valid especially when the displacement field u(r, £) is smooth across space {see below). Combining (10) and (11) gives a further detailed measurement model in the spatial domain: Zm(r,£) + Ym(r,6) Zu) + me, 6) [V- Sur, 0) = 0 (12) a ’ ’ at ! ’ a ’ = In MRI, the measurements are sampled in the spectral domain. If we would use (12) directly, we would have to reconstruct the image data m(r, £) from the measured k-space data M(K, t). The temporal resolution can be increased by directly using the spectral data, removing the need for an image reconstruction step. To do so, we must transform (12) to the spectral domain. First, we must transform the temporal and spatial derivatives. Given a function f(r, t) with corresponding Fourier transform F (K, t), the temporal derivative can be moved inside the Fourier integral: Zr) =2 [14 Fk De? kr dk att Vr at Rd ’ . a 270k = [na 5: F t)ei27KT gk (13) _p-1]8 =F SF 0)

The spatial derivative in the spectral domain can be written as a multiplication with the spatial frequency Kk: Vir t) =V {ea Fk, te?" dk = fa V[F(k, te'275T] dk = [za iZEKF (K, t)e 27" dk = Fl i2nkF (kK, t)]

Furthermore, we can use the convolution theorem of the Fourier transform: fart) grt) =F Fk A) + Gk t)] (15)

Using (13), (14), and (15), the measurement model as given in (12) can be transformed to the spectral domain, with a summation over the spatial dimensions: a . J . J oo SME +i2n Nh kM 0) UK) + 2n Sj MOK) + kj UK 6) =0 (16) Dynamical Model One way to describe the dynamics of systems/objects of which dynamical information is gathered, is as dynamics of systems comprising discrete particles and with n degrees of freedom.

For example, a single particle that is free to move in both the x- and y-direction has two degrees of freedom.

Systems with nonconservative or external forcing terms are excluded below, especially for sake of simplicity, but can be included if the dynamics of the object to be studied requires so.

The (linearized) equations of motion can be described by a second-order differential equation:

da? I Mp (8) ~= q(t) + Kp(0)q(t) = 0 (17)

Here, q(t) € R" are the generalized coordinates for each of the n degrees of freedom, Mp = MJ € R™" and Kp = K} € R™" are respectively the mass and stiffness matrix, which depend on the dynamical parameters 8. None of the quantities in (17) are dependent on r.

The spatial dependencies are implicitly included in the model, as each generalized coordinate corresponds to one part of the moving system, with its own location in space.

It is preferred that a coordinate system is chosen such that q = © is the equilibrium position.

The generalized coordinates are connected to the displacement field u(r, t) through a set of basis functions, as described below.

Note that despite the names of Mp and Kp, this model is not only applicable to mass-spring systems. Any nonlinear system, such as a swinging pendulum, but also an organ such as a heart, can also be described by (17), in particular when the displacement from equilibrium is small. Furthermore, (17) can directly be used in the spectral domain, since neither r nor K appear explicitly.

The preferred dynamical model is premultiplied with (Mp (8))72. This results in a better scaling of the parameters and prevents the optimization from converging towards parameters for which the mass matrix was very small. Using the matrix 05 (8) = (Mp(8)) 1Kp (8), whose eigenvalues are the squared natural frequencies of the system, the dynamical model can be written IO as: a? za) + 0p (@)q() =0 (18) Basis Functions The displacement field u(r, t) is defined at many spatial locations (of the measured volume or of the object), while the dynamical systems only have a limited number of degrees of freedom. Therefore, we can reduce the dimensionality of the optimization problem (5) by using a set of n spatially varying basis functions. We assume that, by using the right choice of basis functions, a relation between the displacement field u(r, t) and the i, generalized coordinate gq; (t) can be defined by the iy, basis function ¢;(r) € R%: u(r, t) = Xi digit) (19) Next, (19) can be converted to the spectral domain since the Fourier transform is a linear transformation: Uk) = IL, Oda) 20) Discretization As described above, the k-space can be sampled at a high temporal resolution. This allows for the approximation of the temporal derivatives in (16) and (18) using finite differences. Combining (16) and (20), the first temporal derivative using finite differences is defined as: 2 a Mti) Mti) at MK tn) ~ vy (21)

0 _ ities) 7 dili) ZU t) ~ Tk, Di) Lalli) (22) It is noted that other finite difference schemes or more generally numerical differentiation methods can alternatively be used. Applying finite differences to the second temporal derivative in (18) gives: a» Gtr) 29) alt) 0 Eat)» (23) If we have measurements from p different timepoints, and n degrees of freedom, we can concatenate all vectors q(t;) fori = 1 ...p to form one vector q € R™ which contains all generalized coordinates at every timepoint. By combining (16), (21), and (22), we can rewrite the measurement model as a linear expression in q. Therefore, we can also write our measurement model as a system of linear equations: Aq b=0 (24) A € CNPXxTP and b € CNP contain the parts of the measurement model that are respectively linear in and independent of the degrees of freedom q, and N is the number of sample points per readout.

Similarly, by combining (18) and (23), the dynamical model can be written using the dynamical model matrix C (8) € R™*"P_ which depends on the dynamical parameters 8: C(8)q=0 (25) Replacing G and F by (24) and (25) respectively, the optimization problem (5) can be rewritten as: min | Aq—b 2+ 211 C(0)q lI (26) Example Dynamical Systems The above shown and derived embodiments of the spectro-dynamical models are shown to be effective by several experiments with the dynamical systems of Fig. 2. Different dynamical systems were used, with increasing degrees of complexity. Fig. 2 shows schematic overviews of four different dynamical systems: (a) simple pendulum, (b) coupled pendula, (c) spherical pendulum, and (d) coupled spherical pendula. All pendula have length I and mass m, my, or mj. They are subject to a gravitational force in the vertical direction, with gravitational acceleration g. The position of each pendulum is described by two orthogonal coordinates x and y, both in the horizontal plane. Systems (b) and (d) also contain a spring, attached to the two pendula at a distance lj, from the pivot point, with spring constant k.

These example dynamical models describing spring systems are implementations of F (see eg. 4 and eq. 17). Dynamical models of organs exist that can be used with the herein described methods. For instance, Maxime Sermesant, in “When Cardiac Biophysics Meets Groupwise Statistics: Complementary Modelling, Approaches for Patient-Specific Medicine”, Signal and Image processing, Université de Nice — Sophia, Antipolis, 2016, the contents of which are incorporated herein by reference in its entirety, describes various dynamical models of the heart in section 3.5.1, in particular as equation 3.2 (page 59). These dynamical models are similar to eq. (17) and therefore the dynamical models herein can be adapted according to eq. 3.2 of Sermesant.

For instance, eq. 17 can be adapted, or can be replaced, by the dynamical model 3.2 of Sermesant. For instance, Y as used by Sermesant is equivalent to the displacement field u in the spatial domain. Equivalents of M, K, and L as used above and below are also provided as M, K, and C respectively. The term Z described above and below is equivalent to any one or more of Fp, Fc, Fg, and Fj given in eq. 3.2 of Sermesant.

Simple Pendulum To start, a simple pendulum with one degree of freedom (Fig. 2(a)) was considered. Using the angle ¢ the pendulum makes with respect to the vertical position, the differential equation of this system is: mL? + mgsing = 0 (27) The dynamics of this system can be linearized around the stable equilibrium ¢ = 0. This small-angle approximation results in a differential equation analogous to the one of a mass connected to a spring. Written in the format of (17), this results in: n =1 = [x Mo) = 0 (28) Kp(8) =1[9]

Since the gravitational acceleration g is fixed to 9.81 m/s’. the only remaining unknown dynamical parameter is the length of the pendulum : 6 =] (29)

Note that the polar coordinate ¢ has been replaced by the Cartesian coordinate x.

Since there is only one degree of freedom, there is only a single basis function, which has a constant value of 1 over the entire spatial domain: oi) =1Vr (30) Coupled Pendula Two simple pendula, as described by (27), can be coupled together with a spring to form a slightly more complicated dynamical system (Fig. 2(b)). This coupled oscillator shows nonlinear dynamics.

However, we can again linearize these equations using the small-angle approximation and rewrite them in terms of the Cartesian coordinates x; and x, of both pendula.

The system then becomes equivalent to two masses connected with three springs, as described by: n =2 x1 9 = [| 2 wo [0 a Ky(8) = me +k} kl} | klit mygl + kl

Here, x4 and x, are respectively the x-coordinates of the first and second pendelum, m, and Mm, are their masses, and k is the spring constant.

Both pendula have the same length I, and the spring is attached at a distance I, from the pivot point.

In total there are six unknown dynamical parameters in (31). However, because of the structure of the matrices, only three parameters can be estimated, Again, the gravitational acceleration g is fixed to 9.81 m/s? such that the length [ can be estimated.

The remaining parameters are captured by 7, and 7,, which are defined as: n=2% 09) Ty = Kik (33) m;l2

This leaves only three unknown parameters to be estimated: 9=[LT,T2] (34)

It is assumed that the maximum displacement of each pendulum is not larger than half the distance between the pivot points, such that each pendulum is confined to a separate half-plane.

Two basis functions are used, one for each pendulum: (Tm) =1 Hr), p(X) = He(r) (35) where H, (I) is the Heaviside step function in the x-direction: 0, forx <0 ‚ Hr) = forx = 0 (36)

Spherical Pendulum If a simple pendulum is able to swing in two directions, it is called a spherical pendulum since its tip traces the surface of a sphere (Fig. 2¢). The dynamics of this system can be described in spherical coordinates by two coupled, nonlinear differential equations.

Again, these equations must be linearized around equilibrium.

For small displacements, this system is analogous to a single mass connected to two identical springs placed along the x- and y-directions: n =2 x 9 =|] Ee 0 (37) Mo(®) = 1] 19 © 20) = 0 gl Again, the length of the pendulum is the only unknown parameter: 6 =[!] (38) Two basis functions were used for the displacement in the x- and y-direction respectively:

TL _ [0 $10) =, vem =] vr 69

Coupled Spherical Pendula Finally, we couple two spherical pendula together by placing a spring between the two swinging rods (Fig. 2{d)). This results in a system with four degrees of freedom, The linearized differential equations are equivalent to the ones for two masses, each connected to two orthogonal springs, with another spring connecting the two masses: n =4 X1 _ | q - Xo V2 m,l2 0 0 0 0 ml? 0 0 Mp0) = ! (40) n(8) 0 m;l2 0 | 0 0 0 mol? migl+kit 0 klè 0 0 migl 0 0 K(@) =|_, ' 2 kl 0 mygl + kl; 0 0 0 0 m,gl The parameters 7, and 7, as described by (32) and (33) are used again, resulting in three unknown parameters: 6 =[l1715] (41) During reconstruction, four basis functions were used; two for the x- and y-displacements in one half of the spatial domain, and two for those in the other half: 1 —Hy(r) i | = , = 42 $10) =|, e= (42) H,(r) H | ry=|.~* , r)= 43 $20) = [7] =| (43)

Again, this assumes that the masses are confined to nonintersecting half-planes.

Data Acquisition for the example dynamic systems A setup was made using LEGO (The Lego Group, Billund, Denmark) that allowed for testing all the described dynamical systems.

Since this construction was made entirely from plastic, it did not generate any signal during the scan that could interfere with the measurements.

It consisted of two rigid arms with a pivot allowing for motion in two directions. The two pendula could additionally be coupled using a plastic spring. Since the linearization of the differential equations for the coupled systems required that the spring’s orientation remains straight, the spring was attached near the top of the pendula to minimize its displacement (see Fig. 2(b), 2(d), and 3).

Figure 3 shows photos of the setup from three different points of view. Two spherical pendula of length I were connected together with a spring (light blue arrow), attached at distance I from the pivot point. The complete setup was made from plastic. At the tip of each pendulum, a gel-filled glass vial was placed (red arrow). Thin cotton strings allowed operation of the pendula from a distance. At the extremity of each arm, a gel-filled glass vial (TOS, Eurospin II test system, Scotland) was placed in a vertical position (see Fig. 3), generating MR signal. Thin cotton strings with negligible weight were attached to each pendulum to allow manual operation of the pendula from outside the MRI scanner bore. The pendula were brought from equilibrium just before data acquisition started, and were able to swing freely during the measurements. By removing the second pendulum, or by not giving any initial displacement in the y-direction, all four systems as described above were measured using this single setup. The experiments were performed with a 1.5T MRI scanner (Ingenia, Philips Healthcare, Best, The Netherlands). A spoiled gradient echo sequence was used with the scan parameters as reported in Table 1. The body coil was used as receive array to get a homogeneous receive sensitivity. Data acquisition was done in a horizontal plane perpendicular to the orientation of the gel-filled tubes, with a Cartesian sampling pattern and a right-left readout direction. Every readout line was repeated Ny times before moving on to the next line, as in Fig. 1(b). The data was thus fully sampled along the k,-direction, while only a single k,, frequency was sampled per readout. Table 1: Data acquisition parameters used for the example dynamical systems Parameter Value ~ 1D Value - 2D Experiments | Experiments Repetition time (TR Echo time (TE Flip angle 9 degrees Field of view 420 mm x 195 mm Slice thickness Matrix size 212 x12 Number of repetitions (Nep) 2400 100 Total number of readouts (p 2400 1200

Reconstruction of displacements of the example systems The noise in the measured (complex-valued) k-space data was smoothed out by applying a Gaussian filter with a standard deviation of 10 ms in the temporal dimension (by convolution) and 1 cm in the readout direction (by windowing the k-space). The generalized coordinates q and the dynamical parameters 8 were estimated by solving the optimization problem as given in (26). A suitable value for the regularization parameter A was chosen by comparing the residuals of the measurement model and the dynamical model using the L-curve approach. Since the optimization is linear in q, variable projection (VARPRO) could be used to solve for q and 0 separately. Thus, the generalized coordinates could be estimated efficiently using a complete orthogonal decomposition as implemented in Matlab (The MathWorks Inc., Natick, MA, USA). The dynamical parameters 8 were optimized using Matlab’s Isgnonlin function.

To validate the estimated displacements, a 1D Fast Fourier Transform (FFT) was performed on each readout line. The maximum intensity of each glass vial was determined and, after smoothing over time, used as a reference position for that tube in the readout direction. Note that in these examples the data is fully sampled in the readout direction, but only one sample point is available in the phase-encoding direction (orthogonal to the readout direction) is available at each repetition interval. This reference signal for the displacements was compared to the reconstructed displacements by calculating the root mean square error (RMSE) for the 1, generalized coordinates g;(£) in the x-direction: RMSE, = [ES Su (0:6) = Gires 2)” (44) The estimated dynamical parameters 8 were validated by taking the FFT of the reference displacements in the temporal dimension. The resulting spectrum gives the contribution of each frequency in the dynamics. The frequencies with the highest contribution should be concentrated around the natural frequencies of the system, which can be found by solving the following eigenvalue problem for w: det(Kp (0) — w?Mp(8)) = 0 (45) The validation signal was zero-padded before taking the FFT, increasing the spectral resolution of the frequency spectrum. This allowed for better visibility of the frequencies corresponding to the peaks in the spectrum.

Undersampling The acquired data used has been fully sampled for the 1D acquisitions, and only undersampled in the phase-encode direction for the 2D experiments. In 3D experiments, an additional dimension is also undersampled. To investigate the behaviour of the reconstructions with even fewer available data, the readout direction was retrospectively undersampled by selecting a small subset of N samples per readout. These few points are preferably selected symmetrically around the center of k-space (except for N = 1). Again, the estimated displacements in the x-direction were validated using the RMSE as given by (44). The number of samples N are preferably at least equal to the number of degrees of freedom n, or else the measurement model could be underdetermined and it would be more difficult or no longer possible to estimate the displacement field.

Results The L-curves of the residuals for both models can be seen in Fig. 4, which shows L-curves for the experiments with (a) a simple pendulum, (b) coupled pendula, (c) a spherical pendulum, and (d) coupled spherical pendula. The horizontal axis is the residual of the measurement model, while the residual of the dynamical model is on the vertical axes. The 61 values for A used to create these curves ranged from 10 to 1.0 x 107, spaced evenly on a logarithmic scale. The red dot indicates the chosen value 4 = 5.0 x 103.

For all four experiments, a regularization parameter of 1 = 5.0 x 107 resulted in a good trade-off between suppressing noise in the measurement model, without introducing too much bias towards the dynamical model.

In Fig. 5, the 1D FFT of the data can be seen. Fig. 5 in particular shows 1D FFT of all readout lines, for the experiments with (a) a simple pendulum, (b) coupled pendula, (c) a spherical pendulum, and (d) coupled spherical pendula. The maximum intensity of each tube was determined and smoothed over time to obtain the red reference lines. Note that for the 2D acquisitions in (¢) and (d), only the data in the x-direction was reconstructed. Furthermore, the contrast changes after every 100 repetitions, as the next k-coordinate is sampled.

For the 1D experiments, which are fully sampled, the oscillations with one or two frequencies are clearly visible. The data of the 2D experiments can only be visualized along the x- direction, as only one k,, frequency was sampled at every point in time. Furthermore, transitions in the image contrast are visible every 100 readouts, as different k, -coordinates were acquired. Still, this 1D reconstruction can be used to generate a reference line for the displacements in the x- direction.

Table 2: Reconstructed dynamical parameters Experiment | RMSE, Im Simple 0.43 03 N.A. N.A. Pendulum Coupled 0.35 284 7.72 7.37 Pendula Spherical 1.3 02 N.A. N.A. Pendulum Coupled 1.1 301 7.72 8.06 Spherical Pendula T= (klj)/ (ml). 1, = (KE) (mpl?) The reconstructed displacements of all degrees of freedom can be seen in Fig. 6, which shows the estimated displacements. The black lines indicate the estimated displacements for the experiments with (a) a simple pendulum, (b) coupled pendula, (c) a spherical pendulum, and (d) coupled spherical pendula. For the 2D acquisitions in (c) and (d), the x- and y-displacements are plotted separately. The red line is the reference line as determined in Fig. 5.

The estimated displacements accurately follow the reference line, with a temporal resolution of one TR (4.4 ms). This is also evident from the RMSE,, which is in the order of 0.1 mm for the (fully sampled) 1D acquisitions, and around 1 mm for the (undersampled) 2D acquisitions (Table 3). Note that the amplitude of the estimated displacements is gradually decreasing, despite the fact that no friction was included in the dynamical model. This shows the robustness of our method against imperfections in the dynamical model, in particular due to the use of A.

For the displacements in the y-direction, no reference line is available. From the theoretical model, a harmonic oscillation with a slowly decreasing amplitude due to friction would be expected in this direction. However, as can be seen in Fig. 6c and 6d, the amplitude of the reconstructed displacements is not constant, nor is it gradually decreasing. Again, this shows the robustness of the method.

The estimated dynamical parameters are given in Table 3. The length of the pendulum is consistently estimated to be around 30 cm. The actual distance between the pivot point and the center of the tube is 32 cm. However, comparing these values is not always useful, as the pendulum is not a point mass. Additionally, since the masses of the two pendula are equal, the other two parameters for the coupled systems (7, and 72) should be approximately equal, which is the case as well. Note that, in these examples, these parameters are only estimated up to a certain scaling, as changing the masses or the position of the spring’ s attachment point has the same effect on the dynamics as changing the spring stiffness.

Table 3: Undersampled reconstructions Sor | Fi " Estimating the displacements in both directions is not possible with one measurement The natural frequencies of the dynamical systems are compared using the estimated dynamical parameters as in (45) to the frequency spectrum of the reference line used for validation.

Figure 7 shows natural frequencies of the estimated systems for the experiments with (a) a simple pendulum, (b) coupled pendula, (c) a spherical pendulum, and {d) coupled spherical pendula. The plotted spectra are the FFTs of the validation data, as determined in Fig. 5. The reference positions are zero-padded in the time domain for a better visualization of the peaks in the frequency spectrum. The black dashed lines indicate the natural frequencies for each system as calculated by (45), using the estimated dynamical parameters 8.

As can be seen in Fig. 7, the location of the natural frequencies correspond well to the peaks in the zero-filled FFT spectram. This indicates that the estimated dynamical parameters are accurately identified, and can be used to predict the dynamics of the systems.

Finally, the RMSE, values and estimated lengths for the retrospectively undersampled data can be seen in Table 3. Even with only 2 sample points per readout, the dynamics could still be determined, albeit with a slightly larger error compared to the reference. Fig. 8 shows estimated displacements for the undersampled datasets of (a) a simple pendulum, and (b) a spherical pendulum. The red line is the reference line as determined in Fig. 5. Note how the estimated displacements with only one or two samples every readout are almost identical to those estimated with all 212 samples in the readout direction. This indicates that our method is robust against high undersampling.

In the above, embodiments of the spectro-dynamic MRI are shown to provide an effective way for the identification of dynamical systems from MRI measurements at a very high temporal resolution. A measurement model has been derived to connect the measured data to the motion fields. By using the measured data directly in k-space, a temporal resolution of a few milliseconds is achieved. To handle highly undersampled data, the dynamical model has been added as a regularization term. These two models were optimized simultaneously to estimate the displacement field, as well as a set of dynamical system parameters. The method is validated on four relatively simple systems, and we have shown that it is possible to accurately estimate both ID and 2D dynamics with only two data points per readout.

The estimated displacements in Fig. 6 have a good agreement with the reference motion fields. Note that the amplitude of the reconstructed displacements gradually decreases. This is caused by the friction of the pendula. Even though no friction force has been included in the dynamical model, this phenomenon could still be reconstructed from the data through the measurement model together with the correct dynamical parameters. This shows the effectiveness of the trade-off coefficient A in (26), as dynamical model imperfections can still be reconstructed, especially when its value has been chosen correctly as described in relation to Fig. 4.

The sampling pattern can be chosen to improve the reconstruction. In the above examples, a Cartesian sampling pattern was used, with one readout line per time point. Thus, the data is intrinsically undersampled in the phase-encode direction, as only a single k „-coordinate is sampled during every readout. As a result, the motion along the readout direction could be estimated more accurately. The accuracy of the estimated displacements was dependent on which k„-coefficient was sampled. For most objects, most of the energy of the signal is located around the center of k- space where the low spatial frequencies are sampled, while the noise is distributed evenly over the entire k-space. Sampling at high spatial frequencies thus results in a low signal-to-noise ratio (SNR). On the other hand, at very low spatial frequencies, only global intensity changes are encoded, which are less sensitive to small displacements. Therefore, most of the useful information about the dynamics is located at those spatial frequencies with a sufficiently high level of detail, while still having sufficient SNR.

When performing a single measurement in the phase-encode direction, the convolution in (16) cannot be performed in this dimension. By choosing basis functions that are independent on the y-coordinate this limitation was bypassed. For other systems, the convolution can be evaluated, for instance by iteratively reconstructing the motion fields together with the missing data, enabling the calculation of the convolution.

In Fig. 6¢ and 6d, the estimated y-displacements during the acquisition of these samples is solely based on the dynamical model, effectively interpolating the dynamics between the closest surrounding points for which k,, # 0. Another option would be to use an interleaved sampling pattern such as the one of Fig. Ic, by sampling not one but a few different phase-encode coefficients, after which this pattern is repeated. This can be seen as a hybrid between Fig. Ia) and 1(b). An interleaved pattern will increase the amount of data sampled per time interval, at the cost of a larger At. With a few interleaved phase-encoding steps, this decrease in temporal resolution is acceptable for most dynamical systems, such as longues and the heart. By sampling several consecutive phase-encode lines and assuming that the displacement field is smooth and continuous in this direction, a way to calculate the convolution in the measurement model is provided.

The spectro-dynamic MRI model is generic and is applicable to 3D measurements. This will increase the effective undersampling factor since there are two phase-encode directions (k,, and k,), but the above shows that the method comprising the measurement and dynamical models are even effective when only 2 samples per readout are available. This robustness against undersampling shows the effectiveness of the method in 3D acquisitions. Optionally, more regularization can be added to the displacement fields, as is done in image registration. Prior knowledge about the dynamical parameters can also be added to regularize the reconstraction.

In the above examples, the error of (24) is minimized using ordinary least squares, assuming that only the vector b is subject to noise. It is an option to use a total least squares method to take possible errors in the matrix A into account during minimization of (26). This makes the optimization more complex, but results in an even better estimation of the dynamics.

The example dynamical systems described above are constructed using nondeformable objects. However, for instance the pendula connected with a spring can be seen together as a deformable object. The third term in the measurement model (16) also provides for obtaining results for nondeformable objects.

In summary, spectro-dynamic MRI enables the analysis of dynamical systems, such as the heart or the motion of joints, in vivo at a high temporal resolution. The measurement model does not depend on the system under investigation and can directly be applied, while the dynamical model is adjusted to model the object under investigation, e.g. deformable organ tissues, with different dynamical parameters.

Example of Velocity Field Approximation Assume we have an object in a certain reference configuration at time to. The location of one specific particle within this object at £, is Tg. As this object moves through space, the location r of this particle changes. The displacement field is a function of both space and time. It can be written as a function of the location of each particle in the reference configuration, called the Lagrangian description: u(ry, t) =r(ry,t) — ry (46)

Alternatively, we use the Eulerian description of the displacement field, in which the current coordinate is used as independent variable: u(r, t) =r —ry(r,t) 47) Notice that the initial coordinate of the particle is now no longer an independent variable, but it has become a function of the current coordinate. Since the magnetization is measured at fixed coordinates, it is more convenient to use the Eulerian description for the displacement field in the measurement model, as is done in Section

2.3. This requires an Eulerian expression of the velocity field as well. The velocity field is the rate of change in the position of each particle. That means we need to keep rp constant when taking the temporal derivative: dar Dr v= = == 48 dtly, =constant Dt (48) |l D. . Lo : . . Here, 57 is the material derivative. When the Lagrangian description as in (46) is used, the velocity field can be readily obtained from the displacement field: D a v(t) = —U(lg,t) = —u(ly,t) (49) Dt at However, for the Eulerian description as in (47), the total derivative with respect to time must be taken: v(r,t) =—u(rt) > Tr Dt > Zur) + Vu) - = 50) - Tat (r, (r, dt 0

J = u(r, t) + vu(r, t) - v(r,t) Solving for v{r, £) gives us an expression for the velocity field in the current configuration: -1 90 ‚ v(r,t) = (I — Vu(r, t)) tur, t) (51) where J is the identity matrix.

When the deformation of the moving object is small, Vu(r, £) is small as well. Note that this assumption does not limit the displacements themselves, only the spatial derivative of the displacement field. In this case, the velocity field can be approximated as the partial temporal derivative of the displacement field in the Eulerian description: v(t) = 2 u(r, t) (52) This assumption allows us to express the measurement model as a linear function of the displacement field, as can be seen in (12). This ends the velocity field approximation example.

Now referring to Fig. 9, the device 700, which is an embodiment of the device for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal {M) measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, comprises one or more processors 710 configured to perform any one or more of the methods described above. The device 700 may comprise a storage medium 720 for storing any of the model, parameters, and/or other data required to perform the method steps. The storage medium 720 may also store executable code that, when executed by the processor 710, executes one or more of the method steps as described above. The device 700 may also comprise a network interface 730 and/or an input/output interface 740 for receiving user input. Although shown as a single device, the device 700 may also be implemented as a network of devices such as a system for parallel computing or a supercomputer or the like. A person of skill in the art would readily recognize that steps of various above-described methods can be performed by programmed computers. Herein, some embodiments are also intended to cover program storage devices, e.g., digital data storage media, which are machine or computer readable and encode machine-executable or computer-executable programs of instructions, wherein said instructions perform some or all of the steps of said above-described methods. The program storage devices may be, e.g., digital memories, magnetic storage media such as a magnetic disks and magnetic tapes, hard drives, or optically readable digital data storage media. The embodiments are also intended to cover computers programmed to perform said steps of the above-described methods. The functions of the various elements shown in the figures, including any functional blocks labelled as “units”, “processors” or “modules”, may be provided through the use of dedicated hardware as well as hardware capable of executing software in association with appropriate software. When provided by a processor, the functions may be provided by a single dedicated processor, by a single shared processor, or by a plurality of individual processors, some of which may be shared. Moreover, explicit use of the term “unit”, “processor” or “controller” should not be construed to refer exclusively to hardware capable of executing software, and may implicitly include, without limitation, digital signal processor (DSP) hardware, network processor, application specific integrated circuit (ASIC), field programmable gate array (FPGA), read only memory (ROM) for storing software, random access memory (RAM), and non volatile storage. Other hardware, conventional and/or custom, may also be included. Similarly, any switches shown in the FIGS. are conceptual only. Their function may be carried out through the operation of program logic, through dedicated logic, through the interaction of program control and dedicated logic, or even manually, the particular technique being selectable by the implementer as more specifically understood from the context.

1t should be appreciated by those skilled in the art that any block diagrams herein represent conceptual views of illustrative circuitry embodying the principles of the invention. Similarly, it will be appreciated that any flow charts, flow diagrams, state transition diagrams, pseudo code, and the like represent various processes which may be substantially represented in computer readable medium and so executed by a computer or processor, whether or not such computer or processor is explicitly shown.

Whilst the principles of the described methods and devices have been set out above in connection with specific embodiments, it is to be understood that this description is merely made by way of example and not as a limitation of the scope of protection which is determined by the appended claims.

The disclosure comprises the following clauses:

1. Computer implemented method for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal (M) measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, the method comprising: providing a measurement model (G) defining a relation between the time dependent displacements of the moving object (U, q) in the spatial frequency domain and the measured SFDMR signal (M) in the spatial frequency domain; providing a dynamical model (F) dependent on a set of dynamical parameters (©) modelling dynamics of the moving object in the spatial frequency domain in relationship to the time dependent displacements of the moving object (U, q) in the spatial frequency domain; and performing optimization with an objective function based on the measurement model and the dynamical model for the time dependent displacements of the moving object (U, q) and the dynamical parameters (©) in order to determine an optimized or final set of the time dependent displacements of the moving object (U, q) and the dynamical parameters (©).

2. Computer implemented method according to clause 1, wherein the dynamical model (F) is a first differential equation defining a relation between the dynamical parameters (©), the displacements of the moving object (U, q) and temporal derivatives of the displacements.

3. Computer implemented method according to clause 1 or 2, wherein the measurement model {G) is a second differential equation defining a relation between the measured SFDMR signal (M), the displacements of the moving object (U, q) and temporal derivatives of the displacements, wherein preferably the relation defined by the second differential equation further comprises temporal derivatives of the SFDMR signal (M).

4. Computer implemented method according to clause 2 or 3, wherein the temporal derivatives of the displacements of the moving object (U, q) are approximated from a plurality of finite differences in the SFDMR signal (M) measured for a plurality of times.

5. Computer implemented method according to any of the preceding clauses, wherein the dynamical parameters (®) comprise one or more of a mass matrix, a stiffness matrix and a damping matrix associated with the moving object.

6. Computer implemented method according to any one of the preceding clauses, wherein the SFDMR signal is measured for a plurality of spatial frequencies.

7. Computer implemented method according to clause 6, wherein the SFDMR signal comprises a plurality of signal values measured at a same spatial frequency (k) of the plurality of spatial frequencies (k) and respective times, wherein a time difference between consecutive ones of the signal values measured at the same spatial frequency (Kk) is between 1 and 10 repetition times (TR), preferably between 1 and 5 repetition times (TR).

8. Computer implemented method according to clause 6 or 7, wherein the displacements of the moving object are represented by a displacement field (U) describing the displacements of the moving object in the spatial frequency domain for each time of the plurality of times and for each of the plurality of spatial frequencies (k); or a generalized coordinate vector (q) comprising generalized coordinates associated with the moving object for each time of the plurality of times.

9. Computer implemented method according to clause 8, in dependence of clause 5, wherein the dynamical model is implemented as a differential equation comprising one or more of: d2 Mp(6) TE q(t), Ly(8)3- q(t). and Kp(8) q(t),

wherein Mp, is the mass matrix, Lp is the damping matrix, Kp is the damping matrix, q is the generalized coordinate vector, and t is the time, wherein preferably the dynamical model is implemented as: d? d Mp(8) 7 0) + Cp(6) 77 q(0) + Kp(@)a(t) +Z = 0, wherein Z is a force term.

10. Computer implemented method according to any of clauses 6-9, and at least clause 3, wherein the measurement model is defined as: ZMK O+, KMD LUDO + 2m %; MK, 6) + ky Uy (Kk, t) = 0. wherein M is the SFDMR signal in the spatial frequency domain, k a vector comprising spatial frequencies for each datapoint of the measured SFDMR signal, t the time, j the spatial dimension, k; the spatial frequency along the j-th dimension, and U, a displacement field describing displacements of the moving object in the spatial frequency domain for spatial dimension j.

11. Computer implemented method according to any of the preceding clauses, wherein the step of performing optimization comprises solving an optimization problem based on the measurement model and the dynamical model, wherein the dynamical model is used as a regularisation term on the measurement model.

12. Computer implemented method according to clause 11, wherein the optimization problem is defined as the minimization, for the displacements of the moving object (U, q) and the dynamical parameters (©), of a sum of the norm of the measurement model and the dynamical model multiplied by a trade-off parameter that changes an amount of regularization by the dynamical model on the measurement model.

13. Computer implemented method according to clause 12 and at least clause 8, wherein the optimization problem is defined as: min heMEK 6), UK OKO IR +2IFUKD, 8,ko) |P, wherein: G is the measurement model, F is the dynamical model, M is the measured SFDMR signal, U is a displacement field describing displacements of the moving object in the spatial frequency domain, k is a vector comprising spatial frequencies for each datapoint of the measured SFDMR signal, t is the time, © is a vector comprising the dynamical parameters,

2 is the trade-off-parameter being a scalar parameter that changes the amount of regularization of the dynamical model on the measurement model, a is the power to which results of the respective norms (||..|) are raised, wherein preferably ais 1 or 2.

14. Computer implemented method according to clause 12 and at least clause 8, wherein the optimization problem is defined as: min I GIM(K, t), q(t), kK, ©) IR + 21 F(q(t),0,k, t) If, wherein: G is the measurement model, F is the dynamical model, M is the measured SFDMR signal, q is the generalized coordinate vector, k is a vector comprising spatial frequencies for each datapoint of the measured SFDMR signal, t is the time, © is a vector comprising the dynamical parameters, ò is the trade-off-parameter being a scalar parameter that changes the amount of regularization of the dynamical model on the measurement model, a is the power to which results of the respective norms are raised, wherein preferably a is 1 or2.

15. Computer implemented method according to any one of the preceding clauses, wherein the SFDMR signal is an undersampled SFDMR signal.

16. Computer implemented method according to any one of the preceding clauses, wherein the SFDMR signal comprises a plurality of signal values measured at a same spatial frequency (k) and respective times, wherein a time difference between consecutive ones of the signal values measured at the same spatial frequency (k) is smaller than a threshold time, wherein preferably the threshold time is between 1 and 100 repetition times (TR), more preferably between 1 and 50 repetition times (TR).

17. Device for determining time dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance, SFDMR, signal (M) measured for a plurality of times and emitted from the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, the device comprising one or more processors configured to: provide a measurement model (G) defining a relation between the time dependent displacements of the moving object (U, q) in the spatial frequency domain and the measured SFDMR signal (M) in the spatial frequency domain;

provide a dynamical model (F) dependent on a set of dynamical parameters (8) modelling dynamics of the moving object in the spatial frequency domain in relationship to the time dependent displacements of the moving object (U, q) in the spatial frequency domain; and perform optimization with an objective function based on the measurement model and the dynamical model for the time dependent displacements of the moving object (U, q) and the dynamical parameters (©) in order to determine an optimized or final set of time dependent displacements of the moving object (U, q) and the dynamical parameters (©).

18. A computer program product comprising computer-executable instructions for performing the method of any one of the clauses 1-16, when the program is run on a computer.

Claims (18)

CONCLUSIONS 1. Computer-implemented method for determining time-dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance (SFDMR) signal (M) measured for a plurality of times and emitted by the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, the method comprising: providing a measurement model (G) defining a relationship between the time dependent IO displacements of the moving object (U, q) in the spatial frequency domain and the measured SFDMR signal (M) in the spatial frequency domain; providing a dynamic model (F) which is dependent on a set of dynamic parameters (©) which model dynamics of the moving object in the spatial frequency domain in relation to the time dependent displacements of the moving object (U, g) in the spatial frequency domain; performing optimization with an objective function, based on the measurement model and the dynamic model, for the time-dependent displacements of the moving object (U, q) and the dynamic parameters (©) in order to obtain an optimized or final set of the time-dependent displacements of the moving object (U, q) and the dynamic parameters (©). A computer-implemented method according to claim 1, wherein the dynamic model {F) is a first differential equation defining a relationship between the dynamic parameters (©), the displacements of the moving object (U, q) and time derivatives of the displacements. A computer-implemented method according to claim 1 or 2, wherein the measurement model (G) is a second differential equation defining a relationship between the measured SFDMR signal (M), the displacements of the moving object (U, q) and time derivatives of the displacements preferably, the relationship defined by the second differential equation further comprising time derivatives of the SFDMR signal (M). The computer-implemented method of claim 2 or 3, wherein the tidal derivatives of the displacements of the moving object (U, q) are approximated from a plurality of finite differences in the SFDMR signal (M) measured for a plurality of times. A computer-implemented method according to any preceding claim, wherein the dynamic parameters (©) comprise one or more of a mass matrix, a stiffness matrix and a damping matrix associated with the moving object. A computer-implemented method according to any preceding claim, wherein the SFDMR signal is measured for a plurality of spatial frequencies. The computer-implemented method of claim 6, wherein the SFDMR-IO signal comprises a plurality of signal values measured at the same spatial frequency (k) of the plurality of spatial frequencies (k) and respective times, wherein a time difference between successive of the signal values measured at the same spatial frequency (k) is between 1 and 10 repetition times ("repetition times", TR), preferably between 1 and 5 repetition times (TR). A computer-implemented method according to claim 6 or 7, wherein the displacements of the moving object are represented by a displacement field (U) describing the displacements of the moving object in the spatial frequency domain for each time of the plurality of times and for each of the plurality of spatial frequencies (k); or a generalized coordinate vector (q) comprising generalized coordinates associated with the moving object for each time of the plurality of times. The computer-implemented method of claim 8 when dependent on claim 5, wherein the dynamic model is implemented as a differential equation comprising one or more of: d2 Mp(8) za), Ly(8) = q(t), and Kp (8)q(t). where Mp is the mass matrix, Lp is the damping matrix, Kp is the stiffness matrix, q is the generalized coordinate vector and t is time, preferably the dynamic model is implemented as: d? d Mp (0) a8) + Cp(@) = a(t) + Kn(0)a(®) +2 = 0, where Z is an expletive. A computer-implemented method according to any one of claims 6 to 9 and at least claim 3, wherein the measurement model is defined as: = Mk, D+; kM) * Zj , t) + 21 ny ; M(kt)* kj LU, £) = 0, where M is the SFDMR signal in the spatial frequency domain, k is a vector containing spatial frequencies for each data point of the measured SFDMR signal, t is the time, j is the spatial dimension is, k; is the spatial frequency along the jth dimension and U; is a displacement field that describes displacements of the moving object in the spatial frequency domain for spatial dimension j. A computer-implemented method according to any preceding claim, wherein the step of performing optimization comprises solving an optimization problem based on the measurement model and the dynamic model, wherein the dynamic model is used as a regularization term on the measurement model. A computer-implemented method according to claim 11, wherein the optimization problem is defined as the minimization, for the displacements of the moving object (U, q) and the dynamic parameters (0), of a sum of the norm of the measurement model and the dynamic model multiplied by a weighting parameter that adjusts an amount of regularization by the dynamic model to the measurement model. A computer-implemented method according to claim 12 and at least claim 8, wherein the optimization problem is defined as: min IGMKOUKOkO IR +All F(UK,L),0,Kkt) 12, where: G is the measurement model, F is the dynamic model , M is the measured SFDMR signal, U is a displacement field describing displacements of the moving object in the spatial frequency domain, k is a vector containing spatial frequencies for each data point of the measured SFDMR signal, t is time, © a is vector which includes the dynamic parameters, 2. is a weighting parameter which is a scalar parameter that adjusts the amount of regularization by the dynamic model to the measurement model, a is the exponent of the result of the respective norms, preferably a is 1 or 2. A computer-implemented method according to claim 12 and at least claim 8, wherein the optimization problem is defined as: min Ih GMK, 0), q(t), Kk, t) I? + A Il F(q(t),0,k,t) II*, where: G is the measurement model, F is the dynamic model, M is the measured SFDMR signal, q is the generalized coordinate vector, k is a vector which includes spatial frequencies for each data point of the measured SFDMR signal, t is the time, © is a vector containing the dynamic parameters, è is a weighting parameter which is a scalar parameter representing the amount of regularization by the dynamic model on the measurement model adjusts, a is the exponent of the result of the respective norms, where preferably a is 1 or 2. A computer implemented method according to any preceding claim, wherein the SFDMR signal is an undersampled SFDMR signal. A computer-implemented method according to any preceding claim, wherein the SFDMR signal comprises a plurality of signal values measured at the same spatial frequency (k) and respective times, wherein a time difference between successive of the signal values measured at the same spatial frequency (k) is less than a threshold time, preferably the threshold time is between 1 and 100 repetition times (TR), more preferably between 1 and 50 repetition times (TR). 17. Apparatus for determining time-dependent displacements of a moving object based on a measured spatial frequency domain magnetic resonance (SFDMR) signal (M) measured for a plurality of times and emitted by the moving object after magnetic resonance excitation of the moving object according to an applied pulse sequence, the device comprising one or more processors configured to: provide a measurement model (G) defining a relationship between the time-dependent displacements of the moving object ( U, q) in the spatial frequency domain and the measured SFDMR signal (M) in the spatial frequency domain; providing a dynamic model (F) which is dependent on a set of dynamic parameters (©) which model dynamics of the moving object in the spatial frequency domain in relation to the time dependent displacements of the moving object (U, q) in the spatial frequency domain; performing optimization with an objective function, based on the measurement model and the dynamic model, for the time-dependent displacements of the moving object (U, q) and the dynamic parameters (©) in order to obtain an optimized or final set of the time-dependent displacements of the moving object (U, q) and the dynamic parameters (©). A computer program product comprising computer executable instructions for performing the method of any one of claims 1 to 16 when the program is run on a computer.