Method for reconstructing MRI images with correction of off-resonance effects

Abstract

A computer-implemented method for reconstructing an MRI image of a body, comprising: a) setting MRI signal samples (s q a) acquiring a k-space image (211) of a magnetic resonance imaging system (200) at different k-space positions along a k-space trajectory, b) acquiring or estimating a B inhomogeneity map ΔB(r), where B is a static magnetic field for polarizing nuclear spins in a body and r represents a position in image space, c) initializing an image (f), d) updating the image by performing a data consistency operation in k-space (211), where the transformation between image space and k-space is performed using a non-uniform discrete pseudo-Fourier transform operator taking into account the B inhomogeneity map obtained from step b) to correct for off-resonance effects, and e) further updating the image by using an image correction neural network (231).

Description

[Technical Field]

[0001] The present invention relates to a method for reconstructing magnetic resonance imaging (MRI) images that aims to enable efficient correction of off-resonance effects, especially in combination with non-Cartesian k-space sampling patterns. [Background technology]

[0002] Although magnetic resonance imaging (MRI) is one of the most powerful imaging techniques used in clinical routine today, the procedure remains time-consuming, especially when acquisition of large and / or high-resolution images containing millions of pixels is desired. For example, a 205 × 205 × 52 mm MRI image was acquired using a T2* sequence with an echo time TE of 6.8 ms in an 11.7 T scanner with 8 transmit coils and 32 receive coils, with an acceptable signal-to-noise ratio (SNR) of 7.6. 3 Acquiring a 3D image of the human brain with a field of view of 100 μm and 200 μm resolution may require an acquisition time of approximately 2 h at a short repetition time (TR ≈ 42 ms), which is clearly unacceptable for clinical purposes.

[0003] This is because MRI signals are transient, and the excitation / acquisition procedure must be repeated for several hours to collect useful measurements in the so-called "k-space" region, the spatial frequency spectrum of the image, which is related to physical space by a multidimensional Fourier transform. Sampling theory dictates that the sampling frequency should be at least twice the highest frequency contained in the signal being sampled (which, in the case of MRI, is the multidimensional Fourier transform of the imaged body); otherwise, aliasing artifacts will occur, known as the Nyquist or Nyquist-Shannon criterion. As a result, using conventional acquisition schemes, the number of k-space samples must be at least equal to the number of pixels in the reconstructed image. Furthermore, SNR requirements impose a minimum acquisition time for each sample.

[0004] Several techniques have been developed to reduce acquisition times while avoiding artifacts ("accelerated MRI"). Among these techniques, compressed sensing (CS) (see, for example, [Lustig 2008]) is noteworthy as it achieves orders of magnitude acceleration compared to parallel imaging, especially during imaging at high matrix sizes (either high resolution and narrow field of view or low resolution and wide field of view).

[0005] Compressive sensing techniques rely on three principles: The reconstructed image must be amenable to a sparse (or compressible) representation. In other words, it must be possible to decompose the image on a predetermined basis (e.g., wavelet-based) such that only a small fraction of the decomposition coefficients are non-zero for strict sparsity or significantly greater than zero for compressibility. Typically, for noisy signals, a coefficient is considered significantly greater than zero if its absolute value is at least equal to the standard deviation of the noise. Alternatively, only a predetermined fraction of the coefficients (those with the largest absolute value) can be kept. For example, only the top 1% of the coefficients are kept, resulting in a compression coefficient of 100. Reconstruction must be performed using nonlinear methods that promote sparsity of the image representation as well as consistency with the acquired k-space samples. To accelerate acquisition, k-space must be undersampled in a non-coherent manner. Non-coherent sampling is typically achieved using a pseudo-random undersampling pattern. While undersampling reduces the number of signal acquisitions and thus provides the required acceleration, the pseudo-randomness ensures that subsampling artifacts are non-coherent, i.e., decorrelated or noise-like, in the sparse representation. This non-coherence is crucial, as it measures the degree of correlation between any pair of elements captured in a sparse (e.g., wavelet) and sensed (e.g., Fourier in CS-MRI) basis.

[0006] In compressed sensing in MRI, non-Cartesian pseudo-random undersampling of k-space is often used because this strategy results in improved incoherence (lower correlation between samples). Preferably, the pseudo-random undersampling is non-uniform, and its density matches the energy distribution of the image acquired in k-space. In clinical applications, this usually means using a variable sampling density that is high near the center of k-space (low spatial frequencies) and falls off radially toward the periphery of k-space (i.e., high spatial frequencies).

[0007] From a purely theoretical perspective, pseudorandom undersampling can be achieved by drawing samples according to a predetermined probability distribution corresponding to the desired sampling density. However, in practice, the short lifetime of MR signals necessitates the acquisition of samples through segmented acquisitions along smooth trajectories, which are defined by the spatial coupling of time-varying magnetic field gradients applied in two and three dimensions (2D / 3D) to the imaged body after excitation of nuclear spins with radio frequency (RF) pulses, depending on the slice versus volume (2D vs. 3D) acquisition.

[0008] Commonly used k-space trajectories are parallel lines (leading to Cartesian sampling), spokes (straight lines radiating from the center of k-space), rosettes, uniform density and variable density spirals. All of these have been applied to compressed sensing in both 2D (slice) and 3D (volume) imaging, for example, by performing only a limited number of signal acquisitions on a Cartesian grid or by randomly sampling the spokes, spirals or rosettes.

[0009] However, better results are obtained by using "non-parametric" trajectories that offer greater incoherence. In particular, the so-called SPARKLING ("Spreading Projection Algorithm for Rapid K-space sampling") technique ([Boyer 2016], [Lazarus 2017], [Lazarus 2019], [Chaithya 2022], WO 2019 / 048565) is based on the projection of a predetermined, usually non-uniform, target sampling distribution onto a set of "allowable" 2D or 3D curves, i.e., onto all curves that represent trajectories that can be obtained without exceeding limits on the magnetic field gradient values ​​and the corresponding slew rates.

[0010] Unfortunately, non-Cartesian MRI sampling schemes are more susceptible to off-resonance artifacts than traditional Cartesian sampling schemes due to ΔB, i.e., inhomogeneities in the static magnetic field B. These inhomogeneities are typically found near the air / tissue interfaces around the sinuses (bucconasal region) or ear canal due to variations in magnetic susceptibility at these interfaces. This is because off-resonance artifacts appear across the readout dimension (vs. non-Cartesian sampling), which is one-dimensional (vs. multidimensional) in Cartesian sampling. Thus, B inhomogeneities are mixed across multiple dimensions in non-Cartesian sampling, making it more difficult to compensate for B inhomogeneities.

[0011] [Chaithya 2023] recently proposed a set of non-Cartesian k-space trajectories called "MORE+GoLF SPARKLING" that minimizes the impact of off-resonance effects during acquisition. Nevertheless, at least in some cases, correction for these effects is still necessary during the image reconstruction step.

[0012] N c is the number of spokes, and N s is the number of samples per spoke, M=N c ×N sConsider the case where N MRI samples are acquired in k-space according to a non-Cartesian multi-spoke sampling scheme Ω: N=N x ×N y ×N z The image to be reconstructed

number

number

number

[0013] Equation (1) can be solved iteratively by proximal gradient descent.

number

[0014] To apply the ΔB0 correction, equation (2) must be modified as follows:

number

number

number

number

[0015] The term Δω0(r n )t m depends on both the k-space and image domains and does not fit into a standard Fourier transform. It has been proposed by [Noll 1991] to decompose this exponential term into a sum of products of variables, each dependent on a single domain.

number

[0016] In this way, terms that depend on the k-space domain can be factored out of the pseudo-Fourier transform. Equation (6) becomes a weighted sum of L direct Fourier transforms.

number

[0017] Coefficient B=(b m,l )∈C M,L and C=(c l,n )∈C L,N can be estimated by solving the matrix factorization problem

number

number

[0018] Another technique, known as Multi-Frequency Interpolation (MFI), divides the off-resonance frequency range into L discrete values ​​Δω 0,l Divide into

number

[0019] Yet another alternative approach, known as Multi-Temporal Interpolation (MTI), divides the time window into L discrete values ​​t l Divide into

number

[0020] Whatever the method used to calculate B and C, equation (8) is a non-uniform discrete pseudo-Fourier transform operator F Ω,Σ (Σ represents interpolation) is defined as F Ω,Σ and its auxiliary operators

number

number

[0021] Use spoke redundancy (i.e., N c The computational burden can be significantly reduced by solving a weighted version of equation (9) using the histogram of the ΔB field map (using the same decomposition for each spoke), typically reducing the image dimensions N = 384 × 384 × 208 voxels to N = 1000 bins (see [Fessler 2005]). Thus, the matrix E is reduced from M × N to N s ×N b and thus correction factors can be obtained in a few seconds for high-resolution 3D volumes.

[0022] However, compared to the case where the ΔB0 correction is not performed, the use of the pseudo-Fourier operator (8) and its adjoint operator increases the computational cost by a factor L, which is typically on the order of 10-20 at 3 Tesla and up to 40 at 7 Tesla to achieve satisfactory results. In fact, the computational cost of image reconstruction in the framework described above is in most cases N logN ×I×Q×L, where "I" is the number of iterations of the proximal gradient algorithm in equations (3) and (4), F The ΔB correction is determined by the number of calls, which can result in hours (4-8 hours) of computation time to reconstruct a high-quality image, compared to minutes (10-15) without ΔB correction. [Prior art documents] [Patent documents]

[0023] [Patent Document 1] International Publication No. 2019 / 048565 Brochure Summary of the Invention [Problem to be solved by the invention]

[0024] The present invention aims to completely or partially overcome the above-mentioned drawbacks of the prior art, and more particularly to reduce the computation time required to perform MRI image reconstruction with ΔB correction, especially in the case of non-Cartesian k-space sampling and / or parallel acquisition. [Means for solving the problem]

[0025] According to the present invention, this objective is achieved by combining model-based image reconstruction (i.e., reconstruction based on Equation (8)) with a limited number of interpolation components L, e.g., L = 5, using an image correction neural network to complete the off-resonance correction. The use of an image correction neural network allows the value of L (and optionally Q and I) to be reduced with little increase in computational complexity and without compromising the effectiveness of the ΔB correction and the quality of the reconstructed image.

[0026] Alternating between model-based image reconstruction in k-space and neural network-based image correction in physical space using so-called "unfolded" architectures dedicated to non-Cartesian acquisitions is known, for example, from [Ramzi 2022], but was based solely on conventional non-uniform Fourier transform reconstruction (i.e., without ΔB0 correction). To the best of our knowledge, this is the first time that such an architecture has been used together with pseudo-Fourier transform reconstruction for ΔB0 correction, and the unexpected result is that such an architecture allows for a significant reduction in the number of interpolation components L.

[0027] The object of the present invention is to provide a computer-implemented method for reconstructing MRI images of a body, comprising the steps of: a) acquiring a set of MRI signal samples at different k-space locations along a k-space trajectory; b) acquiring or internally estimating a separate B0 inhomogeneity map, where B0 is the static magnetic field for polarizing the nuclear spins of the body and r represents a position in image space; c) initializing the image; d) updating the image by performing a data consistency operation in k-space, where the transformation between image space and k-space is performed using a non-uniform discrete pseudo-Fourier transform operator taking into account the B0 non-uniformity map obtained from step b) to correct for off-resonance effects; e) further updating the image by using an image correction neural network; A computer-implemented method comprising:

[0028] Another object of the invention is a computer program product comprising instructions that cause a computer to carry out the steps of such a method when the program is executed by a computer.

[0029] Yet another object of the present invention is to provide a magnet or coil that generates a static magnetic field for polarizing the nuclear spins of the body; gradient coils for generating magnetic field gradients that define a k-space trajectory; a transmit coil coupled to a radio frequency pulse generator for applying radio frequency pulses at the Larmor frequency to the polarized nuclear spins; at least one receive coil coupled to a radio frequency receiver for acquiring MRI signals emitted by the nuclear spins; a processor programmed to reconstruct an MRI image from the acquired MRI signals; The processor is a magnetic resonance imaging device, characterized in that it is programmed to carry out the steps of such a method.

[0030] Additional features and advantages of the present invention will become apparent from the subsequent description taken in conjunction with the accompanying drawings. [Brief explanation of the drawings]

[0031] [Figure 1] 1 is a high-level flowchart of a method in accordance with the present invention. [Figure 2] FIG. 2 is a diagram of an image reconstruction pipeline according to an embodiment of the present invention. [Figure 3A] 10 is an experimental result showing the technical effect of the present invention. [Figure 3B] 10 is an experimental result showing the technical effect of the present invention. [Figure 4] FIG. 1 is a block diagram of an MRI scanner according to an embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0032] The first step, designated as step (a), of the method of FIG. 1 consists in acquiring set MRI signal samples at different k-space positions along a k-space trajectory. The trajectory can be any known, physically realizable 2D or 3D trajectory, preferably performing a non-Cartesian sampling of k-space. According to a preferred embodiment of the present invention, the trajectory can be of the SPARKLING or MORE+GoLF SPARKLING type. As is well known in the art, the k-space trajectory is defined by applying appropriate magnetic field gradients to the body of the imaging subject, polarizing the nuclear spins by applying a (nearly) uniform static magnetic field B0. The polarized nuclear spins are excited by a radio frequency pulse at the Larmor frequency γB0 and re-emit a radio frequency MRI signal. Many different excitation RF pulse sequences are known in the art and can be used to implement the present invention. The excitation pulses and / or MRI signals can be emitted / received using a single RF coil or several coils. If several receive coils are used (parallel MRI), their sensitivity maps S q must be known to perform image reconstruction. Sensitivity maps can be obtained through external calibration or by the MRI signal itself (self-calibration), as considered by [El Gueddari 2018].

[0033] The second step (b) consists in obtaining or estimating a ΔB(r) map from a particular MRI acquisition sequence (in which case step (b) may precede step (a)) or using the MRI signal samples acquired in step (a). See, for example, [Daval-Frerot 2022].

[0034] Step (c) consists in initializing the pixelated or voxelized image to be reconstructed. For example, initialization to zero can be used. The initialized image is denoted f0.

[0035] Steps (d) and (e), which are repeated I≧1 times, preferably I>1 times, form the core of the method of the present invention and are shown in detail in FIG.

[0036] In step (d), a data consistency calculation, preferably also a sample density compensation, is performed on the MRI signals s emitted by all Q receive coils. q These k-space operations, represented by block 211 in FIG.

number

[0037] The weighted combination of the individual receive coil contributions and transformation to physical space is performed in block 221 by applying the following operator to the output of block 211:

number

number

[0038] Next, block 231 applies an image correction neural network, such as UNets or another developed neural network architecture (KIKI-Net, Cascade-Net, PDNet...), to the physical space data. An overview and criteria for suitable neural network architectures is provided by [Ramzi 2020].

[0039] Block 241 applies the operator ∂ ... F Ω,Σ S q (14) The transformation to k-space is performed by applying

[0040] The second iteration then begins at block 212, similar to 211, and so on for I iterations. The pipeline ends at blocks 22I and 23I, the output of which is the final reconstructed image

number

[0041] Overall, this process is similar to the proximal gradient descent method of equations (3) and (4), with the proximal operator replaced by the application of a neural network for image denoising.

[0042] According to a particular embodiment, if the interpolation coefficients of the pseudo-Fourier operator are of the SVI type, the matrices B and C are calculated for 2L elements, and alternatively, the L pairs of coefficients b corresponding to the L largest singular values ​​are calculated. m,l , c l,n It may be advantageous to use a first set of L pairs of coefficients corresponding to the L+1 to L+2 largest singular values, and a second set of L pairs of coefficients corresponding to the L+1 to L+2 largest singular values, since the SVI coefficients correspond to different regions of the off-resonance spectrum. 1 / 2 It is called.

[0043] In the training phase, all neural network blocks 23i (i = 1 to I) are trained one by one. For this purpose, block 25i (251...25I) is trained by combining the concatenated output of the first i-th neural network output and the target image

number

[0044] The technical effect of the present invention is illustrated by Figures 3A and 3B, which show exemplary MRI images of a patient's head obtained by different methods. The images were acquired with the following acquisition parameters: an in-plane 24 cm field of view spanning 12.5 cm, an observation (or readout) time T = 20.48 ms centered at an echo time TE = 20 ms. obs ,Repetition time TR=37ms,Various number of spokes N c , and N s Acquired using a full 3D SPARKLING trajectory with 2048 samples per spoke. This image represents a larger dataset (99 training, 11 validation, and 11 testing patients, both male and female, ages 18-90 years, BMI 16.5-40).

[0045] The lines in Figures 3A and 3B correspond to (from top to bottom) axial, sagittal, and coronal slices (or views) of a 3D image of the same patient's brain, with different columns corresponding to different reconstruction methods.

[0046] The first column of FIG. 3A corresponds to the prior art case, where the image reconstruction is performed using N F = 780, corresponding to N i = 20 and N c By adopting = 20, the image is performed by the proximal gradient descent method (Equations (3) and (4)) without ΔB correction. Very strong artifacts are observed, especially in the buccal-nasal region (top of the first line, left of the second line, bottom of the third line).

[0047] The fifth column of Figure 3A corresponds to another case according to the prior art, where the ΔB correction is performed using a modified proximal gradient descent algorithm using the pseudo-Fourier operation of Equation (8) with an interpolation component of L = 20, also known as the "wavelet" approach. Comparison with the images in the first column shows that the artifacts have essentially disappeared, however, N i =20, N c By choosing N = 20 and L = 20, F= 15600, which corresponds to a very long reconstruction time, e.g., on the order of 8 hours. This image is considered a benchmark for evaluating the performance of alternative reconstruction methods.

[0048] The second, third and fourth columns (from left to right) of Figure 3A correspond to different embodiments of the method of the present invention, applying both a pseudo-Fourier operator and a neural network to correct for the ΔB effect. F The second column N corresponds to =45 i =5, N c For N = 5 and L = 1, it can be seen that the use of a single interpolation term already significantly reduces artifacts compared to the "no correction" case in the first column. i =5, N c =5, L=3→N F =135 (third column) and N i =5, N c =5, L=5→N F At = 225 (fourth column), a further reduction in artifacts is achieved, corresponding to a reduction in computation time of over 60 times compared to the baseline time (i.e., hours to minutes). In the latter case, residual blur is only visible in the immediate vicinity of the sinuses.

[0049] The first, fourth, and fifth (last) columns of Figure 3B are identical to the corresponding columns of Figure 3A. Figure 2A shows the same parameters (and therefore approximately the same computation time) as the fourth column, namely N i =5, N c =5, L=5→N F = 225. The artifacts are much more apparent, confirming that the method of the present invention is much more effective at compensating for ΔB than prior art wavelet approaches.

[0050] The third column in Figure 3B corresponds to image reconstruction using Fourier operators and neural networks but without ΔB correction. The image quality is comparable to that of the second column and much lower than that of the fourth column, confirming the importance of ΔB correction.

[0051] A quantitative assessment of image quality for the second, third, and fourth columns in both Figures 3A and 3B can be performed by calculating the root mean square error (RMSE), the peak signal-to-noise ratio (PSNR, expressed in dB), and the structural similarity index (SSIM, ranging from 0 to 1) compared to the fifth column (i.e., the "optimal" reconstruction). These indices are defined in [Ramzi 2020], and their values ​​for the different cases considered here are reproduced in the figures.

[0052] FIG. 4 shows a highly schematic diagram of an apparatus for carrying out the method according to the present invention. This apparatus is essentially an MRI scanner including a magnet or primary coil PC for generating a static magnetic field B0, a single or preferably multiple radio-frequency coils TC for generating RF pulses and receiving NMR signals (different radio-frequency coils can be used for transmitting and receiving, or the same coil can serve both purposes), and three gradient coils GCx, GCy, and GCz (not shown) for generating magnetic field gradients Gx, Gy, and Gz along the x, y, and z axes, respectively. The gradient coils are arranged around a scanner bore SB into which the body BD to be imaged can be introduced. The term "body" should be interpreted broadly and refers to any living or inanimate object containing atomic nuclei, at least some of which are assumed to have non-zero spin. For example, the body BD can be a human or animal body or part thereof (e.g., a head), or a product such as a "phantom." The apparatus also includes a control unit (CTR) for driving the gradient coil and the one or more radio frequency coils according to predetermined drive signals (a reference DSR for the radio frequency coil and a DSG for the gradient coil) to perform readout sequences and (preferably non-Cartesian) k-space trajectories, respectively. The apparatus also includes a signal processing unit (SPU) for acquiring NMR signals from the one or more radio frequency coils and performing image reconstruction. The control unit and the signal processing unit may be separate devices, or a single device may perform both functions. The one or more devices may include one or more appropriately programmed general-purpose computers or digital signal processors, appropriately configured specialized digital circuits, or both. The signal processing unit also includes receiver circuitry for acquiring the NMR signals, which includes a signal amplifier, a sampler, and an analog-to-digital (ADC) converter.

[0053] The scanner of FIG. 4 differs from the prior art essentially in that the signal processing unit SPU is configured or programmed to perform image reconstruction according to the method of the present invention.

[0054] Non-patent literature [Boyer 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. [Chaithya 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. [Chaithya 2023]Chaithya G R et al“Improving spreading projection algorithm for rapid k-space sampling trajectories through minimized off-resonance effects and gridding of low frequencies.”Magnetic Resonance in Medicine(2023). [Daval-Frerot 2022]Daval-Frerot G.et al.“Iterative static field map estimation for off-resonance correction in non-Cartesian susceptibility weighted imaging”Magnetic Resonance in Medicine(2022). [El Gueddari 2018]El Gueddari L,et al.“Self-calibrating nonlinear reconstruction algorithms for variabledensity sampling and parallel reception MRI”2018 IEEE 10th Sensor Array and Multichannel Signal Processing Workshop(SAM).IEEE;2018:415-419. [Fessler 2005]Fessler J.A.et al.“Toeplitz based iterative image reconstruction for MRI with correction for magnetic field inhomogeneity”.IEEE Trans Signal Process.2005;53:3393-3402. [Lazarus 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. [Lazarus 2019]Lazarus,Carole et al.“SPARKLING:variable-density k-space filling curves for accelerated T 2 * - weighted MRI”Magnetic Resonance in Medicine,Wiley,2019,81(6),pp.3643-3661. [Lustig 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 [Ramzi 2020]Ramzi,Z.et al.Benchmarking MRI Reconstruction Neural Networks on Large Public Datasets”,Appl.Sci.2020,10,1816 [Ramzi 2022]Ramzi Z.et al.“NC-PDNet:a Density-Compensated Unrolled Networkfor 2D and 3D non-Cartesian MRI Reconstruction”IEEE Transactions on Medical Imaging,2022,10.1109 / TMI.2022.3144619.

Claims

1. 1. A computer-implemented method for reconstructing an MRI image of a body (BD), comprising: a) Set MRI signal samples (s ｑ ) at different k-space locations along a k-space trajectory; b) B ０ Inhomogeneity map ΔB ０ (r), ０ is a static magnetic field for polarizing nuclear spins of the body, and r represents a position in image space; c) initializing an image (f); d) updating the image by performing a data consistency operation in k-space (211), wherein the transformation between image space and k-space is performed by using the B 2 data obtained from step b) to correct for off-resonance effects. ０ - an updating step, which is performed using a non-uniform discrete pseudo-Fourier transform operator taking into account the non-uniformity map; e) further updating the image by using an image correction neural network (231); 11. A computer-implemented method comprising:

2. The method of claim 1 , wherein steps d) and e) are repeated multiple times.

3. 3. The method of claim 1, wherein step d) is performed by decomposing the non-uniform discrete pseudo-Fourier transform operator into a sum of non-uniform discrete Fourier transform operators weighted by interpolation coefficients.

4. the non-uniform pseudo-Fourier transform operator F decomposed into a sum of non-uniform discrete Fourier transform operators weighted by interpolation coefficients; Ω teeth, s=F Ω，Σ f where s is the k-space position k ｍ = k (t ｍ ) at each time t ｍ M>1 signal samples s acquired at ｍ where k is the k-space trajectory, m ranges from 1 to M, and f is the position r ｎ The real space image x corresponding to ｎ is an N-component vector representing N>1 pixels or voxels of [Equation 1] The index l ranges from 1 to L, where L>1 is the number of the non-uniform discrete Fourier transform operators, and the coefficient b ｍ，ｌ and c ｌ，ｎ teeth, [Equation 2] But Δω ０ (r ｎ ) = γΔB ０ (r ｎ ) [Equation 3] 4. The method of claim 3, wherein γ is a gyromagnetic ratio of the nuclear spins.

5. The interpolation coefficients are [Equation 4] The method of claim 4, wherein the singular value decomposition of

6. The coefficients b ｍ，ｌ and c ｌ，ｎ The method of claim 5 when dependent on claim 2, wherein different sets of are used in successive iterations of steps d) and e).

7. 7. The method according to claim 1, wherein step a) is performed by parallel acquisition with Q>1 receive coils (RFC), and wherein signal samples emitted by said receive coils are combined with respective sensitivity maps during step d).

8. The data consistency calculation of step d) is [Equation 5] and where Q≧1 is the number of receiver coils acquiring MRI signal samples, and s ｑ is the vector of signal samples acquired by the qth receiver coil, f is the pixel or voxel vector representing the image, α is a positive real value, and F Ω is the non-uniform discrete pseudo-Fourier transform operator, and S ｑ 8. The method according to claim 1, wherein ∑ is the sensitivity map of the q receiving coil, D is a density correction operator, and the superscript H denotes the conjugate transpose of a complex matrix.

9. 9. The method according to any one of claims 1 to 8, wherein the image correction neural network used in step e) is a fully convolutional neural network.

10. The method of claim 9 , wherein the image correction neural network is Unet.

11. 11. The method of claim 9 or 10 when dependent on claim 2, wherein each iteration of step e) uses a different image correction neural network trained on a data set including the outputs of all previous iterations of step d).

12. The method of any one of claims 1 to 11, wherein the k-space trajectory defines a non-Cartesian sampling of k-space.

13. 13. The method of claim 12, wherein the k-space trajectory is a continuous trajectory such that the k-space locations at which samples are acquired define a pseudo-random sampling of the k-space consistent with a predetermined target sampling density.

14. A computer program product comprising instructions that, when said program is executed by a computer (SPU), cause said computer to carry out said steps of the method according to any one of claims 1 to 13.

15. A magnetic resonance imaging apparatus, A static magnetic field (B) for polarizing the nuclear spins of the body (BD) ０ a magnet or coil (PC) for generating a Gradient coils (GCx, GCy, GCz) for generating magnetic field gradients that define a k-space trajectory; a transmit coil (RFC) coupled to a radio frequency pulse generator for applying radio frequency pulses to the polarized nuclear spins at their Larmor frequency; at least one receive coil (RFC) coupled to a radio frequency receiver for acquiring MRI signals emitted by said nuclear spins; a processor (SPU) programmed to reconstruct an MRI image of the body from the acquired MRI signals; Including, A magnetic resonance imaging apparatus, characterized in that the processor is programmed to carry out the steps of the method according to any one of claims 1 to 13.

Images