Highly robust fast imaging method and device combining k-space and image space reconstruction

Abstract

The present invention discloses a highly robust fast imaging method and device combining k-space and image space reconstruction, belonging to the field of magnetic resonance imaging. The method reconstructs variable acceleration undersampled image frames formed in the magnetic resonance parameter measurement of the object under test: first, the k-space automatic calibration method is used to reconstruct the undersampled image frames with lower acceleration factors with automatic calibration signals, and then the automatically optimized sensitivity map and correction factor map are derived therefrom to reconstruct other undersampled image frames with higher acceleration factors. The present invention can combine the traditional sensitivity map with the sensitivity map obtained by subspace model decomposition to achieve automatic optimization, thereby achieving highly robust and accurate reconstruction from highly undersampled data, and the speed is four times or more than that of ordinary methods. As a self-calibrated reconstruction method and without the need for full sampling frames, the method is particularly suitable for accelerating three-dimensional multi-frame magnetic resonance imaging.

Description

Technical Field

[0001] The invention belongs to the field of magnetic resonance imaging and can realize highly robust and fast imaging for magnetic resonance parameter measurement requiring multi-frame image data. Background Art

[0002] Magnetic resonance imaging (MRI) can be used to measure various physiological parameters and has been shown to be useful for quantitatively assessing various pathological information, such as longitudinal relaxation time (T1), transverse relaxation time (T2), apparent diffusion coefficient (ADC), amide proton transfer (APT), etc. Quantitative measurement of physiological parameters requires the acquisition of multiple data sets with pulse sequence parameter modulation, such as echo time (TE), flip angle (FA), inversion time (TI), or saturation frequency. The signal evolution in the multi-part data set is then fitted to a parametric model to achieve pixel-by-pixel parameter estimation. The long scanning time for MRI parameter measurement due to multiple data acquisitions limits its clinical application. Long acquisition time not only reduces patient throughput, but may also affect the effect of parameter quantification due to patient motion.

[0003] In order to solve the limitation of long scanning time, many researchers have developed related methods to accelerate multi-frame imaging. Taking T1 and T2 quantification as an example, there are compressed sensing (CS), REPCOM, LLR-SPIRiT, kt SLR, kt PCA, magnetic resonance fingerprint (MRF), ALOHA and other methods. Compressed sensing takes advantage of the sparsity of images in specific transform domains, such as wavelet transform, and some technologies have been successfully applied to MR parameter mapping. Since time dependence can be modeled as the low rank of spatiotemporal data, low rank is usually used as a mandatory constraint for optimization problems. kt PCA uses matrix decomposition to decompose the training data into a set of principal components, from which parameters are estimated. MRF uses a simulation-based dictionary and model matching algorithm to estimate multiple parameter mappings simultaneously.

[0004] Parallel imaging is the most widely used acceleration method in clinical practice, among which GRAPPA and SENSE are the most widely used. GRAPPA models the missing points in k-space as a linear combination of adjacent acquired data points of all channels; SENSE regards the reconstruction problem as a linear inversion in image space. Theoretically, if the obtained sensitivity map is accurate, the SENSE method is optimal and can obtain an exact solution in the sense of least squares. However, the sensitivity map is usually inaccurate. Recently, a variable acceleration sensitivity encoding method (vSENSE) applied in the field of multi-frame imaging has been proposed to accelerate parameter quantification. By assuming that the sensitivity maps of different frames in multi-frame imaging are consistent, vSENSE can achieve a higher multiple of acceleration. As an improvement on the vSENSE method, the KIPI method realizes automatic calibration, does not require the acquisition of additional sensitivity maps, does not require full sampling frames, and allows more flexible variable acceleration strategies. However, the sensitivity maps obtained by the KIPI method are not always accurate, especially in low-signal areas of the image, such as the scalp area in the inversion recovery (IR) relaxation time mapping sequence. Errors in these areas will cause signal aliasing in the image, resulting in artifacts.

[0005] Therefore, in the field of multi-frame imaging, it is of great significance to propose a method that can automatically optimize the sensitivity map. Summary of the invention

[0006] The purpose of the present invention is to solve the problems of long scanning time of multi-frame imaging technology and inaccurate coil sensitivity map of existing acceleration technology, and to provide a highly robust fast imaging method combining k-space and image-space reconstruction (English name is joint K-space and Image-space Parallel Imaging with high Robustness, hereinafter referred to as RobustKIPI method) and device, which can realize automatic optimization of sensitivity map and accurately reconstruct images from highly under-sampled data.

[0007] The specific technical solutions adopted by the present invention are as follows:

[0008] In a first aspect, the present invention provides a highly robust fast imaging method combining k-space and image space reconstruction, which is used for reconstructing undersampled image frames in multi-frame imaging, wherein the undersampled image frames include a first undersampled image frame with an automatic calibration signal and a second undersampled image frame with an acceleration factor not less than 2, and the acceleration factor of the first undersampled image frame is not higher than the acceleration factor of the second undersampled image frame;

[0009] The imaging method steps are as follows:

[0010] S1: for the first undersampled image frame, a complete image frame is reconstructed using a k-space automatic calibration reconstruction method as a calibration frame; the k-space data of the calibration frame is Fourier transformed to obtain a channel image of each channel, and all channel images are merged to obtain a channel merged image;

[0011] S2: Using the ESPIRiT method on the automatic calibration signal of the first undersampled image frame, obtaining the ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map, and identifying the first support area that completely covers the target imaging area based on the ESPIRiT eigenvalue map; dividing each channel image by the channel merged image to obtain the original coil sensitivity map of each channel; identifying the second support area after removing the noise and the hole area in the first support area based on the channel merged image; the ESPIRiT coil sensitivity map and the original coil sensitivity map are both composed of corresponding amplitude maps and phase maps;

[0012] S3: adding the sensitivity amplitude map of the ESPIRiT coil in the first support area but not in the second support area to the original coil sensitivity amplitude map in the second support area, and obtaining an optimized sensitivity amplitude map after smoothing filtering; interpolating and completing the noise and hole areas based on the original coil sensitivity phase map in the second support area to obtain an optimized sensitivity phase map; combining the optimized sensitivity amplitude map and the optimized sensitivity phase map to finally obtain an optimized sensitivity map;

[0013] S4: using the optimized sensitivity map, performing SENSE reconstruction with an acceleration factor of 1 on the calibration frame to generate a reference image; at the same time, retrospectively undersampling the calibration frame, and the acceleration factor of the retrospective undersampling is the same as the acceleration factor of the second undersampled image frame, and then performing SENSE reconstruction on the retrospectively undersampled data to generate a retrospective reconstructed image with potential aliasing artifacts; obtaining a correction factor map based on the reference image and the retrospectively reconstructed image, and the pixel value of each position in the correction factor map is the quotient of the pixel value of the corresponding position in the reference image and the retrospectively reconstructed image;

[0014] S5: For the remaining second under-sampled image frames except the first under-sampled image frame, perform SENSE reconstruction using the optimized sensitivity map, and multiply the SENSE reconstructed image with the correction factor map to obtain a complete image frame with suppressed artifacts.

[0015] As a preference of the first aspect, both the first under-sampled image frame and the second under-sampled image frame are two-dimensional images; the acceleration factor of the first under-sampled image frame is preferably 2, and the acceleration factor of the second under-sampled image frame is preferably 2 to 4.

[0016] As a preference of the first aspect, both the first under-sampled image frame and the second under-sampled image frame are three-dimensional images; the acceleration factor of the first under-sampled image frame is preferably 4, and the acceleration factor of the second under-sampled image frame is preferably 4 to 12.

[0017] As a preferred embodiment of the first aspect, in S1, the reconstruction method for automatic k-space calibration is preferably the GRAPPA method, by calculating the coil weights and applying them to the under-sampling area to reconstruct a complete image frame; preferably, the GRAPPA method is GRAPPA with Tikhonov regularization.

[0018] As a preferred embodiment of the first aspect, in S1, all coil images are merged by a square root method or an adaptive method.

[0019] As a preferred embodiment of the first aspect, in S2, the ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map are both a set of ESPIRiT coil sensitivity maps and ESPIRiT eigenvalue maps corresponding to the maximum eigenvalue after eigenvalue decomposition of the ESPIRiT method (each channel has one ESPIRiT coil sensitivity map, but all channels correspond to one ESPIRiT eigenvalue map).

[0020] As a preferred embodiment of the first aspect, in S4, the correction factor map needs to be filtered to remove abnormal values.

[0021] As a preferred embodiment of the first aspect, in S5, when performing SENSE reconstruction using the optimized sensitivity map, a truncated singular value regularization method is adopted.

[0022] In a second aspect, the present invention provides a magnetic resonance imaging data processing device, which includes a memory and a processor;

[0023] The memory is used to store computer programs;

[0024] The processor is used to implement the highly robust fast imaging method combining k-space and image space reconstruction as described in any one of the first aspects when executing the computer program.

[0025] In a third aspect, the present invention provides a magnetic resonance imaging device, comprising a magnetic resonance scanner and a control unit, wherein the control unit stores a computer program, and when the computer program is executed, it is used to implement a highly robust fast imaging method combining k-space and image space reconstruction as described in any one of the first aspects; the magnetic resonance scanner is used to acquire the first under-sampled image frame and the second under-sampled image frame data.

[0026] Compared with the prior art, the present invention has the following beneficial effects:

[0027] The prior art requires accurate identification of noise areas and signal areas, and requires a significant extrapolation of the sensitivity map, which usually results in a high number of outliers in the sensitivity map. The present invention can automatically fill in the coil sensitivity map of the area around the hole or the signal through the coil sensitivity map generated by the ESPIRiT method, thereby achieving automatic optimization.

[0028] The prior art is limited by the data signal-to-noise ratio. When the signal-to-noise ratio is poor, the coil sensitivity map is also poor. The present invention introduces the ESPIRiT method to combine the sensitivity map obtained by direct division with the sensitivity map obtained by model eigenvalue decomposition to produce a more robust and accurate coil sensitivity map.

[0029] Compared with the prior art, the present invention improves the image reconstruction effect while ensuring the acceleration effect. The reconstructed image is consistent with the true value and has no obvious artifacts. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 Flowchart for combining the original coil sensitivity amplitude map and the ESPIRiT coil sensitivity amplitude map for RubustKIPI. The first support area A is the sensitivity range to be formed; the second support area B is the remaining area after excluding the noise and hole areas of the original sensitivity map in the first support area A, and the difference between area A and area B is the noise and hole area C. The combined coil sensitivity amplitude map is the combination of the original sensitivity amplitude map in area B and the ESPIRiT coil sensitivity amplitude map in area C.

[0031] Figure 2 Magnitude and phase of three different coil sensitivity maps (channel 7). Left column: original coil sensitivity map, dividing the image of each channel by the RMS image. Middle: coil sensitivity map calculated with ESPIRiT. Right: combining the original and ESPIRiT to get the combined sensitivity map and fitting. White arrows indicate areas of incorrect sensitivity map due to noise or holes.

[0032] Figure 3Reconstructed water phantom source images and corresponding T1 parameter images obtained from undersampled IR TSE data; (A): IR TSE images (TI = 300ms) were reconstructed from 4x undersampled data using GRAPPA, ESPIRiT, KIPI, and RubustKIPI. KIPI uses the original (Raw) sensitivity map and the ESPIRiT sensitivity map, respectively. (B): The error between the full sampled true value and the reconstructed image, with nRMSE at the bottom. (CD): T1 parameter image (C) and reconstruction error (D). The arrows in (AB) indicate artifacts in the source image, and the white arrows in (C) indicate artifacts in the T1 parameter image. The source image error image is magnified 20 times, and the T1 parameter image error map is magnified 10 times. nRMSE: Normalized root mean square error.

[0033] Figure 4 IR-TSE images with TI=50ms reconstructed using GRAPPA, ESPIRiT, KIPI and RubustKIPI (A). The second row (B) shows the difference from the full sampled true value. KIPI uses the original (Raw) sensitivity map and the ESPIRiT sensitivity map respectively. The gray arrows indicate the aliasing artifacts of the GRAPPA and ESPIRiT results. The white arrows indicate relatively large errors, and the source image error image is the result of 20 times magnification. The nRMSE value is shown at the bottom.

[0034] Figure 5 For Figure 4 The T1 parameter images (A) are calculated from the source images shown. The second row (B) shows the difference from the fully sampled true values. KIPI uses the Raw sensitivity map and the ESPIRiT sensitivity map, respectively. The white arrows indicate relatively large artifacts in KIPI using the Raw and ESPIRiT coil sensitivity maps. The T1 parameter image error map is a 10x magnification. The nRMSE values ​​are shown at the bottom.

[0035] Figure 6 : is an algorithm flow chart of the RubustKIPI method of the present invention, wherein the first undersampled image frame that completes the reconstruction is the calibration frame. DETAILED DESCRIPTION

[0036] The present invention is further described and illustrated below in conjunction with the accompanying drawings and specific embodiments. The technical features of each embodiment of the present invention can be combined accordingly without conflicting with each other.

[0037] In the field of magnetic resonance parameter measurement that requires multi-frame imaging, the data set consists of a series of images acquired from the same volume with different pulse sequence parameters, which means that the coil sensitivity maps of different frames are theoretically the same. The traditional SENSE method explicitly uses the sensitivity map from a separate reference scan, while the GRAPPA method implicitly utilizes sensitivity encoding through the k-space correlation estimated by the automatic calibration signal (ACS). However, it is difficult to obtain accurate sensitivity maps in SENSE reconstruction, which limits its maximum acceleration factor in practice. For GRAPAA, all image frames need to repeatedly acquire ACS data, resulting in a relatively long scan time.

[0038] The recently proposed KIPI method utilizes variable acceleration factors (AF) of different CEST frames and combines the advantages of SENSE and GRAPPA. KIPI samples CEST frames with a low acceleration factor (e.g., AF=2), fills the unsampled k-space with GRAPPA, obtains sensitivity maps from the full k-space generated by GRAPPA, and then reconstructs other highly undersampled (AF>2) frames using SENSE and an artifact suppression algorithm. The coil sensitivity map obtained by dividing the image of each channel of KIPI by the combined image of the channels is not always accurate, especially in some sequences where there are holes, such as the scalp area in the inversion recovery sequence; in applications with low signal-to-noise ratio, the sensitivity map obtained by KIPI is also relatively inaccurate.

[0039] The ESPIRiT method performs SVD decomposition on the constructed calibration matrix A and obtains the subspace and null space by truncating the smaller singular values. The ESPIRiT model regards the reconstruction problem as a subspace completion problem, and finally there is a small calibration equation for each pixel position, as follows:

[0040] ψ q csm q =csm q [1]

[0041] The subscript q represents the position, is the correlation matrix here, is the coil sensitivity vector at the pixel position, and Nc is the number of channels. Therefore, the coil sensitivity vector csm q is the correlation matrix ψ q The eigenvalue of is 1, corresponding to the eigenvector. qBy performing eigenvalue decomposition, Nc eigenvalues ​​and corresponding Nc eigenvectors can be obtained. All eigenvalues ​​will not be greater than 1, so there are Nc sets of coil sensitivity maps, but the one with the largest eigenvalue is relatively accurate. By iterating all pixel positions, the ESPIRiT coil sensitivity map decomposed by the model can be obtained. It has the characteristics of smoothness and insensitivity to noise, but the ESPIRiT coil sensitivity map is relatively fuzzy and the phase information is inaccurate.

[0042] The present invention combines the KIPI sensitivity map and the ESPIRiT coil sensitivity map to obtain a more accurate and robust sensitivity map for reconstruction. The process of the highly robust fast imaging method combining k-space and image space reconstruction provided by the present invention is as follows: Figure 6 As shown, the specific implementation process is described below.

[0043] The imaging method reconstructs undersampled image frames acquired by a magnetic resonance machine. It is divided into a first undersampled image frame, which is undersampled with a lower acceleration factor and contains ACS; and a second undersampled image frame, which is undersampled with a higher acceleration factor (AF>2). The acceleration factor of the second undersampled image frame is not lower than that of the first undersampled image frame. There is only one first undersampled image frame, and the rest are second undersampled image frames. The first undersampled image frame can be used to obtain an accurate sensitivity map and correction factor map, and the second undersampled image frame that is undersampled with a higher acceleration factor can be accurately reconstructed, achieving a balance between image quality and acquisition speed.

[0044] It should be noted that the multi-frame image data collected by the magnetic resonance machine can be 2D data or 3D data, and the variable acceleration factor is changed according to the actual situation. If it is a 2D image, the acceleration factor of the first under-sampled image frame is preferably 2, and the acceleration factor of the second under-sampled image frame is preferably 2 to 4; if it is a 3D image, the acceleration factor of the first under-sampled image frame is preferably 4, and the acceleration factor of the second under-sampled image frame is preferably 4 to 12.

[0045] The highly robust fast imaging method combining k-space and image space reconstruction has the following specific steps:

[0046] S1: For the first undersampled image frame, a complete image frame is reconstructed using a k-space automatic calibration reconstruction method as a calibration frame; a channel image of each channel is obtained after Fourier transforming the k-space data of the calibration frame, and all channel images are merged to obtain a channel merged image.

[0047] In the present invention, the reconstruction method of k-space automatic calibration can adopt the GRAPPA method, which calculates the coil weight and applies it to the under-sampling area to reconstruct a complete image frame. Specifically, GRAPPA belongs to the prior art, which calculates the weight through ACS to linearly combine the data of each channel to reconstruct the complete image frame. Further, if the GRAPPA method is used, the present invention recommends using GRAPPA with Tikhonov regularization.

[0048] However, it should be noted that, in this step, although GRAPPA is recommended to realize the reconstruction of the first under-sampled image frame, other self-calibrated parallel imaging methods in k-space such as SAKE, SPIRiT, LORAKS, etc. may also be adopted according to the situation.

[0049] In addition, in this step, all coil images may be combined by using a square root (RSS) method, and of course, may also be combined by using an adaptive combine method.

[0050] S2: Using the ESPIRiT method on the automatic calibration signal of the first undersampled image frame, obtain the ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map, and identify the first support area that completely covers the target imaging area based on the ESPIRiT eigenvalue map; divide each channel image by the channel merged image to obtain the original coil sensitivity map of each channel; and identify the second support area after removing noise and hole areas from the first support area based on the channel merged image.

[0051] It should be noted that the above-mentioned ESPIRiT coil sensitivity map and the original coil sensitivity map are both in complex form and can be decomposed into corresponding amplitude maps and phase maps. The ESPIRiT coil sensitivity map is composed of the ESPIRiT coil sensitivity amplitude map and the ESPIRiT coil sensitivity phase map, and the original coil sensitivity map is composed of the original coil sensitivity amplitude map and the original coil sensitivity phase map.

[0052] Since the present invention does not consider the case where the imaging object is larger than the FOV, the ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map are both a set of coil sensitivity maps and eigenvalue maps corresponding to the maximum eigenvalue after the eigenvalue decomposition of the ESPIRiT method.

[0053] It should be noted that the first support area is equivalent to the range area of ​​the target sensitivity map to be finally formed, and it should completely cover the target imaging area. For example, for brain magnetic resonance imaging, the target imaging area is the brain area. However, when the first support area covers the target imaging area, a margin should be left at the edge to avoid imaging noise errors at the edge. The so-called margin at the edge means that the edge of the first support area is slightly expanded relative to the edge of the target imaging area, generally 5 to 20 pixels. And according to the recognition process of the two support areas, it can be seen that the second support area should be completely contained in the first support area.

[0054] In addition, in this step, the first support area and the second support area are preferably identified by a threshold value method. The so-called threshold value identification method refers to extracting an area greater than the threshold value as the corresponding support area by setting a threshold value. The specific threshold value used to extract the first support area and the second support area can be adjusted according to actual conditions. In the present invention, the following threshold value identification method can be used to realize the identification and extraction of the first support area and the second support area:

[0055] 1) Identification and extraction of the first support area: The Otsu method (maximum inter-class variance method) is used to determine the Otsu area for the channel merged image, and the proportion of this Otsu area in the entire channel merged image is calculated as A1. The pixel values ​​in the ESPIRiT eigenvalue map are sorted from large to small, and the pixel value with percentile A1 is selected as the first threshold. The area in the ESPIRiT eigenvalue map with pixel values ​​greater than the first threshold is identified as the first support area;

[0056] 2) Identification and extraction of the second support area: Find the maximum pixel value Q from the channel merged image, then calculate the second threshold value as Q*B1 according to the preset ratio value B1, and identify the area in the channel merged image that is greater than the second threshold value as the second support area. The above B can generally be taken as 0.1-0.2, so as to identify the noise and hole areas in the sensitivity map as much as possible.

[0057] S3: Add the sensitivity amplitude map of the ESPIRiT coil that is in the first supporting area but not in the second supporting area (i.e., the part of the ESPIRiT coil sensitivity amplitude map that is in the first supporting area but not in the second supporting area) and the original coil sensitivity amplitude map in the second supporting area (i.e., the part of the original coil sensitivity amplitude map that is located in the second supporting area), and obtain an optimized sensitivity amplitude map after smoothing filtering; interpolate and complete the noise and hole areas based on the original coil sensitivity phase map in the second supporting area (i.e., the part of the original coil sensitivity phase map that is located in the second supporting area) to obtain an optimized sensitivity phase map; combine the optimized sensitivity amplitude map and the optimized sensitivity phase map to finally obtain an optimized sensitivity map.

[0058] It should be noted that in the present invention, the noise and hole areas are interpolated and supplemented based on the original coil sensitivity phase map in the second supporting area. Various different forms of interpolation algorithms can be used, such as a local weighted polynomial regression algorithm, as long as it is possible to use valuable pixels to interpolate and supplement null-value pixels.

[0059] In this step, adding the sensitivity amplitude map of the ESPIRiT coil in the first support region but not in the second support region to the original coil sensitivity amplitude map in the second support region refers to adding at corresponding pixel positions. Combining the optimized sensitivity amplitude map and the optimized sensitivity phase map refers to combining into an optimized sensitivity map (complex number), and fitting can be used for further smoothing.

[0060] S4: Using the optimized sensitivity map, perform SENSE reconstruction with an acceleration factor of 1 on the calibration frame to generate a reference image; at the same time, perform retrospective undersampling on the calibration frame, and the acceleration factor of the retrospective undersampling is the same as the acceleration factor of the second undersampled image frame, and then perform SENSE reconstruction on the data after the retrospective undersampling to generate a retrospective reconstructed image with potential aliasing artifacts; obtain a correction factor map based on the reference image and the retrospective reconstructed image, and the pixel value of each position in the correction factor map is the quotient of the pixel values ​​of the corresponding position in the reference image and the retrospective reconstructed image.

[0061] It should be noted that the correction factor graph obtained by direct division may contain outliers, so it is best to filter it through a filter to remove outliers before using it in the next step.

[0062] It should be noted that in the present invention, there are multiple frames of the second undersampled image frame, and the acceleration factors may be different, so it is necessary to generate a corresponding correction factor map for each acceleration factor used. For example, if the acceleration factor of a frame of the second undersampled image frame is X, then in S4, it is necessary to retrospectively undersample the calibration frame based on the acceleration factor X, and then reconstruct the retrospectively reconstructed image through SENSE, and obtain the correction factor map corresponding to the acceleration factor X based on the reference image and the retrospectively reconstructed image; when the acceleration factor of another frame of the second undersampled image frame is Y, then in S4, it is necessary to retrospectively undersample the calibration frame based on the acceleration factor Y, and then reconstruct the retrospectively reconstructed image through SENSE, and obtain the correction factor map corresponding to the acceleration factor Y based on the reference image and the retrospectively reconstructed image.

[0063] S5: For the remaining second under-sampled image frames except the first under-sampled image frame, perform SENSE reconstruction using the optimized sensitivity map, and multiply the SENSE reconstructed image with the corresponding correction factor map to obtain a complete image frame with suppressed artifacts.

[0064] It should be noted that, in this step, when reconstructing SENSE, it is preferred to use the truncated singular value regularization method to discard some smaller singular values.

[0065] It should be noted that in this step, multiplying the image reconstructed by SENSE with the corresponding correction factor map refers to multiplying two pixel values ​​at corresponding positions as the pixel values ​​at corresponding positions in the complete image frame.

[0066] It should be noted that when the SENSE reconstructed image is multiplied with the corresponding correction factor map in this step, the correction factor map corresponding to the acceleration factor should be used. For example, if the SENSE reconstructed image is reconstructed from undersampled data with an acceleration factor of X, the correction factor map corresponding to the acceleration factor X should be multiplied therewith; if the SENSE reconstructed image is reconstructed from undersampled data with an acceleration factor of Y, the correction factor map corresponding to the acceleration factor Y should be multiplied therewith.

[0067] Therefore, the above S1 to S4 constitute the highly robust fast imaging method (RobustKIPI method) combining k-space and image space reconstruction of the present invention.

[0068] Similarly, based on the same inventive concept, another preferred embodiment of the present invention further provides a magnetic resonance imaging data processing device corresponding to the highly robust fast imaging method combining k-space and image space reconstruction provided in the above embodiment, which includes a memory and a processor;

[0069] The memory is used to store computer programs;

[0070] The processor is used to implement the highly robust fast imaging method combining k-space and image space reconstruction as shown in the above S1 to S5 when executing the computer program.

[0071] It is understandable that the above-mentioned memory may include random access memory (RAM) and non-volatile memory (NVM), such as at least one disk memory. At the same time, the storage medium may also be various media that can store program codes, such as a USB flash drive, a mobile hard disk, a magnetic disk or an optical disk. Of course, with the widespread application of cloud servers, the above-mentioned software programs can also be installed on the cloud platform to provide corresponding services, so the computer-readable storage medium is not limited to the form of local hardware.

[0072] It can be understood that the above-mentioned processor can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components.

[0073] It should also be noted that those skilled in the art can clearly understand that for the convenience and simplicity of description, the specific working process of the device described above can refer to the corresponding process in the aforementioned method embodiment, and will not be repeated here. In each embodiment provided in this application, the division of steps or modules in the device and method is only a logical function division, and there may be other division methods in actual implementation, such as multiple modules or steps can be combined or integrated together, and a module or step can also be split.

[0074] In addition, the logic instructions in the above-mentioned memory can be implemented in the form of software functional units and can be stored in a computer-readable memory when sold or used as an independent product. Based on such an understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product, which is stored in a memory and includes a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention.

[0075] Similarly, based on the same inventive concept, another preferred embodiment of the present invention also provides a magnetic resonance imaging device corresponding to the highly robust fast imaging method combining k-space and image space reconstruction provided in the above embodiment, which includes a magnetic resonance scanner and a control unit. Among them:

[0076] The control unit stores a computer program, which, when executed, is used to implement the highly robust fast imaging method combining k-space and image space reconstruction as shown in S1 to S5 above; the magnetic resonance scanner is used to acquire the first under-sampled image frame data and the second under-sampled image frame data.

[0077] The above-mentioned magnetic resonance scanner can be realized by using existing technology, and can adopt mature commercial products, and the specific model is not limited.

[0078] It should be noted that, in addition to storing the program for implementing the RobustKIPI method, the control unit should also have the imaging sequence and other software programs necessary for implementing multi-frame imaging.

[0079] Of course, the above-mentioned control unit can be an independent control unit or a control unit provided by the magnetic resonance scanner, that is, the above-mentioned highly robust and fast imaging method combining k-space and image space reconstruction can be integrated into the control unit of the magnetic resonance imaging device in the form of a data processing program, so that the magnetic resonance scanner can directly output the reconstruction result without adding an additional control unit.

[0080] The following will further demonstrate the technical effect that can be achieved by the aforementioned highly robust fast imaging method (RobustKIPI method) combining k-space and image space reconstruction of the present invention through an embodiment, so that those skilled in the art can better understand the essence of the present invention.

[0081] Example

[0082] This embodiment applies the highly robust fast imaging method (RobustKIPI method) combining k-space and image space reconstruction shown in S1 to S5 above to a specific example. The specific data source of the RobustKIPI method and the specific implementation method of each step are described in detail below:

[0083] 1. MRI experiment

[0084] The water phantom and human experiments were performed on a 3T Siemens MRI scanner (MAGNETOM Prisma) using a 20-channel head receiving coil. The water phantom used a standard uniform water phantom provided by Siemens. The human experiments were approved by the local institutional review board. The sequence used was the inversion recovery (IR) TSE imaging sequence, and the MRI parameter measured was the longitudinal relaxation time (T1).

[0085] For the human body and water phantom experiments, the sequence parameters were as follows: TR = 3000 ms, TE = 7.5 ms, FOV = 220 × 220 mm 2 , resolution = 1.7 × 1.7 mm 2 , slice thickness = 5 mm, acquisition matrix = 128 × 128, turbine factor = 16. Full sampling data were collected at 7 different inversion times, TI = 50, 150, 300, 500, 800, 1300, 2000 ms.

[0086] 2. Image reconstruction and analysis

[0087] All processing and analysis were performed offline using MATLAB (MathWorks, Natick, MA) software written on a PC computer (3.2 GHz).

[0088] For the water model, the seventh frame (TI = 2000ms) was selected as the first undersampled image frame, AF = 2. The other six frames were the second undersampled frames, AF = 4. For the first undersampled image frame, the middle 24 phase encoding lines were retained as ACS, and the second undersampled image frame had no ACS. The full sampled data was retrospectively undersampled in the above manner to obtain data with variable acceleration factors. For human experiments, the above retrospective undersampling method was also adopted.

[0089] The RobustKIPI of this embodiment is divided into the following steps:

[0090] Step 1: Perform GRAPPA reconstruction on the first undersampled image frame and fill in ACS to obtain a calibration frame. Perform Fourier transform on the calibration frame to obtain the channel image of each channel; use the square root (RSS) method to merge all channel images to obtain the merged RSS image. The GRAPPA kernel size is 2×5, that is, 2 acquired phase encoding lines and 5 frequency encoding points, and each point is fitted using 10×20 (number of channels) points.

[0091] The second step is to calculate ESPIRiT, the original sensitivity map, and the first and second support regions. The ESPIRiT method is used on the automatic calibration signal of the first undersampled image frame to obtain the ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map. The kernel size of ESPIRiT is 6×6; the subspace threshold is 0.02, that is, the right singular vectors corresponding to the singular values ​​greater than 2% of the maximum singular value are listed as the basis of the subspace; the first support region is confirmed by implementing a threshold recognition operation on the ESPIRiT eigenvalue map, and the threshold in this case is 0.99. The ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map are both a set of coil sensitivity maps and eigenvalue maps corresponding to the maximum eigenvalue after the eigenvalue decomposition of the ESPIRiT method, and multiple sets are not required.

[0092] Divide each channel image by the RSS image to obtain the original coil sensitivity map of each channel. Perform a threshold operation on the RSS image to identify the second support region. In this example, the threshold is 0.1. The second support region is contained in the first support region.

[0093] The third step is to calculate the optimized sensitivity map. The sensitivity amplitude map of the ESPIRiT coil that is in the first support area but not in the second support area is added to the sensitivity amplitude map of the original coil in the second support area to obtain an optimized sensitivity amplitude map. Since the phase of the ESPIRiT coil sensitivity is significantly different from the original phase, the original coil sensitivity phase map in the second support area is extrapolated and interpolated to obtain an optimized sensitivity phase map; after combining the optimized sensitivity amplitude map and the optimized sensitivity phase map into a complex sensitivity map, fitting is performed to obtain the final optimized sensitivity map. In this example, both extrapolation, interpolation and fitting are performed using a local weighted polynomial regression (LWPR) operation with a cubic weighted kernel. The polynomial used in this embodiment is a second-order polynomial with a window width of 12.

[0094] The fourth step is to calculate the correction factor map for artifact suppression. Using the optimized sensitivity map, perform SENSE reconstruction with an acceleration factor of 1 on the calibration frame to generate a reference image ρ0; at the same time, perform retrospective undersampling on the calibration frame (the acceleration factor of retrospective undersampling is 2 or 4), and then perform SENSE reconstruction on the undersampled data to generate a retrospective reconstructed image ρ1 with potential aliasing artifacts; obtain the correction factor map based on the reference image and the retrospective reconstructed image, and the pixel value at each position in the correction factor map is the quotient of the pixel value at the corresponding position in the reference image and the retrospective reconstructed image. In addition, the correction factor map in this embodiment needs to be filtered by a median filter with a window of 3×3 to remove outliers. The correction factor map C can be expressed as follows:

[0095] C=ρ0 / ρ1

[0096] In this example, since the second under-sampled image frame has two acceleration factors (AF=2 or 4), this process is performed twice to obtain a correction factor map corresponding to AF=2 and a correction factor map corresponding to AF=4.

[0097] Step 5: Reconstruct all second undersampled image frames. Use the optimized sensitivity map to perform SENSE reconstruction on the second undersampled image frame to obtain a reconstructed image with potential artifacts, and use truncated singular value regularization to discard singular values ​​less than 2% of the maximum value. Then perform artifact suppression, multiply the second undersampled image frame with the corresponding correction factor map, and obtain the final reconstructed artifact-suppressed image.

[0098] The above steps 1 to 5 can be found in Figure 6 The algorithm flow chart is shown.

[0099] For comparison, GRAPPA, ESPIRiT, and KIPI reconstructions were performed on the retrospective undersampled data. The GRAPPA and ESPIRiT methods used the ACS calculation weights or ESPIRiT sensitivity maps obtained from the first undersampled image frame to reconstruct all image frames. The KIPI method used two sensitivity maps, one was the original coil sensitivity map obtained by direct division, and the other was the ESPIRiT coil sensitivity map. The RubustKIPI proposed in the present invention used the above-mentioned automatically optimized coil sensitivity map.

[0100] T1 parameter image is fitted by We can get , where y is the reconstructed pixel value, u and v are model adaptation variables. Given 7 sets of data, solve for 3 unknowns including T1.

[0101] The accuracy of the RubustKIPI method was evaluated by comparing the normalized root mean square error (nRMSE) of GRAPPA, ESPIRiT, KIPI, and RubustKIPI with the true value of the full sampling.

[0102] 3. Results Analysis

[0103] Figure 1 The flowchart of RubustKIPI combining the original coil sensitivity amplitude map and the ESPIRiT coil sensitivity amplitude map is shown. The first support region A is the sensitivity range to be formed in the end, which should completely cover the range of the brain area to be imaged and have a slight margin for edge expansion; the second support region B is the area after excluding the noise and hole areas of the original sensitivity map in region A. After comparison, it can be seen that the first support region is larger than the second support region. The area in the first support region but not in the second support region can be regarded as the noise and hole area C. RubustKIPI can retain the accurate original sensitivity map and complete the sensitivity map in the noise and hole area C with the ESPIRiT sensitivity map.

[0104] Figure 2 The amplitude and phase of three different coil sensitivity maps for channel 7 are shown. They are the original coil sensitivity map, the ESPIRiT coil sensitivity map, and the RubustKIPI optimized coil sensitivity map. It can be seen that the original coil sensitivity map has noisy edges and non-smooth jumps (arrows), while the ESPIRiT sensitivity map is very smooth but a bit blurry. The sensitivity map optimized by RubustKIPI retains the details and avoids the interference of noise and holes. Its phase map is highly consistent with the original phase map.

[0105] Figure 3The reconstructed water phantom source image and the corresponding T1 parametric image obtained from the undersampled IR TSE data are shown. In this case, due to the absence of void regions and the high signal-to-noise ratio, the water phantom itself is relatively homogeneous. Both KIPI and RubustKIPI achieve high-quality results, while GRAPPA and ESPiRiT methods produce obvious artifacts in both the source image and the T1 parametric image (arrows).

[0106] Figure 4 IR TSE source images reconstructed for TI=50ms for GRAPPA, ESPIRiT, KIPI and RubustKIPI. KIPI uses the original (Raw) sensitivity map and the ESPIRiT sensitivity map. The results of GRAPPA and ESPIRiT have obvious unfolding artifacts, indicated by the arrows in the second and third columns, which are not present in the results of KIPI. KIPI using the ESPIRiT sensitivity map also shows some artifacts, while RubustKIPI has the smallest error (nRMSE=0.021).

[0107] Figure 5 Shows from Figure 4 The T1 parameter image calculated from the source image. Obvious unfolding artifacts can be seen in the reconstruction results of the GRAPPA and ESPIRiT methods; the KIPI method using the original (Raw) sensitivity map shows slight artifacts due to the holes and noise areas; and the KIPI method using the ESPIRiT sensitivity map also shows certain artifacts because the ESPIRiT sensitivity map itself lacks some details. However, there are no artifacts in the RubustKIPI results, the errors in gray matter and white matter are very small, and the total error is also the smallest (nRMSE = 0.0686). The T1 parameter image generated by the RubustKIPI method is most consistent with the true value.

[0108] The above-described embodiment is only a preferred solution of the present invention, but it is not intended to limit the present invention. A person skilled in the relevant technical field may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, any technical solution obtained by equivalent replacement or equivalent transformation falls within the protection scope of the present invention.

Claims

1. A highly robust fast imaging method combining k-space and image space reconstruction for reconstructing undersampled image frames in multi-frame imaging, characterized in that: The under-sampled image frames include a first under-sampled image frame with an automatic calibration signal and a second under-sampled image frame with an acceleration factor not less than 2, and the acceleration factor of the first under-sampled image frame is not higher than the acceleration factor of the second under-sampled image frame; The imaging method steps are as follows: S1: for the first undersampled image frame, a complete image frame is reconstructed using a k-space automatic calibration reconstruction method as a calibration frame; the k-space data of the calibration frame is Fourier transformed to obtain a channel image of each channel, and all channel images are merged to obtain a channel merged image; S2: using the ESPIRiT method on the automatic calibration signal of the first under-sampled image frame to obtain an ESPIRiT coil sensitivity map and an ESPIRiT eigenvalue map, and identifying a first support area that completely covers the target imaging area based on the ESPIRiT eigenvalue map; Dividing each channel image by the channel merged image to obtain an original coil sensitivity map of each channel; identifying a second support region after removing noise and hole regions from the first support region based on the channel merged image; The ESPIRiT coil sensitivity map and the original coil sensitivity map are both composed of corresponding amplitude maps and phase maps; S3: adding the sensitivity amplitude map of the ESPIRiT coil in the first support area but not in the second support area to the original coil sensitivity amplitude map in the second support area, and obtaining an optimized sensitivity amplitude map after smoothing filtering; interpolating and completing the noise and hole areas based on the original coil sensitivity phase map in the second support area to obtain an optimized sensitivity phase map; combining the optimized sensitivity amplitude map and the optimized sensitivity phase map to finally obtain an optimized sensitivity map; S4: using the optimized sensitivity map, performing SENSE reconstruction with an acceleration factor of 1 on the calibration frame to generate a reference image; at the same time, retrospectively undersampling the calibration frame, and the acceleration factor of the retrospective undersampling is the same as the acceleration factor of the second undersampled image frame, and then performing SENSE reconstruction on the retrospectively undersampled data to generate a retrospective reconstructed image with potential aliasing artifacts; obtaining a correction factor map based on the reference image and the retrospectively reconstructed image, and the pixel value of each position in the correction factor map is the quotient of the pixel value of the corresponding position in the reference image and the retrospectively reconstructed image; S5: For the remaining second under-sampled image frames except the first under-sampled image frame, perform SENSE reconstruction using the optimized sensitivity map, and multiply the SENSE reconstructed image with the correction factor map to obtain a complete image frame with suppressed artifacts.

2. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: The first under-sampled image frame and the second under-sampled image frame are both two-dimensional images; the acceleration factor of the first under-sampled image frame is 2, and the acceleration factor of the second under-sampled image frame is 2 to 4.

3. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: The first under-sampled image frame and the second under-sampled image frame are both three-dimensional images; the acceleration factor of the first under-sampled image frame is 4, and the acceleration factor of the second under-sampled image frame is 4 to 12.

4. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: In S1, the k-space automatic calibration reconstruction method is the GRAPPA method, which calculates coil weights and applies them to under-sampling areas to reconstruct a complete image frame; the GRAPPA method is GRAPPA with Tikhonov regularization.

5. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: In S1, all coil images are merged by a square root method or an adaptive method.

6. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: In S2, the ESPIRiT coil sensitivity map and the ESPIRiT eigenvalue map are both a set of coil sensitivity maps and eigenvalue maps corresponding to the maximum eigenvalue after eigenvalue decomposition of the ESPIRiT method.

7. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: In S4, the correction factor map needs to be filtered to remove abnormal values.

8. The highly robust fast imaging method combining k-space and image space reconstruction as claimed in claim 1, characterized in that: In S5, when performing SENSE reconstruction using the optimized sensitivity map, a truncated singular value regularization method is adopted.

9. A magnetic resonance imaging data processing device, characterized in that: including memory and processor; The memory is used to store computer programs; The processor is used to implement the highly robust fast imaging method combining k-space and image space reconstruction as described in any one of claims 1 to 8 when executing the computer program.

10. A magnetic resonance imaging device, characterized in that: The invention comprises a magnetic resonance scanner and a control unit, wherein a computer program is stored in the control unit, and when the computer program is executed, it is used to implement the highly robust fast imaging method combining k-space and image space reconstruction as described in any one of claims 1 to 8; the magnetic resonance scanner is used to acquire the first under-sampled image frame data and the second under-sampled image frame data.

Images