Image reconstruction method using subspace-based parallel imaging, medium and electronic device

Abstract

An image reconstruction method using subspace-based parallel imaging, a medium, and an electronic device, relating to the field of magnetic resonance imaging. The method comprises the following steps: acquiring k-space calibration data of a magnetic resonance scanner, and constructing a calibration matrix therefrom; performing singular value decomposition on the calibration matrix, and extracting a subspace basis vector, so as to generate a k-space convolution kernel on the basis of the subspace basis vector; and reconstructing undersampled data using the obtained k-space convolution kernel. Compared with conventional methods, the method improves the tolerance to sensitivity map defects and has extremely high robustness.

Description

Image reconstruction methods, media, and electronic devices based on subspace parallel imaging Technical Field

[0001] This invention belongs to the field of magnetic resonance imaging, specifically relating to an image reconstruction method, medium, and electronic device based on subspace parallel imaging. Background Technology

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

[0003] Parallel imaging is the most widely used acceleration method in clinical magnetic resonance imaging (MRI). It can be broadly divided into two categories: image domain reconstruction algorithms based on explicit coil sensitivity distribution, such as SENSE; and k-space methods based on the correlation of adjacent k-space data points in multiple channels, such as GRAPPA and SPIRiT. GRAPPA treats reconstruction as an interpolation problem in k-space. In the GRAPPA algorithm, missing points are represented as a linear combination of all adjacent k-space data points of the coils. GRAPPA does not require additional sensitivity maps and is relatively robust. However, the correlation of k-space data only holds true within a small range, and the reconstruction effect of GRAPPA is poor when the acceleration factor is high. The SENSE method directly uses the physical principle of sensitivity weighting, theoretically allowing for optimal reconstruction. Due to interference from factors such as the dielectric effect of biological tissue, 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, making it difficult to obtain accurate sensitivity maps. Even small errors in the sensitivity map can lead to significant artifacts in the reconstructed image.

[0004] The SENSE-LORAKS method incorporates low-rank constraints into the reconstruction model, which can comprehensively utilize sensitivity map information, phase constraints, and image support constraints, thus improving reconstruction quality. However, it is still affected by sensitivity defects.

[0005] Therefore, it is important to propose a robust method for reconstruction using sensitivity maps. Summary of the Invention

[0006] The purpose of this invention is to solve the problem of the lack of robustness of parallel imaging methods based on explicit sensitivity map reconstruction, and to provide an image reconstruction method and apparatus based on subspace-based parallel imaging (hereinafter referred to as SPAN method), which can convert the sensitivity map into a k-space convolution kernel to improve the robustness of reconstruction.

[0007] The specific technical solution adopted in this invention is 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 comprising the following steps:

[0009] S1: Acquire k-space calibration data from the magnetic resonance scanner;

[0010] S2: In the k-space calibration data, a sliding window is used to extract data point by point. The data blocks corresponding to the sliding window are converted into one-dimensional vectors and arranged in order to form the Hankel calibration matrix.

[0011] S3: Perform singular value decomposition (SVD) on the Hankel calibration matrix. Based on the magnitude of the singular values ​​obtained from the decomposition and the truncation threshold, divide the singular vector matrix into signal subspace basis vectors and zeroing subspace basis vectors.

[0012] S4: Construct k-space convolution kernels based on subspace basis vectors, where: when the subspace basis vectors are selected from the signal subspace basis vectors, the corresponding k-space convolution kernels constructed are self-generated k-space convolution kernels; when the subspace basis vectors are selected from the zeroed subspace basis vectors, the corresponding k-space convolution kernels constructed are zeroed k-space convolution kernels.

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

[0014] As a preferred embodiment of the first aspect above, in S1, the k-space calibration data is an automatic calibration signal (ACS) directly acquired by a magnetic resonance scanner, or the center portion of the coil sensitivity map of the magnetic resonance scanner is taken after Fourier transform.

[0015] As a preferred embodiment of the first aspect above, in S2, the method for constructing the Hankel calibration matrix is ​​as follows: the data block of each sliding window is vectorized, wherein the data of different channels are concatenated in the column direction to form a one-dimensional column vector with the same total number of elements as the data block. Then, the one-dimensional column vector is converted into a one-dimensional row vector through conjugate transpose. Finally, the one-dimensional row vectors corresponding to all data blocks are sequentially filled into the matrix row by row to obtain the Hankel calibration matrix.

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

[0017] As a preferred embodiment of the first aspect above, in S3, the threshold for singular values ​​is 0.02 times the maximum singular value, and 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 vectors, while the singular vectors corresponding to all singular values ​​less than the truncation threshold form the zeroing subspace basis vectors.

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

[0019] S41. Based on the selected subspace basis vector V * The resulting 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, Project the matrix V * V * H According to M 2 Divide the data into blocks, with each block consisting of a row and a unit of M. 2 The rows are treated as a separate block matrix; for each block matrix, each M is... 2 A row vector of length L is reshaped into a weight block of size M×M×L. Each weight block is then filled into a convolutional kernel of size (2M-1)×(2M-1)×L with all initial values ​​being 0, forming 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 ], M 2 The sum of the intermediate convolutional kernels divided 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 * V is the basis vector of the signal subspace. || At that time, the obtained k-space convolution kernel KN is a self-generated k-space convolution kernel G, and the subspace basis vector V * Let V be the basis vector of the null subspace. ⊥ When the k-space convolution kernel KN is obtained, it is the zeroed k-space convolution kernel N.

[0021] As a preferred embodiment of the first aspect above, in step S5, when the k-space convolution kernel is a self-generated k-space convolution kernel G, the constructed reconstruction equation is:

[0022] When the k-space convolution kernel is a zero-set k-space convolution kernel N, the constructed reconstruction equation is:

[0023] Where X is the fully sampled multi-channel k-space data, which consists of the uncollected k-space data Xu and the collected k-space data Xa, i.e., X = Xu + Xa; It is a compound convolution operation, which means performing k-space convolution on each channel and then summing them.

[0024] As a preferred embodiment of the first aspect mentioned above, the algorithm for solving the reconstruction equation is the conjugate gradient descent algorithm, and the solution object is the uncollected k-space data Xu.

[0025] As a preferred embodiment of the first aspect above, in step S4, the coil sensitivity map is further calculated based on the k-space convolution kernel.

[0026] In a second aspect, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the image reconstruction method based on subspace parallel imaging as described in any of the solutions of the first aspect above.

[0027] Thirdly, the present invention provides a computer electronic device, which includes a memory and a processor;

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

[0029] The processor is configured to, when executing the computer program, implement the image reconstruction method based on subspace parallel imaging as described in any of the first aspects above.

[0030] Compared with the prior art, the present invention has the following beneficial effects.

[0031] This invention can transform sensitivity maps into k-space convolution kernels and then reconstruct undersampled data, greatly improving tolerance to sensitivity map defects and exhibiting higher robustness than the traditional SENSE method. Of course, for data with ACS but no sensitivity map, this invention can also directly generate k-space convolution kernels from the ACS; moreover, compared to existing technologies, this invention has higher computational efficiency and shorter computation time in solving for k-space convolution kernels. Attached Figure Description

[0032] Figure 1 is a schematic diagram of the SPAN method. (A) The sensitivity map is Fourier transformed into k-space data. (B) A calibration matrix is ​​constructed from the generated ACS data, and its subspace is obtained by SVD. (C) A k-space self-generating convolution kernel G and a zero-space convolution kernel N are generated from the signal subspace or the null space basis vectors.

[0033] Figure 2 shows the 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. Sensitivity map estimates are obtained by voxel-by-voxel eigenvalue decomposition of the image domain maps. (B): For... The eigenvectors corresponding to the eigenvalue = "1" form a sensitivity map (red box). (D): For The eigenvectors corresponding to the eigenvalues ​​= "0" form a sensitivity map (red box).

[0034] Figure 3 shows a comparison of calibration times for different algorithms. The figures illustrate the time required to compute convolutions using the original SPIRIT algorithm (vertical lines), the SPIRIT algorithm based on Cholesky decomposition (horizontal lines), and SPAN (diagonal lines) on 16, 24, and 32-channel data.

[0035] Figure 4 compares the reconstruction results of sensitivity maps using different signal support region sizes. SENSE and SPAN reconstruct the image with R=4 (B) using sensitivity maps with signal support regions that are too large (C), suitable (D), and too small (E). The error map compared to the fully sampled image (A) is also shown in the figure, with nRMSE displayed at the bottom of the error map. White circles represent signal regions.

[0036] Figure 5 shows a comparison of reconstructed images of partial Fourier undersampled data. SENSE, SENSE-LORAKS, SPAN, and SPAN-LORAKS reconstructed the partial Fourier undersampled data (B) using sensitivity maps of signal support regions with values ​​of too large (C), appropriate (D), and too small (E). Error maps compared to the fully sampled image (A) are also shown in the figure, with nRMSE displayed at the bottom of each error map. White circles indicate signal regions.

[0037] Figure 6 shows a comparison of reconstruction results using sensitivity maps with holes. SENSE, SENSE-LORAKS, SPAN, and SPAN-LORAKS reconstruction results using sensitivity maps with holes (B) on partially Fourier undersampled data (D), including image results (E) and k-space results (G). Errors (F, H) compared to fully sampled reference standards (A, C) are also shown in the figure, with nRMSE shown at the bottom of the image domain error plot. White circles indicate signal regions. Detailed Implementation

[0038] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in various embodiments of the present invention can be combined accordingly without mutual conflict.

[0039] In a preferred embodiment of the present invention, an image reconstruction method based on subspace parallel imaging (hereinafter referred to as the SPAN method) is provided, which includes steps S1 to S5, and an exemplary flowchart is shown in Figure 1. The overall flowchart will be described in detail below with reference to Figure 1.

[0040] S1: Acquire k-space calibration data from the magnetic resonance scanner. The k-space calibration data can be acquired in two ways: For additionally acquired sensitivity maps, as shown in Figure 1(A), this invention performs a Fourier transform on the sensitivity map to obtain k-space data, and selects the central portion of the k-space data as the ACS to be acquired in this step; if the acquired data contains an automatic calibration signal, it can be directly used as the ACS to be acquired in this step.

[0041] S2: In the k-space calibration data, a sliding window is used to extract data point by point. The data blocks corresponding to the sliding window are converted into one-dimensional vectors and arranged in order to form a structured block Hankel calibration matrix A.

[0042] For ease of description, the number of channels in the k-space calibration data is denoted as L, and these L channels are sequentially labeled Ch1 to ChL; the size of the sliding window is denoted as M×M. Therefore, the sliding window can slide across the k-space calibration data from top to bottom and from left to right, with a sliding step of one k-space point. As shown in Figure 1(B), the size of each data block extracted by the sliding window is M×M×L.

[0043] In an embodiment of the present invention, continuing to refer to Figure 1(B), the method for constructing the Hankel calibration matrix A based on all extracted data blocks is as follows: For each M×M×L data block, vectorize it, first arranging the data of each channel into a matrix with M elements. 2 A column vector of size ×1 is created, and then the data from different channels are concatenated along the column direction to form a vector M with the same total number of elements as the data block. 2 An L×1 one-dimensional column vector b, then M 2 An L×1 one-dimensional column vector b is transformed into a 1×M vector using its conjugate transpose. 2 A one-dimensional row vector b of LH Finally, the one-dimensional row vector b corresponding to all data blocks is... H Fill the matrix row by row in order to obtain K×M 2 L is the Hankel calibration matrix A, where K represents the total number of data blocks extracted by the sliding window from the k-space calibration data ACS.

[0044] S3: Perform singular value decomposition (SVD) on the Hankel calibration matrix. Based on the magnitude of the singular values ​​obtained from the decomposition and the truncation threshold, divide the singular vector matrix into signal subspace basis vectors and zeroing subspace basis vectors.

[0045] In embodiments of the present invention, the Hankel matrix typically possesses good low-rank properties, therefore, SVD can be performed on A to partition the subspace: A = USV H [1]

[0046] SVD yields two singular vector matrices: a left singular vector matrix U and a right singular vector matrix V. Each singular value vector in the singular vector matrix corresponds to a singular value. The subspace basis vectors can be partitioned based on the magnitude of the singular values ​​obtained from the decomposition and a set truncation threshold. In this invention, the right singular vector matrix V is preferably used for partitioning the subspace basis vectors. Specifically, a truncation threshold can be set based on the magnitude of the singular values. The singular vectors corresponding to all singular values ​​in the singular vector matrix that are greater than or equal to the truncation threshold form the signal subspace basis vector V. || The singular vectors corresponding to all 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 actual conditions. In the embodiments of the present invention, the truncation threshold is preferably the maximum singular value s after SVD. max 2%, or 0.02s max Therefore, referring to Figure 1(B), this step is equivalent to truncating the threshold by 0.02s. max The singular vector matrix V is divided into signal subspace basis vectors V| || The basis vectors of the zero subspace V are set. ⊥ .

[0047] S4: Construct a k-space convolution kernel KN based on subspace basis vectors, where the selected subspace basis vectors V * It can be the basis vector V of the signal subspace. || It can also be the zeroed subspace basis vector V ⊥ If the basis vectors of the subspace V * Choose the signal subspace basis vector V || Construct a self-generating k-space convolution kernel G, if the subspace basis vectors V * Choose the zeroed subspace basis vector V⊥ Then, a zero-set k-space convolution kernel N is constructed. The basic principle of constructing a k-space convolution kernel KN based on two types of subspace basis vectors is similar. Using subspace basis vector V... * To represent V || and V ⊥ Construct k-space convolution kernel KN(V) || and V ⊥ The methods corresponding to G and N respectively are as follows:

[0048] S41. Based on the selected subspace basis vector V * The resulting size is M 2 L×M 2 The projection matrix V of L * V * H .

[0049] S42, Project the matrix V * V * H According to M 2 Divide the data into blocks, with each block consisting of a row and a unit of M. 2 The rows are treated as a separate block matrix; for each block matrix, each M is... 2 A row vector of length L is reorganized into a weight block of size M×M×L. Each weight block is then filled into a convolutional kernel of size (2M-1)×(2M-1)×L, where all original values ​​are 0, forming an intermediate convolutional kernel (the filled region is replaced by the weight block, and the remaining positions remain 0). 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 ], M 2 The sum of the intermediate convolutional kernels divided 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 * V is the basis vector of the signal subspace. || At that time, the obtained k-space convolution kernel KN is a self-generated k-space convolution kernel G, and the subspace basis vector V * Let V be the basis vector of the null subspace. ⊥ When the k-space convolution kernel KN is obtained, it is the zeroed k-space convolution kernel N.

[0050] In an embodiment of the present invention, referring to Figure 1(C), the signal subspace basis vector V || For example, any k-space data block b satisfies: V || V ||H b = b [2]

[0051] Where b is a vectorized multi-channel k-space data block.

[0052] V || V || H All points in a data block can be generated from the k-space data block itself, but V || V || H It is not a convolution kernel itself. In order to extract from V... || V || H Generate a convolution kernel, and consider a V corresponding to an M×M sliding window. || V || H (M is 3 in Figure 1). When this V || V || H As the sliding motion continues throughout the k-space, each point in the k-space will be reconstructed as M. 2 This is because each point will have M points in an M×M sliding window. 2 Different positions. To illustrate, let's take point r in channel 1 as an example. First, point r is located at the top left corner of the M×M sliding window. Note that V || V || H The inner product of the first row weights and b generates the first data point in b, i.e., point r; and this inner product operation is equivalent to taking V... || V || H The first row of weights is arranged in an M×M×L kernel, which is then subjected to a compound convolution operation with the M×M×L data block represented by b. To conform to the property of the convolution kernel—that the generated points are centered within the block—V... || V || H The first row of weights, forming an M×M×L kernel, combines with the surrounding zeros to form a convolution kernel of size (2M-1)×(2M-1)×L (5×5×L in Figure 1), where point r is located at the center of (2M-1)×(2M-1). This process continues for V. || V || H The weight of the first row to the Mth row 2 Row weights can generate M 2 There are (2M-1)×(2M-1)×L intermediate state convolution kernels. From a least-squares perspective, the optimal convolution kernel for generating channel 1 is this M... 2 The average of the convolutional kernels, i.e., the value at position M 2 The average of each weight value. Then, for each M... 2Row weights can be used to obtain the convolution kernel that generates one channel of data. In the same way, the convolution kernels of all L channels can be obtained. The convolution kernels of all channels constitute a self-generating k-space convolution kernel G, with a size of (2M-1)×(2M-1)×L×L (5×5×L×L in Figure 1).

[0053] On the other hand, with the signal subspace basis vector V || Similarly, in constructing a basis vector V based on a null subspace... ⊥ When setting the zero-space convolution kernel N, V ​​can be used. ⊥ V ⊥ H Replace V || V || H As input to the construction process in Figure 1(C). At this time, due to V ⊥ V ⊥ H With b=0, a zero-space convolution kernel N can be generated similarly.

[0054] S5: Construct a reconstruction equation using the above k-space convolution kernel KN (self-generated k-space convolution kernel G or zeroed k-space convolution kernel N), and reconstruct the undersampled k-space data by solving the equation.

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

[0056] in and g ji It is the weight corresponding to the j-th channel when generating the k-space data of the i-th channel, x i This represents the k-space data of the i-th channel, where L is the number of channels. (Symbol) Representative compound convolution operation:

[0057] In the formula: * represents the standard convolution operation.

[0058] For a zero-spaced k-space convolution kernel N, we also have:

[0059] Therefore, formulas [3] and [5] can be used as models for reconstructing the complete k-space, and they are equivalent; either formula can be used to reconstruct the complete k-space. The fully sampled complete multi-channel k-space data X consists of the uncollected k-space data Xu and the collected k-space data Xa, i.e., X = Xu + Xa.

[0060] Taking formula [3] as an example, this invention transforms the reconstruction equation corresponding to formula [3] into a solution to the following optimization problem to reconstruct undersampled data:

[0061] Taking formula [5] as an example, this invention transforms the reconstruction equation corresponding to formula [5] into a solution of the following optimization problem to reconstruct undersampled data:

[0062] The optimization problem of the above reconstruction equation can be solved using the conjugate gradient descent algorithm, with the uncollected k-space data Xu as the solution object.

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

[0064] In addition to reconstructing the complete fully sampled k-space data X, the SPAN method of this 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 below:

[0065] For a self-generated k-space convolution kernel G, it can be first padded with zeros to the image size, and then the self-generated image domain map can be obtained through inverse Fourier transform. As shown in Figure 2(A), each voxel in the figure corresponds to an L×L correlation matrix. Therefore, from the self-generated image domain map The correlation matrix for each voxel can be obtained from this. For each voxel After performing eigenvalue decomposition, theoretically, the eigenvectors corresponding to eigenvalue 1 can form a set of coil sensitivity maps, as shown in Figure 2(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.

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

[0067] 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.

[0068] Therefore, 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 embodiments, which includes a memory and a processor;

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

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

[0071] Furthermore, the logical instructions in the aforementioned memory can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention.

[0072] Therefore, based on the same inventive concept, the present invention provides a computer-readable storage medium corresponding to the image reconstruction method based on subspace parallel imaging. The storage medium stores a computer program, which, when executed by a processor, can realize the image reconstruction method based on subspace parallel imaging as described above.

[0073] Therefore, based on the same inventive concept, the present invention provides a computer program product, including a computer program / instruction, which, when executed by a processor, can realize the image reconstruction method based on subspace parallel imaging as described above.

[0074] Specifically, in the above three embodiments, the computer program is executed by the processor and can perform the aforementioned steps S1 to S5.

[0075] It is understood that the aforementioned storage media may include random access memory (RAM) or non-volatile memory (NVM), such as at least one disk storage device. Furthermore, the storage media may also be various media capable of storing program code, such as USB flash drives, external hard drives, magnetic disks, or optical discs.

[0076] It is understood that the processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0077] It should also be noted that those skilled in the art will understand that, for the sake of convenience and brevity, the specific working process of the system described above can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here. In the embodiments provided in this application, the division of steps or modules in the system and method is merely a logical functional division, and there may be other division methods in actual implementation. For example, multiple modules or steps may be combined or integrated together, and a module or step may also be split.

[0078] The following example will further demonstrate the technical effects that the SPAN method described in S1 to S5 of the present invention can achieve, so that those skilled in the art can better understand the essence of the present invention.

[0079] Example

[0080] 1. MRI experiment

[0081] The human trial was conducted on a 3T Siemens MRI scanner (MAGNETOM Prisma) using a 64-channel head receiver coil. The human trial was approved by the local institutional review board. The sequence used was a TSE imaging sequence. The sequence parameters were as follows: TR = 3000 ms, TE = 7.5 ms, FOV = 212 × 212 mm. 2 Resolution = 1.1 × 1.1 mm 2Layer thickness = 5mm, acquisition matrix = 192×192, echo train length = 24.

[0082] A reference scan was performed using the GRE sequence to obtain the sensitivity map required for the standard SENSE. The acquisition parameters were 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 Data collection time = 51s.

[0083] 2. Image Reconstruction and Analysis

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

[0085] Since formulas [3] and [5] are completely equivalent, this embodiment will reconstruct the image based on formula [5]. Specifically, in this embodiment, when reconstructing the uncollected k-space data Xu using the SPAN method corresponding to steps S1 to S5 above, the following optimization problem is solved:

[0086] Here, Xa represents the collected data, and Xu + Xa = X. This optimization problem can be solved using the conjugate gradient descent algorithm.

[0087] Additionally, when calibration information comes from additionally acquired sensitivity data, SPAN can incorporate S-type LORAKS constraints to reconstruct data with partial Fourier downsampling: strank(P(Xu+Xa))≤μ [9]

[0088] Where P is the LORAKS operator, which transforms k-space data into the S matrix in the LORAKS model; μ is the expected rank parameter of the LORAKS matrix. This embodiment uses an iterative max-min algorithm to solve this optimization problem.

[0089] Furthermore, to analyze the kernel computation efficiency of the SPAN method of this invention, it is compared with the SPIRiT method. The 24 central k-space lines were selected as the ACS from the full sampled data. For the SPAN method, the sliding window size for generating the Hankel calibration matrix was set to 7×7, and the singular value cutoff was 2% of the maximum singular value. The tolerance of the conjugate gradient descent algorithm was set to 1e-3. For the SPIRiT method, the same sliding window size and tolerance were used. Based on previous research, two methods were used to compute the SPIRiT kernel, one with a computational complexity of O(L...). 4 The original calibration algorithm has two forms: one with a computational complexity of O(L) and the other with a computational complexity of O(L).3 A calibration algorithm based on Cholesky decomposition was proposed. The k-space data was compressed into 16, 24, and 32 channels using a channel compression algorithm to compare the calibration times of different methods.

[0090] To analyze the robustness of the SPAN method to sensitivity maps, retrospective uniform downsampling without preserving ACS was performed on the fully sampled k-space data. Sensitivity maps used for SPAN and SENSE reconstruction were calculated from the SENSE reference scan using locally weighted multinomial regression, registration, and interpolation. To test the robustness of SPAN to the signal support region, sensitivity maps with three different signal support region sizes were generated, and then undersampled data with an acceleration factor R = 4 was reconstructed using SENSE and SPAN, respectively. Additionally, partial Fourier downsampling experiments were performed on sensitivity maps with three different signal support region sizes and a sensitivity map with holes, reconstructing the undersampled data using SENSE, SENSE-LORAKS, SPAN, and SPAN-LORAKS, with a total acceleration factor R = 4.5 and a partial Fourier factor PF = 2 / 3. μ in SPAN-LORAKS was set to the number of basis vectors in the signal subspace.

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

[0092] 3. Results Analysis

[0093] Figure 2 shows the coil sensitivity plot calculated from the k-space convolution kernel using compressed four-channel data. The image domain plot is obtained by applying zero-padding and inverse Fourier transform to the self-generated k-space convolution kernel G. (A) The dimensions are 192×192×4×4. For ease of display, [the dimensions will be...]. The various 192×192 plots are displayed sequentially, arranged in a 4×4 pattern; therefore, each voxel has a 4×4 correlation matrix. After voxel-wise eigenvalue decomposition (B), the eigenvectors corresponding to eigenvalues ​​= "1" form a sensitivity map (red box). Similarly, the image domain map obtained by applying zero-padding and inverse Fourier transform to the zeroed k-space convolution kernel N is obtained. (C). Regarding After performing voxel-by-voxel eigenvalue decomposition (D), the eigenvectors corresponding to eigenvalues ​​= "0" form a sensitivity map (red box).

[0094] Figure 3 shows the computation time of convolution kernels for the original SPIRIT (vertical line), the Cholesky-based SPIRIT algorithm (horizontal line), and SPAN (diagonal line) on 16, 24, and 32-channel data. Consistent with previous studies, the Cholesky-based SPIRIT algorithm has a shorter computation time than the original SPIRIT algorithm. Notably, although the calibration computation complexity of the SPAN method and the Cholesky-based SPIRIT algorithm is of the same order of magnitude, the computation time of SPAN is significantly shorter than that of the Cholesky-based SPIRIT algorithm.

[0095] 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 completely cover the image signal. Clearly, regardless of 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 (containing the image signal as much as possible without noise) can the SENSE method produce relatively good reconstruction results (D). However, the SPAN method shows stable results in all cases and has smaller errors.

[0096] Figure 5 shows that the SPAN method can be well combined with LORAKS constraints to reconstruct data with partial Fourier downsampling. When the signal support region is too large or too small, the results produced by the SENSE and SENSE-LORAKS methods have obvious artifacts (C, E); when the signal support region is appropriate, SENSE-LORAKS shows better results (D). In contrast, the SPAN method has a smaller error than SENSE in all cases; the SPAN-LORAKS method produces stable results in all cases and has a smaller error than SENSE-LORAKS in the optimal case.

[0097] Figure 6 shows the reconstruction results of SENSE(-LORAKS) and SPAN(-LORAKS) when a hole is introduced into the sensitivity map. SENSE and SENSE-LORAKS fail to reconstruct the signal value at the hole, and the signal leaks to other locations. 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 of part of the Fourier region completely, resulting in a clearer image with smaller errors.

[0098] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. An image reconstruction method based on subspace parallel imaging, used to reconstruct undersampled k-space data acquired by a magnetic resonance scanner, characterized in that, The method includes the following steps: S1: Acquire k-space calibration data from the magnetic resonance scanner; S2: In the k-space calibration data, a sliding window is used to extract data point by point. The data blocks corresponding to the sliding window are converted into one-dimensional vectors and arranged in order to form the Hankel calibration matrix. S3: Perform singular value decomposition on the Hankel calibration matrix, and divide the singular vector matrix into signal subspace basis vectors and zeroing subspace basis vectors according to the magnitude of the singular values ​​obtained by decomposition and the truncation threshold. S4: Construct k-space convolution kernels based on subspace basis vectors, where: when the subspace basis vectors are selected from the signal subspace basis vectors, the corresponding k-space convolution kernels constructed are self-generated k-space convolution kernels; when the subspace basis vectors are selected from the zeroed subspace basis vectors, the corresponding k-space convolution kernels constructed are zeroed k-space convolution kernels. S5: Construct a reconstruction equation using the k-space convolution kernel, and reconstruct the undersampled k-space data by solving the equation.

2. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In step S1, the k-space calibration data is an automatic calibration signal directly acquired by a magnetic resonance scanner, or the center part of the magnetic resonance scanner's coil sensitivity map is obtained after Fourier transform.

3. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In S2, the method for constructing the Hankel calibration matrix is ​​as follows: the data block of each sliding window is vectorized, and the data of different channels are concatenated in the column direction to form a one-dimensional column vector with the same total number of elements as the data block. Then, the one-dimensional column vector is converted into a one-dimensional row vector through conjugate transpose. Finally, the one-dimensional row vectors corresponding to all data blocks are filled into the matrix row by row in order to obtain the Hankel calibration matrix.

4. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In S3, the singular vector matrix is ​​a right singular vector matrix.

5. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In S3, the threshold for singular values ​​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 basis vectors of the signal subspace, and the singular vectors corresponding to all singular values ​​less than the truncation threshold form the basis vectors of the zeroed subspace.

6. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In S4, based on the selected subspace basis vector V * The method for constructing the k-space convolution kernel is as follows: S41. Based on the selected subspace basis vector V * The resulting 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; S42, Project the matrix V * V * H According to M 2 Divide the data into blocks, with each block consisting of a row and a unit of M. 2 The row is treated as a separate block matrix; For each block matrix, take each M... 2 A row vector of length L is reshaped into a weight block of size M×M×L. Each weight block is then filled into a convolutional kernel of size (2M-1)×(2M-1)×L with all initial values ​​being 0, forming 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 ], M 2 The sum of the intermediate convolutional kernels divided by M 2 This yields 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; wherein the subspace basis vector V * V is the basis vector of the signal subspace. || At that time, the obtained k-space convolution kernel KN is a self-generated k-space convolution kernel G, and the subspace basis vector V * Let V be the basis vector of the null subspace. ⊥ When the k-space convolution kernel KN is obtained, it is the zeroed k-space convolution kernel N.

7. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In step S5, when the k-space convolution kernel is a self-generated k-space convolution kernel G, the constructed reconstruction equation is: When the k-space convolution kernel is a zero-set 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 uncollected k-space data Xu and the collected k-space data Xa, i.e., X = Xu + Xa; This is a compound convolution operation, which means performing k-space convolution on each channel and then adding them together; Preferably, the algorithm for solving the reconstruction equation is the conjugate gradient descent algorithm, and the solution object is the uncollected k-space data Xu.

8. The image reconstruction method based on subspace parallel imaging as described in claim 1, characterized in that, In step S4, the coil sensitivity map is further calculated based on the k-space convolution kernel.

9. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, which, when executed by a processor, implements the image reconstruction method based on subspace parallel imaging as described in any one of claims 1 to 8.

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

Images