Joint k-space and image-space reconstruction imaging method and device

Abstract

To provide a joint k-space and image-space reconstruction imaging method and a device, relating to the field of magnetic resonance imaging.SOLUTION: By means of a k-space parallel imaging method, under-sampled image frames having auto-calibration signals, obtained from multi-frame imaging required for the measurement of magnetic resonance parameters of a measurement object, are reconstructed, and then other under-sampled image frame data without auto-calibration signals are reconstructed in image space. The present invention generates an accurate sensitivity map required for the SENSE method from under-sampled data without the need of additional acquisition. The method generates images almost identical to those obtained by the traditional GRAPPA method and without obvious artifacts, but at a speed of four times or more than that of the traditional GRAPPA method. The method does not require all sampling frames as in the reconstruction method of auto-calibration, and is particularly suitable for accelerating three-dimensional (3D) multi-frame magnetic resonance imaging.SELECTED DRAWING: Figure 8

Description

[Technical Field]

[0001] This application claims priority from a Chinese patent application having a filing date of April 8, 2021 and application number 202110379060.9, the entire contents of which are incorporated herein by reference.

[0002] The present invention is in the field of magnetic resonance imaging and can provide high speed acquisition for measurements of magnetic resonance parameters that require multi-frame image data. [Background technology]

[0003] Magnetic resonance imaging (MRI) may be used to measure tissue parameters including nuclear spin density, longitudinal relaxation time (T1), transverse relaxation time (T2), apparent diffusion coefficient (ADC), cerebral blood flow velocity, magnetization transition velocity, etc. The signal sources of these parameters are the nuclei of endogenous compounds in vivo (e.g., in water). 1 H), and may be exogenous substances or tracers. Direct quantization of tissue MRI parameters can provide a variety of pathological and clinical information of additional value in neurological, musculoskeletal, hepatic, and myocardial imaging.

[0004] Imaging techniques that measure these parameters include relaxation time imaging, diffusion imaging, perfusion imaging, functional magnetic resonance imaging (fMRI), and chemical exchange saturation transfer (CEST) imaging, which generally require the acquisition of multiframe image datasets with pulse sequence parameter modulation for the same imaging object. These modulatable pulse sequence parameters include echo time (TE), flip angle, or gradient strength. For example, diffusion imaging requires the acquisition of data at different gradient fields for the same imaging object. Estimation of target parameters can be achieved by fitting signal changes in the multiframe image datasets to a mathematical model. The need for multiple acquisitions results in long scan times, which limits its clinical application. Furthermore, the technique itself is relatively sensitive to motion, which is likely to cause additional motion due to longer scan times. Therefore, shortening data acquisition times and accelerating the imaging process not only provides convenience to patients but also helps improve image quality.

[0005] Researchers both in China and abroad have proposed many methods to accelerate parameter measurement in multi-frame imaging. Examples of CEST techniques include Keyhole-CEST, SLAM, kw ROSA, and CS-CEST. Keyhole-CEST combines acquired low-resolution images with high-resolution reference images to reduce the data collection time, but at the expense of high-frequency image information. SLAM methods generate regional CEST measurements directly from arbitrarily shaped regions of interest using a priori positioning knowledge acquired from reconnaissance scans. This method has a much higher signal-to-noise ratio than traditional single-voxel and multi-voxel methods, but it loses information about tissue heterogeneity within the region. Among model-based methods, kw ROSA decomposes the Z-spectrum signal based on a subspace and combines it with a measurement model to directly estimate the asymmetric Z-spectrum of interest from complete or incomplete measurements, but its acceleration capabilities are limited. CS-CEST methods based on compressed sensing exploit the sparsity in the transform domain of images to accelerate the image, but have yet to achieve full random undersampling in practical applications.

[0006] Parallel imaging methods are widely used clinically due to their practicality and robustness. Currently, they are basically divided into two types: the self-correcting k-space method represented by GRAPPA, and the image-space method based on explicit coil sensitivity represented by SENSE. GRAPPA is a self-correcting k-space parallel imaging algorithm that fits missing points by linear combination of multi-channel data. The correction equation for GRAPPA is shown in Equation 1 below:

number

[0007] SENSE treats parallel reconstruction as an inversion problem of a linear equation in image space. If the sensitivity map is perfectly accurate, SENSE can find the optimal solution in the least-squares sense. However, in practice, accurate sensitivity maps are difficult to obtain, and small errors in the sensitivity map can cause severe wraparound artifacts in the reconstructed image. Recently, a variable acceleration sensitivity coding (vSENSE) method has been proposed to accelerate the acquisition of multiframe image data and parameter measurement. In multiframe imaging, images from different frames are acquired from the same imaging object and theoretically have the same coil sensitivity. Based on this assumption, the vSENSE method obtains an adjusted sensitivity map from one image frame by full sampling or undersampling with a lower factor, and then reconstructs another image frame by undersampling with a higher acceleration factor. Although vSENSE has higher acceleration capabilities, it has three obvious drawbacks. First, obtaining the initial sensitivity map from an additional reference scan rather than the automatic correction method may cause inconsistencies. Second, in the latest 3D vSENSE method, since a full sampling image cannot be obtained, the SENSE image frame accelerated by 2 times is considered accurate, and the reconstructed image still has potential artifacts. Third, the highest acceleration factor in the 3D vSENSE method is 8, and more artifacts already appear, so the acceleration effect is still limited.

[0008] Therefore, in the field of MR parameter measurement requiring multi-frame imaging, it is extremely important to provide a method that allows automatic corrective reconstruction and has high robustness by allowing for higher acceleration factors. Summary of the Invention [Problem to be solved by the invention]

[0009] The objective of the present invention is to solve the problems of the prior art in which automatic correction is impossible in multi-frame imaging reconstruction, and the reconstructed images have potential artifacts, resulting in a greatly limited acceleration effect, and to provide a high-speed imaging method and apparatus (known as joint K-space and Image-space Parallel Imaging, hereinafter also abbreviated as KIPI method) that combines k-space and image-space for reconstruction, thereby recovering accurate images from highly undersampled k-space data. [Means for solving the problem]

[0010] The specific technical solutions used in the present invention are as follows:

[0011] In a first aspect, there is provided an imaging method for combined k-space and image space reconstruction according to the present invention, for reconstructing undersampled image frames, the undersampled image frames including a first undersampled image frame having an auto-correction signal and a second undersampled image frame having no auto-correction signal and an acceleration factor of 2 or more, wherein the acceleration factor of the first undersampled image frame is less than or equal to the acceleration factor of the second undersampled image frame; The reconstruction method includes: Step S1: reconstructing the first undersampled image frame by a parallel imaging method with automatic correction in k-space to obtain a complete image frame as a correction frame, Fourier transforming the k-space data of the correction frame to obtain a coil image of each channel, and combining the coil images of all channels to obtain a channel combination image; Step S2: Dividing the coil image of each channel by the channel-combined image to obtain a coil sensitivity map of each channel, identifying a region of support from the channel-combined image, and performing smoothing to remove noise from the region of support in the sensitivity map, and extrapolating the non-support region in the sensitivity map to obtain an optimized sensitivity map; S3: performing SENSE reconstruction with an acceleration factor of 1 on the correction frame using the optimized sensitivity map to generate a reference image; performing backward undersampling on the correction frame to make the acceleration factor of the backward undersampling the same as that of the second undersampled image frame; and performing SENSE reconstruction on the undersampled data to generate a backward reconstructed image with potential aliasing artifacts; obtaining a correction factor map based on the reference image and the backward reconstructed image; and calculating the pixel value of each position in the correction factor map as a quotient of the pixel values of the corresponding positions in the reference image and the backward reconstructed image; and step S4 of performing SENSE reconstruction on the remaining second undersampled image frame other than the first undersampled image frame using the optimized sensitivity map, and multiplying the SENSE reconstructed image by the correction coefficient map to obtain a complete image frame that suppresses artifacts.

[0012] In a preferred embodiment of the first aspect, in step S1, the self-correcting parallel imaging method in k-space is the GRAPPA method, which calculates and applies coil weights to the undersampled region to reconstruct and acquire a complete image frame.

[0013] Furthermore, the GRAPPA method is preferably GRAPPA with Tikhonov regularization.

[0014] In a preferred example of the first aspect, the first undersampled image frame and the second undersampled image frame are both two-dimensional images, and the acceleration factor of the first undersampled image frame is preferably 2, and the acceleration factor of the second undersampled image frame is preferably 2 to 4.

[0015] In a preferred example of the first aspect, the first undersampled image frame and the second undersampled image frame are both three-dimensional images, the acceleration factor of the first undersampled image frame is preferably 2×2, and the acceleration factor of the second undersampled image frame is preferably N×M, where 2≦N≦4, 2≦M≦4, and the total acceleration factor is 12 or less.

[0016] In a preferred example of the first aspect, in step S1, all coil images are combined using the square root method or the self-adaptive method.

[0017] As a preferred example of the first aspect, in step S2, noise removal is performed by smoothing the support region in the sensitivity map using a fitting method.

[0018] As a preferred example of the first aspect, in step S3, the correction coefficient map needs to be filtered by a filter so as to remove abnormal values.

[0019] As a preferred example of the first aspect, in step S4, it is preferable to use a cutoff singular value regularization method when performing SENSE reconstruction using the optimized sensitivity map.

[0020] In a second aspect, there is provided a magnetic resonance imaging apparatus according to the present invention, comprising a magnetic resonance scanner and a control unit, the control unit storing a computer program which, when executed, is adapted to realize the method for reconstructing a multi-frame image according to any one of the first aspects, and the magnetic resonance scanner for acquiring the first undersampled image frame data and the second undersampled image frame data. [Effects of the Invention]

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

[0022] The present invention does not require ACS data for every frame, significantly improving the effect of pure acceleration. The present invention can generate accurate sensitivity maps using undersampled image frames with an acceleration factor of 2 or more. Because the first undersampled image frame is first reconstructed using a self-correcting parallel imaging method in k-space (e.g., GRAPPA), the method can achieve self-correction without requiring additionally acquired sensitivity maps. Because the present invention takes into account the robustness of the self-correcting parallel imaging method in k-space, it can further accelerate the acquisition speed compared to the conventional vSENSE method.

[0023] Since it is not necessary to use all sampled image frames, this method is particularly applicable to 3D multi-frame imaging. When acquiring source images required for parameter measurement, the present invention allows the forward acceleration factor in the phase coding and hierarchical coding directions to be increased by up to 12 times, and produces source images and final target parameter images that are consistent with the truth results. Meanwhile, the present invention ensures that images reconstructed from highly undersampled data have no obvious aliasing artifacts. [Brief explanation of the drawings]

[0024] [Figure 1]Figure 1 shows the reconstruction results of an example, comparing +6 ppm phantom images reconstructed by GRAPPA and KIPI. Figure 1a shows the root-sum-square (RSS) reconstruction result of the full sampled k-space. Figure 1b shows the phantom source image reconstructed by GRAPPA when the acceleration factor (AF) is 4. Figure 1c shows the phantom source image reconstructed by KIPI when AF is 4. Figure 1d shows the Z spectrum and full k-space results (solid lines) acquired by GRAPPA (dashed line) and KIPI (double line) in the region of interest of the selected image a. Figures 1e and 1f show the difference maps of images b, c, and a, respectively. Images b and c use the same variable acceleration factor undersampling data. The first +3.5 ppm is selected as the first undersampling image frame with AF=2, S0 and 3.5 ppm (excluding the first undersampling image frame) are selected as AF=2, and the remaining 46 frames are set to AF=4. [Figure 2] Figure 2 shows amide proton transfer weight (APTw) parameter images of an example water model. Figure 2a shows an APTw image calculated from all sampling data. Figure 2b shows an APTw image obtained by reconstructing a second undersampled image frame using weights obtained from the ACS data of a first undersampled image frame using GRAPPA, which shows obvious artifacts. Figure 2c shows an APTw image obtained by reconstructing a second undersampled image frame using a sensitivity map derived from the correction frame using KIPI. Figures 2b and 2c use the same variable acceleration factor undersampling data, with the first +3.5 ppm image frame selected as the first undersampled image frame with AF=2, S0 and ±3.5 ppm (excluding the first undersampled image frame) selected as AF=2, and the remaining 46 frames selected as AF=4. [Figure 3]Figure 3 shows +6 ppm brain images reconstructed by GRAPPA and KIPI in an embodiment. Figure 3a shows the RSS reconstruction of the full sampling k-space, Figure 3b shows the source image of a healthy volunteer reconstructed by GRAPPA when AF=4, Figure 3c shows the source image of a healthy volunteer reconstructed by KIPI when AF=4, Figure 3d shows the Z spectrum and full k-space results (solid lines) acquired by GRAPPA (dashed line) and KIPI (double line) in the region of interest of the selected image a, and Figures 3e and 3f show the difference maps of images b, c, and a, respectively. Images b and c of Figure 3 use the same variable acceleration factor undersampling data, with the first +3.5 ppm selected as the first undersampling image frame with AF=2, S0 and ±3.5 ppm (excluding the first undersampling image frame) selected as AF=2, and the remaining 46 frames selected as AF=4. [Figure 4] Figure 4 shows 2D APTw parameter images of a healthy volunteer in accordance with an embodiment of the present invention. Figure 4a shows an APTw image calculated from all sampling data. Figure 4b shows an APTw image reconstructed by GRAPPA when AF=4, using weights obtained from the corrected frame ACS data to reconstruct the remaining frames. Figure 4c shows an APTw image reconstructed by KIPI when AF=4, using sensitivity maps derived from the corrected frames in the SENSE step. Images 4b and 4c use the same variable acceleration factor undersampling data, with the first +3.5 ppm image frame selected as the first undersampling frame with AF=2, S0 and ±3.5 ppm (excluding the first undersampling image frame) selected as AF=2, and the remaining 46 frames selected as AF=4. [Figure 5]Figure 5 shows -4 ppm 3D-CEST source images of a healthy volunteer (each set contains five maps selected from 72 maps) in an example. Figure 5a shows a -4 ppm source image obtained by reconstruction using conventional GRAPPA with AF=2×1. Figure 5b shows a source image of the healthy volunteer reconstructed using GRAPPA from undersampled data with a variable acceleration factor. Figure 5c shows a source image of the healthy volunteer reconstructed using KIPI from undersampled data with a variable acceleration factor. Figure 5d and Figure 5e show difference maps of the reconstruction results using GRAPPA and KIPI with respect to the true value a, respectively. The reconstruction errors (RNMSE) are 0.026 and 0.014, respectively. The arrows indicate aliasing artifacts reconstructed using GRAPPA. For the variable acceleration undersampling data, the +3.5 ppm frame has ACS data, and AF=2×2 is the first undersampling image frame. The remaining 6 frames have no ACS data, and AF=2×4 is the second undersampling image frame. [Figure 6] Figure 6 shows APTw parameter images of a healthy volunteer in an example (each set contains five maps selected from 72 maps). Figure 6a shows an APTw image obtained by reconstructing all frames using GRAPPA with AF=2×1 and ACS data. Figure 6b shows an APTw image of a healthy volunteer reconstructed using GRAPPA from variable acceleration factor undersampled data. Figure 6c shows an APTw image of a healthy volunteer reconstructed using KIPI from variable acceleration factor undersampled data. The arrows indicate artifacts in the variable acceleration factor GRAPPA APTw image. The +3.5 ppm frame contains ACS data, with AF=2×2 representing the first undersampled image frame. The remaining six frames do not contain ACS data, with AF=2×4 representing the second undersampled image frame. [Figure 7]Figure 7 shows APTw images of a healthy volunteer in an example (each set contains five maps selected from 72 maps). Figure 7a shows an APTw image obtained by reconstructing all frames using GRAPPA with AF=2×1 and ACS data. Figure 7b shows an APTw image of a healthy volunteer reconstructed using GRAPPA from variable acceleration factor undersampled data. Figure 7c shows an APTw image of a healthy volunteer reconstructed using KIPI from variable acceleration factor undersampled data. The arrows indicate artifacts in the variable acceleration factor GRAPPA APTw image. The +3.5 ppm frame contains ACS data, with AF=2×2 representing the first undersampled image frame. The remaining six frames do not contain ACS data, with AF=4×3 representing the second undersampled image frame. [Figure 8] FIG. 8 is a schematic flow chart of the KIPI method of the present invention, in which the reconstructed first undersampled image frame is the corrected frame. DETAILED DESCRIPTION OF THE INVENTION

[0025] The present invention will be further described in detail and illustrated below by specific embodiments with reference to the drawings. Unless conflicting with each other, the technical features of each embodiment of the present invention may be combined correspondingly.

[0026] In the SENSE reconstruction method, the i (1≦i≦N)th coil image m i is known, the sensitivity map SE of the i-th coil i is m i can be obtained by dividing by the coil combination image ρ.

number

[0027] Therefore, for backward SENSE reconstruction with an acceleration factor (AF) of 2 (without loss of generality), the following equation holds:

number

[0028] The least squares solution to equation [3] is shown below in equation 4:

number

[0029] In practical situations, the channel images are generally unknown. Note that the GRAPPA reconstruction results are i Substituting as above, equations [2-4] indeed explain that applying sensitivity maps derived from GRAPPA reconstruction results to SENSE can produce results that are consistent with GRAPPA reconstruction, regardless of the acceleration factor of SENSE.i Whether or not is completely accurate, both equations [3] and [4] hold if the sensitivities used meet the definition in equation [2]. Because different image frames in multi-frame imaging techniques have the same sensitivities, the sensitivity map derived from a GRAPPA reconstructed frame can be applied to other frames.

[0030] Based on the above principle, during sampling, the present invention selects one image frame in multi-frame imaging and undersamples it with a lower acceleration factor (e.g., AF=2 for 2D imaging and AF=2×2 for 3D imaging) to preserve the central autocorrection signal (ACS) and records it as the first undersampled image frame. The remaining frames are undersampled with a higher acceleration factor to eliminate the ACS and record it as the second undersampled image frame. The present invention performs GRAPPA reconstruction on the first undersampled image frame with the lower acceleration factor to obtain a corrected frame. Then, it calculates the corresponding sensitivity map according to Equation [2], and then applies this sensitivity map to perform SENSE reconstruction on the second undersampled image frame. In a general sense, the accuracy of SENSE reconstruction is related only to the sensitivity map used, not to the specific image contrast. The above derivation proves that applying SENSE reconstruction to the corrected frame can obtain the same image quality as GRAPPA with a low acceleration factor. Since the corrected frame and another second undersampled image frame have the same sensitivity (ideally accurate but not obtainable) and the imaging object is the same, if SENSE reconstruction is performed on another frame using the above derived sensitivity map (actual but obtainable), an image with quality close to that of the corrected frame can also be obtained. The process of the imaging method for reconstructing k-space and image space in combination according to the present invention is shown in Figure 8, and the specific implementation process will be described in detail below.

[0031] The multi-frame imaging reconstruction method is primarily for reconstructing undersampled image frames acquired by a magnetic resonance scanner, where the undersampled image frames should include a first undersampled image frame with an ACS and a second undersampled image frame without an ACS. The acceleration factor AF of the second undersampled image frame should be greater than or equal to 2, and the acceleration factor of the first undersampled image frame should be less than or equal to the acceleration factor of the second undersampled image frame, thereby shortening the time required to acquire the second undersampled image frame. While only one first undersampled image frame is required, multiple second undersampled image frames are required. Therefore, if a correction frame can be acquired using the first undersampled image frame and then a second undersampled image frame can be reconstructed for each remaining frame, the acquisition speed can be significantly improved.

[0032] The undersampled image frames acquired by the magnetic resonance scanner may be 2D images or 3D images. The specific acceleration factors of the undersampled image frames need to be adjusted according to actual conditions. If the undersampled image frames are 2D images, the acceleration factor of the first undersampled image frame is preferably 2, and the acceleration factor of the second undersampled image frame is preferably 2 to 4. If the undersampled image frames are 3D images, the acceleration factor of the first undersampled image frame is preferably 2×2, and the acceleration factor of the second undersampled image frame is preferably N×M, where 2≦N≦4, 2≦M≦4, and the total acceleration factor is 12 or less.

[0033] The imaging method for combining and reconstructing the k-space and the image space specifically includes the following steps:

[0034] In step S1, the first undersampled image frame is reconstructed by the automatic correction parallel imaging method in k-space (GRAPPA) to obtain a complete image frame as a correction frame. Then, the k-space data of the correction frame is Fourier transformed to obtain the coil image of each channel, and the coil images of all channels are combined to obtain a channel-combined image.

[0035] The GRAPPA method belongs to the prior art, which calculates and applies coil weights to undersampled regions to reconstruct and obtain a complete image frame. Furthermore, in the case of the GRAPPA method, the present invention recommends adopting GRAPPA with Tikhonov regularization.

[0036] In this step, we recommend that the reconstruction of the first undersampled image frame be achieved by GRAPPA, but other parallel imaging methods with automatic correction in k-space, such as SPIRiT or CAIPIRINHA, may also be used in some cases.

[0037] Also, in this step, the combination of all coil images can be realized by the square root (RSS) method, and of course, can also be realized by the adaptive combine method.

[0038] In step S2, the coil image of each channel is divided by the channel-combined image to obtain a coil sensitivity map of each channel, and the support regions in the sensitivity map are identified from the channel-combined image, and noise is removed by smoothing the support regions in the sensitivity map, and the non-support regions in the sensitivity map are extrapolated to obtain an optimized sensitivity map.

[0039] In this step, noise removal by smoothing can be performed on the support region in the sensitivity map using a fitting method. The fitting method includes two methods: local fitting and global fitting. Of course, noise removal by smoothing the support region can also be performed using a filtering method.

[0040] In step S3, SENSE reconstruction with an acceleration factor of 1 is performed on the corrected frame using the optimized sensitivity map to generate a reference image. At the same time, backward undersampling is performed on the corrected frame to make the acceleration factor of the backward undersampling the same as that of the second undersampled image frame, and SENSE reconstruction is then performed on the undersampled data to generate a backward reconstructed image with potential aliasing artifacts. Next, a correction coefficient map is obtained based on the reference image and the backward reconstructed image, and the pixel value at each position in the correction coefficient map is calculated as the quotient of the pixel values at corresponding positions in the reference image and the backward reconstructed image, i.e., the quotient of the reference image and the backward reconstructed image is calculated for each pixel to obtain the correction coefficient map.

[0041] It should be noted that the directly obtained correction coefficient map may contain outliers, so it is desirable to filter out the outliers before performing the next step.

[0042] In step S4, SENSE reconstruction can be performed on the second undersampled image frame of each remaining frame other than the first undersampled image frame using the optimized sensitivity map, and by multiplying the image after SENSE reconstruction by the correction coefficient map, a complete image frame with suppressed artifacts can be obtained.

[0043] In this step, when the image after SENSE reconstruction is multiplied by the correction coefficient map, multiplying two images means multiplying the pixel values at the same position in the two images point by point to obtain the pixel value at the corresponding position in the complete image frame.

[0044] It should be noted that the second undersampled image frame of the present invention includes multiple frames, and different second undersampled image frames may have different acceleration factors. When performing SENSE reconstruction on each second undersampled image frame, the correction factor map used must also be obtained based on the same acceleration factor. Specifically, if the acceleration factor of one second undersampled image frame is X, then in step S3, the correction frame must be backward undersampled based on the acceleration factor X, and then a backward reconstruction image must be reconstructed using SENSE, and the correction factor map must be obtained based on the reference image and the backward reconstruction image. If the acceleration factor of another second undersampled image frame is Y, then in step S3, the correction frame must be backward undersampled based on the acceleration factor Y, and then a backward reconstruction image must be reconstructed using SENSE, and the correction factor map must be obtained based on the reference image and the backward reconstruction image.

[0045] In this step, when performing SENSE reconstruction with the optimized sensitivity maps, it is preferable to use a cutoff singular value regularization method to discard some of the smaller singular values.

[0046] As a result, the above steps S1 to S4 constitute an imaging method (KIPI) for combining k-space and image space for reconstruction according to the present invention. In practical application, the above KIPI method can be integrated into the control unit of a magnetic resonance imaging apparatus to form a magnetic resonance imaging apparatus capable of automatically performing multi-frame image reconstruction. The magnetic resonance imaging apparatus includes a magnetic resonance scanner and a control unit, and the control unit stores a computer program that can realize the above KIPI method when the computer program is executed. Then, undersampling image data (first undersampling image frame and second undersampling image frame) required for the KIPI method are acquired in advance by the magnetic resonance scanner.

[0047] The above magnetic resonance scanner may be realized by conventional technology and belongs to mature commercial products, so a detailed description is omitted.

[0048] The control unit may be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), or the like, or may be a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic device, a discrete gate or transistor logic device, or a discrete hardware component.

[0049] In addition to storing the program for implementing the KIPI method, the control unit should also have imaging sequences and other software programs required for implementing multi-frame imaging.

[0050] The imaging method and apparatus for combined k-space and image space reconstruction according to the present invention may be applied to any multi-frame imaging technique that requires multiple acquisitions of the same imaging object in magnetic resonance images, including, but not limited to, relaxation time (T1, T2) imaging, diffusion imaging, perfusion imaging, functional imaging (fMRI), magnetic resonance spectroscopy (MRS) imaging, and chemical exchange saturation transfer (CEST) imaging.

[0051] In order to allow those skilled in the art to better understand the essence of the present invention, the technical effects that can be achieved by the above KIPI method of the present invention will be further described below in an embodiment based on CEST imaging.

[0052] Example 1. MRI experiment All phantom and human experiments were performed on a 3 Tesla (T) Siemens scanner (MAGNETOM Prisma, Siemens Healthcare, Erlangen, Germany) using a 64-channel receive head coil. The phantom consisted of a flask filled with 2% agarose gel and two test tubes. One test tube was filled with 10% bovine serum albumin (BSA) dissolved in phosphate-buffered saline (PBS), and the other was filled with 5% bovine serum albumin also dissolved in PBS. The human studies were approved by the local institutional review board.

[0053] The sequence used is a CEST imaging sequence, and the MRI parameter measured is the height of the amide proton transfer effect.

[0054] For the phantom, the CEST scan was performed using a 1.0-s duration, 2-μT saturation pulse followed by a fat-suppressed axial 2D multi-spin echo (TSE) sequence with acquisition parameters of TE = 6.7 ms, TR = 3 s, FA = 90, and FOV = 212 × 186 mm. 2 , resolution = 2.2 × 2.2 mm 2 The slice thickness was 5 mm, the acquisition matrix size was 96 × 96, and the turbo coefficient was 42. A total of 51 frames were acquired at different frequency offsets, including an unsaturated frame S0 and saturated frames between 6 and -6 ppm, with a step width of 0.5 ppm. The signal mean average (NSA) of each saturated frame was 2. The 2D human study used the same parameters as the phantom study.

[0055] For the 3D human brain experiment, data collection was performed using the SPACE-CEST sequence, with the following execution parameters: TE = 17 ms, TR = 3 s, FOV = 212 × 212 × 201 mm. 3 , resolution = 2.8 × 2.8 × 2.8 mm 3The acquisition matrix size is 76 × 76 × 72, turbo factor is 140, NSA is 1.2, and GRAPPA acceleration factor is 2 × 1 (in the phase-coding and hierarchical-coding directions, respectively). ACS uses embedded acquisition, with a matrix size of 24 × 76 × 72. For amide proton transfer weighting (APTw) imaging, a total of seven CEST saturation offset frames are acquired, including (S0), ±3, ±3.5, and ±4 ppm.

[0056] For B0 correction, a GRE sequence with the same field of view, orientation, and resolution as the CEST sequence was used in the 2D and 3D experiments, with a TR of 30 ms. The GRE sequence was performed in double-echo mode, with TEs of 4.92 ms and 9.84 ms, respectively.

[0057] 2. Image reconstruction and analysis All processing and analysis was performed offline using MATLAB (MathWorks, Natick, MA) software, both written on a PC computer (3.2 GHz).

[0058] In the 2D experiment, the first +3.5 ppm frame was selected as the first undersampled image frame, with AF = 2. The second undersampled image frame was selected for undersampling at S0 and ±3.5 ppm (excluding the first undersampled image frame), with AF = 2, and for the other 46 frames, with AF = 4. In the 3D experiment, the +3.5 ppm frame was selected as the first undersampled image frame, with AF = 2 × 2. The remaining six frames were the second undersampled image frames, with AF = 2 × 4 (first 3D experiment) or AF = 4 × 3 (second 3D experiment). Reducing the dimensions of the ACS matrix from 24 × 76 × 72 to 24 × 76 × 24 for the first undersampled image frame meant preserving only a portion of the original ACS. Note that in either the 2D or 3D experiment, ACS data was preserved only in the first undersampled image frame.

[0059] KIPI in this embodiment is divided into the following four steps.

[0060] In step 1, the first undersampled image frame is reconstructed using GRAPPA. This is then reconstructed using GRAPPA with Tikhonov regularization, and the ACS data is inserted to create a corrected frame. The k-space data of the reconstructed corrected frame is Fourier transformed to obtain coil images for each channel. Then, the coil images for all channels are combined using root-square root (RSS) reconstruction. The combined channel image is referred to as the RSS image. The 2D GRAPPA kernel is 4x5, which indicates four phase-coding lines and five frequency-coding points are collected, meaning one missing point is fitted using 20 points in each channel. Similarly, the 3D GRAPPA kernel is 4x5x4 (phase-coding, frequency-coding, and layer-coding directions, respectively).

[0061] In step 2, the coil sensitivity map is then calculated. The coil image for each channel is divided by the RSS image, and the original sensitivity distribution map for each channel is calculated from the reconstructed corrected frame image. The sensitivity map calculated in this manner has the same geometric parameters as the machine-scanned image, so registration is not required. Therefore, support regions are identified from the RSS image by thresholding (the threshold is set close to 0.1), and morphological imaging methods are used to fill cavities in the sensitivity map and smooth the region boundaries. In this embodiment, local weighted projection regression (LWPR) fitting with a cubic weighting kernel is used to smooth noise removal for support regions in the sensitivity map, and non-support regions in the sensitivity map are extrapolated to obtain an optimized sensitivity map. The polynomial used in this embodiment is a quadratic polynomial, with a window width of 12 for support regions and 24 for non-support regions.

[0062] In step 3, a correction coefficient map for suppressing artifacts is calculated. SENSE reconstruction with an acceleration factor of 1 is performed on the correction frame image reconstructed using the optimized sensitivity map to generate a reference image ρ0. Next, backward undersampling is performed on the correction frame to set the acceleration factor of the backward undersampling to the same as that of the second undersampled image frame that needs to be reconstructed subsequently. SENSE reconstruction is then performed on the undersampled data to generate a backward reconstructed image ρ1 with potential aliasing artifacts. The correction coefficient map is defined as the aliasing image divided point by point by the backward reconstructed image, i.e., the pixel value at each position in the correction coefficient map is the quotient of the pixel values at corresponding positions in the reference image and the backward reconstructed image. Furthermore, to remove outliers, the correction coefficient map in this embodiment needs to be filtered using a 3x3 window median filter. The calculation formula for the correction coefficient map C is shown in Equation 5 below:

number

[0063] Finally, in step 4, all remaining second undersampled image frames are reconstructed. In this step, SENSE reconstruction is performed on the second undersampled image frames using the cutoff singular value regularization method to discard singular values 2% smaller than the maximum value. Next, the images generated by the SENSE method are multiplied point-by-point by the correction coefficient map using the artifact suppression method to further reduce the error and obtain complete image frames with artifact suppression.

[0064] It should be noted that since the second undersampled image frame in the 2D experiment of this embodiment has two acceleration factors, i.e., AF=2 and AF=4, in the third step, backward undersampling should be performed according to AF=2 and AF=4, respectively, to generate two different correction factor maps, and the corresponding correction factor maps should also be used for the SENSE reconstruction of the second undersampled image frame in the fourth step.

[0065] For comparison, we perform variable-acceleration GRAPPA reconstruction on backward undersampled data. In such cases, the GRAPPA weights obtained from the ACS of the first undersampled image frame are applied to all other frames. We evaluate the accuracy of the KIPI method by comparing the normalized root mean square error (RNMSE) of GRAPPA, KIPI, and the truth value. When comparing 2D images, the fully sampled image is taken as the truth value; for 3D images, the regular 2x1 GRAPPA image is taken as the truth value because the total sampling time is too long.

[0066] The APTw parameter image is calculated as follows: First, the source image is registered to the first undersampled frame (3.5 ppm). Next, the phase difference between the GRE images acquired with two different TEs is calculated as a B0 map. Then, corrected +3.5-ppm and -3.5-ppm signal values are generated for each voxel based on the calculated B0 map. Finally, the corrected 3.5-ppm and +3.5-ppm images are subtracted to obtain the APTw parameter image.

[0067] 3. Analysis of results As can be seen in Figure 1, the errors of the GRAPPA method (Figure 1b, e) are significantly larger (arrows) than those of the KIPI method (Figure 1c, f). KIPI produces high-quality images (reconstruction error RNMSE = 0.008). Furthermore, with variable acceleration factor undersampling, the enclosed z-spectrum obtained by standard GRAPPA (dashed line in Figure 1d) has a large error compared to the true value (solid line in Figure 1d). However, the results obtained by KIPI (double line in Figure 1d) are nearly identical to the true value.

[0068] Figure 2 shows the APTw parameter image calculated from the source image (Figure 1). The APTw image generated by the KIPI method (Figure 2c) is almost identical to the true value (fully sampled, Figure 2a), with only subtle differences. However, when using AF=4 GRAPPA, the result is a large area of low signal with obvious phase-coding direction aliasing features. KIPI and GRAPPA use the same variable acceleration factor undersampling data, with the first +3.5 ppm selected as the first undersampling image frame with AF=2, S0 and ±3.5 ppm (excluding the first undersampling frame) selected as AF=2, and the remaining 46 frames with AF=4.

[0069] Figure 3 shows an image of a healthy human brain under the same experimental conditions. Similar to the phantom study, the results verify that KIPI's reconstruction is more accurate than GRAPPA's. On the one hand, the agreement between the source image generated by KIPI (Figure 3c) and the true value (Figure 3a) is better than that of the source image reconstructed from the same data by GRAPPA (Figure 3b, with reconstruction errors of 0.018 and 0.029, respectively). On the other hand, the z-spectrum generated by the KIPI method (Figure 3d, double line) is nearly indistinguishable from the full k-space spectrum (Figure 3d, solid line), while the error caused by GRAPPA is significant (Figure 3d, dashed line).

[0070] Figure 4 shows the APTw parameter image generated from the source image (Figure 3). Due to the inaccuracy of the z spectrum, the GRAPPA result (Figure 4b) shows obvious artifacts at the corresponding positions enclosed by the solid lines in Figure 3a. However, compared to the real value (Figure 4a), the image quality obtained by the KIPI method (Figure 4c) is almost unaffected. The lack of enhancement of the APTw image observed here is characteristic of healthy subjects. KIPI and GRAPPA use the same variable acceleration factor undersampling data. The first +3.5 ppm frame is selected as the first undersampling image frame with AF=2, S0 and ±3.5 ppm (excluding the first undersampling image) are selected as AF=2, and the remaining 46 frames are set to AF=4.

[0071] Figure 5 shows a -4 ppm source image acquired by conventional GRAPPA with AF = 2 × 1 (phase-coding and hierarchical-coding directions), as well as the results of the same undersampled data reconstructed by KIPI and conventional GRAPPA. Due to scan time limitations, a fully sampled 3D CEST acquisition cannot be obtained, and therefore the conventional 2 × 1 GRAPPA scan (Figure 5a) is considered to be the true value. Even with the same variable acceleration factor undersampled data, the source image generated by KIPI (Figure 5c) is more true than the GRAPPA (Figure 5b) reconstruction. Aliasing artifacts are also evident in the GRAPPA images (Figures 5b and 5d, arrows).

[0072] Figure 6 shows the APTw parameter image generated by applying B0 correction and image registration to the source image from Figure 5. The APTw image reconstructed by the conventional GRAPPA method (Figure 6b) shows many artifacts, primarily aliasing artifacts in the layer direction. In comparison, the APTw image generated by the KIPI method (Figure 6c) is almost identical to the true image (Figure 6a). The variable acceleration factor undersampling data has ACS data in the +3.5 ppm frame, AF=2×2, and no ACS data in the remaining 6 frames, AF=2×4.

[0073] Figure 7 shows the APTw parameter images of a healthy volunteer when using higher acceleration factors. The +3.5 ppm frame is selected as the first undersampled image frame, with AF = 2 × 2, and the remaining six frames with AF = 4 × 3. Unlike Figure 6(b), the folding artifacts of GRAPPA here appear mainly in the phase-coding direction. Similarly, there are obvious artifacts in the APTw image reconstructed by GRAPPA, but almost none in the KIPI results. A single source image frame with KIPI can reach an acceleration factor of at most 12. In short, the pure effective acceleration factor can reach 8.

[0074] The above-described embodiments are merely preferred solutions of the present invention and are not intended to limit the present invention. Those skilled in the art may make various modifications and variations without departing from the spirit and scope of the present invention. Therefore, any technical solutions obtained by equivalent substitution or equivalent transformation are included in the protection scope of the present invention.

Claims

1. 1. A combined k-space and image space reconstruction imaging method for use in reconstructing undersampled image frames in multi-frame imaging, comprising: The undersampled image frames include a first undersampled image frame having an auto-correction signal and a second undersampled image frame having no auto-correction signal and an acceleration factor of 2 or more, wherein the acceleration factor of the first undersampled image frame is less than or equal to the acceleration factor of the second undersampled image frame; The reconstruction method includes the following steps S1 to S4: In step S1, the first undersampled image frame is reconstructed by a parallel imaging method with automatic correction in k-space to obtain a complete image frame as a correction frame, and the k-space data of the correction frame is Fourier transformed to obtain a coil image of each channel, and the coil images of all channels are combined to obtain a channel combination image; In step S2, the coil image of each channel is divided by the channel-combined image to obtain a coil sensitivity map of each channel, a support region is identified from the channel-combined image, and noise is removed by smoothing for the support region in the sensitivity map, and a non-support region in the sensitivity map is extrapolated to obtain an optimized sensitivity map; In step S3, using the optimized sensitivity map, perform SENSE reconstruction with an acceleration factor of 1 on the correction frame to generate a reference image, and perform backward undersampling on the correction frame to make the acceleration factor of the backward undersampling the same as that of the second undersampling image frame, and further perform SENSE reconstruction on the undersampling data to generate a backward reconstructed image with potential aliasing artifacts, obtain a correction factor map based on the reference image and the backward reconstructed image, and calculate the pixel value of each position in the correction factor map as a quotient of the pixel values of the corresponding positions in the reference image and the backward reconstructed image; In step S4, SENSE reconstruction is performed on the remaining second undersampled image frames other than the first undersampled image frame using the optimized sensitivity map, and the image after SENSE reconstruction is multiplied by the correction coefficient map to obtain a complete image frame that suppresses artifacts. An imaging method for reconstructing a k-space and an image space by combining the k-space and the image space.

2. In step S1, the self-correcting parallel imaging method in k-space is the GRAPPA method, and coil weights are calculated and applied to the undersampled region to reconstruct and acquire a complete image frame.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

3. The GRAPPA method is GRAPPA with Tikhonov regularization 3. The imaging method for reconstructing a k-space and an image space in combination according to claim 2.

4. The first undersampled image frame and the second undersampled image frame are both two-dimensional images, and the acceleration factor of the first undersampled image frame is preferably 2, and the acceleration factor of the second undersampled image frame is preferably 2 to 4.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

5. The first undersampled image frame and the second undersampled image frame are both three-dimensional images, and the acceleration factor of the first undersampled image frame is preferably 2×2, and the acceleration factor of the second undersampled image frame is preferably N×M, where 2≦N≦4, 2≦M≦4, and the total acceleration factor is 12 or less.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

6. In step S1, all coil images are combined using the square root method or the self-adaptive method.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

7. In step S2, noise is removed by smoothing the support region in the sensitivity map using a fitting method.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

8. In step S3, the correction coefficient map needs to be filtered to remove outliers.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

9. In step S4, when performing SENSE reconstruction using the optimized sensitivity map, a cutoff singular value regularization method is used.

2. The imaging method for reconstructing k-space and image space in combination according to claim 1.

10. A magnetic resonance imaging apparatus, The imaging method according to any one of claims 1 to 9 is carried out, comprising a magnetic resonance scanner and a control unit, the control unit storing a computer program, the computer program being executed when the magnetic resonance scanner is adapted to acquire the first undersampled image frame data and the second undersampled image frame data. A magnetic resonance imaging apparatus characterized by:

Images