Image reconstruction method based on subspace parallel imaging, medium and electronic equipment

Abstract

The invention discloses an image reconstruction method based on subspace parallel imaging, a medium and electronic equipment, and belongs to the field of magnetic resonance imaging. The method is used for reconstructing under-sampling data, collected by a magnetic resonance scanner, of a measured object and comprises the following steps that k-space calibration data of the magnetic resonance scanner are obtained, and a calibration matrix is constructed from the k-space calibration data; singular value decomposition is carried out on the calibration matrix, a subspace base vector is extracted, and a k-space convolution kernel is generated based on the subspace base vector; and reconstructing undersampled data by using the obtained k-space convolution kernel. Compared with a traditional method, the method has the advantages that the tolerance to sensitivity graph defects is improved, and the robustness is extremely high.

Description

Technical Field

[0001] The present invention belongs to the field of magnetic resonance imaging, and particularly relates to an image reconstruction method, a medium, and an electronic device based on subspace parallel imaging. Background Art

[0002] Magnetic Resonance Imaging (MRI) is based on the principle of nuclear magnetic resonance. It excites magnetic atomic nuclei to generate signals through radio frequency pulses and uses gradient magnetic fields for spatial encoding, thereby non-invasively obtaining multi-parameter images with high soft tissue contrast. MRI has no ionizing radiation and is highly sensitive to early pathological changes. It is widely used in clinical practice and has become an important means for diagnosing the nervous system, joints, and tumors. However, due to the need for multiple phase encodings, the acquisition time of magnetic resonance is relatively long. And quantitative measurement of physiological parameters requires the acquisition of multi-frame data sets with modulated pulse sequence parameters, such as echo time (TE), flip angle (FA), inversion time (TI), or saturation frequency. This further prolongs the scanning time of MRI.

[0003] Parallel imaging is the most widely used acceleration method in clinical magnetic resonance scanning. It can be roughly divided into two categories: image domain reconstruction algorithms based on explicit coil sensitivity distributions, such as SENSE; k-space methods based on the correlation of multi-channel adjacent k-space data points, such as GRAPPA and SPIRiT. GRAPPA regards reconstruction as an interpolation problem in k-space. In the GRAPPA algorithm, missing points are represented as a linear combination of data at adjacent positions in the k-space of all coils. GRAPPA does not require an additionally acquired sensitivity map and is relatively robust. However, the correlation of k-space data only holds within a small range. When the acceleration factor is relatively high, the reconstruction effect of GRAPPA is poor. The SENSE method directly uses the physical principle of sensitivity weighting and theoretically allows for optimal reconstruction. Due to factors such as the dielectric effect of biological tissues, the inhomogeneity of the main magnetic field, and patient movement, the coil sensitivity maps obtained through pre-scanning often have spatial matching errors and signal-to-noise ratio attenuation, and it is difficult to obtain an accurate sensitivity map. And a tiny error in the sensitivity map will result in significant artifacts in the reconstructed image.

[0004] The SENSE-LORAKS method adds a low-rank constraint to the reconstruction model, can comprehensively utilize sensitivity map information, phase constraints, and image support constraints, and improves the reconstruction quality, but is still affected by sensitivity defects.

[0005] Therefore, it is very important to propose a method that can robustly use the sensitivity map for reconstruction. Summary of the Invention

[0006] The object of the present invention is to solve the problem that the parallel imaging method based on explicit sensitivity map reconstruction is not robust, and to provide an image reconstruction method (abbreviated as SPAN method for Subspace-based parallelimaging hereinafter) and device based on subspace parallel imaging, which can convert the sensitivity map into a k-space convolution kernel and improve the robustness of reconstruction.

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

[0008] In a first aspect, the present invention provides an image reconstruction method based on subspace parallel imaging for reconstructing undersampled k-space data acquired by a magnetic resonance scanner. The method includes the following steps:

[0009] S1: Obtain the k-space calibration data of the magnetic resonance scanner;

[0010] S2: Use a sliding window to slide point by point in the k-space calibration data for data extraction, convert the data block corresponding to the sliding window into a one-dimensional vector, and arrange and combine them in order to form a Hankel calibration matrix;

[0011] S3: Perform singular value decomposition (SVD) on the Hankel calibration matrix, and divide the singular vector matrix into a signal subspace basis vector and a null subspace basis vector according to the magnitudes of the singular values obtained from the decomposition and the truncation threshold;

[0012] S4: Construct a k-space convolution kernel based on the subspace basis vector. Among them, when the subspace basis vector selects the signal subspace basis vector, the corresponding k-space convolution kernel constructed is a self-generated k-space convolution kernel; when the subspace basis vector selects the null subspace basis vector, the corresponding k-space convolution kernel constructed is a null k-space convolution kernel;

[0013] S5: Use the k-space convolution kernel to construct a reconstruction equation, and reconstruct the undersampled k-space data by solving the equation.

[0014] As a preference of the above first aspect, in S1, the k-space calibration data is an auto-calibration signal (ACS) directly acquired by the magnetic resonance scanner, or the central part is taken after the coil sensitivity map of the magnetic resonance scanner is Fourier-transformed into k-space data.

[0015] As a preference of the above first aspect, in S2, the construction method of the Hankel calibration matrix is as follows: vectorize the data block of each sliding window, where the data of different channels are concatenated in the column direction to form a one-dimensional column vector with the same number of elements as the data block, then convert the one-dimensional column vector into a one-dimensional row vector through conjugate transpose, and finally fill the one-dimensional row vectors corresponding to all data blocks into the matrix row by row in order to obtain the Hankel calibration matrix.

[0016] As a preferred embodiment of the first aspect, in S3, the singular vector matrix is a right singular vector matrix.

[0017] As a preferred embodiment of the above-mentioned first aspect, in S3, the threshold of the singular value is 0.02 times the maximum singular value, the singular vectors corresponding to all singular values greater than or equal to the truncation threshold in the singular vector matrix constitute the signal subspace basis vectors, and the singular vectors corresponding to all singular values less than the truncation threshold constitute the zero subspace basis vectors.

[0018] As a preferred embodiment of the first aspect, in S4, according to the selected subspace basis vector V * , the method of constructing the k-space convolution kernel is as follows:

[0019] S41, according to the selected subspace basis vector V * , and the size is M 2 L×M 2 The projection matrix V of L * V * H , M×M is the size of the sliding window, and L is the number of channels of the k-space calibration data.

[0020] S42, the projection matrix V * V * H According to M 2 The row is divided into blocks as a unit, each M 2 row as a separate block matrix; for each block matrix, each M 2 The L-length row vector is reshaped into a M×M×L-sized weight block, and each weight block is filled into a convolution kernel of size (2M-1)×(2M-1)×L and the original value is 0 to form an intermediate convolution kernel, where the i-th element of the weight block corresponding to the i-th row should be placed at the (2M-1)×(2M-1) window center of the intermediate convolution kernel, i∈[1,M 2 ], M 2 The intermediate convolution kernels are added and divided by M 2 , get the convolution kernel corresponding to the current block matrix; the convolution kernels corresponding to all L block matrices form a k-space convolution kernel KN of size (2M-1)×(2M-1)×L×L; where the subspace basis vector V * is the signal subspace basis vector V || When the k-space convolution kernel KN is obtained, it is the self-generated k-space convolution kernel G, and the subspace basis vector V * is the zero subspace basis vector V ⊥ When , the obtained k-space convolution kernel KN is the zero-set k-space convolution kernel N.

[0021] Preferably, in the above first aspect, in S5, when the k-space convolution kernel is the self-generated k-space convolution kernel G, the constructed reconstruction equation is:

[0022]

[0023] When the k-space convolution kernel is the zeroed k-space convolution kernel N, the constructed reconstruction equation is:

[0024]

[0025] Where X is the fully sampled multi-channel k-space data, which is composed of the unacquired k-space data Xu and the acquired k-space data Xa, that is, X = Xu + Xa; is a composite convolution operation, representing performing k-space convolution on each channel and then adding them.

[0026] Preferably, in the above first aspect, the algorithm for solving the reconstruction equation is the conjugate gradient descent algorithm, and the object to be solved is the unacquired k-space data Xu.

[0027] Preferably, in the above first aspect, in S4, a coil sensitivity map is further calculated based on the k-space convolution kernel.

[0028] In a second aspect, the present invention provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, an image reconstruction method based on subspace parallel imaging as described in any one of the above first aspect solutions is implemented.

[0029] In a third aspect, the present invention provides a computer electronic device, which includes a memory and a processor;

[0030] The memory is used to store a computer program;

[0031] The processor is used to implement an image reconstruction method based on subspace parallel imaging as described in any one of the above first aspect solutions when executing the computer program.

[0032] The present invention has the following beneficial effects compared with the prior art.

[0033] The present invention can transform the sensitivity map into a k-space convolution kernel, and then reconstruct the undersampled data, greatly improving the tolerance to defects in the sensitivity map and having higher robustness than the traditional SENSE method. Of course, for data with ACS but without a sensitivity map, the present invention can also directly generate a k-space convolution kernel from the ACS; moreover, compared with the prior art, the present invention has high computational efficiency in solving the k-space convolution kernel and short computational time. Description of the Drawings

[0034] Figure 1 Schematic diagram of the SPAN method. (A) Fourier transform the sensitivity map into k-space data. (B) Construct a calibration matrix from the generated ACS data and perform SVD on it to obtain subspaces. (C) Generate the k-space self-generated convolution kernel G and the nulling convolution kernel N from the signal subspace or null space basis vectors.

[0035] Figure 2 Coil sensitivity maps calculated from the k-space convolution kernel using compressed four-channel data. (A, C): Image domain maps obtained by applying zero-padding and inverse Fourier transform to the k-space convolution kernels G and N. The sensitivity map estimates can be obtained after performing voxel-by-voxel eigenvalue decomposition on the image domain maps. (B): The eigenvectors corresponding to the eigenvalue = "1" form the sensitivity map (red box). (D): The eigenvectors corresponding to the eigenvalue = "0" form the sensitivity map (red box).

[0036] Figure 3 Calibration time comparison of different algorithms. The time for calculating the convolution using the original SPIRiT (vertical lines), the SPIRiT algorithm based on Cholesky decomposition (horizontal lines), and SPAN (diagonal lines) on 16-, 24-, and 32-channel data.

[0037] Figure 4 Reconstruction result comparison of sensitivity maps using different signal support region sizes. SENSE and SPAN reconstruct the image with R = 4 (B) using sensitivity maps with overly large (C), appropriate (D), and overly small (E) signal support regions. The error maps compared to the fully sampled image (A) are also shown in the figure, and the nRMSE is shown at the bottom of the error maps. The white circles indicate the signal regions.

[0038] Figure 5 Reconstruction image comparison of partial Fourier undersampled data. SENSE, SENSE-LORAKS, SPAN, and SPAN-LORAKS reconstruct the partial Fourier undersampled data (B) using sensitivity maps with overly large (C), appropriate (D), and overly small (E) signal support regions. The error maps compared to the fully sampled image (A) are also shown in the figure, and the nRMSE is shown at the bottom of each error map. The white circles indicate the signal regions.

[0039] Figure 6Results comparison for reconstruction using a sensitivity map with holes. Reconstruction results of SENSE, SENSE-LORAKS, SPAN, and SPAN-LORAKS using a sensitivity map with holes (B) for partial Fourier undersampled data (D), including image results (E) and k-space results (G). Errors (F, H) compared to the fully sampled reference standard (A, C) are also shown in the figure, and the nRMSE is shown at the bottom of the image domain error map. The white circles indicate the signal regions. Detailed implementation manners

[0040] To make the above objects, features, and advantages of the present invention more obvious and understandable, the following detailed description of the specific implementation manners of the present invention will be given in conjunction with the accompanying drawings. Many specific details are set forth in the following description to fully understand the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without departing from the connotation of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed below. The technical features in various embodiments of the present invention can be combined correspondingly without conflict.

[0041] In a preferred embodiment of the present invention, an image reconstruction method based on subspace parallel imaging (denoted as the SPAN method) is provided, which includes steps S1 to S5, and the exemplary process is as Figure 1 shown. The following will specifically describe the overall process in conjunction with Figure 1 the accompanying drawings.

[0042] S1: Obtain the k-space calibration data of the magnetic resonance scanner. Among them, there are two ways to obtain the k-space calibration data: for the additionally acquired sensitivity map, as Figure 1 (A) shows, the present invention performs Fourier transform on the sensitivity map to obtain k-space data, and selects the central part of the k-space data as the ACS required in this step; if there is an auto-calibration signal in the acquired data, it can be directly used as the ACS required in this step.

[0043] S2: Use a sliding window to slide point by point in the k-space calibration data for data extraction, convert the data block corresponding to the sliding window into a one-dimensional vector, and arrange and combine them in order to form a structured block Hankel calibration matrix A.

[0044] For ease of description, the number of channels of the k-space calibration data is denoted as L, and the L channels are sequentially denoted as Ch1 to ChL; at the same time, the size of the sliding window is denoted as M×M. Thus, the sliding window can slide in the k-space calibration data in the order from top to bottom and from left to right, and the sliding step of the sliding window is 1 k-space point. As Figure 1 (B) shows, the size of each data block extracted by the sliding window is M×M×L.

[0045] In an embodiment of the present invention, continuing to refer to Figure 1 (B) as shown, the method for constructing the Hankel calibration matrix A based on all the extracted data blocks is as follows: Vectorize each M×M×L data block. First, arrange the data of each channel as a column vector with the number of elements being M 2 ×1, and then splice the data of different channels in the column direction to form a one-dimensional column vector b with the same number of elements as the total number of elements of the data block, which is M 2 L×1. Then, convert the one-dimensional column vector b of M 2 L×1 into a one-dimensional row vector b of 1×M 2 L through conjugate transpose H , and finally fill the one-dimensional row vectors b corresponding to all the data blocks H into the matrix row by row in sequence to obtain the Hankel calibration matrix A of K×M 2 L, where K represents the total number of data blocks extracted by the sliding window from the calibration data ACS in the k space.

[0046] S3: Perform singular value decomposition (SVD) on the Hankel calibration matrix, and divide the singular vector matrix into a signal subspace basis vector and a zeroed subspace basis vector according to the magnitudes of the decomposed singular values and the truncation threshold.

[0047] In an embodiment of the present invention, the Hankel matrix usually has good low-rank properties, so SVD can be performed on A for subspace division:

[0048] A = USV H [1]

[0049] SVD can obtain two singular vector matrices, namely the left singular vector matrix U and the right singular vector matrix V. Each singular value vector in the singular vector matrix corresponds to a singular value, and the subspace basis vectors can be divided according to the magnitudes of the decomposed singular values and the set truncation threshold. In the present invention, the right singular vector matrix V is preferably used for the division of the subspace basis vectors. Specifically, a truncation threshold can be set according to the magnitudes of the singular values. The singular vectors corresponding to all the singular values greater than or equal to the truncation threshold in the singular vector matrix form the signal subspace basis vector V || , and the singular vectors corresponding to all the singular values less than the truncation threshold form the zeroed subspace basis vector V ⊥ . The specific value of the truncation threshold can be optimized according to the actual situation. In an embodiment of the present invention, the truncation threshold is preferably 2% of the largest singular value s max after SVD, that is, 0.02s max . Thus, continuing to refer to Figure 1 (B) as shown, this step is equivalent to passing through the truncation threshold 0.02s maxThe singular vector matrix V is partitioned into the signal subspace basis vectors V || and the zeroed subspace basis vectors V ⊥ .

[0050] S4: Construct the k-space convolution kernel KN based on the subspace basis vectors, where the selected subspace basis vectors V * can be the signal subspace basis vectors V || , or can be the zeroed subspace basis vectors V ⊥ . If the subspace basis vectors V * select the signal subspace basis vectors V || , construct the self-generated k-space convolution kernel G. If the subspace basis vectors V * select the zeroed subspace basis vectors V ⊥ , then construct the zeroed k-space convolution kernel N. The basic principle of constructing the k-space convolution kernel KN based on the two types of subspace basis vectors is similar. Using the subspace basis vectors V * to generically represent V || and V ⊥ , the method for constructing the k-space convolution kernel KN(V || and V ⊥ corresponding to G and N respectively) is as follows:

[0051] S41. According to the selected subspace basis vectors V * , obtain the projection matrix V 2 V 2 of size M * L×M * H .

[0052] S42. Partition the projection matrix V * V * H by taking M 2 rows as a unit. Each M 2 rows forms a separate block matrix; for each block matrix, reorganize each M 2 L-long row vector into a weight block of size M×M×L, and fill each weight block into a convolution kernel of size (2M - 1)×(2M - 1)×L with all original values being 0 to form an intermediate state convolution kernel (the filled area is replaced by the weight block, and the remaining positions remain 0). Among them, the i-th element of the weight block corresponding to the i-th row should be placed at the center of the (2M - 1)×(2M - 1) window of the intermediate state convolution kernel, i ∈ [1, M 2 . Add the M 2 intermediate state convolution kernels and then divide by M 2, the convolution kernel corresponding to the current block matrix is obtained; the convolution kernels corresponding to all L block matrices form a k-space convolution kernel KN of size (2M - 1)×(2M - 1)×L×L; wherein the subspace basis vector V * is the signal subspace basis vector V || When, the obtained k-space convolution kernel KN is the self-generated k-space convolution kernel G, and the subspace basis vector X * is the zeroed subspace basis vector V ⊥ When, the obtained k-space convolution kernel KN is the zeroed k-space convolution kernel N.

[0053] In an embodiment of the present invention, as shown in Figure 1 (C), taking the signal subspace basis vector V || as an example, any k-space data block b satisfies:

[0054] V || V || H b = b[2]

[0055] where b is the vectorized multi-channel k-space data block.

[0056] V || V || H can generate all points in the data block from the k-space data block itself, but V || V || H itself is not a convolution kernel. To generate a convolution kernel from V || V || H , consider a V corresponding to an M×M sliding window || V || H ( Figure 1 where M is 3). When this V || V || H slides throughout the k-space, each point in the k-space will be reconstructed M 2 times because each point will have M 2 different positions in the M×M sliding window. For illustration, take the point r in channel 1 as an example. First, the point r is located at the upper left corner of the M×M sliding window. Note that, the inner product of the first row weights of V || V || H with b can generate the first data point in b, that is, the point r; and this inner product operation is equivalent to multiplying V || V || HThe weight arrangement of the first row is a kernel of M×M×L that performs a composite convolution operation with the M×M×L data block represented by b. To conform to the property of the convolution kernel - the generated points are at the center of the block, V || V || H The M×M×L kernel formed by the weights of the first row of V is combined with the surrounding zeros to form a convolution kernel of size (2M - 1)×(2M - 1)×L ( Figure 1 is 5×5×L in [])), where the position of point r is at the center of (2M - 1)×(2M - 1). And so on, V || V || H The weights from the first row to the M 2 th row of V can generate M 2 intermediate-state convolution kernels of size (2M - 1)×(2M - 1)×L. From the perspective of least squares, the optimal convolution kernel for generating channel 1 is the average of these M 2 convolution kernels, that is, the average of the M 2 weight values at the corresponding positions. Then, every M 2 rows of weights can obtain the convolution kernel corresponding to generating the data of one channel. According to the same method, the convolution kernels of all L channels can be obtained, and the convolution kernels of all channels together form the self-generated k-space convolution kernel G, with a size of (2M - 1)×(2M - 1)×L×L ( Figure 1 is 5×5×L×L in []).

[0057] On the other hand, similar to the signal subspace basis vector V || when constructing the nulling k-space convolution kernel N based on the nulling subspace basis vector V ⊥ , V ⊥ V ⊥ H can be used to replace V || V || H as Figure 1 (C) the input of the construction process. At this time, since V ⊥ V ⊥ H b = 0, the nulling k-space convolution kernel N can be generated similarly.

[0058] S5: Use the above k-space convolution kernel KN (self-generated k-space convolution kernel G or nulling k-space convolution kernel N) to construct a reconstruction equation, and reconstruct the undersampled k-space data by solving the equation.

[0059] In the present invention, the complete multi-channel k-space data X satisfies:

[0060]

[0061] where and is the weight corresponding to the j-th channel when generating the k-space data of the i-th channel, x i is the k-space data of the i-th channel, and L is the number of channels. The symbol represents a composite convolution operation:

[0062]

[0063] In the formula: * is the standard convolution operation.

[0064] For the zeroed k-space convolution kernel N, there is also:

[0065]

[0066] Therefore, Formula [3] and Formula [5] can be used as models for reconstructing the complete k-space, and they are equivalent. Any one of the formulas can be used to reconstruct the complete k-space. The fully sampled complete multi-channel k-space data X is composed of the unacquired k-space data Xu and the acquired k-space data Xa, that is, X = Xu + Xa.

[0067] Taking Formula [3] as an example, the present invention converts the reconstruction equation corresponding to Formula [3] into reconstructing the undersampled data by solving the following optimization problem:

[0068]

[0069] Taking Formula [5] as an example, the present invention converts the reconstruction equation corresponding to Formula [5] into reconstructing the undersampled data by solving the following optimization problem:

[0070]

[0071] The optimization problem of the above reconstruction equation can be solved using the conjugate gradient descent algorithm, and the object to be solved is the unacquired k-space data Xu.

[0072] The above S1 to S5 constitute the entire process of the SPAN method. The SPAN method can solve for the unacquired k-space data Xu, and then reconstruct the complete fully sampled k-space data X.

[0073] In addition, in addition to reconstructing the complete fully sampled k-space data X, the SPAN method of the present invention can further calculate the coil sensitivity map based on the k-space convolution kernel (self-generated k-space convolution kernel G or zeroed k-space convolution kernel N). The calculation methods of the coil sensitivity map are described separately below:

[0074] For the self-generated k-space convolution kernel G, it can be first zero-padded to the image size, and then the self-generated image domain map can be obtained through inverse Fourier transform such asFigure 2 As shown in (A), each voxel in the figure corresponds to an L×L correlation matrix Therefore, from the self-generated image domain map the correlation matrix of each voxel can be obtained For each voxel's After performing eigenvalue decomposition, theoretically, the eigenvectors corresponding to eigenvalue 1 can form a set of coil sensitivity maps, as Figure 2 shown in (B). However, eigenvalue 1 is a theoretical value. In practical applications, a threshold close to 1 (such as 0.9) can be set, and the eigenvectors corresponding to each eigenvalue greater than this threshold can form a set of coil sensitivity maps

[0075] Correspondingly, for the zeroed k-space convolution kernel N, it can also be padded with zeros to the image size first, and then the zeroed image domain map can be obtained through inverse Fourier transform As Figure 2 shown in (C), each voxel in the figure corresponds to an L×L correlation matrix Therefore, from the zeroed image domain map the correlation matrix of each voxel can be obtained For each voxel's After performing eigenvalue decomposition, theoretically, the eigenvectors corresponding to eigenvalue 0 can form a set of coil sensitivity maps, as Figure 2 shown in (D). However, eigenvalue 0 is a theoretical value. In practical applications, a threshold close to 0 (such as 0.1) can be set, and the eigenvectors corresponding to each eigenvalue less than this threshold can form a set of coil sensitivity maps

[0076] It should be noted that the method steps shown in S1 to S5 above can essentially be implemented in the form of a computer program

[0077] Thus, based on the same inventive concept, the present invention also provides a computer electronic device corresponding to the image reconstruction method based on subspace parallel imaging provided in the above embodiment, which includes a memory and a processor

[0078] The memory is used to store a computer program

[0079] The processor is used to implement the image reconstruction method based on subspace parallel imaging as described above when executing the computer program

[0080] In addition, when the logical instructions in the above-mentioned memory are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on such understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention.

[0081] Therefore, based on the same inventive concept, the present invention provides a computer-readable storage medium corresponding to an image reconstruction method based on subspace parallel imaging. A computer program is stored on the storage medium, and when the computer program is executed by a processor, it can implement the image reconstruction method based on subspace parallel imaging as described above.

[0082] Therefore, based on the same inventive concept, the present invention provides a computer program product, including a computer program / instructions. When the computer program / instructions are executed by a processor, they can implement the image reconstruction method based on subspace parallel imaging as described above.

[0083] Specifically, when the computer program in the above three embodiments is executed by a processor, the steps of S1 to S5 described above can be executed.

[0084] It can be understood that the above storage medium may include a random access memory (RAM), and may also include a non-volatile memory (NVM), such as at least one disk memory. At the same time, the storage medium may also be various media such as a USB flash drive, a mobile hard disk, a magnetic disk, or an optical disc that can store program codes.

[0085] It can be understood that the above-mentioned processor may be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it may 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.

[0086] It should be further noted that those skilled in the art can clearly understand that for the convenience and conciseness of description, the specific working process of the system described above can refer to the corresponding process in the foregoing method embodiments, and will not be elaborated herein. In the embodiments provided in the present application, the division of steps or modules in the system and method is only a logical function division, and there may be other division methods in actual implementation. For example, multiple modules or steps can be combined or integrated together, and a module or step can also be split.

[0087] Next, an embodiment will be used to further demonstrate the technical effects that can be achieved by the SPAN method described in S1 to S5 of the present invention, so that those skilled in the art can better understand the essence of the present invention.

[0088] Embodiment

[0089] 1. MRI Experiment

[0090] The human experiment was carried out on a 3T Siemens magnetic resonance scanner (MAGNETOM Prisma), using a 64-channel head receive coil. The human experiment was approved by the local institutional review board. The sequence used was the TSE imaging sequence. The sequence parameters are as follows: TR = 3000 ms, TE = 7.5 ms, FOV = 212×212 mm 2 , Resolution = 1.1×1.1 mm 2 , Slice thickness = 5 mm, Acquisition matrix = 192×192, Echo train length = 24.

[0091] The GRE sequence was used for reference scanning to obtain the sensitivity map required for standard SENSE. The acquisition parameters are as follows: Flip angle = 1°; TR = 4 ms; TE = 1.9 ms; FOV = 450×300×300 mm 3 ; Resolution = 4.7×4.7×3 mm 3 ; Acquisition time = 51 s.

[0092] 2. Image Reconstruction and Analysis

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

[0094] Since Equation [3] and Equation [5] are completely equivalent, in this embodiment, the image will be reconstructed based on Equation [5]. Specifically, when the SPAN method corresponding to the above S1 to S5 steps reconstructs the unacquired k-space data Xu in this embodiment, it is achieved by solving the following optimization problem:

[0095]

[0096] Where Xa represents the collected data, and Xu + Xa = N. This optimization problem can be solved by the conjugate gradient descent algorithm.

[0097] In addition, when the calibration information comes from the additionally collected sensitivity, SPAN can combine the S-type LORAKS constraint to reconstruct the data with partial Fourier downsampling:

[0098]

[0099] Where P is the LORAKS operator that transforms the k-space data into the S matrix in the LORAKS model; μ is the expected rank parameter of the LORAKS matrix. This embodiment uses the iterative maximum-minimization algorithm to solve this optimization problem.

[0100] In addition, in order to analyze the convolution kernel calculation efficiency of the SPAN method of the present invention, it is compared with the SPIRiT method. 24 central k-space lines are selected from the fully sampled data as the ACS. For the SPAN method, the sliding window size for generating the Hankel calibration matrix is set to 7×7, and the truncation value of the singular values is 2% of the largest singular value. The tolerance of the conjugate gradient descent algorithm is set to 1e-3. For the SPIRiT method, the same sliding window size and tolerance are used. According to previous studies, two methods are used to calculate the convolution kernel of SPIRiT. One is the original calibration algorithm with a computational complexity of O(L 4 ) and the other is the calibration algorithm based on Cholesky decomposition with a computational complexity of O(L 3 ). The channel compression algorithm is used to compress the k-space data into 16-, 24-, and 32-channel data to compare the calibration times of different methods.

[0101] To analyze the robustness of the SPAN method to the sensitivity map, retrospective uniform downsampling is performed on the fully sampled k-space data without retaining the ACS. The sensitivity maps for SPAN and SENSE reconstructions are calculated from the SENSE reference scan through local weighted polynomial regression, registration, and interpolation. To test the robustness of SPAN to the signal support region, sensitivity maps with three different signal support region sizes are generated, and then the undersampled data with an acceleration factor R = 4 are reconstructed using SENSE and SPAN respectively. In addition, partial Fourier downsampling experiments are performed on the sensitivity maps with three different signal support region sizes and a sensitivity map with a hole. The undersampled data are reconstructed using SENSE, SENSE-LORAKS, SPAN, and SPAN-LORAKS, and the total acceleration factor R = 4.5, where the partial Fourier factor PF = 2 / 3. The μ of SPAN-LORAKS is set to the number of basis vectors of the signal subspace.

[0102] In this embodiment, the SPAN method uses adaptive channel merging. The normalized root mean square error (nRMSE) is calculated as ‖x - y‖2 / ‖x‖2, where x is the fully sampled reference standard image and y is the reconstructed image.

[0103] 3. Result Analysis

[0104] Figure 2 Shows the coil sensitivity maps calculated from the k-space convolution kernel using compressed four-channel data. The image domain maps obtained by applying zero-padding and inverse Fourier transform to the self-generated k-space convolution kernel G has a size of 192×192×4×4. For convenience of display, the individual 192×192 maps are displayed sequentially and arranged in a 4×4 pattern. Thus, each voxel has a 4×4 correlation matrix. After performing per-voxel eigenvalue decomposition on (B), the eigenvectors corresponding to the eigenvalue = "1" form the sensitivity map (red frame). Similarly, the image domain map obtained by applying zero-padding and inverse Fourier transform to the zeroed k-space convolution kernel N (C). After performing per-voxel eigenvalue decomposition on (D), the eigenvectors corresponding to the eigenvalue = "0" form the sensitivity map (red frame).

[0105] Figure 3 Shows the convolution kernel calculation times of the original SPIRiT (vertical lines), the SPIRiT algorithm based on Cholesky decomposition (horizontal lines), and SPAN (diagonal lines) on 16, 24, and 32-channel data. Consistent with previous studies, the calculation time of the SPIRiT algorithm based on Cholesky decomposition is less than that of the original SPIRiT algorithm. It is worth noting that although the calibration calculation complexities of the SPAN method and the SPIRiT algorithm based on Cholesky decomposition are of the same order of magnitude, the calculation time of SPAN is significantly less than that of the SPIRiT algorithm based on Cholesky decomposition.

[0106] Figure 4 Shows the reconstruction results of SENSE and SPAN when using sensitivity maps with three different signal support region sizes. When the signal support region is too large, the sensitivity map will contain noise; while when the signal support region is too small, the sensitivity map will not fully cover the image signal. Obviously, whether the sensitivity map is too large or too small, the SENSE method will result in strong artifacts (C, E); only when the signal support region is appropriate (as much as possible to contain the image signal but not contain noise), the SENSE method can produce relatively good reconstruction results (D). However, the SPAN method shows stable results in all cases and has smaller errors.

[0107] Figure 5 It shows that the SPAN method can well combine the LORAKS constraint to reconstruct data with partial Fourier downsampling. When the signal support region is too large or too small, obvious artifacts are generated in the results of the SENSE and SENSE-LORAKS methods (C, E); when the signal support region is appropriate, SENSE-LORAKS shows better results (D). In contrast, the error of the SPAN method is less than that of SENSE in all cases; the SPAN-LORAKS method produces stable results in all cases, and the error is less than that of SENSE-LORAKS in the optimal case.

[0108] Figure 6 It shows the reconstruction results of SENSE(-LORAKS) and SPAN(-LORAKS) when a hole is introduced in the sensitivity map. SENSE and SENSE-LORAKS cannot reconstruct the signal value at the hole, and the signal here will leak to other positions. The SPAN method is almost unaffected by the hole and can correctly reconstruct the signal value at the hole; while SPAN-LORAKS further restores the k-space data in the partial Fourier region, and the image is clearer and the error is smaller.

[0109] The above-described embodiments are only a preferred solution of the present invention, but they are not intended to limit the present invention. Those of ordinary skill in the relevant technical field can still make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all technical solutions obtained by adopting the equivalent replacement or equivalent transformation method fall within the protection scope of the present invention.

Claims

1. An image reconstruction method based on subspace parallel imaging for reconstructing undersampled k-space data acquired by a magnetic resonance scanner, characterized in that, The method includes the following steps: S1: Obtain the k-space calibration data of the magnetic resonance scanner; S2: Use a sliding window to slide point by point in the k-space calibration data for data extraction, convert the data block corresponding to the sliding window into a one-dimensional vector, and arrange and combine them in order to form a Hankel calibration matrix; S3: Perform singular value decomposition on the Hankel calibration matrix. According to the magnitudes of the singular values obtained from the decomposition and the truncation threshold, divide the singular vector matrix into a signal subspace basis vector and a null subspace basis vector; S4: Construct a k-space convolution kernel based on the subspace basis vector. Among them: when the subspace basis vector selects the signal subspace basis vector, the corresponding k-space convolution kernel constructed is a self-generated k-space convolution kernel; when the subspace basis vector selects the null subspace basis vector, the corresponding k-space convolution kernel constructed is a null k-space convolution kernel; S5: Use the k-space convolution kernel to construct a reconstruction equation, and reconstruct the undersampled k-space data by solving the equation.

2. The image reconstruction method based on subspace parallel imaging according to claim 1, wherein In S1, the k-space calibration data is an auto-calibration signal directly collected by the magnetic resonance scanner, or the coil sensitivity map of the magnetic resonance scanner is Fourier-transformed into k-space data and then the central part is taken.

3. The image reconstruction method based on subspace parallel imaging according to claim 1, wherein In S2, the construction method of the Hankel calibration matrix is as follows: Vectorize the data block of each sliding window, where the data of different channels are concatenated in the column direction to form a one-dimensional column vector with the same number of elements as the data block, then convert the one-dimensional column vector into a one-dimensional row vector through conjugate transpose, and finally fill all the one-dimensional row vectors corresponding to the data blocks into the matrix row by row in order to obtain the Hankel calibration matrix.

4. The image reconstruction method based on subspace parallel imaging according to claim 1, wherein In S3, the singular vector matrix is a right singular vector matrix.

5. The image reconstruction method based on subspace parallel imaging according to claim 1, characterized in that, In S3, the threshold of the singular value is 0.02 times the maximum singular value. The singular vectors corresponding to all singular values greater than or equal to the truncation threshold in the singular vector matrix form the signal subspace basis vector, and the singular vectors corresponding to all singular values less than the truncation threshold form the null subspace basis vector.

6. The image reconstruction method based on subspace parallel imaging according to claim 1, wherein In S4, according to the selected subspace basis vector V * , the method for constructing the k-space convolution kernel is as follows: S41. According to the selected subspace basis vectors V * , a projection matrix V of size M 2 L×M 2 L is obtained * V * H , where M×M is the size of the sliding window and L is the number of channels of the k-space calibration data; S42. Divide the projection matrix V * V * H into blocks with M 2 rows as a unit, and each M 2 rows forms a separate block matrix; For each block matrix, each M 2 -length row vector is reshaped into a weight block of size M×M×L, and each weight block is filled into a convolutional kernel of size (2M - 1)×(2M - 1)×L with all original values being 0 to form an intermediate convolutional kernel. The i-th element of the weight block corresponding to the i-th row should be placed at the center of the (2M - 1)×(2M - 1) window of the intermediate convolutional kernel, where i ∈ [1, M 2 , and M 2 intermediate convolutional kernels are added and then divided by M 2 to obtain the convolutional kernel corresponding to the current block matrix; The convolution kernels corresponding to all L block matrices form a k-space convolution kernel KN of size (2M - 1) × (2M - 1) × L × L; where the subspace basis vector V * is the signal subspace basis vector V || When, the obtained k-space convolution kernel KN is the self-generated k-space convolution kernel G, and the subspace basis vector V * is the zeroed subspace basis vector V ⊥ When, the obtained k-space convolution kernel KN is the zeroed k-space convolution kernel N.

7. The image reconstruction method based on subspace parallel imaging according to claim 1, characterized in that, In S5, when the k-space convolution kernel is a self-generated k-space convolution kernel G, the constructed reconstruction equation is: In S5, when the k-space convolution kernel is a null k-space convolution kernel N, the constructed reconstruction equation is: Where X is the fully sampled multi-channel k-space data, which consists of the unacquired k-space data Xu and the acquired k-space data Xa, i.e., X = Xu + Xa; is a composite convolution operation, representing performing k-space convolution on each channel and then adding them up; Preferably, the algorithm for solving the reconstruction equation is the conjugate gradient descent algorithm, and the object to be solved is the unacquired k-space data Xu.

8. The image reconstruction method based on subspace parallel imaging according to claim 1, wherein In S4, further calculate the coil sensitivity map based on the k-space convolution kernel.

9. A computer-readable storage medium, characterized in that, A computer program is stored on the storage medium. When the computer program is executed by a processor, it implements the image reconstruction method based on subspace parallel imaging according to any one of claims 1 to 8.

10. A computer electronic device, characterized in that, Comprising a memory and a processor; The memory is used to store a computer program; The processor is used to implement the image reconstruction method based on subspace parallel imaging according to any one of claims 1 to 8 when executing the computer program.

Images