Method and apparatus for performing accelerated magnetic resonance imaging with reduced oversampling at low spatial frequencies

Abstract

The invention provides a method of performing magnetic resonance imaging on a body using a magnetic resonance imaging scanner, the method comprising applying to the body a time-varying magnetic field gradient (Gx, Gy, Gz) defining a continuous trajectory (ST) in k-space, the trajectory complying with a set of constraints comprising constraints on a maximum amplitude and a maximum slew rate of the time-varying magnetic field gradient, sampling points (KS) belonging to the trajectory define pseudo-random sampling of k-space approximate to a predetermined target sampling density, the trajectory in k-space minimizing a cost function under the set of constraints, the cost function being defined by a difference between a first term, referred to as an attraction term, and a second term, referred to as a rejection term, the first term promotes the distribution of sampling points in the k-space to follow a predetermined target sampling density, the second term promotes the separation of the sampling points in the k-space, and the exclusion term represents the sum of contributions corresponding to each pair of sampling points; wherein the trajectory (ST, MT) in the k-space defines a Cartesian sample of a central region (S) of the k-space.

Description

TECHNICAL FIELD

[0001] The present invention relates to a method and apparatus for minimizing the adverse effects of oversampling of the central region of k-space (i.e. low spatial frequencies) while performing accelerated magnetic resonance imaging. BACKGROUND

[0002] Magnetic resonance imaging (MRI) is one of the most powerful imaging techniques routinely used in clinical today, but it is still a time-consuming process, especially when large and / or high-resolution images containing millions of pixels need to be acquired. For example, acquiring a three-dimensional image of the human brain with a field of view of 205x205x52 mm3, a resolution of 200 pm using a T2 * weighted sequence with an acceptable signal-to-noise ratio (SNR) of 7.6 can require an acquisition time of about three hours, which is clearly unacceptable for clinical purposes.

[0003] This is because MRI images are obtained by sampling the so-called "k-space" domain (the spatial frequency spectrum of the image), which is related to physical space by a multi-dimensional Fourier transform. Sampling theory teaches that the sampling frequency should be at least twice the highest frequency contained in the signal to be sampled (in the case of MRI, this is the multi-dimensional Fourier transform of the body to be imaged), otherwise aliasing artifacts will occur; this is known as the Nyquist or Nyquist-Shannon criterion. Therefore, using a traditional acquisition scheme, the number of k-space samples must be at least equal to the number of pixels of the image to be reconstructed. Moreover, SNR requirements set a minimum acquisition time for each sample.

[0004] Several techniques have been developed to reduce the acquisition time while avoiding artifacts ("accelerated MRI").

[0005] Some of these techniques, such as simultaneous multi-slice imaging (SMS) and parallel MRI, require the use of dedicated hardware containing multiple receiving coils to acquire the magnetic resonance signal. Therefore, their implementation is expensive. Moreover, they provide limited acceleration, since the image quality rapidly decreases with the acceleration factor. Even if the two techniques are used simultaneously, the combined acceleration factor does not exceed 8 in practice.

[0006] Other techniques are compatible with the use of a single receiving coil (although multiple coils can also be used), such as partial Fourier imaging, exploiting the redundancy in the k-space information, non-Cartesian k-space filling such as radial or spiral, and compressed sensing (CS). While partial Fourier techniques provide only very limited acceleration factors (typically lower than two), compressed sensing allows higher acceleration, especially when imaging with high matrix size (high resolution and small field of view or low resolution and large field of view), see [Haldar, 2014].

[0007] ​A review of CS MRI can be found in [Lustig et al., 2008].

[0008] The CS technique is based on three principles: - The image to be reconstructed must allow a sparse (or compressible) representation. In other words, it must be possible to decompose it on a predetermined basis (e.g. a wavelet basis) such that, for a strict sparsity, only a small fraction of the decomposition coefficients are non-zero; or for compressibility, only a small fraction of the decomposition coefficients are significantly larger than zero. Typically, in the case of a noisy signal, a coefficient is considered significantly larger than zero if its absolute value is at least equal to the noise standard deviation. Alternatively, only a predetermined fraction of the coefficients can be kept - those with the largest absolute values. For example, only the first 1% of the coefficients can be kept, resulting in a compression factor of 100.

[0009] - The reconstruction must be performed using a non-linear method to improve the sparsity of the image representation and to be consistent with the acquired samples.

[0010] - In order to accelerate the acquisition, the k-space must be undersampled in a pseudo-random pattern. The undersampling reduces the number of signal acquisitions and thus provides the required acceleration, while the pseudo-randomness ensures that, in the sparse representation, the subsampling artifacts are incoherent, i.e. de-correlated or noise-like. This incoherence is extremely important, as it measures the degree of correlation between any pair of elements taken in the sparse basis (e.g. wavelets) and the sensing basis (e.g. Fourier transform in CS-MRI).

[0011] CS usually uses non-Cartesian pseudo-random sampling of the k-space, as this strategy provides an improved incoherence (lower correlation between samples). Preferably, the pseudo-random sampling is non-uniform. This is because, in a sparse basis like the wavelet transform, the image decomposition consists of approximation coefficients and detail coefficients. The approximation coefficients encode the low-frequency information and are thus not sparse, while the detail coefficients encode the high-frequency content and are sparse. The sparsity can be imposed to recover the detail coefficients, while the information related to the approximation coefficients (low-frequency content) must be collected in the data. Therefore, we must use a variable density sampling approach, which acquires more measurements in the low frequencies than in the high frequencies. Another view is that variable density sampling (VDS) is a theoretical means to minimize the number of samples collected, thus guaranteeing accurate image reconstruction with high probability. In clinical applications, this usually means using a variable sampling density that is highest near the center of the k-space (low spatial frequencies) and decreases towards high spatial frequencies.

[0012] From a purely theoretical point of view, pseudo-random sampling could be achieved by drawing the sampling points according to a predefined probability distribution, which corresponds to the desired sampling density. In practice, however, the short lifetime of the MR signal requires that the samples are acquired by piecewise acquisition along a smooth trajectory, which is defined by time-varying magnetic field gradients applied to the body to be imaged after the nuclei spins of the body to be imaged have been excited by means of radio frequency (RF) pulses.

[0013] Let denote the magnetic field gradients applied to the body to be imaged. These magnetic field gradients define a trajectory in k-space, which is represented by (1) The sampling is performed by acquiring the MRI signal generated by the excited nuclei spins at predetermined times, which correspond to points along the trajectory.

[0014] Due to hardware limitations and physiological limitations in clinical applications, the gradient field amplitudes and their slew rates cannot exceed respective limits Gmaxand Smax. Therefore, only sufficiently regular trajectories are allowed.

[0015] When only the nuclei spins within a thin slice of the body are excited, these trajectories can be two-dimensional in 2D MRI; in 3D MRI techniques, the trajectories can be three-dimensional, where the excitation involves the entire body or a thick slab thereof. In the following, for the sake of simplicity, only the case of 2D trajectories is considered; however, the present invention is also applicable to 3D and even 4D (i.e. dynamic) MRI.

[0016] An important property of MRI is that the signal decays exponentially after the application of an excitation radio frequency pulse, typically vanishing. This limits the duration of the signal acquisition, and thus the length of each individual k-space trajectory. Therefore, several excitation RF pulses are required, each followed by MRI signal acquisition along a respective base trajectory (or "shot") to perform a complete k-space sampling. The repetition time TRof these excitation RF pulses, which imposes an upper limit on the duration of the signal acquisition, also depends on the imaging technique used, in particular the type of contrast sought (T1, T2, T2 * ……).

[0017] Commonly used k-space trajectories are parallel lines (resulting in Cartesian sampling), spokes (straight lines diverging radially from the center of k-space), rosettes, uniform and variable density spirals. All these trajectories have been applied to compressed sensing, e.g. by performing a limited number of signal acquisitions on a Cartesian grid, or by randomly sampling spokes, spirals or rosettes.

[0018] However, better results can be obtained by using "nonparametric" trajectories that provide greater incoherence.

[0019] In particular, the so-called SPARKLING technique (“Extended Projection Algorithm for Fast K-space Sampling”) ([Boyer et al., 2016], [Lazarus et al., 2017], [Lazarus et al., 2019], [Chaithya et al., 2020], WO2019 / 048565) is based on projecting a predetermined, typically non-uniform, target sampling distribution (TSD) onto a set of “acceptable” 2D or 3D curves, i.e., onto all curves representing trajectories obtainable without exceeding the magnetic field gradient value and the corresponding conversion rate limit.

[0020] [Chaithya et al., 2022] disclosed an improvement to SPARKLING, called MORE SPARKING.

[0021] A drawback of some non-Cartesian MRI sampling schemes (such as Sparkling) is that they involve significant oversampling of the central portion of the k-space, at sampling rates far exceeding the Nyquist sampling rate—sometimes averaging up to 20 k-space samples per Nyquist voxel. Indeed, while high spatial frequency components (the “details” in wavelet decomposition) can be undersampled and then recovered via sparsity constraints—which is indeed central to CS-MRI—low spatial frequencies (the “approximations” in wavelet decomposition) must be fully sampled to ensure perfect recovery. However, in most non-Cartesian CS-MRI implementations, including Sparkling, each shot is constrained to pass through the center of the k-space at echo time TE (echo time constraint, see [Chauffert et al., 2016]), particularly for good target contrast and to allow correction of motion artifacts, resulting in oversampling, which is disadvantageous for at least two reasons. On the one hand, these additional k-space samples could alternatively be used to sample higher frequencies. On the other hand, since k-space data is obtained by integration on different k-space trajectories, such repeated sampling of the same Nyquist voxel through trajectories of different paths may lead to the accumulation of B0 artifacts.

[0022] Partially, the solution to this problem involves updating the TSD by incorporating TE constraints and reducing the density around the center of the k-space to limit the number of samples per Cartesian k-space voxel. However, this does not completely alleviate the problem, as there are still a few samples at the center of the k-space (one per fire), resulting in dense sampling of the surrounding region.

[0023] Even with the removal of the echo-time constraint (which is possible for sufficiently low scan times, albeit at the cost of reduced target contrast), oversampling can still occur if the firing mechanism gets stuck in a local minimum. Furthermore, removing the echo-time constraint can lead to center sampling in the k-space at different times, thus compromising contrast and increasing B0 artifacts. Summary of the Invention

[0024] The objective of this invention is to overcome, in whole or in part, the aforementioned deficiencies of the prior art. More specifically, this invention aims to provide a method for accelerating MRI by ensuring that the central portion of k-space is sampled according to or close to the Nyquist criterion, thereby avoiding both oversampling and undersampling.

[0025] According to the present invention, this objective is achieved by modifying the SPARKLING trajectory by introducing appropriate affine constraints to ensure fundamental Cartesian sampling of the central region of the k-space. This method is called GoLF-SPARKLING, as it stands for Low-Frequency Meshization (GoLF) SPARKLING.

[0026] One object of the present invention is a method for performing magnetic resonance imaging on a body using a magnetic resonance imaging scanner, the magnetic resonance imaging scanner including a scanner aperture, a primary coil, a radio frequency coil, a gradient coil, and a signal processing unit, the method comprising the following steps: a. Position the body in the scanner aperture, wherein a static and substantially uniform magnetic field oriented along a direction called the longitudinal direction is established by the primary coil, and this magnetic field is called the longitudinal field; b. Using all or part of the radio frequency coil, at least one radio frequency pulse is emitted toward the body, the radio frequency pulse being adapted to excite nuclear spins within the body; c. After each of the aforementioned radio frequency pulses, a time-varying magnetic field gradient defining a trajectory (ST) in k-space is applied to the body using the gradient coil, and simultaneously, samples of the magnetic resonance signal emitted by the excited nuclear spin are acquired using all or part of the radio frequency coil, each sample corresponding to a point (KS) in k-space belonging to the trajectory; and d. Using the signal processing unit, a nonlinear reconstruction algorithm is applied to the acquired sample to reconstruct the magnetic resonance image of the body; The trajectory in the k-space is a continuous trajectory that conforms to a set of constraints, including constraints on the maximum amplitude and maximum conversion rate of the time-varying magnetic field gradient. Thus, the point in the k-space corresponding to the sample, referred to as a sampling point, defines a pseudo-random sampling of the k-space, approximating a predetermined target sampling density. The trajectory in the k-space minimizes a cost function under the constraints described above. This cost function is defined by the difference between a first term, referred to as an attraction term, and a second term, referred to as a repulsion term. The first term promotes the distribution of sampling points in the k-space to follow the predetermined target sampling density, and the second term promotes separation between sampling points in the k-space. The repulsion term is expressed as the sum of contributions corresponding to each sampling point. The characteristic feature is that the trajectory in the k-space defines a Cartesian sample of the central region of the k-space.

[0027] Another object of the present invention is a computer-implemented method for calculating trajectories in k-space for magnetic resonance imaging sample acquisition, the method comprising the following steps: i. Determine or receive the number N of fires constituting the trajectory from user input. c ≥1 and the number N of sampling points along each of the firing points s >1; ii. Determine or receive a predetermined target sampling density from user input; iii. Determine or receive a set of constraints from user input, including constraints on the magnitude of the discrete-time derivative of the trajectory; iv. The trajectory is computed by minimizing a cost function under said set of constraints, the cost function being defined by the difference between a first term called the attraction term and a second term called the repulsion term, the first term promoting the distribution of sampling points in the k-space to follow the predetermined target sampling density, the second term promoting separation between sampling points in the k-space, and the repulsion term being expressed as the sum of contributions corresponding to each pair of sampling points; The characteristic feature is that the trajectory in the k-space defines a Cartesian sample of the central region of the k-space.

[0028] Another object of the present invention is a computer program comprising instructions that, when executed by a computer, cause the computer to perform the method of calculating a trajectory in k-space for magnetic resonance imaging sample acquisition.

[0029] Furthermore, another object of the present invention is a set of driving signals for gradient coils of a magnetic resonance imaging scanner, which, when applied to the gradient coils, drive the gradient coils to generate a time-varying magnetic field gradient. This time-varying magnetic field gradient defines a trajectory in k-space that passes through a plurality of points called sampling points, which define pseudo-random sampling in k-space, approximating a predetermined target sampling density; wherein the trajectory: - It is continuous and conforms to a set of constraints, which include constraints on the maximum amplitude and maximum conversion rate of the time-varying magnetic field gradient; - Minimize the cost function under the aforementioned set of constraints, which is defined by the difference between a first term called the attraction term and a second term called the repulsion term. The first term promotes the distribution of sampling points in the k-space to follow the predetermined target sampling density, and the second term promotes the separation between sampling points in the k-space. The repulsion term is expressed as the sum of the contributions corresponding to each pair of sampling points. The characteristic feature is that the trajectory in the k-space defines a Cartesian sample of the central region of the k-space.

[0030] Another object of the present invention is a magnetic resonance imaging scanner, comprising: - Scanner aperture, the body can be positioned within the scanner aperture; - A primary coil, configured to establish a static and substantially uniform magnetic field in the scanner aperture along a direction known as the longitudinal direction, referred to as the longitudinal field; - A radio frequency coil configured to emit at least one radio frequency pulse toward the body, the radio frequency pulse being adapted to excite nuclear spins within the body, and the radio frequency coil being configured to acquire samples of magnetic resonance signals emitted by the excited nuclear spins; - A gradient coil configured to apply a time-varying magnetic field gradient to the body, the time-varying magnetic field gradient defining a trajectory in k-space; - A signal processing unit that applies a nonlinear reconstruction algorithm to samples of the acquired magnetic resonance signals to reconstruct a magnetic resonance image of the body; and - A controller configured to generate a drive signal for the RF coil and a drive signal for the gradient coil; The controller is configured to generate the drive signal for the gradient coil to produce the time-varying magnetic field gradient, such that the trajectory in k-space passes through a plurality of points called sampling points, which define a pseudo-random sampling of the k-space approximating a predetermined target sampling density. The trajectory is continuous and conforms to a set of constraints, including constraints on the maximum amplitude and maximum conversion rate of the time-varying magnetic field gradient. The trajectory in k-space also minimizes a cost function under the constraints, defined by the difference between a first term called an attraction term and a second term called a repulsion term. The first term promotes the distribution of sampling points in k-space to follow the predetermined target sampling density, and the second term promotes separation between sampling points in k-space. The repulsion term is expressed as the sum of contributions corresponding to each pair of sampling points. The controller is also configured to generate drive signals for the radio frequency coils, such that they acquire samples of the magnetic resonance signals emitted by the excited nuclear spin, the samples corresponding to the sampling points of the trajectory in k-space; The characteristic feature is that the trajectory in the k-space defines a Cartesian sample of the central region of the k-space.

[0031] Another object of the present invention is a non-transient computer-readable information storage support—for example, magnetic storage support such as hard disks or magnetic tapes; optical storage support such as CDs or DVDs; or semiconductor non-volatile memory such as flash memory or solid-state drives (SSHDs)—to store the computer program and / or the set of drive signals. Attached Figure Description

[0032] Additional features and advantages of the invention will become apparent from the following description taken in conjunction with the accompanying drawings, which show: - Figure 1 A timing diagram for implementing a modified GRE (gradient echo) pulse sequence of a SPARKLING trajectory according to the prior art or a GoLFSPARKLING trajectory according to an embodiment of the present invention; - Figure 2 Example 2D Sparkling trajectory based on existing technology; - Figure 3A and 3B The figure shows an oversampling plot of low spatial frequencies obtained using a SPARKLING trajectory according to the prior art; - Figure 4 Examples of Golf-Sparkling trajectories according to various embodiments of the present invention; - Figure 5 and 6 The experimental results obtained on the phantom demonstrate the technical effects of the present invention; - Figure 7 and 8 The experimental results obtained on healthy volunteers also demonstrate the technical effects of the present invention; and - Figure 9 A block diagram of an MRI scanner according to an embodiment of the present invention. Detailed Implementation

[0033] Figure 1 The pulse sequence is shown—more specifically, for T2. * A time-series diagram of a gradient echo (GRE) sequence for weighted anatomical imaging, suitable for implementing a SPARKLING k-space trajectory or a GoLF-SPARKING trajectory according to existing techniques.

[0034] Suppose the body to be imaged (e.g., a patient's head) is placed in a static and uniform magnetic field B0 oriented along the "longitudinal" z-direction. This magnetic field induces partial alignment of the nuclear spins of the body. Depending on the value of the magnetic field B0 and the value of the nuclear magnetomotor, the aligned nuclear spins can be excited (i.e., flipped) by a radio frequency pulse RFP at an appropriate frequency (Ramohr frequency). In the 2D sequence, as shown in the figure above, a magnetic field gradient G is also applied. z That is, a magnetic field along the z-direction, the amplitude of which varies linearly along the same direction. This has the effect of making the Larmor frequency variable along z; therefore, only the nuclear spins within the thin slices of the bulk resonate with the RF pulse and can be excited. As is known in the field of MRI, this "slice-selective" gradient pulse is followed by a brief G... z The negative abrupt change in the magnetic field gradient (“refocusing gradient”) compensates for the dispersion of nuclear spin orientation in the my plane perpendicular to the z-direction.

[0035] Due to the use of the layer selection gradient G z Therefore, only 2D images of the selected body slices are acquired, which requires 2D k-space representation. a k i Sampling is performed on a plane; in the following text, the term "k-space" will be used to specify the three-dimensional k-space. x k y k z The space and the two-dimensional plane within it. In an alternative embodiment, instead of using slicing to select gradients, it is necessary to perform 3D k... x k y k z Sampling is performed on the spatial domain.

[0036] k-space (more precisely, 2D k-space) x k y The trajectory in the plane is obtained by playing G after the RF excitation pulse ends. x and G yDefined by gradient. It's important to emphasize that the applied magnetic field is always oriented along the z-direction, but its magnitude shows a linear variation along the x and y directions. First, apply G... x and G y Pulse, to reach the k to be sampled x k y A suitable point on the boundary of the planar region. Then play "arbitrary" G. x and G y Waveform, defined by its overall radial orientation and direction k x k y The winding trajectory at the center of the plane.

[0037] In a “3D” implementation without layer selection gradients, a G-like approach will be applied. x and G y G z Waveform, to reach k x k y k z Find suitable points on the boundary of the 3D domain in space and define a three-dimensional trajectory in k-space.

[0038] While the gradient waveform is being played, samples of MRI signals emitted by the stimulated nucleus are acquired by one or more radio frequency coils connected to suitable acquisition circuitry, including an analog-to-digital converter (ADC). The acquisition period is determined by... Figure 1 The OP symbol in the attached figure indicates the duration T of the acquisition cycle. obs Limited by MRI signal attenuation. At the end of the acquisition cycle, a final G-force is applied. x Gradient pulses (“scratching gradients”) are applied to eliminate any remaining transverse magnetization. The repetition time TR corresponds to the interval from exactly before the radio frequency pulse RFP to the end of the scratching gradient pulse.

[0039] The sequence is repeated several times with different gradient waveforms, each defining its own k-space trajectory. These k-space trajectories together provide the required k-space sampling. The set of excitation RF pulses and associated gradient waveforms is called a "shot"; each shot corresponds to a basic trajectory. Figure 2 Displays 2D "multi-fire" SPARKLING trajectory MT, including N c =60 basic trajectories or firing STs with radial orientation, for k x k y Non-uniform sampling is performed on the plane.

[0040] In fact, the magnetic field gradient G x G y G zThe gradient changes in discrete time points, separated by intervals of duration Δt (“gradient grating time”). Sampling is also performed at fixed intervals of duration dt (“ADC dwell time”). The ADC dwell time δt is preferably less than or at most equal to the gradient grating time (δt ≤ Δt) to collect several samples between two consecutive gradient steps. However, reducing the ADC dwell time beyond a certain limit will reduce the SNR to an unacceptable level. Therefore, for each specific embodiment of the invention, an optimal value for δt can be found.

[0041] Let k i (1) is the location of the starting point of the trajectory associated with the i-th firing in k-space. This corresponds to the first sample of the MRI signal acquired at this point. Other sampling points correspond to k-space locations given by the following formula: (2) Where m∈[2:M] is an integer index, M is the total number of samples collected along the trajectory, and q and r are the modulus and remainder of the Euclidean division of the collection time, respectively: (3) If δt = Δt, then r = 0, and the number of ADC samples matches the number of gradient time steps. If δt < Δt, then the number of ADC samples is greater than the number of gradient time steps.

[0042] “SPARKLING” relies on an optimization-based approach that involves projecting a target sampling distribution onto a set of discrete forward measurements, particularly with the support of a smooth trajectory [Lazarus et al., 2017; Boyer et al., 2016]. Mathematically, the problem can be transformed into a nonconvex variational optimization problem under possible nonconvex constraints [Boyer et al., 2016]. These constraints are typically represented by the maximum acceptable values ​​of gradient magnitude and transformation rate, but other affine constraints can also be used—for example, forcing the trajectory to pass through specific points in k-space at defined time points (e.g., for echo time definitions).

[0043] According to (Boyer et al., 2016; Lazarus et al., 2017; Lazarus et al., 2019; and Chaithya et al., 2020), the SPARKLING trajectory is computed by minimizing a cost function defined by the difference between two terms under a set of constraints: an "attraction term" that causes the distribution of sampled points in k-space to follow a predetermined (usually non-uniform) target sampling density; and a "repulsion term" that causes separation between sampled points in k-space. This can be written as: (4) in (5) As an attraction, (6) It is an exclusion term, and - N is the number of sampling points along the trajectory. c ≥1 represents the number of shots or basic trajectories that constitute the entire trajectory, N s >1 represents the number of sampling points per firing; - Q Nc It is a set of curves in k-space, including N c Each firing or segment comprises N fires or segments. s There are points that satisfy the constraint; - K= It is N c =k i (t) The k-space sampling mode composed of firing parameters; - Ω is the sampling region in k-space, normalized to [-1, 1] 3 (For 3D trajectories); - x represents a general point in the sampling region Ω; - ρ is the predetermined target sampling density; - It is a norm, preferably the L2 norm.

[0044] This set of constraints can be defined as follows: (7) in: - These are the first and second discrete-time derivatives of the trajectory in k-space, respectively; - For any three-dimensional vector c, ; - in; and , where j∈{1, …c i}, c i It represents the number of affine constraints on the i-th k-space firing, and this constraint is modeled. More precisely, A i se - α and β are normalized parameters representing the hardware constraints of the gradient magnitude and gradient transformation rate; and - This represents Ns points in Ω.

[0045] Typically, the trajectory design process begins with initializing the trajectory (e.g., by N...). c Starting with radial spoke formation, the optimization problem (4) is then solved under constraint (7) using the projective gradient descent algorithm, as discussed in [Chaithya et al., 2020].

[0046] Figure 2 Each individual shot of the trajectory passes through the center of the k-space, resulting in sampling of the central (low-frequency) region of that space at a rate much higher than the Nyquist rate. For example, as... Figure 3A As shown, N c =4096, N s At δ = 2048 and δt = 2µs, each Nyquist voxel has almost 20 k-space samples.

[0047] A Nyquist voxel is a voxel defined on a Cartesian grid in k-space with a size Δk = 1 / FOV, where FOV is the field of view in meters (e.g., 20 cm = 0.2 m). Collecting samples of each Nyquist voxel allows for linear reconstruction in the central region of k-space without aliasing.

[0048] As mentioned above, this low spatial frequency oversampling is suboptimal because these additional k-space samples can be alternatively used to sample other higher frequencies. Furthermore, since the k-space data is acquired as an integral over different k-space trajectories with varying signal amplitudes, resampling the same Nyquist voxel with different trajectories can lead to the accumulation of B0 artifacts. This can be... Figure 3B As seen in the image, the mean and standard deviation of the k-space data are displayed.

[0049] According to the present invention (“GoLF-SPARKLING”), the problem of low spatial frequency oversampling is solved by using modified SPARKLING class trajectories that define Cartesian sampling in the central region of k-space.

[0050] exist Figure 4 In the exemplary embodiment of the present invention shown, the radius is r s A sphere S is defined at the center of k-space. The entire k-space trajectory (part A in the figure) is divided into a non-Cartesian part (e.g., a sparkling trajectory) outside the sphere (part B) and a Cartesian part (part C). In the example shown, N c =256 and N x =N y =Nz =Ñ=64. The internal region of space k inside sphere S is formed by the region along k. x Sampling is done using straight lines with directional orientation, and these lines are based on k. y k z Cartesian grid arrangement in a plane. Figure 4 Parts D and E respectively show the sphere S according to k y k z and k x k z A cross-section of a plane. For each dimension, the Cartesian grid has a pixel size of 2 / ε in k-space.

[0051] Since each k-space firing forms a Cartesian line, the area of ​​the Cartesian sampling circle in the central tangent layer is given by the following formula. (8) From this we can obtain (9) Based on the Nyquist criterion, to obtain an aliasing-free reconstructed image, at least each Sampling is performed. However, since k-space data is sampled at every dwell time δt < Δt, the scanner can actually play the k-space trajectory to sample at least once every Δt. In summary, the k-space velocity at the center of k-space can be introduced as a dimensionless parameter v≥0, which is the velocity of the gradient grating time Δt in the readout direction k. x The Nyquist voxel step size Δk performed on x The number. Specifically, the step size taken by the k-space trajectory at each Δt is... The goal of v is to efficiently utilize the gradient hardware to traverse the center of the k-space at the most feasible speed, while maintaining the Nyquist criterion after sampling at the ADC at each δt. Taking a higher v value can improve k-space coverage, but it will reduce the signal-to-noise ratio because the MR signal accumulates over shorter time intervals (see [reference]). Figure 6 (And the corresponding discussion below). Therefore, a trade-off is necessary. For example, for T2 of 3T. * For imaging, v≤1 is preferred.

[0052] Without loss of generality, suppose the i-th k-space is fired k i The Cartesian part as a k-space line Starting at and in End of section (see) Figure 4 The number of samples at the center of the k-space (i.e., within the sphere s) of the firing (part E) is given by the following formula. (10) The k-space locations of these Nyquist points are given by the following equation. (11) This led to Constraint, the constraint in and Between The index is applied to the i-th spatial firing k. i Where "L" and "H" represent "low" and "high" respectively. The constraint is then expressed as: (12) (13) The affine constraint set of equations (12) and (13) allows Cartesian sampling over the central region of the k-space S. However, in order to force the Cartesian sampling at the Nyquist rate, a suitable target sampling density (TSD) must be chosen, which differs between the central region and the outside of the k-space.

[0053] We consider a normalized sampling domain Ω = [-1, 1]. 3 According to an exemplary embodiment, the target sampling density can be parameterized using a cutoff value C and an attenuation value D, as follows: (14) If N=N c xN S If the total number of sampling points is , then by The number of samples in the central region S of the defined k-space is (15) Furthermore, the number of Nyquist points in the central region of k-space is given by the following formula. (16) If the trajectory velocity defined above is given by v, then in the central region of k-space... From the previous equation, we get: (17) Non-Descartes (N) nc The number of Nyquist samples compared to Cartesian samples (N) c The ratio of the number of Nyquist samples to the number of samples is given by the reciprocal of the corresponding volume ratio: (18) From (18), and considering that the density required for Cartesian sampling is ,available (19) However, in practice, Poisson disk sampling cannot be implemented due to limitations in the velocity and acceleration of the trajectory. To prevent any k-space holes, it is preferable to sample the non-Cartesian region with a period of 1.5 times the Nyquist period, therefore taking... (20) Finally, due to It is a distribution that must be normalized, therefore its sum must equal 1: (twenty one) By substituting (19) and (20) into (21), D can be solved iteratively.

[0054] Figure 5 A comparison between the left portion (MRI image of the test pattern obtained using the prior art 3D SPARKLING trajectory - NIST phantom) and its right portion (MRI image of the same pattern obtained using the GoLF-SPARKING trajectory) shows that the present invention can achieve significantly better flat contrast.

[0055] GoLF-SPARKLING can be advantageously combined with MORE-SPARKLING techniques, as described in European Patent Application 22305592.2 filed by the applicant on April 21, 2022, and [Chaithya et al., 2022]. This technique modifies the SPARKLING exclusion term by introducing a weight W in each individual term, which varies with the time interval of sampling points along the trajectory in k-space. And increase: (twenty two) Preferably, W is Exponential function: (twenty three) Where τ≥1 is a scalar parameter.

[0056] The weight W ensures that the resulting k-space trajectory is temporally uniform, meaning that regions close to each other in k-space are sampled at similar times. The value of the parameter τ must be chosen (e.g., through grid search) to force temporally smooth k-space sampling while avoiding unwanted k-space holes. Combining MORE-SPARKLING with Golf-SPARKLING (MORE+GolfSPARKLING) does not present any particular difficulty, as the first modifies the cost function of the SPARKLING trajectory design algorithm, and the second modifies the affine constraints and TSD to which the trajectory is projected.

[0057] Figure 6 Experimental results for the MORE+GoLF-SPARKLING trajectory at different values ​​of dimensionless k-space velocity are presented: v=0.6 (Fig. i), 1 (ii), 2 (iii), and 3 (iv). The NIST phantom can help optimize the tuning of v because it embodies a resolution inset, which can be used to quantify image quality. For comparison, results are also shown for Cartesian sampling (A), "pure" SPARKLING with τ=0 (B), MORE SPARKLING with τ=1 (C), and GoLF SPARKLING with v=1 (D) using the conventional "GRAPPAp4" k-space trajectory. For the MORE-SPARKLING trajectory, in the context of... Some artifacts were observed near the regions where the value of v changed (marked by arrows in Figure B). These artifacts were significantly reduced using GoLF (Figure D). Furthermore, in GoLF+MORE-SPARKLING (Figures I–iv), it can be seen in the magnified resolution insets that using v>1 results in a significant loss of detail. This careful analysis suggests that v≤1 is more suitable for 3T. Imaging. Within this range, for a wider k-space coverage, the preferred choice is v=1, that is, to spend the shortest time at the center of the k-space at the highest possible speed and collect more samples in the high-frequency region, thereby producing sharper edges.

[0058] Figure 7 In vivo results are shown for “pure” SPARKLING (C) with τ=0, MORE SPARKLING (D) with τ=1, GoLF SPARKLING (E) with v=1, and MORE+GoLF SPARKLING (F) with τ=1 and v=1. In all cases, an acceleration factor of AF=15 was used, resulting in a scan time of 3 minutes and 21 seconds. Figure A corresponds to the conventional “GRAPPA p4” k-space trajectory with Cartesian sampling, and Figure B corresponds to the ΔB0 field map. As expected, MORE-SPARKLING showed greater robustness to ΔB0 inhomogeneities than classic SPARKLING and exhibited minimized artifacts. GoLF+MORE-SPARKLING provided significantly improved image quality with significantly reduced noise levels and contrast closer to the Cartesian reference, paving the way for improved clinical applications of SPARKLING trajectories. However, the GoLF+MORE-SPARKLING trajectory is slightly less robust to ΔB0 inhomogeneity than the MORE-SPARKLING trajectory, and... Figure 7 Some residual artifacts can be seen in F.

[0059] Figure 8 The in vivo results for the MORE-SPARKLING (column A) and GoLF+MORE-SPARKLING (column B) trajectories at acceleration factors AF=10 (row 1), 15 (row 2), 20 (row 3), and 40 (row 5) are shown. The aim of this experiment was to evaluate how image quality changes over scan time. Overall, the GoLF feature provides less noise and more detail compared to the MORE-SPARKLING trajectory alone. In fact, for the GoLF+MORE-SPARKLING trajectory, image quality is maintained up to AF=15, while for the MORE-SPARKLING trajectory alone, a decrease in image quality is observed at AF=15. Furthermore, a direct diagonal comparison can be made between the two methods: the image quality of the GoLF+MORE-SPARKLING trajectory at AF=15, 20, 30, and 40 is comparable to that of the MORE-SPARKLING trajectory at AF=10, 15, 20, and 30, respectively.

[0060] Figure 9 The apparatus for carrying out the method of the present invention is shown very schematically. The apparatus is essentially an MRI scanner, comprising a primary coil PC for generating a longitudinal field B0, a single or preferably multiple radio frequency coils TC for generating RF pulses and receiving MRI signals (different RF coils may be used for both transmission and reception, or the same coil may be used for both purposes), and magnetic field gradients G for generating magnetic field gradients G along the x, y, and z axes, respectively. x G y G z Three gradient coils, GCx, GCy, and GCz (not shown), are arranged around the scanner aperture SB, into which the body BD to be imaged can be introduced. The term "body" should be interpreted broadly to refer to any biological or non-biological object containing atomic nuclei, wherein at least some of the nuclei should have non-zero spin. For example, the body BD can be a human or animal body or a portion thereof (e.g., a head), or an article such as a "phantom." The device also includes a control unit CTR for driving the gradient coils and one or more of the radio frequency coils according to predetermined drive signals (a reference DSR for the radio frequency coils and a DSG for the gradient coils) to implement a readout sequence (see, for example, [link to relevant documentation]). Figure 1The device also includes a k-space trajectory of the Golf-Sparkling type. It further includes a signal processing unit (SPU) for acquiring MRI signals from one or more radio frequency coils and performing image reconstruction. The control unit and the signal processing unit can be separate devices or a single device capable of implementing both functions. The device or devices may include one or more appropriately programmed general-purpose computers or digital signal processors, appropriately configured dedicated digital circuitry, or both. The signal processing unit also includes receiving circuitry for acquiring MRI signals, which includes a signal amplifier, a sampler, and an analog-to-digital (ADC) converter.

[0061] Figure 9 The main difference between this scanner and existing technologies is that the control unit CTR is configured or programmed to drive gradient coils in such a way that they define the Golf-Sparkling k-space trajectory. For example, the drive signal DRG that defines the Golf-Sparkling trajectory can be stored in the memory of the control unit.

[0062] The present invention has been described with reference to specific examples, but is not limited thereto. For example, the predetermined sampling density may be different from the predetermined sampling density represented by equations (14)-(21).

[0063] This invention is applied to 3D imaging, but also to 2D multi-layer k-space sampling and 4D imaging, where acceleration can be greatly increased.

[0064] Trajectory optimization can be initialized from any classic k-space filling support (including Cartesian lines, helices, radial lines, etc.), not just radial spokes as in the exemplary embodiment.

[0065] The shape of the central region of the sampled Cartesian k-space does not have to be spherical or elliptical, even if these constitute the preferred embodiment. It can have other shapes, such as cubes or cuboids.

[0066] The method of this invention is applicable to any type of segmented or single-shot MR readout scheme, from GRE (see detailed examples above), Turbo FLASH (also known as MPRAGE for brain applications) to spin echo (SE) or fast spin echo (TSE), TrueFisp (also known as bSSFP), and EPI.

[0067] It is also applicable to any type of MR sequence weighted T1, T2, T3 *MRI can be performed using p (spin density), fMRI (functional MRI), or BOLD MRI, and preparations including but not limited to ASL (arterial spin labeling), DWI (diffusion-weighted imaging) and all its variants (DTI, DSI, IVIM, kurtosis, NODDI), any type of MTR (magnetic transmissibility) (including CEST), and quantitative MRI (including simultaneous multiparameter techniques such as quantitative sensitivity mapping (QSM), MRF (MR fingerprinting), static or dynamic MR angiography). This includes more exotic MRI applications such as MR thermography or electromagnetic property computed tomography (EPT). It can also be applied to heteronuclear imaging, such as sodium or phosphorus.

[0068] It is compatible with parallel imaging using coil phased arrays and simultaneous multi-slice techniques.

[0069] Depending on the implementation, Cartesian sampling can be performed in one or more regions of the k-space, rather than its central region.

[0070] References [Andersson et al., 2003] JLR Andersson, S. Skare, J. Ashburner, “How to correct magnetic susceptibility distortion in spin echo echo plane images: an application in diffusion tensor imaging.” Neuroimaging, 20(2):870-888, 2003.

[0071] [Boyer et al., 2016] Boyer, Claire et al., “On the generation of magnetic resonance imaging sampling schemes.” Journal of SIAM Imaging Sciences, 9.4 (2016): 2039-2072.

[0072] [Chaithya et al., 2020] Chaithya GR et al., “For High Resolution T2” * "Optimized Full 3D Sparkling Trajectory for Weighted Magnetic Resonance Imaging". 2020 arXiv preprint arXiv:2108.02991.

[0073] [Chaithya et al., 2022] Chaithya GR et al., “MORE-SPARKLING: Non-Cartesian Trajectories with Minimized Off-Resonance Effects”, Proceedings of the ISMRM-ESMRMB 2022 Joint Annual Meeting and ISMRT Annual Meeting, London, UK, 7-12 May 2022.

[0074] [Chauffert et al., 2016] Chauffert, Nicolas et al., “Projection algorithm for gradient waveform design in magnetic resonance imaging”, IEEE Transactions on Medical Imaging, IEEE Transactions on Medical Imaging, 2016, 35(9), pp. 2026-2039.

[0075] [Daval-Frerot et al., 2021] Daval-Frérot, Guillaume, “Estimation of Partial Resonance Correction for Non-Cartesian SWI Using Internal Field Maps”, International Society for Magnetic Resonance in Medicine, May 2021, online, USA.

[0076] [Haldar, 2014] Haldar, Justin P. “Low-rank modeling (LORAKS) for local k-space neighborhoods in constrained MRI”, IEEE Transactions on Medical Imaging, 33.3 (2014): 668-680.

[0077] [Lazarus et al., 2017] Lazarus, Carole et al., “SPARKLING: A novel non-Cartesian sampling scheme for accelerating 2D anatomical imaging at 7T using compressed sensing,” 25th Annual Meeting of the International Society for Magnetic Resonance Imaging, Honolulu, USA, April 2017.

[0078] [Lazarus 2019] Lazarus, Carole et al., “SPARKLING: for accelerating T2” * "Variable density k-space filling curves of weighted MRI", Medical Magnetic Resonance, Wiley, 2019, 81(6), pp. 3643-3661.

[0079] [Lustig et al., 2007] Lustig M., Donoho D., Pauly JM “Sparse MRI: Application of Compressed Sensing in Fast Magnetic Resonance Imaging”, Magn Reson Med 2007;58: 1182-1195.

[0080] [Sutton et al., 2003] Sutton, Bradley P et al., “Iterative image reconstruction of MRI under field inhomogeneity.” IEEE Transactions on Medical Imaging, 2003; 22(2): 178–88.

Claims

1. A method of performing magnetic resonance imaging of a body using a magnetic resonance imaging scanner, the magnetic resonance imaging scanner comprising a scanner bore (SB), a primary coil (PC), a radio frequency coil (RFC), gradient coils (GC x , GC y , GC z ) and a signal processing unit (SPU), the method comprising the steps of: a. positioning said body (BD) in said scanner bore (SB), wherein a static and substantially uniform magnetic field (B0) oriented along a direction referred to as longitudinal direction (z) is established by said primary coil (PC), said magnetic field being referred to as longitudinal field; b. emitting at least one radio frequency pulse (RFP) to said body using all or part of said radio frequency coil (RFC), said radio frequency pulse being adapted to excite nuclear spins within said body; c. after the or each radio frequency pulse, applying time-varying magnetic field gradients (G x , G y , G z ) defining a trajectory (ST) in k-space to the body using the gradient coils (GC x , GC y , GC z ) and simultaneously acquiring samples of magnetic resonance signals emitted by the excited nuclear spins using all or part of the radio frequency coils, each of the samples corresponding to a point (KS) of the k-space belonging to the trajectory. and d. applying a non-linear reconstruction algorithm to said samples acquired using said signal processing unit (SPU) to reconstruct a magnetic resonance image of said body; wherein said trajectories (ST, MT) in k-space are continuous trajectories complying with a set of constraints including constraints on maximum amplitude and maximum slew rate of said time-varying magnetic field gradients, so that said points (KS) in k-space, referred to as sampling points, corresponding to said samples define a pseudo-random sampling of said k-space, approximating a predetermined target sampling density, said trajectories in k-space minimizing a cost function under said set of constraints, said cost function being defined by a difference between a first term, referred to as an attractive term, promoting a distribution of sampling points in k-space following said predetermined target sampling density, and a second term, referred to as a repulsive term, promoting separation between sampling points in k-space, said repulsive term being expressed as a sum of contributions corresponding to each pair of sampling points; characterized in that said trajectories (ST, MT) in k-space define a Cartesian sampling of a central region (S) of said k-space.

2. The method of claim 1, wherein, said set of constraints including an affine constraint imposing said Cartesian sampling of said central region (S) of said k-space.

3. The method according to any of the preceding claims, wherein, said central region (S) of said k-space has a spherical or elliptical shape.

4. The method according to any of the preceding claims, wherein, said predetermined target sampling density is such that said trajectories are adapted to sample said central region of said k-space at a rate comprised between the Nyquist rate and 1.5 times the Nyquist rate.

5. The method according to any one of the preceding claims, wherein, said attractive term is defined by and said repulsive term is defined by wherein: - K is a two or three dimensional vector representing coordinates of a sampling point in k-space; - is the number of sampling points along the trajectory, each trajectory consisting of N c ≥1 shots, each shot containing N s >1 sampling points; - Q Nc is a set of curves in k-space comprising N c segments, each segment comprising N s points and complying with said constraints; - Ω is a sampling region in k-space; - x represents a generic point of said sampling region Ω; - p is said predetermined target sampling density. - is a norm; - t i , t j is the time of arrival of the two sample points of the same shot identified by the integer indices i and j; and - W is a monotonically increasing function of the weight of the the weight of the patient.

6. The method according to any one of the preceding claims, wherein, said method further comprising a preliminary step of calculating said trajectories (ST) in k-space.

7. A computer-implemented method of calculating trajectories in k-space for samples acquisition in magnetic resonance imaging, said method comprising the steps of: i. determining or receiving from user input the number N of strokes that make up the trajectory c ≥ 1 and the number N of sampling points along each of the strokes s > 1; ii. determining or receiving a predetermined target sampling density from user input; iii. determining or receiving a set of constraints from user input, said set of constraints including constraints on amplitudes of discrete time derivatives of said trajectories; iv. calculating said trajectories in k-space under said set of constraints. iv. said trajectory is computed by minimizing a cost function under said set of constraints, said cost function being defined by the difference between a first term, called an attractive term, which promotes the distribution of the sampling points in the k-space to follow said predetermined target sampling density, and a second term, called a repulsive term, which promotes the separation between sampling points in the k-space, said repulsive term being expressed as a sum of contributions corresponding to each pair of sampling points; characterized in that said trajectory (ST, MT) in the k-space defines a Cartesian sampling of a central region (S) of the k-space.

8. The method of claim 7, wherein, said set of constraints comprises an affine constraint imposing said Cartesian sampling of said central region (S) of the k-space.

9. The method of claim 8, wherein, said central region (S) of the k-space has a spherical or elliptical shape.

10. The method of any one of claims 7 to 9, wherein, said predetermined target sampling density is such that said trajectory is adapted to sample said central region of the k-space at a rate comprised between the Nyquist rate and 1.5 times the Nyquist rate.

11. The method of any one of claims 7 to 10, wherein, said attractive term is defined by the following equation and said repulsive term is defined by the following equation wherein: - K is a two or three dimensional vector representing the coordinates of a sampling point in the k-space; - is the number of sampling points along the trajectory; - Q Nc is a set of curves in k-space comprising N c segments, each segment comprising N s points and complying with said constraints; - Ω is a normalized sampling region in the k-space; - x represents a generic point of said sampling region Ω; - p is said predetermined target sampling density. - is a norm; - t i , t j is the time of arrival of the two sample points of the same shot identified by the integer indices i and j; and - W is a monotonically increasing function of the weight of the the weight of the patient.

12. A computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method according to any one of claims 7 to 11.

13. A set of drive signals (DSG) for a gradient coil of a magnetic resonance imaging scanner, which, when applied to the gradient coil, drive the gradient coil to produce time-varying magnetic field gradients (G x , G y , G z ), which define a trajectory (ST, MT) in k-space, which trajectory passes through a plurality of points (KS) referred to as sampling points, which define a pseudo-random sampling of the k-space, approximating a predetermined target sampling density; wherein the trajectory: - is continuous and complies with a set of constraints comprising constraints on the maximum amplitude and maximum slew rate of the time-varying magnetic field gradient; - a cost function is minimized under said set of constraints, said cost function being defined by the difference between a first term, called an attractive term, which promotes the distribution of the sampling points in the k-space to follow said predetermined target sampling density, and a second term, called a repulsive term, which promotes the separation between sampling points in the k-space, said repulsive term being expressed as a sum of contributions corresponding to each pair of sampling points; characterized in that said trajectory (ST, MT) in the k-space defines a Cartesian sampling of a central region (S) of the k-space.

14. A magnetic resonance imaging scanner comprising: - a scanner bore (SB) within which a body can be positioned; - a primary coil (PC) configured to establish in the scanner bore a static and substantially uniform magnetic field (B0) oriented along a direction called longitudinal direction (z), said magnetic field being called longitudinal field; - a radio frequency coil (RFC) configured to emit at least one radio frequency pulse (RFP) to the body, said radio frequency pulse being adapted to excite nuclear spins within the body, and said radio frequency coil being configured to acquire samples of magnetic resonance signals emitted by the excited nuclear spins; - a gradient coil (GC x , GC y , GC z ) configured to apply a time-varying magnetic field gradient (G x , G y , G z ) to the body, the time-varying magnetic field gradient defining a trajectory (ST, MT) in k-space; - a signal processing unit (SPU) applying a non-linear reconstruction algorithm to the acquired samples of magnetic resonance signals to reconstruct a magnetic resonance image of the body; and - a display configured to display the reconstructed magnetic resonance image. - a controller (CTR) configured to generate drive signals (DSR) for the radio frequency coils and drive signals (DSG) for the gradient coils; wherein the controller is configured to generate the drive signals (DSG) for the gradient coils to produce the time-varying magnetic field gradients (G x , G y , G z ) such that the trajectory (ST) in k-space passes through a plurality of points (KS) in k-space, referred to as sampling points, defining a pseudo-random sampling of the k-space, approximating a predetermined target sampling density, the trajectory being continuous and complying with a set of constraints including constraints on maximum amplitudes and maximum slew rates of the time-varying magnetic field gradients, the trajectory in k-space also minimizing a cost function under the set of constraints, the cost function being defined by a difference between a first term, referred to as an attractive term, promoting a distribution of sampling points in the k-space to follow the predetermined target sampling density, and a second term, referred to as a repulsive term, promoting a separation between sampling points in k-space, the repulsive term being expressed as a sum of contributions corresponding to each pair of sampling points; and wherein the controller is further configured to generate drive signals (DSR) for the radio frequency coils such that they acquire the sample of magnetic resonance signals emitted by the excited nuclear spins, the sample corresponding to the sampling points (KS) of the trajectory in k-space; and characterized in that the trajectory (ST, MT) in k-space defines a Cartesian sampling of a central region (S) of the k-space.

Images