A Fast Online Dynamic Magnetic Resonance Reconstruction Method Based on Subspace Tracking

Abstract

This invention discloses a rapid online dynamic magnetic resonance imaging (MRI) reconstruction method based on subspace tracking, relating to the field of medical image processing technology. The method includes: Step 1, obtaining a coil sensitivity map from pre-scanned full-sample data and constructing an encoding operator using an undersampling template; Step 2, decomposing the image into two components and sequentially solving for the component values ​​corresponding to each frame, thereby rapidly acquiring the image of the current frame. This method can be applied to interventional procedures guided by MRI images. It achieves rapid scanning through undersampling and rapid reconstruction of undersampling data using the proposed online reconstruction algorithm. Physicians can promptly determine the position of catheters, guidewires, or puncture needles during intervention based on the rapidly reconstructed MRI image.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, specifically to a fast online dynamic magnetic resonance reconstruction method based on subspace tracking. Background Technology

[0002] To achieve the goal of rapidly acquiring the current magnetic resonance imaging frame, it is necessary to simultaneously achieve rapid scanning to acquire k-space data and rapid reconstruction of k-space data to obtain the image. In the rapid scanning stage, compressed sensing theory is typically used to overcome the limitations of the Nyquist sampling theorem, performing incoherent undersampling of the data. Because the amount of data acquired is reduced, the scanning time is accelerated. However, the undersampling strategy can introduce aliasing artifacts into the image. Therefore, it is necessary to utilize prior information about the image, such as its temporal or spatial sparseness or transform sparseness, to formulate corresponding optimization equations to solve for the image and remove undersampling artifacts to obtain a high-quality image. A mainstream existing method is the low-rank plus sparse decomposition method, i.e., the L+S method. It utilizes the similarity of the background between image frames and the sparsity of dynamic information, dividing the spatiotemporal matrix composed of all frames into a low-rank component and a sparse component for separate solutions. First, the aliased images of all frames are vectorized, and then concatenated into a two-dimensional spatiotemporal matrix. Then, singular value decomposition is performed, and a threshold is set to filter out smaller singular values ​​to obtain the low-rank component L. Next, a soft thresholding method is used to obtain the sparse components. Finally, the residual is subtracted through a data consistency operation to obtain the image matrix M for the current iteration. k This iterative process is repeated until convergence is achieved, yielding a reconstructed image. While this method guarantees reconstruction quality, the iterative process results in a long reconstruction time. Furthermore, it's an offline reconstruction method, requiring all k-space data to be acquired before reconstruction can begin. This clearly doesn't meet the needs of online reconstruction, which is a process of scanning and rebuilding simultaneously, relying solely on information from previous frames. Moreover, the image quality of low-rank and sparse methods is highly dependent on the number of frames; with fewer frames, image quality is poorer, making them even less suitable for online scenarios.

[0003] Subsequently, Ongie et al., based on subspace tracking theory, applied the Grassmannian Rank-One Update Subspace Estimation (GROUSE) algorithm to dynamic magnetic resonance reconstruction and introduced the RUFFed GROUSE algorithm. They first proposed an online batch processing model, processing a small batch of data at a time. Each batch of images was modeled as a sum of a low-rank component and a sparse component, where the low-rank component is the product of a basis matrix and a weight vector. The sparse component S was then solved first. j Then, the least squares method is used to update the weight vector w corresponding to each image in each batch.t Obtain the weight matrix W j Then, the basis matrix U for the next batch is solved using the Grassmann manifold descent algorithm. j For the initial value U0 of the basis matrix, the common low-frequency component of all k-space data is selected as input, and the Grouse algorithm is used to solve it. Although this method reduces the computational space complexity, in some online reconstruction processes, we cannot obtain the k-space data of all frames during the calculation process, and the estimation of the initial value is crucial to the image reconstruction quality. Due to the initial value estimation, this method also has drawbacks when applied to online reconstruction.

[0004] Recently, Babu et al. also proposed an online subspace tracking method, namely the online altGDmin method, the algorithm flowchart of which is shown below. Figure 3 As shown. Assume each magnetic resonance imaging frame can be decomposed into a mean component z1 and a low-rank component U. k b k and an unstructured residual signal component e k During online reconstruction, the first batch of k-space data is selected to estimate z1 and the basis matrix U. k Then keep z1 and U k The value remains unchanged, and b is updated every frame. k and e k While this approach aligns with the pattern of online reconstruction, the mean component z1 and the basis matrix U... k The estimation relies solely on data from the first few frames. This online method can result in inconsistent image quality as motion changes occur because the fixed mean component cannot be updated and adapted to variations in motion. Summary of the Invention

[0005] The purpose of this invention is to provide a fast online dynamic magnetic resonance reconstruction method based on subspace tracking. After obtaining the undersampled k-space data of the current frame using magnetic resonance machine scanning, the subspace tracking algorithm is used to quickly reconstruct the k-space data. This achieves low latency in acquisition and reconstruction, and the image of the current frame can be obtained quickly.

[0006] The technical problem solved by this invention is: how to reconstruct magnetic resonance images online while ensuring the reconstruction quality with low latency.

[0007] This invention can be achieved through the following technical solution: a fast online dynamic magnetic resonance reconstruction method based on subspace tracking, comprising the following steps:

[0008] Step 1: Obtain the encoding operator E t

[0009] k-space data is acquired using multi-coil acquisition, and a relationship d between each frame of k-space data and the corresponding image is constructed. t =E t x t +ε t If the number of coils is set to c, then d t ∈￡ mc×1 Let x represent the k-space data at time t. t ∈￡ n×1 ε represents the vectorized image at the corresponding time point. t This represents the noise at the corresponding moment.

[0010] Furthermore, in the case of multi-coil undersampling, the encoding operator E in the relation is... t The acquisition process involves scanning the full sampled data before online reconstruction, obtaining the coil sensitivity map using the reconstructed time-averaged image, and acquiring the encoding operator E based on the undersampled template for each frame. t ;

[0011] Step 2: Establishing the Online Model

[0012] x in step one t It is decomposed into two components, namely the subspace component s. t and error component e t As shown in the following formula: x t =s t +e t =u t w t +e t ;where u t Let w denote an orthogonal basis for a subspace. t This represents the weight value; and u is solved sequentially for each frame. t w t and e t This value allows for the rapid acquisition of the image of the current frame.

[0013] A further improvement of this invention lies in: solving for u in each frame t w t and e t The process of valuing includes:

[0014] (1) Use the L+S algorithm to accurately reconstruct the data of the first five frames and obtain the initial value u0 of the orthogonal basis of the subspace;

[0015] (2) On the Grassmann manifold based on u t-1 Use the Grouse algorithm to obtain the current orthogonal basis u. t Then, the weight value w is solved using the least squares algorithm. t This completes the subspace component s t Solving for;

[0016] (3) Utilizing subspace components s t The solution result for the error component e t Optimization is performed by updating the error component e of the current frame using the CGLS algorithm. t ;

[0017] (5) Based on the currently obtained u t and e t Repeat steps (3) and (4) to continue the reconstruction of the next frame.

[0018] A further improvement of the present invention is that the step of obtaining the initial value u0 of the orthogonal basis of the subspace includes:

[0019] S1: Before performing online reconstruction, scan and acquire five frames of undersampled K-space data.

[0020] S2: Reconstruct the first five frames of data using the L+S algorithm to obtain an image with undersampled artifacts removed.

[0021] S3: Vectorize the five reconstructed images and form a two-dimensional spatiotemporal matrix M. Perform singular value decomposition on M and use the obtained left singular vector as the initial value u0 of the estimated orthogonal basis.

[0022] A further improvement of the present invention is that the update of the GROUSE algorithm at each time step is based on the natural incremental gradient descent method. When solving the optimization equation of the subspace components by the GROUSE algorithm, ut and wt are updated alternately to obtain the subspace components of the current frame.

[0023] A further improvement of the present invention lies in: assuming u t Zhang Cheng's subspace lies within the Grassmann manifold of the one-dimensional subspace, and the optimization equations for the subspace components are as follows: stspan(u t )∈G(n,1).

[0024] A further improvement of the present invention is that: based on the solved subspace components s t and relation d t =E t x t +ε t The error component e is obtained. t Optimization equation The error components are set as all-zero matrices, and the CGLS algorithm is used to update the error components e in each frame. t .

[0025] A further improvement of the present invention lies in: in solving the error component e tA warm-start strategy is introduced, which uses the error component value obtained in the previous frame as the initial value for the current calculation, thereby improving image quality.

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

[0027] 1. This method models a single frame of image, which is more suitable for online reconstruction scenarios compared to the L+S algorithm that models based on all frames. The algorithm uses the Grouse algorithm to calculate subspace components and obtains error components through CGLS. Compared to other dMRI reconstruction algorithms, this method achieves higher image reconstruction quality.

[0028] 2. This invention eliminates the need for iterative solving when performing the Grouse algorithm, thus significantly reducing image reconstruction time compared to the iterative L+S algorithm. For subspace initialization estimation, the L+S algorithm is used to estimate the initial values ​​of the subspace orthogonal basis. Experiments show that the L+S method can obtain good initial value estimates using only information from the first few frames, thereby improving image quality, especially in the first few frames. This method can achieve low reconstruction time while maintaining image quality, and reconstruction is based on a single frame image. It can reconstruct the current frame image immediately after scanning a frame of undersampled k-space data, making it applicable to some MRI-guided interventional surgeries. Attached Figure Description

[0029] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings.

[0030] Figure 1 This is a schematic diagram of the method execution flow of the present invention;

[0031] Figure 2 This is a schematic diagram of the online reconstruction mode of the present invention;

[0032] Figure 3 This is a schematic diagram of the existing altGDmin online reconstruction process. Detailed Implementation

[0033] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided.

[0034] Please see Figure 1-2 As shown, a fast online dynamic magnetic resonance reconstruction method based on subspace tracking includes the following steps:

[0035] Step 1: Obtain the encoding operator E t

[0036] In dynamic magnetic resonance imaging, k-space data is acquired using multiple coils. The relationship between each frame of k-space data and the corresponding image is as follows:

[0037] d t =E t x t +ε t (1)

[0038] Let the number of coils be c, then d t ∈￡ mc×1 Let x represent the k-space data at time t. t ∈￡ n×1 ε represents the vectorized image at the corresponding time point. t This represents the noise at the corresponding moment.

[0039] In the case of multi-coil undersampling, E t Represents the encoding operator, and Where S c ∈￡ n×n It is a diagonal matrix, and its diagonal elements are the values ​​of the sensitivity spectrum of the c-th coil; F∈£ n×n It is the Fourier transform operator, θ t This is the undersampling matrix.

[0040] Before performing online reconstruction, the full sampled data is scanned, and the coil sensitivity map is obtained using the reconstructed time-averaged image. Furthermore, the encoding operator E is obtained based on the undersampled template for each frame. t Among them, the undersampling template is randomly generated before the formal scanning begins, and the undersampling method of k-space data is determined based on the undersampling template.

[0041] Step 2: Establishing the Online Model

[0042] x in step one t It is decomposed into two components, namely the subspace component s. t and error component e t As shown in the following formula: x t =s t +e t =u t w t +e t (2)

[0043] Wherein, subspace component s t ∈￡ n×1 This represents the main information of the t-th frame image, with the error component e. t ∈￡ n×1 Represents detailed information, and the error component e t This can be viewed as a fine-tuning based on the subspace tracking and reconstruction results. t ∈￡n×1 For an orthogonal basis of a subspace, w t ∈£ represents the weight value. Therefore, we only need to calculate u for each frame. t w t and e t This means that the image of the current frame can be obtained.

[0044] Solve for u in each frame t w t and e t The specific process is as follows:

[0045] (1) One-dimensional variable density Cartesian sampling method is used to collect k-space data, and the data of the first five frames are reconstructed by L+S (Low Rank Plus Sparsity). The initial value u0 is obtained by singular value decomposition of the final result. It should be noted that in the improved initial value strategy, the L+S algorithm is used to accurately reconstruct the data of the first five frames, so that the initial value u0 of the orthogonal basis of the obtained subspace is more accurate, so as to further improve the reconstruction quality of the subsequent frames.

[0046] Specifically, the initial values ​​of the orthogonal basis play a crucial role in the reconstruction quality. Since it is an online reconstruction, only limited data is available for initial value estimation. Therefore, based on the limited data, the initial values ​​of the subspace components are obtained through accurate reconstruction using the static RPCA method, i.e., the L+S algorithm, as shown in the following equation:

[0047]

[0048] First, input the spatiotemporal matrix consisting of the first 5 frames of the acquired images, i.e., M∈£. n×5 The L+S algorithm is used to solve the reconstruction equation:

[0049] The specific steps are as follows: For equation (3), solve it using the iterative soft threshold method.

[0050] S1: First, define the soft threshold operator as Λ λ (x) = x / |x|gmax(|x|-λ,0), where λ is a real threshold;

[0051] S2: The solution begins with an image M containing undersampling artifacts as input. First, perform singular value decomposition (SVD) on the image to obtain the left singular vector U and the right singular vector V. H And the singular value matrix Σ, i.e., M = UΣV H The diagonal elements of the singular value matrix Σ are the singular values ​​arranged in descending order. After soft thresholding, only the larger singular values ​​are retained in the diagonal elements. The updated singular value matrix is ​​then multiplied by U and V on the left and right, respectively.H Obtaining the low-rank matrix L, this step is defined as the Singular Value Thresholding (SVT) operation. λ (M)=UΛ λ (Σ)V H .

[0052] S3: After obtaining matrix L, perform a sparse transformation on the values ​​in ML and then perform a soft thresholding. The resulting matrix is ​​then inversely transformed to obtain a sparse matrix S.

[0053] S4: Finally, perform a data consistency operation to remove artifacts, resulting in M ​​= L + SE. * (E(L+S)-d). Based on this value of M, repeat the above steps to iteratively update L, S, and M until the error between two iterations is less than a set value, at which point the process stops.

[0054] S5: After obtaining M, perform singular value decomposition on it, and select the left singular vector corresponding to the largest singular value as the estimated value of u0.

[0055] (2) On the Grassmann manifold based on u t-1 Use the Grouse algorithm to obtain the current orthogonal basis u. t Then, the weight value w is solved using the least squares algorithm. t This completes the subspace component s t Finally, the error component e is solved using the Conjugate Gradient Least Squares (CGLS) method. t This allows for the estimation of the current frame;

[0056] Specifically, first, subspace components s are processed. t Solving for:

[0057] During optimization, the subspace components are first solved based on the current k-space data, as shown in equation (3). Since u t Spanning a one-dimensional subspace, and given that the Grassman manifold is a set of subspaces of a given dimension, therefore assuming u t If Zhang Cheng's subspace lies within the Grassmann manifold of the one-dimensional subspace, then the optimization equation can be obtained as follows:

[0058]

[0059] Solving the above equation using the GROUSE algorithm alternately applies u. t and w t Update the solution to the subspace components of the current frame;

[0060] First, we fix the u from the previous frame.t-1 At this point, w can be solved using equation (4):

[0061]

[0062] The value of w is obtained using the least squares algorithm: w = (E t u) -1 d t .

[0063] Gradient descent is performed on the Grassmann manifold to solve for the basis matrix u. t :

[0064] First, define the cost function.

[0065] The GROUSE algorithm updates at each time step based on the natural incremental gradient descent method. First, it calculates the cost function at time t with respect to u. t Partial derivatives:

[0066] in The obtained derivative multiplied by the projection matrix Q = I - uu T Project the gradient onto the Grassmann manifold:

[0067] in right Perform singular value decomposition.

[0068] Let p = uw, which is the predicted value, and the geodesic along the Grassmann manifold and The singular values ​​and singular vectors are related. At this point, the step size is set to η, and u is updated along the geodesic to obtain the value at this moment. t :

[0069]

[0070] Choose the step size η in each iteration t satisfy In obtaining u t Then, using the formula w = (E t u) -1 d t Get the current weight value w t Therefore, it can be concluded that in solving w t and u t In this case, we only need to first obtain w using the least squares method, and then calculate p. t and Perform a rank-one update on the Grassmann manifold.

[0071] Secondly, in the subspace component st After the solution is completed, the error component e is then analyzed. t Optimization is performed by subtracting the obtained subspace component st, and then the optimization equation is obtained.

[0072] First, the error components are set as an all-zero matrix. Then, CGLS is used to update the e of each frame. t Meanwhile, a warm-start strategy is introduced, which uses the error component value obtained in the previous frame as the initial value for this calculation. This method can improve image quality.

[0073] The specific steps are as follows: First, calculate the residual value r. t =d t -E t s t -E t e t-1 and search direction Then calculate the step size. The final update yields the current value e. t =e t-1 +α t p t .

[0074] (3) Based on the currently obtained u t and e t Repeat step (2) to continue the reconstruction of the next frame.

[0075] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A fast online dynamic magnetic resonance reconstruction method based on subspace tracking, characterized in that: Includes the following steps: Step 1: Obtain the encoding operator The relationship between each frame of k-space data and the corresponding image is as follows: If the number of coils is set to c, then This represents the k-space data at time t. This represents the vectorized image at the corresponding time point. This represents the noise at the corresponding moment. Before the undersampling scan begins, the encoding operators in the relation need to be configured. The acquisition process includes obtaining coil sensitivity maps and undersampled templates; specifically: scanning the full sampled data, using the reconstructed time-averaged image to obtain coil sensitivity maps, and randomly generating undersampled templates. By combining the coil sensitivity maps and the undersampled templates for each frame, the encoding operator can be determined. ; Step 2: Establishing the Online Model In step one Decomposed into two components, namely the subspace component. and error components As shown in the following formula: ;in Describes an orthogonal basis for a subspace. Represent the weight values; and solve for each frame sequentially. , and This value allows for the rapid acquisition of the image of the current frame.

2. The fast online dynamic magnetic resonance reconstruction method based on subspace tracking according to claim 1, characterized in that, Solve for each frame , and The process of valuing includes: (1) Use the L+S algorithm to accurately reconstruct the data of the first five frames and obtain the initial values ​​of the orthogonal basis of the subspace. ; (2) Based on Grassmann manifolds Use the Grouse algorithm to obtain the current orthogonal basis. Then, the weight values ​​are solved using the least squares algorithm. This completes the subspace component. Solving for; (3) Utilizing subspace components The solution results for the error components Optimization is performed by updating the error components of the current frame using the CGLS algorithm. ; (4) Based on the currently obtained and Repeat steps (2) and (3) to continue the reconstruction of the next frame.

3. The fast online dynamic magnetic resonance reconstruction method based on subspace tracking according to claim 2, characterized in that, Obtain the initial values ​​of the orthogonal basis of the subspace. The steps include: S1: Before performing online reconstruction, scan and acquire five frames of undersampled K-space data; S2: Reconstruct the first five frames of data using the L+S algorithm to obtain an image with undersampled artifacts removed; S3: Vectorize the five reconstructed images and form a two-dimensional spatiotemporal matrix M. Perform singular value decomposition on M and use the obtained left singular vectors as the initial values ​​for the estimated orthogonal basis. .

4. The fast online dynamic magnetic resonance reconstruction method based on subspace tracking according to claim 2, characterized in that, The update of the GROUSE algorithm at each time step is based on the natural incremental gradient descent method. When solving the optimization equations for the subspace components using the GROUSE algorithm, the update is alternately performed on... and Update the subspace components of the current frame.

5. The fast online dynamic magnetic resonance reconstruction method based on subspace tracking according to claim 4, characterized in that, Assumption Zhang Cheng's subspace lies within the Grassmann manifold of the one-dimensional subspace, and the optimization equations for the subspace components are as follows: .

6. The fast online dynamic magnetic resonance reconstruction method based on subspace tracking according to claim 2, characterized in that, Based on the solved subspace components and relational expressions The error components are obtained. Optimization equation The error components are set as all-zero matrices, and the CGLS algorithm is used to update the error components for each frame. .

7. The fast online dynamic magnetic resonance reconstruction method based on subspace tracking according to claim 6, characterized in that, In solving error components A warm-start strategy is introduced, which uses the error component value obtained in the previous frame as the initial value for the current calculation, thereby improving image quality.

Images