Technique for processing magnetic resonance imaging data

Abstract

The invention relates to a method and to a device for processing magnetic resonance imaging, MRI, data. In particular, the present invention relates to a technique for determining spin phase fluctuations between separate shots of acquisitions of a multi-shot MRI data collection scan. The method comprises providing (S100) MRI data acquired by a multi-shot MRI data collection scan of an MRI scanner, wherein the MRI data comprises overlapped k-space regions commonly sampled in separate shots of acquisitions. The method further comprises determining (S200) spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan from calibration regions, wherein the overlapped k- space regions in the MRI data between separate shots are used as the calibration regions based on which spatially-varying spin phase fluctuations in separate shots are extracted. Advantageously, a robust, navigator-free, computational efficient multi-shot method without SNR penalty may be provided.

Description

[0001] Technique for processing magnetic resonance imaging data

[0002] Field of the invention

[0003] The invention relates to a method and to a device for processing magnetic resonance imaging, MRI, data. In particular, the present invention relates to a technique for determining spin phase fluctuations between separate shots of acquisitions of a multi-shot MRI data collection scan.

[0004] Technical background

[0005] In the present specification, reference is made to the prior art references [1] -

[0049] as listed at the end of this specification. The references illustrate technical background of the invention and related prior art techniques.

[0006] As essential tools for neuroimaging, functional [1-3] and diffusion weighted [4] MRI scans (i.e., fMRI, DW-MRI) are commonly acquired using the single-shot echo-planar imaging (EPI) [5] sequence, enabling rapid imaging by capturing a segment of k-space after one spin excitation. However, the spatial resolution achievable with this method is constrained by the finite width of the k- space acquisition window, due to the short T2 or T2* relaxation times of transverse magnetization. Additionally, the long EPI echo train can lead to geometric distortion in the presence of magnetic field inhomogeneity.

[0007] To surpass these limits by single-shot coverage, multiple shots of scans that samples distinct k-space bands can be combined into one image with higher resolution or less distortion. Nevertheless, the presence of uncontrollable shot-to-shot phase variations [6] poses a significant challenge, leading to ghosts artifacts in image from combined k-space bands with inconsistent phase. These phase fluctuations can arise from factors such as motion (e.g., respiration, heartbeat, mechanical vibration), system imperfections (e.g., varying eddy currents between shots), and are especially pronounced in higher magnetic field strengths [7] and during the presence of diffusion sensitizing gradients, impeding the widely adoption of this high-resolution strategy.

[0008] For decades, numerous methods for multi-shot acquisitions have been explored to mitigate the notorious shot-dependent phase errors, aiming for robust sampling of distinct k-space regions with phase consistency. One established prior art approach leverages the navigator echoes additional to the imaging echoes to explicitly estimate the specific phase evolution or map during each shot. These navigators - ranging from one- [8-11], two-[12-18] to three-dimensional - can be integrated into the sequence of line-by-line scans or multi-shot EPI using mosaic segmentation and interleaved sampling schemes. However, additional navigator acquisitions prolong the scan time, may contain inconsistent information compared to the imaging echoes, and can suffer from signal decay due to short transverse relaxation time in higher magnetic field. While self-navigating multi-shot MRI is possible, it is limited to certain non-conventional sampling patterns, such as spiral (e.g., SNAIL)

[0019] , concentric rectangular stripes (e.g., PROPELLER

[0020] ), or "Butterfly" navigators

[0021] that redundantly samples k-space center along the imaging trajectory during each acquisition.

[0009] Alternatively, post-processing techniques have been investigated to remove the ghosts' artifacts originating from the shot-dependent phase errors in moderately undersampled EPI interleaves. In each interleave, with sufficient data allowing for a reasonable parallel imaging reconstruction, the phase fluctuation maps can be estimated by iterative reconstruction

[0022] , or explicitly obtained by taking total variations of the SENSE reconstruction for each interleave

[0023] , Moreover, exploiting low-rank penalty terms formulated in either k-space

[0024] ,

[0025] or local image-space

[0026] during iterative reconstruction can mitigate ghosts due to inter-shot phase errors. Unfortunately, given highly undersampled data, these approaches could end up with residue image artifacts or substantially increased computational complexity.

[0010] Finally, the smoothly varying inter-shot phase variations can be eliminated in several magnitudebased low-resolution scans with subpixel modulation in shifts, which are used to reconstruct a high- resolution image [27-31], These methods import the essential elements of the well-established super-resolution optical microscopy [32-42], However, suppressing partial voxel signals in each low-resolution scan inevitably cause SNR loss because of energy conservation, as another crucial factor in functional and particularly diffusion MRI due to its inherently low SNR nature.

[0011] Objective of the invention

[0012] In view of the above-mentioned shortcomings of the prior art, it is an objective of the invention to provide an improved technique for dealing with shot-to-shot phase variations in multi-shot MRI scans. In particular, it is an object of the invention to provide a more efficient, preferably navigator-free technique, for dealing with shot-to-shot phase variations in multi-shot MRI scans.

[0013] Summary of the invention These objectives are solved by a method and / or a device comprising the features of the independent claims. Advantageous embodiments and applications of the invention are defined in the dependent claims.

[0014] According to a first general aspect of the invention, the above objective is solved by a method of processing magnetic resonance imaging, MRI, data. The method comprises the step of providing MRI data acquired by a multi-shot MRI data collection scan of an MRI scanner, wherein the provided MRI data comprises overlapped k-space regions commonly sampled in separate shots of acquisitions. In other words, the k-space data regions in separate shots of acquisitions partially overlap. The method further comprises determining spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan from predetermined k-space regions (hereinafter referred to as calibration regions). The overlapped k-space regions in the MRI data between separate shots are used as the calibration regions based on which spatially-varying spin phase fluctuations in separate shots are extracted. Advantageously, a robust, navigator-free, computational efficient multi-shot method without SNR penalty may be provided.

[0015] In other words, the MRI data are acquired using a multi-shot technique, wherein a segment of k- space data is acquired after each RF excitation. Namely, data from partial k-space regions are acquired in separate shots and concatenated to produce a single image. A k-space region corresponds to a "location" in k-space grid. The provided MRI data corresponds to MRI data acquired in Radio Frequency, RF, receiver coils of the MRI scanner. According to prior art approaches, usually only a few k-space pixels along readout dimension are overlapped to avoid sampling "holes", but according to the present invention, overlapped k-space regions commonly sampled in separate shots of acquisitions are exploited to extract spatially varying spins phase fluctuations or spin phase fluctuation maps in separate shots. The MRI data in respectively overlapping k-space regions are not identical, so that the MRI data in overlapping k-space regions can be used to extract the spin phase fluctuation between them. In this specification, the term spin phase fluctuations is also referred to simply as phase fluctuations and the term spin phase fluctuations map is also referred to simply as phase fluctuations map.

[0016] According to preferred embodiments, the spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan are determined from calibration regions preferably based on an approach using k-space interpolation or based on an approach using eigen- value decompositions. Advantageous embodiments of such approaches using k-space interpolation or eigenvalue decompositions will be described further below.

[0017] According to a preferred embodiment, the MRI data based on which the spin phase fluctuations between shots are determined do not comprise any 2D or 3D navigator data that are additionally acquired in each shot of acquisition to monitor the spin phase fluctuations in each shot. Advantageously, time-consuming navigators can be avoided and a faster and more robust MRI technique is provided. For example, about 30ms - 50 ms per shot of acquisition can be saved compared to a 2D navigator-based multi-shot MRI technique.

[0018] According to a further aspect, the spin phase fluctuations may be determined by estimating, from the MRI data in the overlapped k-space region, a kernel convolved with k-space signals, wherein the kernel corresponds to the Fourier transform of a relative image-space spin phase fluctuation map during each shot, and being shift-invariant across all k-space locations and radio frequency (RF) receivers channels. Depending on the definition of Fourier transform and convolution in a specific context, the kernel may be a flipped and complex conjugate version of the Fourier transform of the image space spin phase fluctuation map.

[0019] In other words, an image-space spin phase error map during each shot may be mathematically described as a convolution kernel in k-space. The kernel may be used to reconstruct an object's image from the provided MRI data while correcting for shot-dependent phase errors. A key insight according to this aspect is that these convolution kernels are shift-invariant across the entire k-space, and thus not limited to the k-space central region (k-space center) as usually used for parallel imaging reference scans due to higher signal-to-noise-ratio, SNR. Therefore, the kernel for the relative image-space spin phase fluctuation maps may be estimated from overlapped k-space regions commonly sampled in separate shots of acquisitions in order to obtain the phase fluctuations, e.g. for correcting for the finally reconstructed high-resolution images. Such phase fluctuation correction methods based on shift-invariant kernel extraction in subspace can successfully improve the resolution and reduce geometric distortion compared to single-shot EPI, without 2D navigators. In this context, the term image-space spin phase fluctuation map should be understood as a two or three- dimensional distribution of the spatially varying spin phase in each shot of the MRI data collection scan.

[0020] According to a further aspect, the overlapped k-space regions may comprise k-space regions centred around a k-space location different than the k-space centre. Preferably, the overlapped k- space regions may comprise at least one k-space region centred around a k-space location different than the k-space centre. Advantageously, the provided MRI data or the MRI data collection scans providing the MRI data with the overlapped k-space regions are not confined to the k-space centre.

[0021] According to a further preferred embodiment, the step of extracting the spatially-varying spin phase fluctuations in separate shots comprises establishing a linear interpolation relationship. The linear interpolation relationship interpolates a k-space data point at a location ks in a k-space data region from at least one patch of k-space data centred around the location ks in at least one other overlapped k-space data region acquired by different shots of the MRI data collection scan. Thus, the k-space data point and the patch of k-space data are from different shots of k-space data but within an overlapped k-space region at k-space locations commonly sampled in separate shots of acquisitions. This interpolation may also be more broadly described as a mapping of patch of k- space data in one shot of acquisition to a k-space data point in another shot of acquisition using acquisition data from overlapped (overlapping) k-space regions commonly sampled in separate shots of acquisitions. In this context, the term patch (k-space patch) should be understood as a k- space subset, also referred to as a k-space neighbourhood or k-space stamp in the prior art. The linear interpolation relationship may be established as a linear equations system.

[0022] Establishing the linear interpolation relationship has the advantage that information regarding shot-dependent phase errors contained in the MRI data is extracted from a relatively small k-space regions at any k-space locations, i.e. the overlapped k-space regions in k-space segments acquired in different shots of acquisitions, using a mathematical technique that can efficiently be solved to extract the spin phase fluctuation between different shots of acquisitions.

[0023] The linear interpolation relationship may be solved by direct inversion to obtain or estimate localized cardinal functions in k-space. The direct inversion may comprise regularizations to obtain the localized cardinal functions in k-space, e.g. for correcting spatially-varying spin phase fluctuations between separate shots of acquisitions.

[0024] The linear interpolation relationship may also be solved by eigenvalue decompositions with thresholding in eigenvector space to obtain image-space spin phase fluctuations maps, e.g., for correcting spatially-varying spin phase fluctuations between separate shots of acquisitions.

[0025] The above two approaches to solve the linear interpolation relationship have the advantage that the mathematical frameworks previously used to analyse MRI parallel imaging, e.g., reproducing kernel Hilbert space (RKHS) and subspace methods (ESPIRiT), serve as theoretical backgrounds with mathematical rigor to explain the invented computational efficient methods for extracting spatially- varying spin phase fluctuations of multi-shot MRI scans.

[0026] According to a further aspect, the step of establishing a linear interpolation relationship comprises determining a calibration matrix by selecting subsets of MRI data from the calibration region. By way of example, the calibration matrix may be determined by taking a patch of MRI data using a mathematical sliding window (e.g. the sliding window may define the range a portion of data are selected from, at a moving location at a time) as a row of the calibration matrix, iteratively sliding through all k-space locations and radio frequency, RF, receivers' channels with the same patch size in the calibration region. The calibration matrix has the advantage that spatially-varying spin phase fluctuations in separate shots can be extracted efficiently, e.g. using an approach to directly estimate the k-space interpolation relationships or using an eigenvalue approach, as further described below.

[0027] According to a yet another aspect, each localized cardinal function that is estimated (obtained) from two or multiple overlapped k-space data regions may correspond to a kernel of coefficients corresponding to approximations of the Fourier transform of the relative image-space spin phase between separate shots of acquisitions, interpolating a patch of k-space data in one of the two or multiple overlapped k-space data regions, to a k-space data point in another one of the two or multiple overlapping k-space data regions. This aspect uses that fact that a k-space neighborhood from one acquisition shot (from data in one shot, within a overlapped k-space region) can be interpolated to a data point in another acquisition shot (from data in another shot, within the overlapped k-space region) with the cardinal function representing k-space modulation effects by image phase fluctuations. The coefficients of the kernel may correspond to the flipped and complex conjugate of the approximation of the Fourier transform of image-space phase fluctuations, depending on the definitions of Fourier transform and convolution in a specific context.

[0028] The localized cardinal functions as the unknowns to be solved may be estimated by inverting the linear equations system describing the linear interpolation relationship between the k-space data patches and the points in separate shots within the calibration region. The spin phase fluctuation maps in the corresponding separate shots of acquisitions may be obtained by, optionally flipping and taking complex conjugate operation, zero-filling and inverse Fourier-transforming the estimated cardinal functions into image-space. Additionally or alternatively, the spin phase fluctuation maps in the corresponding separate shots of acquisitions may be used to correct the respective images as the inverse Fourier transform of the k-space data acquired in separate shots, iteratively performed for all shots.

[0029] According to a yet another aspect, when obtaining localized cardinal functions, the method may further comprise the step of correcting k-space data in all shots of acquisitions by, optionally first flipping and taking complex conjugate operation of the estimated cardinal functions, convolving each of these cardinal functions with all the k-space data acquired in the respective shot including but not limited to the overlapped k-space region, to interpolate k-space data without effects of image-space phase fluctuations relative to the shot where a data point is interpolated to, iteratively performed for all shots.

[0030] As mentioned above, a further preferred embodiment uses an eigenvalue approach based on the determined calibration matrix to estimate the spatially-varying spin phase fluctuations in separate shots. According to this embodiment, extracting the spatially-varying spin phase fluctuations in separate shots of acquisitions may comprise performing eigenvalue decompositions on the calibration matrix to determine image-space eigenvector maps, wherein the determined image-space eigenvector maps comprises a primary eigenvector map associated with the largest eigenvalue, the phases of the primary eigenvector map correspond to the spin phase fluctuation maps. Thus, the eigenvector maps associated with the largest eigenvalue are referred to here as primary eigenvector map.

[0031] The primary eigenvector maps may be taken complex conjugate and multiplied with the image corresponding to a partial k-space sampling in each shot for direct corrections of phase fluctuations. Or a so-called forward model may be used considering all acquired k-space data without phase fluctuations corrections and estimated eigenvector maps (in a linear system) to compute an image. This embodiment has the advantage that the spatially-varying spin phase fluctuations in separate shots can be extracted in a manner through signal space that is more robust to noise.

[0032] According to a further aspect of this preferred embodiment, the step of performing eigenvalue decompositions may comprise the following steps:

[0033] - decomposing the calibration matrix using a first singular value decomposition, SVD;

[0034] - discarding resulting singular vectors in noise space corresponding to small singular values below a thresholding value; - transforming the remaining singular vectors in the signal space into image-space, e.g. by re-shap- ing the signal space singular vectors into k-space filters, optionally flipping and taking complex conjugate, zero-filling, and taking a inverse Fourier transform; and

[0035] - performing a second SVD on a matrix formed from data at each image pixel location across shot and eigen-index dimensions, to obtain singular vectors at each pixel location, containing signal fluctuation information, the singular vectors at all pixel locations correspond to image-space eigenvector maps determined from the eigenvalue decompositions of the calibration matrix.

[0036] According to a further embodiment, an amplitude of the image-space eigenvector maps may be used to correct possible amplitude fluctuations between separate shots of acquisitions. Advantageously, the proposed method may also be used for correcting amplitude fluctuations. In addition or alternatively, a secondary eigenvector obtained from image-space eigenvector maps may be used as the secondary amplitude and / or phase fluctuation maps for correcting k-space data or images in separate shots jointly with the primary ones.

[0037] According to a further aspect, the MRI data may be acquired by an MRI data collection scan using an MRI pulse sequence of multi-shot readout-mosaic-segmented Echo Planar Imaging, EPI. According to this aspect, the overlapped k-space data regions are intersections between overlapped EPI readout segment bands, preferably comprising e.g., 8 to 30 overlapped pixels along readout dimension between overlapped k-space bands in a typical 0.5mm in-plane high-resolution readout segmented EPI scan in a commercial MRI scanner.

[0038] According to a further variant of this aspect, a readout gradient polarity may be flipped between shots of acquisitions for overlapped EPI readout segments. This can avoid spatially-varying T2 or T2* signal decay variations along different paths of sampling trajectory in the overlapped k-space pixels, where only the convolution effects by image-space phase fluctuations caused by shot-de- pendent motions or system imperfections remains for the kernel extractions.

[0039] According to a further aspect, the provided MRI data may be are acquired by a functional MRI data collection scan and / or a diffusion weighted MRI data collection scan. According to yet another aspect, the provided MRI data may be acquired by one of the following MRI pulse sequences: a multishot Cartesian EPI data collection scan or a multi-shot non-Cartesian data collection scan. Thus, the proposed method is not limited to EPI MRI imaging but also applicable to multi-shot spiral trajectory (so called spiral EPI) instead of a Cartesian trajectory. According to a second general aspect of the invention, a method of magnetic resonance imaging, MRI, an object is provided. The method comprises the steps of providing an MRI scanner and acquiring MRI data acquired by a multi-shot MRI data collection scan of the MRI scanner, wherein the MRI data comprises overlapped k-space regions commonly sampled in separate shots of acquisitions. The method further comprises determining spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan from calibration regions, wherein the overlapped k-space regions in the MRI data between separate shots are used as the calibration regions based on which spatially-varying spin phase fluctuations in separate shots are extracted. In particular, extracting the spatially-varying spin phase fluctuations in separate shots based on the overlapped k-space regions may be performed according to a method according to the first general aspect.

[0040] A further general subject of the invention is a computer program product, comprising a sequence of machine instructions, which causes a computer performing a method according to any one of the above mentioned aspects, when executing the sequence of machine instructions.

[0041] Another general subject of the invention is a medium on which an embodiment of the just mentioned computer program product is stored. The medium (e.g., non-transitory storage medium) may be magnetic (e.g., a floppy disk or a hard drive) or optical (e.g., a compact disk read only memory, or "CD ROM"), and may be read only or random access. If the computer program product is transmitted from a website, server, or other remote source using a physical cable, digital subscriber line (DSL), or wireless technologies then the physical cable, DSL, or wireless technologies such as infrared, radio, and microwave are included in the definition of medium.

[0042] A further general subject of the invention relates to a medium on which an embodiment of the above mentioned computer program product is stored and which is processable by the computer. A further general subject of the invention relates to a computer on which a computer program product according as mentioned above is stored and which is processable by the computer. Therefore, the computer may be configured to perform arithmetical, logical, and input / output operations, according to the machine instructions encoded in the computer program product. The computer may be configured to only process the imaging data acquired by any MRI scanner, or, according to another aspect, the computer is formed as a control device for an MRI scanning unit, i.e. the computer is directly coupled to an MRI scanning unit. Another general subject of the invention is an MRI scanner configured to perform a method according to any of the above mentioned aspects. The MRI scanner may comprise a control unit configured to perform a method according to any of the above mentioned aspects.

[0043] Brief description of the drawings

[0044] Further details and advantages of the invention are described in the following with reference to the attached drawings, which show in:

[0045] Figure 1: a flowchart illustrating a method of processing MRI data according to an embodiment of the invention;

[0046] Figure 2: a sequence illustration of a multi-shot readout mosaic segmented EPI without 2D navigator acquisitions with an option of diffusion weighting according to an embodiment of the invention;

[0047] Figure 3: a flowchart illustrating a method of processing MRI data according to a further embodiment of the invention;

[0048] Figure 4: an illustration of a linear system to estimate inter-shot phase fluctuation maps according to the embodiment of Fig. 3;

[0049] Figure 5: a flowchart illustrating a method of processing MRI data according to a further embodiment of the invention;

[0050] Figure 6A and 6B: an illustration of an eigenvalue approach to estimate EPI shot-to-shot variations according to the embodiment of Fig. 5;

[0051] Figure 7 an illustration of inverse filtering;

[0052] Figure 8A a comparison of uncorrected and corrected 7-shot readout segmented EPI images;

[0053] Figure 8B a comparison of uncorrected and corrected 7-shot readout segmented EPI images; and

[0054] Figure 9 an MRI device according to an embodiment of the invention.

[0055] Detailed description of preferred embodiments of the invention

[0056] The embodiments shown in the figures correspond, at least in part, so that similar or identical parts are labeled with the same reference numerals and reference is also made to the description of the other embodiments or figures for their explanation in order to avoid repetitions. Furthermore, for reasons of clarity, not all (identical) components appearing more than once have always been referenced separately. Figure 1 shows a flowchart illustrating a method of processing MRI data according to an embodiment of the invention.

[0057] The method comprises the step S100 of providing MRI data acquired by a multi-shot MRI data collection scan of an MRI scanner, wherein the MRI data comprises overlapped k-space regions commonly sampled in separate shots of acquisitions. The MRI scanner may be a conventional MRI scanner configured to perform a multi-shot MRI data collection.

[0058] The embodiments described here are based on MRI data acquired by an MRI data collection scan using a MRI pulse sequence of multi-shot readout-mosaic-segmented EPI, wherein the overlapped k-space data regions are intersections between overlapped EPI readout segment bands. The overlapped k-space bands acquired in separate shots, by way of example, comprising 8 to 30 overlapped pixels along readout dimension within overlapped k-space regions. If the phase encoding axis and RF dimension were counted as well, there would be e.g., about 8x200x32 overlapped pixels, or can be about 15x200x32 overlapped pixels. It is noted, however, that the method is not limited to EPI but it also applicable to other trajectory, such as a multi-shot non-Cartesian data collection scan with a small commonly sampled k-space regions between separate shots.

[0059] As explained in the introductory section, the spin phase errors need to be corrected when conducting multi-shot MRI. The method thus provides a technique for determining the spin phase fluctuations between separate shots of acquisitions.

[0060] The preferred embodiments of the invention are based on a navigator-free method for determining spin phase fluctuations between separate shots of acquisitions. Thus, contrary to many prior art approaches, the provided MRI data do not comprise any 2D or 3D navigator data that are additionally acquired in each shot of acquisition to monitor the spatially-varying spin phase fluctuations in each shot. In other words, the data collection scan does not comprise any 2D or 3D navigator acquisition before or after the imaging echoes within one shot of signal excitation and acquisition.

[0061] This is illustrated in Fig. 2 that shows a sequence illustration of a multi-shot readout mosaic segmented EPI without 2D navigator acquisitions with an option of diffusion weighting according to an embodiment of the invention. Fig. 2 shows an EPI echo train that samples a segment of k-space with low k-space coverage along the readout dimension, after RF excitation (here, spin echo type) and optionally, diffusion weighted encoding. Optional diffusion gradients 7 are shown. The phase blip areas 20 of the readout gradient in different shots are variable, to start signal readout at a corresponding k-space location (i.e., readout dimension) with a specific readout gradient polarity. Here, the readout gradient polarity is preferably flipped between shots of acquisitions for overlapped EPI segments. This can avoid spatially-varying T2 or T2* signal decay variations along different paths of sampling trajectory in the overlapped k-space pixels, with only the convolution effects by image-space phase fluctuations caused by shot-dependent microscopic motions of system imperfections remained for the extraction.

[0062] Fig. 2 further shows that before or after the imaging echoes, there is no additional 2D navigator acquisition before the next shot starts, the beginning of which is illustrated by the vertical dashed line.

[0063] By contrast, in a conventional multi-shot 2D EPI with 2D navigator acquisition (not shown in Fig. 2), after the first EPI echo train and the corresponding imaging echo acquisitions, there would be another refocusing pulse and another EPI echo train that is applied to sample 2D navigator echoes, acquiring a small central k-space band to estimate the phase fluctuation map in the current shot before the next shot start. As can be seen in Fig. 2, such a block for an additional 2D navigator acquisition is eliminated (i.e. not present), making the sequence in every shot substantially shorter. This can speed up long EPI scans with large slice number, diffusion weighted directions, and repetitions. As a result, time-consuming navigators can be avoided and a faster MRI technique is provided, e.g. about 30ms - 50 ms per shot of acquisition can be saved compared to a 2D navigator navigator-based multi-shot MRI technique. Additionally, the 2D navigator may fail due to inconsistency with imaging echoes (e.g., motions, T2 or T2* decay in navigator echoes). Thus, estimating phase fluctuations directly from imaging echoes avoid this problem, providing a simpler, more stable and more robust solution.

[0064] Returning to Fig. 1, in step S200, the spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan are determined from calibration regions in the provided MRI data.

[0065] An image-space phase error map during each shot can be mathematically described as a convolution kernel in k-space, similar to the RF receivers' sensitivity maps. These RF sensitivity maps, through decades of developments for parallel imaging [45-47], can be robustly estimated by extracting their corresponding convolution kernels from the low-resolution reference scans simultaneously acquired in all RF receivers, in the signal subspace spanned by the sensitivity maps (e.g., GRAPPA, ESPIRiT). Here, one key insight is that, these convolution kernels are shift-invariant across the entire k-space, not limited to the k-space central region as usually used for parallel imaging reference scans due to higher SNR. Therefore, the kernel for the relative phase maps can be estimated from overlapped k- space regions commonly sampled in separate shots of acquisitions. Thus, the overlapped k-space regions need not unnecessarily be confined to the k-space center. By contrast, the overlapped k- space regions comprise k-space regions centred around a k-space location different than the k- space centre.

[0066] Consequently, the overlapped k-space regions in the MRI data between separate shots (e.g., overlapped pixels between overlapped bands) are used as the calibration regions based on which spatially-varying spin phase fluctuations in separate shots are extracted. The obtained phase fluctuations may then be used for correcting for the finally reconstructed high-resolution images.

[0067] Thus, the method of the invention uses the overlapped k-space regions in the MRI data between separate shots as the calibration regions.

[0068] The proposed method may use mathematical frameworks and techniques of so-called MRI subspace methods [43,44] to estimate the relative signal variations between shots of acquisition. These mathematical frameworks and techniques applied here for extracting the spatially-varying spin phase fluctuations from the overlapped k-space regions are explained below.

[0069] From the recent sampling theory of parallel imaging

[0044] and nonlinear gradients

[0048] based on reproducing kernel Hilbert space, RKHS, the acquired time-domain (e.g., k-space) signals with distinct spatially (e.g., amplitude, phase) and temporally (e.g., field switching) encoded functions can be interpolated from one to another, with the interpolation weights called the cardinal function. This general term becomes the well-known "GRAPPA weights

[0047] " in parallel imaging, and the noise amplification factor due to a dynamic sampling "stamp" in image acceleration by nonlinear gradient modulations

[0048] ,

[0070] As noted above, according to the embodiments described herein, phase fluctuations estimation based on the shift-invariant kernel extraction is demonstrated using the MRI pulse sequence of readout segmented EPI, although the mathematical treatments can be applicable to other multishot MRI data collection scenarios, as long as the k-space bands acquired in separate shots overlap. In multi-shot EPI, within a region of k-space where two EPI bands acquired in separate shots intersect, the cardinal function can be seen as a kernel of coefficients representing the approximations of the Fourier transform of the relative image phase between shots, interpolating a neighborhood of data in one EPI band to a k-space point in another. Therefore, by establishing linear mapping relationships within their commonly sampled k-space locations as the calibration region, the shotdependent phase fluctuations can be estimated, e.g. by inverting the linear system, optionally with regularization to be more robust to noise.

[0071] Mathematically, the acquired time-domain MRI samples within intersection of EPI bands obtained in separate shots is the Fourier transform of the object with a shot-dependent relative phase map. Thus, f11(t) = T?m, (1), where fl!(t) is the time-domain (i.e., t) acquired k-space signals within the overlapping regions in overlapped k-space bands acquired from shot number 1 to N (N > 2), T is the Fourier transform operator, P = J’2- J’NP isavertically concatenated stack of diagonal matrix with diagonal elements the phase map values in each shot, m is the object image possibly modulated by RF receiver sensitivity detected by all available RF receivers. The maximum EPI shot number is larger or equal than the number of shot N corresponding to an overlapped k-space region. As shown in Fig. 4A, two overlapped readout EPI segments are shown, and thus, N = 2 without loss of generality.

[0072] In the RKHS perspective, as Equation (17)

[0044] or Equation (9)

[0048] in the references, a k-space neighborhood from one acquisition can be interpolated to a data point in another segment acquired in another acquisition shot with the cardinal function representing k-space modulation effects by image phase fluctuations, as illustrated in Fig. 4A: fk(t) = Zt,i fi ( with t,

[0073] Here, t and t represents the acquisition time in EPI segments separately obtained in different shots, i and k is the shot index for an EPI band contributing to the k-space overlapped region, C is the calibration region including the k-space overlapped pixels in all RF receivers, i (t) and fj^(t) are the MRI acquired signals in a localized neighborhood and a point shifted within the calibration region C, respectively. Note the notation k doesn't represent k-space location as in many other literatures. To estimate the cardinal function, Equation (2) can be formulated in linear algebra similar to Equation (4) in the reference

[0043] :

[0074] The calibration matrix A is constructed by sliding a (mathematical) window throughout the calibration region C, taking each block in ^(t) as a row in the matrix denoted by operator Rt.

[0075] In the following, based on the mathematical frameworks described above, two preferred embodiments are described for extracting the spatially-varying spin phase fluctuations from overlapped regions in overlapped k-space bands acquired in separate shots, using the calibration matrix A.

[0076] Figs. 3 and 4 describe a direct estimation of the k-space interpolation relationships, e.g., similar to a GRAPPA-type operation, whereas Figs. 5 and 6 describe an eigenvalue approach, e.g., similar to an ESPIRiT-type operation, connecting the method described in this specification and the two well- known RF parallel imaging techniques with mathematical similarity. It is noted that the calibration matrix illustrated in Fig. 4 is also used in the eigenvalue approach.

[0077] Both embodiments have in common that no additional 2D navigators are required to monitor the shot-dependent phase fluctuations, which substantially reduces pulse sequence duration, and eliminates possibility of navigator failures. Further, according to both embodiments, spin phase fluctuations are determined by estimating, from the MRI data in the overlapped k-space region, a kernel convolved with k-space signals, wherein the kernel corresponds to the Fourier transform of a relative image-space spin phase fluctuation map during each shot, and being shift-invariant across all k- space locations and radio frequency (RF) receivers' channels. Applied to the example of readout segmented EPI, the shot-dependent phase fluctuations are estimated from the overlapped k-space regions between overlapped bands. This is iterated over all overlapped k-space regions for different pairs of EPI bands. By way of example, only two EPI readout segments may overlap. Alternatively, more than two EPI bands may overlap. For example, the width of the kernel (or the localized cardinal function) determines the highest spatial resolution of the spin phase fluctuations estimated. The larger the width, the higher spatial resolution of spin phase fluctuations can be estimated, more shifts (e.g., larger overlapped region) may be needed to invert the linear equations system.

[0078] Fig. 3 shows a flowchart illustrating a method for according to an embodiment of the invention wherein a direct estimation of the k-space interpolation relationships is used to extract the spa- tially-varying spin phase fluctuations from overlapped regions in overlapped k-space bands acquired in separate shots.

[0079] In step S100 in Fig. 3, also as described in Fig. 1, MRI data acquired by a multi-shot MRI data collection scan of an MRI scanner is provided, wherein the MRI data comprises overlapped k-space regions 3 commonly sampled in separate shots of acquisitions.

[0080] Steps S210 to 250 in Fig. 3 describe an approach to directly estimate the k-space interpolation relationships in order to extract the spatially-varying spin phase fluctuations from overlapped regions in overlapped k-space bands acquired in separate shots.

[0081] In step 210 in Fig. 3, within the overlapped k-space region, at least one patch (k-space subset) of k- space data is selected, the patch being centered around the location ksin at least one other overlapped k-space data region acquired by different shots of the MRI data collection scan.

[0082] In step S220 in Fig. 3, within the overlapped k-space regions, a k-space data point at a location ksin a k-space data region is selected.

[0083] By way of example, Fig. 4A illustrates schematically the k-space interpolation in the calibration region, from a patch acquired in one shot to a single point in another shot. Fig. 4A illustrates schematically provided MRI data 2 that comprises overlapped k-space regions 3 commonly sampled in separate shots of acquisitions. In particular, the k-space interpolation between data in overlapped k- space regions 3, here as overlapped EPI bands (cf. Equation (2)), is shown, wherein two readout EPI segments (also referred to as EPI bands) are overlapped, and thus, N = 2 in the example of Fig. 4A.

[0084] In step S225, an interation of steps S210 and S220 is performed over all k-space locations and all RF receivers in an overlapped k-space region from respective shots.

[0085] In step S230 in Fig. 3, a linear interpolation relationship (linear system) (cf. Equation (3)) is then established interpolating the selected k-space data point selected in step S220 from the at least one patch of k-space data selected in step S210. Step 230 is performed for all pairs of overlapped k-space regions centered at different k-space locations. Establishing the linear interpolation relationship comprises determining a calibration matrix A (cf. Equation (3)) by selecting subsets of MRI data 2 from the overlapped k-space data regions 3. As described above, the calibration matrix A is determined by taking a patch of MRI data using a mathematical sliding window with the same size of the patch as a row in the calibration matrix A, iteratively sliding through all k-space locations and RF receiver channels 6 in the calibration region.

[0086] Fig. 4B is a graphical illustration of the entire linear system 4 (also referred to as linear equation system) describing the k-space interpolation relationships between two overlapped EPI segments 3, using two localized cardinal functions 5 with the finite width defined by the size of the sliding window.

[0087] Each localized cardinal function 5 estimated from two or multiple overlapped k-space data regions corresponds to a kernel of coefficients corresponding to approximations of the Fourier transform of the relative image-space spin phase fluctuations between separate shots of acquisitions, interpolating a patch of k-space data in one of the two or multiple overlapped k-space data regions, to a k-space data point in another one of the two or multiple overlapping k-space data regions.

[0088] To implement it, only one cardinal function 5 needs to be solved to estimate the relative phase fluctuations between shots. However, the mapping relations for the two cardinal functions 5 are included in this linear system illustration for mathematical completeness. In the following, the cardinal functions 5 are also referred to as u^'[t] cf. Equation 3. When solving the linear equation system 4, certain portions of the calibration matrix can be cropped out, which is also indicated by the zeros on the right side of the linear equation system 4, because they do not contribute to interpolating to a point in their own shot.

[0089] In step S240 in Fig. 3, this linear system is solved to obtain the (localized) cardinal functions u^'[t]. The linear system may be solved by direct inversion, to obtain the localized cardinal functions in k- space for correcting spatially-varying spin phase fluctuations between separate shots of acquisitions. In other words, the localized cardinal functions as the unknowns to be solved are estimated by inverting the linear equations system describing the linear interpolation relationship between the k-space data patches and the points in separate shots within the calibration region with optional regularizations (e.g., truncated singular value decompositions, TSVD).

[0090] In step S250 in Fig. 3, the image-space spin phase fluctuations are then estimate from the cardinal functions 5. For example, spin phase fluctuation maps in the corresponding separate shots of acquisitions may be obtained by flipping, taking complex conjugate operation, zero-filling and inverse Fourier-transforming the estimated cardinal functions into image-space. It is noted that whether a flipping and taking the complex conjugation operation is needed in a particular situation depends on how the Fourier transform and the convolution is defined.

[0091] The spin phase fluctuation maps in the corresponding separate shots of acquisitions can be used to correct the respective images as the inverse Fourier transform of the partial k-space data acquired in separate shots, iteratively performed for all shots. By optionally using truncated singular value decomposition (i.e., TSVD), the noise space of this linear system can be discarded to produce more accurate phase maps, which is particularly helpful for diffusion weighted multi-shot EPI due to its low SNR property.

[0092] In Fig. 3, the steps S210 to S240 are further illustrated using a more detailed example:

[0093] In Fig. 3, according to step S210, within an overlapped k-space region, a k-space data point is selected at a sampling time location from one shot of acquisition data, denoted as (tx,ty,tz, all s). For convenience only, the k-space data point may be selected as the center location of the patch. The notations tx,ty,tz correspond to the time-domain (i.e., t) sampled k-space locations in different spatial axis, s is the index of a RF receiver channel.

[0094] According to step S220 in Fig. 3, within the overlapped k-space region, a k-space data point is selected in a patch at a sampling time location from one shot, denoted as (tx-tl:tx+tl',ty-t2:ty+t2',tz- t3:tz+t3', for all s).

[0095] According to step S230 in Fig. 3, a linear mapping relationship is established between the patch and point in separate shots of acquisition. Theoretically, if more than two bands are overlapped, one point is selected in a particular shot, and several patches are selected from others.

[0096] The linear equation can be written in general formalism of the inverse problem: y = Ex, where y is a m x 1 vector with each row containing a point, E is A m x n matrix with each row containing a k- space patch. An iteration is performed over all k-space locations and all available RF receivers, wherein data are put in different rows in y and E, to make the linear system more well-conditioned. Thus, x can be calculated by inversion of the obtained linear system, with optional regularization (e.g., by the use of truncated singular value decomposition for regularization).

[0097] According to step S240, the cardinal functions, that are mapping from one of several shots to a particular shot (i.e., the reference shot) within one overlapped region, can be obtained. By Iterating over different overlapped regions, cardinal functions for interpolating overlapped k-space data for different pairs of overlapped k-space bands can be acquired.

[0098] With the cardinal functions all estimated by direct inversion of linear system for k-space interpolation, inter-shot phase fluctuations can subsequently be determined in step S250 in Fig. 3 according to two different approaches.

[0099] According to a first approach, if more than two bands are overlapped, one point is selected in a particular shot, and several patches are selected from others. The phase fluctuation maps can be used to correct the corresponding k-space bands in all shots.

[0100] A first option of the first approach comprises multiplying the complex conjugate of the estimated phase maps to the image corresponding to different partial k-space bands to compensate the relative phase fluctuations, followed by summing up all corrected bands. An inverse filter can then be applied to the summed-up k-space data to remove effects by k-space modulation transfer function (MTF) due to k-space bands overlapping. By performing a Fourier transform, a high-resolution image without artifacts due to shot-dependent phase fluctuations can be obtained.

[0101] A second option of the first approach comprises taking the estimated phase maps into a forward model written as y = Ex, where E contains the MRI encoding matrix (i.e., depending on the trajectory and optional RF receiver sensitivity maps) and the estimated phase map information, y contains the acquired uncorrected k-space data in shots. This linear system can then be solved to also obtain a high-resolution image x, without inverse filtering.

[0102] According to a second approach, the obtained cardinal functions are, optionally flipped and taken complex conjugate, convolved with the entire k-space bands (i.e., not limited to the calibration region) acquired in respective shots, to interpolate partial k-space data in each shot removing relative phase fluctuations compared to the data in the selected reference shot. This is performed iteratively for all overlapped bands properly, so that partial k-space data in all shots of acquisitions are corrected.

[0103] The corrected k-space bands (e.g., EPI segments) are then summed up. An inverse filtering is applied to undo a k-space modulation transfer function, MTF, due to k-space band overlapping. A high- resolution image is obtained without artifacts due to image-space phase fluctuations. According to a further embodiment of the invention, the cardinal functions can also be estimated directly from the calibration matrix A alone, similar to calibrating RF receivers' sensitivity from low- resolution fully-sampled central k-space region in ESPIRiT for parallel imaging

[0043] , Figs. 5 and 6 describe such an eigenvalue approach, e.g. as an ESPIRiT-type operation.

[0104] According to this embodiment, and as illustrated in Fig. 5, the method comprises again steps S100 to S230 as described earlier. However, instead of a direct estimation of the cardinal function, this embodiment uses an eigenvalue approach. According to this approach, the spatially-varying spin phase fluctuations in separate shots of acquisitions comprises performing eigenvalue decompositions on the calibration matrix to determine image-space eigenvector maps, wherein the determined image-space eigenvector maps comprises a primary eigenvector map associated with the largest eigenvalue, the phases of the primary eigenvector map correspond to the spin phase fluctuation maps. Thus, the eigenvector map associated with the largest eigenvalue are referred to here as primary eigenvector map.

[0105] The primary eigenvector map may be used to directly correct the image corresponding to a partial k-space band (i.e., a portion of k-space data) in each shot, or a so-called forward model may be used considering all acquired data and estimated eigenvector maps (in a linear system) to compute an image. This embodiment has the advantage that the spatially-varying spin phase fluctuations in separate shots can be extracted in a manner that is more robust to noise due to signal truncation in noise space.

[0106] By way of example, performing eigenvalue decompositions may comprise step S260 of decomposing the calibration matrix using a first singular value decomposition, SVD. This is followed by step S270, wherein resulting singular vectors in noise space corresponding to small singular values below a thresholding value are discarded. Then, in step S280, the remaining singular vectors in the signal space are transformed into image-space, e.g. by reshaping the signal space singular vectors into k- space filters, optionally flipping and taking complex conjugate, zero-filling, taking inverse Fourier transform.

[0107] Finally, in step S290, a second singular value decomposition, SVD, is performed on a matrix formed from data at each image pixel location across shot and eigen-index dimensions, to obtain singular vectors at each pixel location, containing signal fluctuation information. The singular vectors at all pixel locations correspond to image-space eigenvector maps determined from the eigenvalue decompositions of the calibration matrix. To show that the cardinal functions can also be estimated directly from the calibration matrix A alone using an eigenvalue approach, the null-space of the calibration matrix A can be indicated by rewriting Equation (3) as:

[0108] 0 = Auk'[t] - f{ t]

[0109] = A(u«[t] - etk), (4) where e^kis a vector with "1" in appropriate position that chooses the data of fk[t] in the calibration region, and "0" elsewhere. Therefore, shot-dependent phase fluctuation maps can also be extracted elegantly from signal space 7|| of the calibration matrix A alone. Similar to Equation (6-18) in ESPIRiT for parallel imaging

[0043] , the calibration matrix A can be decomposed by a singular value decomposition, SVD, operation:

[0110] A = USVH, (5)

[0111] The singular vectors V contains signal space V|| associated with singular values larger than a thresholding value and noise space V±associated with singular values smaller than a thresholding value. Then, by applying operation Wf11in Equation (6) which truncates A onto the signal space 7||, through removing singular vectors with below-thresholding singular values in noise space, the solution for estimating the image-space phase fluctuations appears to lie in the subspace spanned by the shot-dependent phase fluctuation maps: where M represents RtHRtwhich is equal to the number of samples in each k-space neighborhood selected by Rt. Expanding the k-space acquisition by Equation (1), point-wise operation at the image-space pixel q can be obtained: ’- ’J’m = J>m, (7) Therefore, in each shot of acquisition, the phase fluctuation map can be obtained by eigenvalue decomposition of all Qq(i.e., Qq= GqGq) along pixel q, and taking the eigenvectors corresponding to eigenvalue equal to one:

[0112] This can be implemented by reshaping the signal space matrix V|| into k-space filters, zero-filled and Fourier transform into image space. At each image pixel location, another SVD operation is performed on the formulated matrix Gqand finally, the resulted primary eigenvector is obtained as the relative signal fluctuation values between separate shots of acquisitions. The phases of the estimated maps are used to correct the phase errors.

[0113] Figs. 6A and 6B graphically illustrate the implementations of the derived Equation (1-9) to further illustrate the eigenvalue approach to estimate EPI shot-to-shot variations by computing the eigenvector of the image-space operator Qqpixelwise.

[0114] As shown in the upper left portion of Fig. 6A, the k-space data from the overlapped k-space regions of two overlapped EPI segment acquired in separate shots is used to formulate a calibration matrix A. By singular value decomposition, the signal subspace V|| is selected (cf. right portion of Fig. 6A), reshaped to k-space patches, and then zero-filled inverse Fourier transform to image space (as illustrated in the lower left portion of Fig. 6A).

[0115] Fig. 6B then illustrates that at each image pixel location in different shot numbers and eigenvector index, a matrix Gqcan be obtained. By performing another SVD operation on this matrix, and taking the primary singular vector corresponding to the largest singular value, the relative signal variations between shots can be obtained at each pixel. The phase of the estimated variations is taken, aligned between different pairs of overlapped segments, and used to correct all EPI segments in the image reconstruction.

[0116] Similar to ESPIRiT in parallel imaging, in non-ideal case (e.g., folder-over artifacts due to smaller FOV than the object), multiple sets (e.g., Mq) of phase maps 8 including from the non-primary eigenvectors in the final estimation step can be obtained according to Equation (10), and also used to reconstruct ghost-free image with a "relaxed" signal model. Therefore, M (e.g., usually M = 1 or 2) sets of eigenvector maps can be used for reconstruction. With the estimated relative phase maps 8 across multiple shots, the final image without ghost artifacts due to shot-to-shot variations can be reconstructed by taking the phase maps 8 in the forward model:

[0117] An alternative simpler manner is to multiply the complex conjugate phase maps 8 (i.e., conjugate operation *) to each uncorrected image m corresponding to a partial k-space band, sum over k- space signals from all bands, and apply an inverse filter L to remove k-space modulation transfer function, MTF, due to overlapped k-space region:

[0118] The k-space MTF and the inverse filter for a MRI sequence of 7-shot readout segmented EPI are shown in Figure 7. The upper portion of Fig. 7 shows the k-space MTF given the EPI band overlapping and the Kaiser filtering on each readout segment. The lower portion of Fig. 7 show the inverse filter calculated based on the MTF to remove the k-space modulation in the finally reconstructed image. As mentioned above, the conjugate phase maps 8 can be multiplied to the images corresponding to readout segments acquired in separate shots. A high-resolution image can be obtained by summing up corrected multi-shot images. In this way, an inverse filter is required to remove the MTF effects due to k-space overlapping. Additionally, to reduce the Gibbs ringing propagation which can lead to errors by multiplying with the estimated smoothly varying phase maps, low-resolution k-space readout segments can be low-pass filtered (e.g., Kaiser coefficient, 3.0). This k-space apodization in low-resolution acquisitions can also be eliminated by the inverse filter process.

[0119] The Eigenvalue approach as illustrated in Figs. 6A and 6B are further illustrated using a more detailed example:

[0120] First, a multi-shot MRI scan is performed, from which k-space data in all scans can be combined into a high-resolution image. The multi-shot MRI scan is performed without a 2D or 3D navigator scan in each acquisition to monitor the phase fluctuations in each shot. Again, by way of example, only, multi-slice readout-segmented EPI scans are used, with overlapped k-space regions from overlapped EPI bands. However, any multi-shot scans could be used, as long as they are conducted in a way that a certain area of overlapped k-space data is commonly sampled in separate shots of acquisitions, such as multi-shot readout segmented mosaic EPI or multi-shot spiral with a small densely sampled central k-space regions in all shots.

[0121] The shot-dependent phase fluctuations are estimated from the overlapped k-space regions between bands by iterating over all overlapped k-space regions for different pairs of overlapping EPI bands. For example, in the present example, only two EPI readout segments overlap. In general, this technique described herein also applies to cases where many EPI bands overlap.

[0122] For each overlapped k-space region as a calibration region, a small sliding window is used to acquire k-space data from a patch at a time location ts in a RF receiver channel. This patch can be generally denoted as (tx-tl:tx+tl,ty-t2:ty+t2,tz-t3:tz+t3, all s). Again, the notations tx,ty,tz correspond to the time-domain (i.e., t) sampled k-space locations in different spatial axis, s is the index of a RF receiver channel.

[0123] As a starting point, a calibration matrix A is initialed with all zero values. The patch data within the calibration region at a k-space location (tx,ty,tz) and at a RF receiver channel (i.e., index s) from all separately acquired bands is put as a single row of the matrix A. For the patch covering data from all separately acquired bands, an iteration is performed over different k-space locations (tx,ty,tz) and all available RF receiver channels (there can be one or many such channels). At each different k-space location and RF channel, the patches covering data in multiple shots are taken as a distinct row of the matrix A.

[0124] In a following step, a singular value decomposition, SVD, is performed on the calibration matrix A: A = USVH.

[0125] Among the diagonal component in resulting eigenvalue matrix S, indices with eigenvalues above a threshold are selected and used to select a certain number of rows in V (cf. equations (5) and (6) above. Therefore, V|| with only a limited number of rows in V (or limited number of columns in VH) is acquired. V|| is then reshaped into k-space patches, in different shots and eigenvector index. In other words, the resulting singular vectors in noise space corresponding to small singular values below a thresholding value are discarded and the the remaining singular vectors in the signal space are transformed into image-space, by reshaping the signal space singular vectors into k-space filters, and zero-filling and Fourier transforming into image space. At each image pixel location, data along shot and eigen index dimensions Gq (shots, eigen index) is taken to perform another singular value decomposition to obtain: Gq = U'S'V'H.

[0126] The singular vectors in U' are filled into each spatial pixel locations. A map of eigenvectors of imagespace operator at each image pixel is obtained (x,y,z, shots, eigen index). Similarly, the eigenvalues in S' can be used to form a map (x,y,z, shots, eigen index), indicating the image support (nonzero object, and zero background). It is noted that the eigenvectors at each pixel location can be phase aligned, by multiplication with the conjugate phase of the eigenvector in shot=l.

[0127] The phase of eigenvector maps at each pixel locations are the estimated spin phase fluctuation maps. Usually the primary eigenvector (with the biggest singular value in practical implementation) is used to produce a spatial map for different shots (x,y,z, shots, 1). A second set of eigenvector may also be used to produce a map in (x,y,z,shots,2).

[0128] The phase fluctuation maps can be used to correct the corresponding k-space bands in all shots. This can be achieved by multiplying the complex conjugate of the estimated phase maps to the image corresponding to different partial k-space bands and by summing up all corrected bands. Then, an inverse filter may be applied to the summed-up k-space. An inverse Fourier transform is then used to obtain the high-resolution image. By Fourier transforming this image-space phase fluctuation maps into k-space filters (kernels) which can be used to convolve with partial k-space bands in each shot, equivalent correction can also be done in k-space, similar to an option in the direct interpolation.

[0129] Alternatively, this can be achieved by taking the estimated phase maps into a forward model y=Ex, where E contains the conventional EPI encoding matrix and the estimated phase map information, y contains the acquired k-space data in shots. Solving this linear system can also obtain a high- resolution image x, without inverse filtering.

[0130] Figs 8A and 8B show in-vivo scans without and with diffusion weight, respectively, with only 8 pixels along readout dimension overlapped between separately acquired k-space bands. By way of example only, the in-vivo multi-shot EPI scans can be obtained by a 3T or 7T human MRI scanner, with and without diffusion weighting. Fig. 8A shows a comparison of an uncorrected and a corrected 7-shot readout segmented EPI image. The scans shown in Fig. 8A have been conducted without diffusion weighting. The left image of Fig. 8A shows an uncorrected image, obtained by a 2D multi-slice 7-shot readout segmented EPI with 0.7mm in-plane resolution. The right image of Fig. 8A shows a corrected version of the left image. The image on the right has been corrected using the techniques as described in this specification to demonstrate the high reconstruction quality invulnerable to shot-dependent phase errors which were estimated based on overlapped k-space regions commonly sampled in separate shots of acquisitions. Compared to uncorrected image, the corrected multi-shot reconstruction achieves significant resolution enhancement without ghost artefacts.

[0131] Similarly, Fig. 8B shows a comparison of uncorrected and corrected 7-shot readout segmented EPI images. The scans shown in Fig. 8B have been conducted with diffusion weighting (b= 1500 s / mmA2). The left image of Fig. 8B shows an uncorrected image (e.g., manifests as local blurring effects according to the characteristics of the readout segmented EPI sequence), obtained by a 2D multi-slice 7-shot readout segmented EPI with 0.7mm in-plane resolution. The right image of Fig. 8B shows a corrected version of the left image. The right image has been corrected using the techniques as described in this specification to demonstrate the reconstruction quality with shot-de- pendent phase errors estimated based on overlapped k-space regions commonly sampled in separate shots of acquisitions.

[0132] The reconstruction of the phase corrections as described in this specification represents only one component in the entire MRI reconstruction pipelines. By way of example, from raw RF receivers' data, the reconstruction steps can additionally include re-gridding, correction of gradient delay (using so-called ID navigator which is different than the 2D / 3D navigator mentioned above), parallel imaging reconstruction, and partial Fourier reconstruction.

[0133] Fig. 9 schematically shows an example of an MRI scanner 10 comprising a control device 11 and an MRI scanning unit 12 for imaging an object 1, e.g. a patient. The MRI scanning unit 12 may be a conventional MRI scanning unit 12 configured to perform a multi-shot MRI data collection. The control device 11 is configured to control the MRI scanning unit 12, i.e. sending machine in-struc- tions to the MRI scanning unit 12 scanner, and configured to perform a method according to the invention. Thus, the control device 11 is configure to control the MRI scanning unit so that MRI data 2 acquired by a multi-shot MRI data collection scan of the MRI scanning unit 12 are provided, wherein the MRI data 2 comprises overlapped k-space regions commonly sampled in separate shots of acquisitions. The control device 11 is further configured to determine spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan from calibration regions, wherein the overlapped k-space regions in the MRI data 2 between separate shots are used as the calibration regions based on which spatially-varying spin phase fluctuations in separate shots are extracted.

[0134] The connection between the scanning unit 12 and the control device 12 can be a wired, wireless, or any other type of data communication line, which allows for a transfer of information between the scanning unit 12 and the control device 12. However, it is also possible that the MRI data 2 acquired by the scanning unit 12 is transferred via a medium 6 or data network to another distant computer and processed there.

[0135] The features of the invention disclosed in the above description, the drawing and the claims can be of significance both individually as well as in combination or sub-combination for the realisation of the invention in its various embodiments.

[0136] References

[0137] 1. Ogawa S, Tank DW, Menon R, et al. Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging. Proc Natl Acad Sci.

[0138] 1992;89(13):5951-5955. doi:10.1073 / pnas.89.13.5951

[0139] 2. Kwong KK, Belliveau JW, Chesler DA, et al. Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc Natl Acad Sci. 1992;89(12):5675- 5679. doi:10.1073 / pnas.89.12.5675

[0140] 3. Bandettini PA, Wong EC, Hinks RS, Tikofsky RS, Hyde JS. Time course EPI of human brain function during task activation. Magn Reson Med. 1992;25(2):390-397. doi:10.1002 / mrm.1910250220

[0141] 4. Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophys J. 1994;66(l):259-267. doi:10.1016 / S0006-3495(94)80775-l

[0142] 5. Mansfield P. Multi-planar image formation using NMR spin echoes. J Phys C Solid State Phys. 1977;10(3):L55-L58. doi:10.1088 / 0022-3719 / 10 / 3 / 004

[0143] 6. Norris DG. Implications of bulk motion for diffusion-weighted imaging experiments: Effects, mechanisms, and solutions. J Magn Reson Imaging. 2001;13(4):486-495. doi:10.1002 / jmri,1072 7. Van de Moortele PF, Pfeuffer J, Glover GH, Ugurbil K, Hu X. Respiration-inducedBO fluctuations and their spatial distribution in the human brain at 7 Tesla. Magn Reson Med.

[0144] 2002;47(5):888-895. doi:10.1002 / mrm.10145

[0145] 8. Ehman RL, Felmlee JP. Adaptive technique for high-definition MR imaging of moving structures. Radiology. 1989;173(l):255-263. doi:10.1148 / radiology,173.1.2781017

[0146] 9. Anderson AW, Gore JC. Analysis and correction of motion artifacts in diffusion weighted imaging. Magn Reson Med. 1994;32(3):379-387. doi:10.1002 / mrm.1910320313

[0147] 10. Ordidge RJ, Helpern JA, Qing ZX, Knight RA, Nagesh V. Correction of motional artifacts in diffusion-weighted MR images using navigator echoes. Magn Reson Imaging. 1994;12(3):455- 460. doi:10.1016 / 0730-725X(94)92539-9

[0148] 11. De Crespigny AJ, Marks MP, Enzmann DR, Moseley ME. Navigated Diffusion Imaging of Normal and Ischemic Human Brain. Magn Reson Med. 1995;33(5):720-728. doi:10.1002 / mrm.1910330518

[0149] 12. Butts K, de Crespigny A, Pauly JM, Moseley M. Diffusion-weighted interleaved echo-planar imaging with a pair of orthogonal navigator echoes. Magn Reson Med. 1996;35(5):763-770. doi:10.1002 / mrm.1910350518

[0150] 13. Butts K, Pauly J, Crespigny AD, Moseley M. Isotropic diffusion-weighted and spiral-navigated interleaved EPI for routine imaging of acute stroke. Magn Reson Med. 1997;38(5):741-749. doi:10.1002 / mrm.1910380510

[0151] 14. Atkinson D, Porter DA, Hill DLG, Calamante F, Connelly A. Sampling and reconstruction effects due to motion in diffusion-weighted interleaved echo planar imaging. Magn Reson Med. 2000;44(l):101-109. doi:10.1002 / 1522-2594(200007)44:l<101::AID- MRM15>3.0.CO;2-S

[0152] 15. Atkinson D, Counsell S, Hajnal JV, Batchelor PG, Hill DLG, Larkman DJ. Nonlinear phase correction of navigated multi-coil diffusion images. Magn Reson Med. 2006;56(5):1135-1139. doi:10.1002 / mrm.21046

[0153] 16. Holdsworth SJ, Skare S, Newbould RD, Guzmann R, Blevins NH, Bammer R. Readout-Segmented EPI for Rapid High Resolution Diffusion Imaging at 3T. Ear J Radiol. 2008;65(l):36- 46. doi:10.1016 / j.ejrad.2007.09.016

[0154] 17. Porter DA, Heidemann RM. High resolution diffusion-weighted imaging using readout-segmented echo-planar imaging, parallel imaging and a two-dimensional navigator-based reacquisition. Magn Reson Med. 2009;62(2):468-475. doi:10.1002 / mrm.22024

[0155] 18. Heidemann RM, Porter DA, Anwander A, et al. Diffusion imaging in humans at 7T using read- out-segmented EPI and GRAPPA. Magn Reson Med. 2010;64(l):9-14. doi:10.1002 / mrm.22480 19. Liu C, Bammer R, Kim D hyun, Moseley ME. Self-navigated interleaved spiral (SNAILS): Application to high-resolution diffusion tensor imaging. Magn Reson Med. 2004;52(6):1388-1396. doi:10.1002 / mrm.20288

[0156] 20. Pipe JG. Motion correction with PROPELLER MRI: Application to head motion and free- breathing cardiac imaging. Magn Reson Med. 1999;42(5):963-969. doi:10.1002 / (SICI)1522- 2594(199911)42:5<963::AID-MRM17>3.0.CO;2-L

[0157] 21. Cheng JY, Alley MT, Cunningham CH, Vasanawala SS, Pauly JM, Lustig M. Nonrigid motion correction in 3D using autofocusing withlocalized linear translations. Magn Reson Med. 2012;68(6):1785-1797. doi:10.1002 / mrm.24189

[0158] 22. Uecker M, Karaus A, Frahm J. Inverse reconstruction method for segmented multishot diffusion-weighted MRI with multiple coils. Magn Reson Med. 2009;62(5):1342-1348. doi:10.1002 / mrm.22126

[0159] 23. Chen N kuei, Guidon A, Chang HC, Song AW. A robust multi-shot scan strategy for high-reso- lution diffusion weighted MRI enabled by multiplexed sensitivity-encoding (MUSE). NeuroImage. 2013;72:41-47. doi:10.1016 / j. neuroimage.2013.01.038

[0160] 24. Mani M, Jacob M, Kelley D, Magnotta V. Multi-shot sensitivity-encoded diffusion data recovery using structured low-rank matrix completion (MUSSELS). Magn Reson Med. 2017;78(2):494-507. doi:10.1002 / mrm.26382

[0161] 25. Chen X, Wu W, Chiew M. Motion compensated structured low-rank reconstruction for 3D multi-shot EPI. Magn Reson Med. 2024; 91: 2443-2458. doi:10.1002 / mrm.30019

[0162] 26. Hu Y, Levine EG, Tian Q, et al. Motion-robust reconstruction of multishot diffusion-weighted images without phase estimation through locally low-rank regularization. Magn Reson Med. 2019;81(2):1181-1190. doi:10.1002 / mrm.27488

[0163] 27. Hennel F, Pruessmann KP. MRI with phaseless encoding: MRI with Phaseless Encoding. Magn Reson Med. 2017;78(3):1029-1037. doi:10.1002 / mrm.26497

[0164] 28. Hennel F, Tian R, Engel M, Pruessmann KP. In-plane "superresolution" MRI with phaseless sub-pixel encoding: Hennel et al. Magn Reson Med. 2018;80(6):2384-2392. doi:10.1002 / mrm.27209

[0165] 29. Tian R, Hennel F, Pruessmann KP. Low-distortion diffusion tensor MRI with improved phaseless encoding. J Magn Reson. 2019;309:106602. doi:https: / / doi.org / 10.1016 / j.jmr.2019.106602

[0166] 30. Lally PJ, Matthews PM, Bangerter NK. Unbalanced SSFP for super-resolution in MRI. Magn Reson Med. 2021;85(5):2477-2489. doi:10.1002 / mrm.28593

[0167] 31. Tian R, Hennel F, Bianchi S, Pruessmann KP. Superresolution MRI with a Structured-Illumination Approach. Phys Rev Appl. 2023;19(3):034074. doi:10.1103 / PhysRevApplied.19.034074 32. Hell SW, Kroug M. Ground-state-depletion fluorscence microscopy: A concept for breaking the diffraction resolution limit. Appl Phys B Lasers Opt. 1995;60(5):495-497. doi:10.1007 / BF01081333

[0168] 33. Heintzmann R, Cremer CG. Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating. In: Optical Biopsies and Microscopic Techniques III. Vol 3568. International Society for Optics and Photonics; 1999:185-197. doi:10.1117 / 12.336833

[0169] 34. Gustafsson MGL. Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy. J Microsc. 2000;198(2):82-87. doi:10.1046 / j,1365- 2818.2000.00710.x

[0170] 35. Frohn JT, Knapp HF, Stemmer A. True optical resolution beyond the Rayleigh limit achieved by standing wave illumination. Proc Natl Acad Sci. 2000;97(13):7232. doi:10.1073 / pnas.130181797

[0171] 36. Heintzmann R, Cremer C, Thomas J. Saturated patterned excitation microscopy— a concept for optical resolution improvement. J Opt Soc Am A. 2002;19(8).

[0172] 37. Heintzmann R. Saturated patterned excitation microscopy with two-dimensional excitation patterns. Micron. 2003;34(6-7):283-291. doi:10.1016 / S0968-4328(03)00053-2

[0173] 38. Gustafsson MGL. Nonlinear structured-illumination microscopy: Wide-field fluorescence imaging with theoretically unlimited resolution. Proc Natl Acad Sci. 2005;102(37):13081- 13086. doi:10.1073 / pnas.0406877102

[0174] 39. Betzig E, Patterson GH, Sougrat R, et al. Imaging Intracellular Fluorescent Proteins at Nanometer Resolution. Science. 2006;313(5793):1642-1645. doi:10.1126 / science.1127344

[0175] 40. Rust MJ, Bates M, Zhuang X. Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM). Nat Methods. 2006;3(10):793-796. doi:10.1038 / nmeth929

[0176] 41. Bretschneider S, Eggeling C, Hell SW. Breaking the Diffraction Barrier in Fluorescence Microscopy by Optical Shelving. Phys Rev Lett. 2007;98(21):218103. doi:10.1103 / PhysRevLett.98.218103

[0177] 42. Chmyrov A, Keller J, Grotjohann T, et al. Nanoscopy with more than 100,000 "doughnuts." Nat Methods. 2013;10(8):737-740. doi:10.1038 / nmeth.2556

[0178] 43. Uecker M, Lai P, Murphy MJ, et al. ESPIRiT-an eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA. Magn Reson Med. 2014;71(3):990-1001. doi:10.1002 / mrm.24751

[0179] 44. Athalye V, Lustig M, Martin Uecker. Parallel magnetic resonance imaging as approximation in a reproducing kernel Hilbert space. Inverse Probl. 2015;31(4):045008. doi:10.1088 / 0266- 5611 / 31 / 4 / 045008 Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): Fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997;38(4):591-603. doi:10.1002 / mrm.1910380414 Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: Sensitivity encoding for fast MRI. Magn Reson Med. 1999;42(5):952-962. doi:10.1002 / (SICI)1522-

[0180] 2594(199911)42:5<952::AID-MRM16>3.0.CO;2-S Griswold MA, Jakob PM, Heidemann RM, et al. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med. 2002;47(6):1202-1210. doi:10.1002 / mrm.10171 Tian R, Uecker M, Davids M, et al. Accelerated 2D Cartesian MRI with an 8-channel local BO coil array combined with parallel imaging. Magn Reson Med. 2024;91(2):443-465. doi:10.1002 / mrm.29799 Layton KJ, Kroboth S, Jia F, et al. Pulseq: A rapid and hardware-independent pulse sequence prototyping framework. Magn Reson Med. 2017;77(4):1544-1552. doi:10.1002 / mrm.26235

Claims

Claims1. Method of processing magnetic resonance imaging, MRI, data, comprising the steps of: providing (S100) MRI data acquired by a multi-shot MRI data collection scan of an MRI scanner, wherein the MRI data (2) comprises overlapped k-space regions (3) commonly sampled in separate shots of acquisitions; and determining (S200) spin phase fluctuations between separate shots of acquisitions of the multi-shot MRI data collection scan from calibration regions, wherein the overlapped k-space regions (3) in the MRI data (2) between separate shots are used as the calibration regions based on which spatially-varying spin phase fluctuations in separate shots are extracted.

2. The method according to claim 1, wherein the spin phase fluctuations are determined by estimating, from the MRI data (2) in the overlapped k-space regions (3), a kernel convolved with k-space signals, wherein the kernel corresponds to the Fourier transform of a relative image-space spin phase fluctuation map (8) during each shot, and being shift-invariant across all k-space locations and radio frequency (RF) receivers channels (6).

3. The method according to claim 1 or 2, wherein the overlapped k-space regions (3) comprises k-space regions centred around a k-space location different than the k-space centre.

4. The method according to an one of the preceding claims, wherein extracting the spatially- varying spin phase fluctuations in separate shots comprises: establishing (S230) a linear interpolation relationship (4) interpolating a k-space data point at a location ksin a k-space data region (3), from at least one patch of k-space data centered around the location ksin at least one other overlapped k-space data region (3) acquired by different shots of the MRI data collection scan, the linear interpolation relationship being solved a) by direct inversion, preferably with regularizations, to obtain localized cardinal functions (5) in k-space, or b) by eigenvalue decompositions with thresholding in eigenvector space, to obtain image-space spin phase fluctuations map (8) for correcting spatially-varying spin phase fluctuations between separate shots of acquisitions.

5. The method according to claim 4, wherein establishing a linear interpolation relationship (4) comprises determining a calibration matrix (A) by selecting subsets of MRI data from the calibration region.

6. The method according to claim 5, wherein the calibration matrix (A) is determined by taking a patch of MRI data using a mathematical sliding window with the same size of the patch as a row in the calibration matrix (A), iteratively sliding through all k-space locations in the calibration region and RF receivers channels (6).

7. The method according to any one of the claims 4 to 6, wherein each localized cardinal function (5) estimated from two or multiple overlapped k- space data regions (3) corresponds to a kernel of coefficients corresponding to approximations of the Fourier transform of the relative image-space spin phase between separate shots of acquisitions, interpolating a patch of k-space data in one of the two or multiple overlapped k-space data regions, to a k-space data point in another one of the two or multiple overlapping k-space data regions (3).

8. The method according to any one of claims 4 to 7, wherein the localized cardinal functions (5) as the unknowns to be solved are estimated by inverting the linear equations system (4) describing the linear interpolation relationship between the k-space data patches and the points in separate shots within the calibration region, and wherein spin phase fluctuation maps in the corresponding separate shots of acquisitions a) are obtained by optionally flipping, taking complex conjugate operation, zero-filling and inverse Fourier-transforming the estimated cardinal functions into image-space; and / or b) are used to correct the respective images as the inverse Fourier transform of the k-space data acquired in separate shots, iteratively performed for all shots.

9. The method according to any one of claims 4 to 7, further comprising correcting k-space data in all shots of acquisitions by optionally first flipping and taking complex conjugate operation of the estimated cardinal functions, convolving each of these cardinal functions with all the k-space data acquired in the respective shot including but not limited to the overlapped k-space region, iteratively performed for all shots.

10. The method according to claim 5 or claim 6, wherein extracting the spatially-varying spin phase fluctuations in separate shots of acquisitions comprises: performing eigenvalue decompositions on the calibration matrix to determine image-space eigenvector maps, wherein the determined image-space eigenvector maps comprises a primary eigenvector map associated with the largest eigenvalue, the phases of the primary eigenvector map correspond to the spin phase fluctuation maps.

11. The method according to claim 10, wherein the step of performing eigenvalue decompositions comprises: decomposing (S260) the calibration matrix using a first singular value decomposition, SVD; discarding (S270) resulting singular vectors in noise space corresponding to small singular values below a thresholding value; transforming (S280) the remaining singular vectors in the signal space into image-space, by reshaping the signal space singular vectors into k-space filters, optionally flipping and taking complex conjugate, zero-filling, taking inverse Fourier transform; and performing (S290) a second SVD on a matrix formed from data at each image pixel location across shot and eigen-index dimensions, to obtain singular vectors at each pixel location, containing signal fluctuation information, the singular vectors at all pixel locations correspond to image-space eigenvector maps determined from the eigenvalue decompositions of the calibration matrix.

12. The method according to claim 10 or 11, wherein an amplitude of the image-space eigenvector maps is used to correct possible amplitude fluctuations between separate shots of acquisitions; and / or wherein a secondary eigenvector obtained from image-space eigenvector maps is used as the secondary amplitude and / or phase fluctuation maps for correcting k-space data or images in separate shots jointly with the primary ones.

13. The method according to any one of the preceding claims, wherein the MRI data (2) based on which the spin phase fluctuations between shots are determined do not comprise any 2D or 3D navigator data that are additionally acquired in each shot of acquisition to monitor the spin phase fluctuations in each shot.

14. The method according to any one of the preceding claims, wherein the MRI data (2) are acquired by a MRI data collection scan using a MRI pulse se-quence of multi-shot readout-mosaic-segmented Echo Planar Imaging, EPI, and wherein the overlapped k-space data regions are intersections between overlapped EPI readout segment bands, preferably comprising 8 to 30 overlapped pixels along readout dimension between overlapped k- space bands.

15. The method according to claim 14, wherein a readout gradient polarity is flipped between shots of acquisitions for overlapped EPI readout segments.

16. The method according to any one of the preceding claims, wherein the provided MRI data are acquired by one of the following MRI data collection scans: a functional MRI data collection scan, or a diffusion weighted MRI data collection scan; and / or and wherein the provided MRI data are acquired by one of the following MRI pulse sequences: a multi-shot Cartesian EPI data collection scan or a multi-shot non-Cartesian data collection scan.

17. Method of magnetic resonance imaging, MRI, an object, comprising the steps of: providing an MRI scanner (10); acquiring MRI data (2) acquired by a multi-shot MRI data collection scan of the MRI scanner (10), wherein the MRI data (2) comprises overlapped k-space regions commonly sampled in separate shots of acquisitions; and determining spin phase fluctuations between separate shots of acquisitions of the multishot MRI data collection scan from calibration regions, wherein the overlapped k-space regions in the MRI data between separate shots are used as the calibration regions based on which spatially- varying spin phase fluctuations in separate shots are extracted.

18. Computer program product comprising a sequence of machine instructions which causes a computer performing a method according to any one of the preceding claims, when executing the sequence of machine instructions.

19. Medium on which a computer program product according to claim 18 is stored.

20. Computer on which a computer program product according to claim 18 is stored and which is processable by the computer.

21. Computer according to claim 20, characterized in that, it is formed as a control device (11) for an MRI scanner (10).

22. Magnetic Resonance Imaging, MRI, scanner (10) comprising a control device (11) which is configured to perform a method according to any one of claims 1-16 or a computer according to claim 20 or 21.

Images