METHOD AND APPARATUS FOR PERFORMING ACCELARATED MAGNETIC RESONANCE IMAGING WITH REDUCED OFF-RESONANCE EFFECT

Abstract

A method of performing magnetic resonance imaging of a body using a magnetic resonance imaging scanner the method includes applying to the body a time-varying magnetic field gradient (G.sub.x, G.sub.y, G.sub.z) defining a continuous trajectory (ST) in k-space complying with a set of constraints including constraints on maximum amplitude and maximum slew rate of the time-varying magnetic field gradient, such that sampling points (KS) belonging to the trajectory define a pseudo-random sampling of the k-space, approximating a predetermined target sampling density, the trajectory in k-space minimizing, subject to the set of constraints, a cost function defined by the difference between a first term, called attraction term, promoting consistency of the distribution of sampling points in k-space with the predetermined target sampling density, and a second term, called repulsion term, promoting separation in k-space between sampling points, the repulsion term being expressed as a sum of contributions corresponding to respective pairs of sampling points; wherein each the contribution is weighted by a weight which increase with temporal separation of the sampling points along the trajectory in k-space.

Claims

1. A method of performing magnetic resonance imaging of a body using a magnetic resonance imaging scanner comprising a scanner bore (SB), a primary coil (PC), radio frequency coils (RFC), gradient coils (GC.sub.x, GC.sub.y, GC.sub.z) and a signal-processing unit (SPU), the method comprising the steps of: a. positioning the body (BD) in a scanner bore (SB) where a static and substantially uniform magnetic field (B.sub.0), called longitudinal field, oriented along a direction (z), called longitudinal direction, is established by the primary coil (PC); b. using all or part of the radio-frequency coils (RFC) to transmit to said body at least one radio-frequency pulse (RFP) adapted to excite nuclear spins inside said body; c. after said or each said radio-frequency pulse, using the gradient coils (GC.sub.x, GC.sub.y, GC.sub.z) to apply to said body a time-varying magnetic field gradient (G.sub.x, G.sub.y, G.sub.z) defining a trajectory (ST) in k-space and simultaneously using all or part of the radio-frequency coils to acquire samples of a magnetic resonance signal emitted by the excited nuclear spin, each sample corresponding to a point (KS) of the k-space belonging to said trajectory; and d. using the signal processing unit (SPU) to apply a nonlinear reconstruction algorithm to the acquired samples to reconstruct a magnetic resonance image of said body; wherein said trajectory (ST, MT) in k-space is a continuous trajectory complying with a set of constraints including constraints on maximum amplitude and maximum slew rate of said time-varying magnetic field gradient, such that the points (KS) of the k-space corresponding to the samples, called sampling points, define a pseudo-random sampling of the k-space, matching a predetermined target sampling density, said trajectory in k-space minimizing, subject to said set of constraints, a cost function defined by the difference between a first term, called attraction term, promoting consistency of the distribution of sampling points in k-space with said predetermined target sampling density, and a second term, called repulsion term, promoting separation in k-space between sampling points, said repulsion term being expressed as a sum of contributions corresponding to respective pairs of sampling points; wherein each said contribution is weighted by a weight which increase with temporal separation of the sampling points along said trajectory in k-space.

2. The method of claim 1, wherein said attraction term is defined by $F_p^a\left(K,\Pi \right)=\underset{K\in Q^{\mathrm{Nc}}}{\mathrm{argmin}}\frac{1}{p}\underset{i=1}{\overset{p}{.Math.}}{\int }_{\Omega }.Math.x-K\left[i\right]\Pi \left(x\right)\mathrm{dx}.Math.$ and said repulsion term is defined by $F_p^r\left(K\right)=\frac{1}{2p^2}\underset{1\le i,j\le p}{.Math.}W\left(.Math.t_i-t_j.Math.\right).Math.K\left[i\right]-K\left[j\right].Math.$ where: K is a two or three dimensional vector representing the coordinates of a sampling point in k-space; p=N.sub.c.Math.N.sub.s is the number of sampling points along the trajectory, each trajectory being constituted of N.sub.c≥1 shots each comprising N.sub.s>1 sampling points; Q.sup.Nc is a set of curves in k-space comprising N.sub.c segments of N.sub.s points each and complying with said constraints; Ω is a sampling region in k-space; x represents a generic point of the sampling region Ω; Π is said predetermined target sampling density; ∥.Math.∥ is a norm; t.sub.i, t.sub.j are the times at which two sampling points of a same shot, identified by integer indices i and j, are reached; and W is a monotonously increasing function of |t.sub.i−t.sub.j| expressing said weight.

3. The method of claim 2, wherein norm ∥.Math.∥ is the L.sub.2 norm.

4. The method of claim 2, wherein W is an exponential function of |t.sub.i−t.sub.j|.

5. The method of claim 1, wherein said set of constraints further comprises at least one linear constraint on the trajectory.

6. The method of claim 1, further comprising a preliminary step of computing said trajectory (ST) in k-space.

7. A computer-implemented method of computing a trajectory in k-space to be used for sample acquisition in magnetic resonance imaging, the method comprising the steps of: i. determining, or receiving from a user input, a number N.sub.c≥1 of shots constituting the trajectory and a number N.sub.s>1 of sampling points along each said shot; ii. determining, or receiving from a user input, a predetermined target sampling density; iii. determining, or receiving from a user input, a set of constraints including constraints on amplitudes of discrete time derivatives of the trajectory; iv. computing said trajectory by minimizing, subject to said set of constraints, a cost function defined by the difference between a first term, called attraction term, promoting consistency of the distribution of sampling points in k-space with said predetermined target sampling density, and a second term, called repulsion term, promoting separation in k-space between sampling points, said repulsion term being expressed as a sum of contributions corresponding to respective pairs of sampling points; wherein each said contribution is weighted by a weight which increases with temporal separation of the sampling points along said trajectory in k-space.

8. The method of claim 7, wherein said attraction term is defined by $F_p^a\left(K,\Pi \right)=\underset{K\in Q^{\mathrm{Nc}}}{\mathrm{argmin}}\frac{1}{p}\underset{i=1}{\overset{p}{.Math.}}{\int }_{\Omega }.Math.x-K\left[i\right]\Pi \left(x\right)\mathrm{dx}.Math.$ and said repulsion term is defined by $F_p^r\left(K\right)=\frac{1}{2p^2}\underset{1\le i,j\le p}{.Math.}W\left(.Math.t_i-t_j.Math.\right).Math.K\left[i\right]-K\left[j\right].Math.$ where: K is a two or three dimensional vector representing the coordinates of a sampling point in k-space; p=N.sub.c.Math.N.sub.s is the number of sampling points along the trajectory; Q.sup.Nc is a set of curves in k-space comprising N.sub.c segments of N.sub.s points each and complying with said constraints; Ω is a normalized sampling region in k-space; x represents a generic point of the sampling region Ω; Π is said predetermined target sampling density; ∥.Math.∥ is a norm; t.sub.i, t.sub.j are the times at which two sampling points of a same shot, identified by integer indices i and j, are reached; and W is a monotonously increasing function of |t.sub.i−t.sub.j| expressing said weight.

9. The method of claim 8, wherein norm ∥.Math.∥ is the L.sub.2 norm.

10. The method of claim 8, wherein W is an exponential function of |t.sub.i−t.sub.j|.

11. The method of claim 10, wherein step iv. is performed through multi-resolution iterations with a decreasing decimation level R, W being a decreasing function of R.

12. The method of claim 7, wherein said set of constraints further comprises at least one linear constraint on the trajectory.

13. A computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method of claim 7.

14. A set of driving signals (DSG) for gradient coils of a magnetic resonance imaging scanner which, when applied to said gradient coils, drive them to generate a time-varying magnetic field gradient (G.sub.x, G.sub.y, G.sub.z) defining a trajectory (ST, MT) in k-space passing through a plurality of points (KS) called sampling points defining a pseudo-random sampling of the k-space, matching a predetermined target sampling density; wherein said trajectory: is continuous and complies with a set of constraints including constraints on maximum amplitude and maximum slew rate of said time-varying magnetic field gradient; minimizes, subject to said set of constraints, a cost function defined by the difference between a first term, called attraction term, promoting consistency of the distribution of sampling points in k-space with said predetermined target sampling density, and a second term, called repulsion term, promoting separation in k-space between sampling points, said repulsion term being expressed as a sum of contributions corresponding to respective pairs of sampling points; wherein each said contribution is weighted by a weight which increases with temporal separation of the sampling points along said trajectory in k-space.

15. A magnetic resonance imaging scanner comprising: a scanner bore (SB) inside which a body can be positioned; a primary coil (PC) configured to establish in the scanner bore a static and substantially uniform magnetic field (Bo), called longitudinal field, oriented along a direction (z), called longitudinal direction; radio frequency coils (RFC) configured to transmit to said body at least one radio-frequency pulse (RFP) adapted for exciting nuclear spins inside said body and to acquire samples of a magnetic resonance signal emitted by the excited nuclear spin; gradient coils (GC.sub.x, GC.sub.y, GC.sub.z) configured to apply to said body a time-varying magnetic field gradient (G.sub.x, G.sub.y, G.sub.z) defining a trajectory (ST, MT) in k-space; a signal processing unit (SPU) to apply a nonlinear reconstruction algorithm to the acquired samples of the magnetic resonance signal for reconstructing a magnetic resonance image of said body; and a controller (CTR) configured to generate driving signals (DSR) for the radio frequency coils and (DSG) for the gradient coils; wherein the controller is configured to generate said driving signals for (DSG) for the gradient coils to generate a said time-varying magnetic field gradient (G.sub.x, G.sub.y, G.sub.z) such that said a trajectory (ST) in k-space passes through a plurality of points (KS) called sampling points defining a pseudo-random sampling of the k-space, matching a predetermined target sampling density, said trajectory being continuous and compliant with a set of constraints including constraints on maximum amplitude and maximum slew rate of said time-varying magnetic field gradient, said trajectory in k-space also minimizing, subject to said set of constraints, a cost function defined by the difference between a first term, called attraction term, promoting consistency of the distribution of sampling points in k-space with said predetermined target sampling density, and a second term, called repulsion term, promoting separation in k-space between sampling points, said repulsion term being expressed as a sum of contributions corresponding to respective pairs of sampling points; and wherein the controller is also configured to generate driving signals (DSR) for the radio frequency coils such that they acquire said samples of a magnetic resonance signal emitted by the excited nuclear spin in correspondence to the sampling points (KS) of the trajectory in k-space; wherein each said contribution to the repulsion term of the cost function is weighted by a weight which increase with temporal separation of the sampling points along said trajectory in k-space.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Additional features and advantages of the present invention will become apparent from the subsequent description, taken in conjunction with the accompanying drawings, which show:

[0056] FIG. 1, the chronogram of a modified GRE (Gradient Recalled Echo) pulse sequence to play implementing a “SPARKLING” trajectory according to the prior art or a “MORE SPARKLING” trajectory according to an embodiment of the invention.

[0057] FIG. 2, an exemplary 2D “SPARKLING” trajectory according to the prior art;

[0058] FIG. 3, an exemplary ΔB.sub.0 map of human brain;

[0059] FIG. 4 illustrates off-resonance artifacts in MRI image reconstruction in the case of Cartesian k-space sampling according to the prior art;

[0060] FIG. 5 illustrates off-resonance artifacts in MRI image reconstruction in the case of SPARKLING k-space sampling according to the prior art;

[0061] FIG. 6 compares a 3D SPARKLING k-space trajectory, according to the prior art (upper panel) to a 3D MORE-SPARKLING trajectory according to an embodiment of the invention (lower panel), showing the impact of temporal weighting;

[0062] FIGS. 7 and 8 illustrate the technical effect of the invention, namely the reduction in off-resonance artifacts obtained using 3D MORE-SPARKLING k-space trajectories instead of 3D SPARKLING ones;

[0063] FIG. 9 compares different 3D MORE SPARKLING trajectories with varying temporal weights; and

[0064] FIG. 10 is a block diagram of a MRI scanner according to an embodiment of the invention.

DETAILED DESCRIPTION

[0065] FIG. 1 shows the chronogram of a pulse sequence—more particularly a Gradient Recalled Echo (GRE) sequence for T2* weighted anatomical imaging suitable for implementing either a SPARKLING k-space trajectory according to the prior art or a MORE-SPARKLING one.

[0066] It is supposed that a body to be imaged (e.g. the head of a patient) is immerged in a static and homogeneous magnetic field B.sub.0 oriented along a “longitudinal” direction z. This magnetic field induces a partial alignment of the spins of the atomic nuclei of the body. The aligned nuclear spins may be excited (i.e. flipped) by a radio-frequency pulse RFP at a suitable frequency (Larmor frequency), depending on the values of the magnetic field B.sub.0 and of the magnetic momentum of the nuclei. In 2D sequences, as the one illustrated on the figure, a magnetic field gradient G.sub.z—i.e. a magnetic field oriented along the z direction, whose magnitude varies linearly along this same direction—is also applied. This has the effect of making the Larmor frequency variable along z; as a consequence, only the nuclear spins within a thin slice of the body are resonant with the RF pulse and can be excited. As known in the art of MRI, this “slice selection” gradient pulse is followed by a short negative blip of the G.sub.z magnetic field gradient (“refocusing gradient”) which compensates for a dispersion of the nuclear spin orientations in the xy plane, perpendicular to the z direction.

[0067] Due to the use of a slice selection gradient G.sub.z, only a 2D image of the selected slice of the body is acquired, which requires sampling of a 2D k.sub.xk.sub.y plane of the k-space; in the following, the expression “k-space” will be used to designate both the three-dimensional k.sub.xk.sub.yk.sub.z space and a two-dimensional plane within it. In alternative embodiments, no slice selection gradient is used, and a domain of the 3D k.sub.xk.sub.yk.sub.z space will have to be sampled.

[0068] A trajectory in the k-space (more precisely, in the 2D k.sub.xk.sub.y plane) is defined by playing G.sub.x and G.sub.y gradients after the end of the RF excitation pulse. It is important to underline that the applied magnetic field is always oriented along the z direction, but its magnitude shows linear variations along the x and y directions. First of all, G.sub.x and G.sub.y pulses are applied to reach a suitable point on the boundary of the region of the k.sub.xk.sub.y plane to be sampled. Then “arbitrary” G.sub.x and G.sub.y waveforms are played, defining a winding trajectory with an overall radial orientation, directed toward the center of the k.sub.xk.sub.y plane.

[0069] In “3D” implementations without slice selection gradient, a G.sub.z waveform similar to G.sub.x and G.sub.y will be applied to reach a suitable point of the boundary of the 3D domain in the k.sub.xk.sub.yk.sub.z space and to define a three-dimensional trajectory in k-space.

[0070] While the gradient waveforms are played, samples of the NMR signal emitted by the excited nuclei are acquired by one or more radio-frequency coils connected to a suitable acquisition circuit including an analog-to-digital (ADC) converter. The acquisition period, whose duration T.sub.obs is limited by the decay of the NMR signal, is designated by reference OP on FIG. 1. After the end of the acquisition period, a final G.sub.x gradient pulse (“spoiling gradient”) is applied to cancel residual transversal magnetization. The repetition time TR corresponds to the interval beginning just before the radio-frequency pulse RFP and ending at the end of the spoiling gradient pulse.

[0071] This sequence is repeated several times with different gradient waveforms defining respective k-space trajectories which, together, provide the required k-space sampling. The ensemble constituted by an excitation RF pulse and the associated gradient waveforms is called a “shot”; each shot corresponds to an elementary trajectory. FIG. 2 shows a 2D “multishot” SPARKLING trajectory MT comprising N.sub.c=60 radially-oriented elementary trajectories or shots ST, inhomogeneously sampling the k.sub.xk.sub.y plane.

[0072] In practice, the magnetic field gradients undergo stepwise changes at discrete time-points separated by intervals of duration Δt (“gradient raster time”). Sampling is also performed at regular intervals of duration dt (“ADC dwell time”). The ADC dwell time dt is preferably lower than, or at most equal to, the gradient raster tile (dt≤Δt) so as to allow collecting several samples between two consecutive gradient steps. At the same time, reducing the ADC dwell time beyond a certain limit decreases the SNR to an unacceptable level. Therefore, for each specific embodiment of the invention there is an optimal value for dt which can be found.

[0073] Let k.sub.i(1) the position, in the k-space, of the starting point of the trajectory associated to the i.sup.th shot. A first sample of the NMR signal is acquired in correspondence to this point. The other sampling points correspond to k-space positions given by:

[00002]$\begin{array}{cc}{k}_i\left(m\right)=k_i\left(1\right)+\gamma \left(\underset{j=1}{\overset{q}{.Math.}}{G}_i\left(j\right)*\Delta t+G_i\left(q+1\right)*r\right)& \left(2\right)\end{array}$

where m∈[2: M] is an integer index, M being the overall number of samples acquired along the trajectory, and q and are respectively the modulus and the rest of the Euclidean division of the acquisition time by:

t.sub.ADC,m=(m−1)*dt=q*Δt+r   (3)

[0074] If dt=Δt, then r=0 and the number of ADC samples matches the number of gradient time steps. If dt<Δt, then the number of ADC samples is larger than the number of gradient time steps.

[0075] “SPARKLING” relies on an optimization-based method that consists of projecting a target sampling distribution over a set of discrete pushforward measures, in particular supported by smooth trajectories [Lazarus et al, 2017; Boyer et al, 2016]. Mathematically, this problem can be cast as a non-convex variational optimization problem under possibly non-convex constraints [Boyer et al, 2016]. The constraints are typical expressed by maximal acceptable values for the gradient amplitude and slew rate, but additional affine constraints may also be used—e.g. imposing that the trajectory passes through a particular point of the k-space at a defined time point (e.g, for echo time definition).

[0076] According to ([Boyer et al, 2016], [Lazarus et al. 2017], [Lazarus et al. 2019] and [Chaithya et al. 2022], a SPARKLING trajectory is computed by minimizing, subject to a set of constraints, a cost function defined by the difference between two terms: an “attraction term”, which promotes consistency of the distribution of sampling points in k-space with a predetermined (usually non-uniform) target sampling density (TSD), Π, and a “repulsion term”, promoting separation in k-space between sampling points. This can be written as:

[00003]$\begin{array}{cc}\hat{K}=\underset{K\in Q^{\mathrm{Nc}}}{\mathrm{argmin}}{F}_p\left(K,\Pi \right)=F_p^a\left(K,\Pi \right)-F_p^r\left(K\right)& \left(4\right)\end{array}$ [0077] Where

[00004]$\begin{array}{cc}{F}_p^a\left(K,\Pi \right)=\underset{K\in Q^{\mathrm{Nc}}}{\mathrm{argmin}}\frac{1}{p}\underset{i=1}{\overset{p}{.Math.}}{\int }_{\Omega }.Math.x-K\left[i\right]\Pi \left(x\right)\mathrm{dx}.Math.& \left(5\right)\end{array}$

is the attraction term which ensures the sampling patterns K matches the prescribed TSD Π,

[00005]$\begin{array}{cc}{F}_p^r\left(K\right)=\frac{1}{2p^2}\underset{1\le i,j\le p}{.Math.}.Math.K\left[i\right]-K\left[j\right].Math.& \left(6\right)\end{array}$

is the repulsion term and: [0078] K is a two or three dimensional vector representing the coordinates of a sampling point in k-space; [0079] p=N.sub.c.Math.N.sub.s is the number of sampling points along the trajectory, N.sub.c≥1 being the number of shots, or elementary trajectories, constituting the whole trajectory and N.sub.s>1 the number of sampling points of each shot; [0080] Q.sup.Nc is a set of curves in k-space comprising N.sub.c shots, or segments, each comprising N.sub.s points and complying with the constraints; [0081] Ω is a sampling region in k-space, normalized to [−1, 1].sup.3 (for 3D trajectories); [0082] x represents a generic point of the sampling region Ω; [0083] ∥.Math.∥ is a norm, preferably the L.sub.2 norm.

[0084] The set of constraints can be defined as follows

[00006]$\begin{array}{cc}\left\{\begin{array}{c}A\circ K=V\\ \forall i=1.Math.N_c,k_i\in {\Omega }^{Ns}\\ {.Math.k_i.Math.}_{\infty }\le 1,{.Math.{\dot{k}}_i.Math.}_{2,\infty }\le \alpha ,{.Math.{\stackrel{.Math.}{k}}_i.Math.}_{2,\infty }\le \beta \end{array}\right\}& \left(7\right)\end{array}$

where: [0085] A∈R.sup.Nc×Ns×3; V∈R.sup.Nc×Ns×3 [0086] A∘K is the Hadamard product between A and K; [0087] A∘K=V models linear constraints on the trajectories, for instance the Echo Time constraint, which ensure that teach elementary trajectory (shot) passes through the center of the k-space at the echo time TE (see [Chauffert et al. 2016]); [0088] {dot over (k)}.sub.i, {umlaut over (k)}.sub.i are, respectively, the first and second time derivative of the trajectory in k-space;

[0089] For an arbitrary three-dimensional vector c,

[00007]${.Math.c.Math.}_{2,\infty }=\underset{0\le n\le \mathrm{Ns}-1}{\sup }\sqrt{{.Math.c_x\left[n\right].Math.}^2+{.Math.c_y\left[n\right].Math.}^2+{.Math.c_Z\left[n\right].Math.}^2;}$ [0090] α, β are normalized parameters expressing hardware constraints on gradient amplitude and gradient slew rate; and [0091] Ω.sup.N.sup.s represents Ns points in Ω.

[0092] Usually, the trajectory design process starts with an initialization trajectory, formed e.g. by N.sub.c radial spokes, then the optimization problem (4) is solved under constraints (7) by the projected gradient descent algorithm, as discussed in [Chaithya et al. 2020].

[0093] In practice, the optimization is performed through multi-resolution iterations which start by spreading N.sub.R.sub.max=p/2.sup.R.sup.max samples at a maximum decimation level R.sup.max (e.g. R.sup.max=5, i.e. only one in 2.sup.R.sup.max=2.sup.5=32 samples is considered) and iterates through a decreasing dyadic process, i.e. R.sup.max=5 (N.sub.0=p/32) to R=0 (N.sub.final=p). This is performed to ensure that the optimization is carried out with faster convergence when the k-space trajectory is optimized coarsely first, (R=R.sup.max) and is then progressively refined a convergence is approached (R=0, meaning no decimation).

[0094] Ideally, for MRI—whatever the sampling strategy of the k-space—the longitudinal magnetic field Bo should be as uniform as possible, and the primary coil of a MRI scanner, responsible for generating this field, is specifically designed to this aim. The body to be imaged, however, is inhomogeneous, and its non-uniform magnetic susceptibility in turns induces spatial inhomogeneities ΔB.sub.0 of the magnetic field B.sub.0. In the case of a human or animal body, these inhomogeneities are particularly strong in the vicinity of air-filled cavities. This is illustrated on FIG. 3 for the case of a human head, where paranasal sinuses are mainly responsible for magnetic field inhomogeneities in the brain.

[0095] Ideally, for a perfectly homogeneous B.sub.0 field, the MRI signal s at time t is given by

s(t)=∫f(r)e.sup.−i2π(k(t).Math.r)dr   (8)

where f(r) is the spatial distribution of the sample magnetization (r representing position) and k(t) is the sampling trajectory in k-space. Equation (8) is inverted to reconstruct f(r) from s(t). When spatial inhomogeneities ΔB.sub.0 of the magnetic field B.sub.0 are taken into account, (8) becomes

s(t)=∫f(r)e.sup.−iω(r)te.sup.−i2π(k(t).Math.r)dr   (9)

where ω(r)=γΔB.sub.0(r) (γ being the gyromagnetic ratio) is the position-dependent variation of the Larmor frequency due to the magnetic field inhomogeneity. If the additional term e.sup.−iω(r)t is not taken into account (which requires the knowledge of ΔB.sub.0(r) and is computationally costly and/or requires an additional acquisition which would make the overall SPARKLING not competitive at acquisition), it introduces artifacts in the reconstructed image. These artifacts also depend on the sampling trajectory in k-space k(r).

[0096] FIG. 4 shows the impact of ΔB.sub.0(r) on a test pattern (left panel: test pattern; right panel: reconstructed image) when a Cartesian line-by-line sampling trajectory in k-space is used. It can be seen that almost no information is lost, even if the patter is slightly deformed.

[0097] FIG. 5 shows the impact of ΔB.sub.0(r) on a test pattern of FIG. 4 (left panel) when a 2D SPARKLING sampling trajectory in k-space is used. As it can be seen in the right panel of the figure, details are lost, spurious ones are created and in some regions a substantial loss of signal is even observed. These artifacts are much more detrimental than the simple geometrical distortion observed with Cartesian sampling.

[0098] The inventors have realized that this detrimental effect is partly because SPARKLING samples the k-space in a temporally inhomogeneous way. Otherwise stated, nearby locations of the k-space may be sampled at very different times, for instance at the beginning and at the end of the trajectory. This is important because the signal amplitude varies during the acquisition time due to T.sub.2 decay. When this decay is taken into account, equation (9) becomes

s(t)=∫f(r)e.sup.−(α(r)+ω(r))te.sup.−i2π(k(t).Math.r)dr   (10)

where α(r) is the inverse of the (spatially dependent) T.sub.2 decay time constant. With this, having k-space samples of s(t) taken in the same region of k-space in a temporally inhomogeneous way would result in varying values of s(t) for the same region of k-space. During reconstruction, these values are averaged out resulting in amplification of the off-resonance artifacts.

[0099] To reduce the impact of the off-resonance artifacts on image reconstruction, the present invention propose to modify the SPARKLING trajectories to make them temporally homogeneous—i.e. to ensure that regions of the k-space that are close to each other are sampled at close times—while keeping their significant advantages in terms of acceleration of the signal acquisition. More particularly, this is achieved by modifying the expression of the repulsion term of the cost function (equation (6)) by the introduction of a weighting function which is a growing function of the temporal separation between sampling points. The repulsion term becomes then

[00008]$\begin{array}{cc}{F}_p^r\left(K\right)=\frac{1}{2p^2}\underset{1\le i,j\le p}{.Math.}W\left(.Math.t_i-t_j.Math.\right).Math.K\left[i\right]-K\left[j\right].Math.& \left(11\right)\end{array}$

where t.sub.i, t.sub.j are the times at which two sampling points of a same shot, identified by integer indices i and j, are reached; and W is a monotonously increasing function of |t.sub.i−t.sub.j|, preferably of value greater than one, expressing said weight.

[0100] In a first embodiment of the invention, W(|t.sub.i−t.sub.j|) is an exponential function:

[00009]$\begin{array}{cc}W\left(.Math.t_i-t_j.Math.\right)=\exp \left(\frac{.Math.t_i-t_j.Math.\tau }{p}\right)& \left(12\right)\end{array}$

with τ≥0 is a scalar repulsion weighting parameter. It is important to note that an excessively strong temporal weighting of the repulsion (i.e. too large a value of τ) results in k-space holes, which is not desirable for good reconstructed image quality. To prevent this, τ needs to be grid-searched appropriately to enforce temporally smooth k-space sampling, while avoiding undesirable k-space holes. This is done by obtaining SPARKLING trajectories for varying values of and by choosing the value which is optimal with respect to robustness to off-resonance effects, while not having k-space holes which leads to poorer reconstruction performance.

[0101] The impact of the temporal weighting on a k-space trajectory is illustrated on FIG. 6, where the shade of gray represent the time at which a specific point of k-space is “visited” by a trajectory. The upper panel represent a 3D SPARKLING trajectory according to the prior art; the temporal inhomogeneity is clearly apparent in the form of strong density contrasts. The lower panel of the figure represents a 3D SPARKLING trajectory according to an embodiment of the invention. It can be seen that the grey level varies smoothly, except on an “equatorial” plane.

[0102] In a second, and preferred, embodiment of the invention,

[00010]$\begin{array}{cc}W\left(.Math.t_i-t_j.Math.\right)=\exp \left(\frac{.Math.t_i-t_j.Math.\tau }{2^{R^{\max }-R}}\right)& \left(13\right)\end{array}$

[0103] The purpose of the weighting

[00011]$\frac{1}{2^{R^{\max }-R}}$

is to shape the amount of temporal repulsion added as a function of the current decimation level, in such a way that the temporal repulsion is stronger at the initial stages of the algorithm and is then substantially reduced near convergence, at finer resolution levels (lower R), in order to prevent the appearance of unwanted k-space holes. Other expressions for W achieving the same result can be envisioned.

[0104] FIGS. 7 and 8 illustrate the technical effect of the invention using a test pattern and the brain of a healthy human volunteer, respectively.

[0105] In FIG. 7, lines 1 and 2 correspond to two different magnetic field configurations, represented by the ΔB.sub.0 maps of column E. The images in column A has been acquired using the SPARKLING method. Significant off-resonance artifacts are visible with strong signal losses. Columns (C)-(E) show the results obtained using the inventive MORE-SPARKLING method, at increasing levels of τ from 1000 to 2000. As the value of τ increases, the impact of the off-resonance is strongly reduced and most of the lost signal in (A) is recovered.

[0106] In FIG. 8, panel B shows axial, sagittal and coronal sections of a 3D ΔB.sub.0 map, where arrows identify the regions of highest inhomogeneity. Panel A shows axial, sagittal and coronal sections of a reference MRI image acquired using a Cartesian sampling scheme, with scan time of 15 minutes and 13 seconds. In this and in the following panels, arrows highlight artifacts, which are identified by the arrows' label: “b” for blurring, “c” for contrast loss and “a” for loss of anatomical details. It can be seen that, even with Cartesian acquisition, magnetic field inhomogeneity leads to localized blurring and the loss of some anatomical details. Panel C corresponds to SPARKLING acquisition with a scan time of 3 minutes and 22 seconds (acceleration by a factor of about 4.5 compared to Cartesian sampling): very strong blurring can be identified. Panel D corresponds to MORE SPARKLING (with τ=1000, temporal weighting function defined by equation 12): it can be seen that the impact of ABo is substantially reduced compared to the case of prior-art SPARKLING, while maintaining the significant acceleration factor compared to the reference Cartesian sampling.

[0107] FIG. 9 illustrates SPARKLING/MORE SPARKLING trajectories generated for different values of the temporal weight parameter τ (N.sub.c=3639; N.sub.s=2048; weighting function W expressed by equation 13). Panel A correspond to the prior art SPARKLING case, and temporal inhomogeneity is clearly visible. For increasing values of τ (panel B: τ=0.2; panel C: τ=0.4; panel D: τ=0.8; panel B: τ=1.0; panel B: τ=1.2), temporal homogeneity increases, but gaps in the k-space begin to appear (highligted b, temporal homogeneity increases, but gaps in the k-space (highligted by arrows) begin to appear. A more detailed analysis shows that increasing τ beyond a value of 0.8 provides minimal benefits in terms of reduction of off-resonance effects, therefore τ≈0.8 is an optimal choice (the exact optimal value depends on the specific setting considered, and on the form of weighting function W).

[0108] FIG. 11 illustrates, very schematically, an apparatus for carrying out a method according to the invention. This apparatus is basically a MRI scanner comprising a primary coil PC for generating the longitudinal field B.sub.0, a single or preferably a plurality of radio-frequency coils TC for generating RF pulses and receiving the NMR signal (different radio-frequency coils may be used for transmission and reception, or the same coils may serve for both purposes), three gradient coils GCx, GCy and GCz for generating magnetic field gradients G.sub.x, G.sub.y, G.sub.z along the x, y and z axes respectively (not represented). The gradient coils are disposed around a scanner bore SB in which a body BD to be imaged can be introduced. The term “body” should be construed broadly, designating any living or non-living object containing atomic nuclei at least some of which are supposed to have nonzero spin. For instance, body BD may be a human or animal body or a portion thereof (e.g. a head) or an article of manufacture such as a “phantom”. The apparatus also comprises a control unit CTR for driving said gradient coils in and said radio-frequency coil or coils according to a predetermined driving signals (reference DSR for the radio-frequency coils and DSG for the gradient coils) to implement a readout sequence (see e.g. FIG. 1) and a k-space trajectory of the MORE SPARKLING type. The apparatus also comprises a signal processing unit SPU for acquiring NMR signals from the radio-frequency coil or coils and performing image reconstruction. The control unit and the signal processing unit may be separate devices or a single device may implement both functions. Said device or devices may include one or more suitably programmed general purpose computers or digital signal processors, suitably configured specialized digital circuits or both. The signal processing unit also includes a receiving circuit for acquiring NMR signals, comprising a signal amplifier, a sampler and an analog to digital (ADC) converter.

[0109] The scanner of FIG. 9 differs from the prior art essentially in that the controlling unit CTR is configured or programmed for driving the gradient coils in such a way that they define a MORE SPARKLING k-space trajectory. For instance, a driving signal DRG defining a MORE SPARKLING trajectory may be stored in a memory of said controlling unit.

[0110] The invention has been described with reference to specific examples, but is not limited to them. For instance, the weighting function may have a form different from that of equation (12).

[0111] The invention applies to 3D imaging, but also to 2D multislice k-space sampling and to 4D imaging where acceleration can be largely increased.

[0112] Trajectory optimization can be initialized starting from any classic k-space filling support (including Cartesian lines, spirals, radial, . . . ), not only radial spokes as in the exemplary embodiments.

[0113] The inventive method may be adapted to any types of MR readout scheme segmented or single-shot, from GRE (cf. the detailed examples above), Turbo FLASH (also called MPRAGE for brain applications) to Spin Echo (SE) or Turbo Spin Echo (TSE), TrueFisp (also called bSSFP), EPI.

[0114] It may also be adapted to any types of MR sequence weighting T1, T2, T2*, ρ (Spin Density), fMRI (functional MRI) or BOLD MRI, and preparation including none-exhaustively, ASL (Arterial Spin Labelling), DWI (Diffusion weighting imaging) with all its variants (DTI, DSI, IVIM, Kurtosis, NODDI), MTR (Magnetization Transfer Ratio) of any type including CEST, quantitative MRI including simultaneous multiparametric techniques such as Quantitative Susceptibility Mapping (QSM), MRF (MR Fingerprinting), MR Angiography both static or dynamic. This includes more exotic MRI applications such as MR thermometry or Electromagnetic Property Tomography (EPT). It may also be applied to heteronuclear imaging such as Sodium or Phosphorus.

[0115] It is compatible with parallel imaging using coil phased array and Simultaneous Multislice technique.

REFERENCES

[0116] [Andersson et al. 2003] J. L. R. Andersson, S. Skare, J. Ashburner “How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging”. Neurolmage, 20(2):870-888, 2003. [0117] [Boyer et al, 2016] Boyer, Claire, et al. “On the generation of sampling schemes for magnetic resonance imaging.” SIAM Journal on Imaging Sciences 9.4 (2016): 2039-2072. [0118] [Chaithya et al. 2020] Chaithya G R et al. “Optimizing full 3D SPARKLING trajectories for high-resolution T2*-weighted MagneticResonance Imaging. 2020 arXiv preprint arXiv:2108.02991 & IEEE Transactions on Medical Imaging, 2022. [0119] [Chauffert et al. 2016] Chauffert, Nicolas et al. “A projection algorithm for gradient waveforms design” in Magnetic Resonance Imaging. IEEE Transactions on Medical Imaging, Institute of Electrical and Electronics Engineers, 2016, 35 (9), pp.2026-2039. [0120] [Daval-Frerot et al. 2021] Daval-Frérot, Guillaume “Off-resonance correction of non-Cartesian SWI using internal field map estimation” International Society for Magnetic Resonance in Medicine, May 2021, Online, United States. [0121] [Haldar, 2014] Haldar, Justin P. “Low-Rank Modeling of Local k-Space Neighborhoods (LORAKS) for Constrained MRI”, IEEE Transactiojn of Medical Imaging 33.3 (2014): 668-680. [0122] [Lazarus et al, 2017] Lazarus, Carole, et al. “SPARKLING: Novel Non-Cartesian Sampling Schemes for Accelerated 2D Anatomical Imaging at 7T Using Compressed Sensing.” 25th annual meeting of the International Society for Magnetic Resonance Imaging, Apr 2017, Honolulu, United States. [0123] [Lazarus 2019] Lazarus, Carole et al. “SPARKLING: variable-density k-space filling curves for accelerated T2*-weighted MRI” Magnetic Resonance in Medicine, Wiley, 2019, 81 (6), pp.3643-3661. [0124] [Lustig et al, 2007] Lustig M., Donoho D., Pauly J. M. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn Reson Med 2007; 58: 1182-1195 [0125] [Sutton et al. 2003] Sutton, Bradley P, et al.“Iterative Image Reconstruction for MRI in the Presence of Field Inhomogeneities”. IEEE Transactions on Medical Imaging 2003;22 (2): 178-88.