Free-breathing fast three-dimensional cardiac magnetic resonance imaging method and system

Abstract

The application discloses a free breathing fast three-dimensional heart magnetic resonance imaging method and system, and belongs to the three-dimensional imaging field, and comprises the following steps: acquiring undersampled K-space data of a heart under different breathing states; constructing a motion solving network and an image reconstruction network based on implicit neural representation; according to the sampling K-space data under different breathing states, the motion solving network and the image reconstruction network are iteratively updated and iteratively reconstructed through a loss function comprising a motion compensation data fidelity term, a motion field smoothness term, an image smoothness term and an image edge sharpness term, so as to obtain a breathing motion field estimation result and a heart nuclear magnetic resonance motion compensation reconstruction image; by means of implicit neural representation, the breathing motion field and the three-dimensional heart magnetic resonance image are directly reconstructed from the undersampled K-space data, and in combination with various regular term constraints, the imaging time is reduced, the image signal-to-noise ratio is improved, the image reconstruction quality can be guaranteed, and no additional training data is needed.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of three-dimensional cardiac magnetic resonance imaging, and particularly relates to a free-breathing fast three-dimensional cardiac magnetic resonance imaging method and system, and particularly relates to a free-breathing fast three-dimensional cardiac magnetic resonance imaging method based on implicit neural representation. BACKGROUND

[0002] Three-dimensional cardiac magnetic resonance imaging can obtain high-resolution three-dimensional images of the heart, and further realize non-invasive assessment of the whole heart anatomical structure and function, thereby providing a basis for diagnosis and treatment of heart diseases. However, compared with conventional two-dimensional imaging, three-dimensional cardiac magnetic resonance imaging faces some key challenges that hinder its wide clinical application. First, due to the slow acquisition speed of magnetic resonance and the large sampling matrix of three-dimensional imaging, three-dimensional high-resolution cardiac magnetic resonance scanning requires a long scanning time, which is challenging for patients who cannot maintain stillness for a long time. High-oversampling acquisition is required in clinical practice to shorten the scanning time. Second, respiratory motion during three-dimensional cardiac scanning often leads to significant respiratory artifacts and reduces image quality. To deal with respiratory motion in free-breathing magnetic resonance scanning, diaphragm navigation technology and self-navigation technology have been proposed. Among them, the self-navigation technology uses data collected at different respiratory states to achieve a scanning efficiency close to 100% and a predictable scanning time. Since the collected data is scattered in different motion states, image registration is generally required in the image reconstruction process to estimate the respiratory motion field between different motion states, which is then used for subsequent motion compensation image reconstruction. However, due to high-oversampling acquisition, the images at different respiratory states often contain significant oversampling artifacts, which reduces the registration accuracy and affects the motion compensation reconstruction quality.

[0003] The methods for motion compensation image reconstruction of under-sampled three-dimensional magnetic resonance imaging data containing respiratory motion proposed in the field at the present stage mainly include two categories. The first category is the traditional iterative method, including NR-SENSE, NR-PROST and the like. This category of method first solves the three-dimensional non-rigid respiratory motion field between different respiratory states by using a free deformation non-rigid registration algorithm based on the smoothness constraint of the motion field; then the solved motion field is substituted into the loss function, and the motion compensation reconstruction image is obtained by minimizing the loss function through the iterative algorithm. This kind of method performs well on low-fold under-sampled data, and as the under-sampling multiple increases, the motion estimation accuracy decreases, and the reconstructed image appears obvious artifacts and noise amplification. The second category is the motion compensation reconstruction method based on deep learning, such as MoCo-MoDL. The network model constructed by this category of method includes a motion estimation network and an image reconstruction network, which first estimates the three-dimensional non-rigid respiratory motion field from the different respiratory state images containing artifacts by using the motion estimation network, and then substitutes the non-rigid respiratory motion field into the image reconstruction network to output the motion compensation reconstruction result. This kind of method can realize higher under-sampling multiple than the traditional iterative method by training on a large amount of magnetic resonance data. However, in most actual situations, it is difficult to collect a large-scale three-dimensional cardiac magnetic resonance training data set; and when the test set and the training set are inconsistent, the trained neural network may have generalization problems, resulting in a significant decrease in reconstruction quality. In order to avoid using large-scale training data, some self-supervised image reconstruction methods are proposed, such as DIP, DIP-CS and the like. However, this kind of method can only deal with the under-sampled data reconstruction problem without motion, that is, it can only reconstruct the under-sampled K-space data under the same respiratory motion state, and cannot utilize the data collected under all respiratory states for motion compensation image reconstruction, which has low data utilization and is difficult to realize high-fold acceleration. Therefore, a scheme for high-fold under-sampling motion compensation reconstruction without training data is needed to improve the clinical usability of three-dimensional cardiac magnetic resonance imaging. SUMMARY

[0004] To solve the above technical problems, the present application provides a free breathing fast three-dimensional cardiac magnetic resonance imaging method and system to solve the problems existing in the prior art.

[0005] To achieve the above purpose, the present application provides a free breathing fast three-dimensional cardiac magnetic resonance imaging method and system, comprising:

[0006] Acquire under-sampled K-space data of the heart under different respiratory states, wherein the under-sampled K-space data is under-sampled Fourier space cardiac magnetic resonance data;

[0007] Construct a motion solving network and an image reconstruction network based on implicit neural representation;

[0008] Construct a loss function for joint motion estimation and motion compensation image reconstruction;

[0009] According to the undersampled K-space data under different respiratory states, the motion solving network and the image reconstruction network are iteratively updated and iteratively reconstructed through a loss function, to obtain a final respiratory motion field estimation result and a motion compensation reconstructed image.

[0010] The loss function includes different loss terms of data fidelity, data smoothing and image detail constraint.

[0011] Optionally, the acquisition process of the undersampled K-space data comprises:

[0012] The magnetic resonance imaging data of a free-breathing heart is acquired through a self-navigation imaging sequence; and for the magnetic resonance imaging data, the data of the center point of the K-space in the magnetic resonance imaging data within the acquisition window of each cardiac cycle is acquired as navigation data;

[0013] A respiratory curve is generated according to the navigation data, the respiratory curve is divided according to the respiratory amplitude of the respiratory curve, to obtain a respiratory phase, wherein each respiratory phase corresponds to a respiratory state from the end of inspiration to the end of expiration;

[0014] The undersampled K-space data under different respiratory states is obtained by integrating the respiratory phase with the magnetic resonance imaging data.

[0015] Optionally, the loss function includes a motion compensation data fidelity term, a motion field smoothing regular term, an image smoothing regular term and an image edge sharpness regular term; wherein the motion compensation data fidelity term is a loss term between the real undersampled K-space data and the K-space data synthesized by the motion compensation reconstructed image, the motion field smoothing regular term is a smoothing measurement loss term of the respiratory motion field, the image smoothing regular term is a smoothing measurement loss term of the motion compensation reconstructed image, and the image edge sharpness regular term is an edge sharpness loss term of the motion compensation reconstructed image.

[0016] Optionally, the loss function is a weighted sum result of the motion compensation data fidelity term, the motion field smoothing regular term, the image smoothing regular term and the image edge sharpness regular term.

[0017] Optionally, the motion solving network adopts a multi-layer perception machine, wherein the motion solving network includes 3 layers of hidden layers, the hidden neurons of each hidden layer are 32, 32 and 3, the first layer of hidden layers and the second layer of hidden layers adopt an Elu activation function, the length of the output layer is 3, and the input is the encoding information of the respiratory motion field, wherein the encoding information of the respiratory motion field is obtained by performing hash encoding on the spatial coordinates of the voxel points of the respiratory motion field and the respiratory state, and the output result is the respiratory motion field estimation result, wherein the respiratory motion field estimation result is the displacement amount in different directions of each voxel point.

[0018] Optionally, the image reconstruction network employs a multilayer perceptron, wherein the image reconstruction network includes 6 hidden layers, the input is the encoded information of the reconstructed image, which is obtained by hashing the spatial coordinates of the voxel points of the reconstructed image; the output is the real and imaginary parts of the voxel signals of the voxel points of the reconstructed image.

[0019] Optionally, the process of iteratively updating and reconstructing the motion solving network and the image reconstruction network includes:

[0020] The loss function is optimized by a semi-quadratic splitting method, and by introducing auxiliary variables, the optimization steps of the loss function for the motion solving network and the image reconstruction network are split into a first sub-problem and a second sub-problem.

[0021] For the first subproblem, the spatial coordinates of voxel points in the reconstructed image and motion field are obtained. The motion solving network and the image reconstruction network are used to reconstruct the motion field based on the spatial coordinates of voxel points in the reconstructed image and motion field, thus obtaining the breathing motion field estimation result and the motion-compensated reconstructed image for the first subproblem. The motion solving network and the image reconstruction network are then optimized based on the breathing motion field estimation result and the motion-compensated reconstructed image for the first subproblem.

[0022] For the second subproblem, the motion-compensated reconstructed image of the first subproblem is denoised and sharpened to obtain the processed image of the second subproblem; the motion solving network and the image reconstruction network are optimized based on the processed image of the second subproblem.

[0023] The optimization steps of the first and second subproblems are executed iteratively and alternately to obtain the final respiratory motion field estimation results and motion compensation reconstruction images.

[0024] Optionally, a loss function is used to optimize the motion solving network and the image reconstruction network for the first subproblem. for:

[0025]

[0026] Where z represents an auxiliary variable, This indicates the respiratory phase, or the respiratory state. Indicates the first K-space data collected from each respiratory state. Indicates the first K-space sampling mask corresponding to each breathing state Indicates Fourier transform, Indicates coil sensitivity. The first term represents the estimation of the motion solver network. The respiratory motion field corresponding to each breathing state denotes a motion compensated reconstructed image estimated by an image reconstruction network, denotes a finite difference operator, , , denotes a weight coefficient corresponding to different terms.

[0027] Optionally, the loss function for optimizing the motion solving network and the image reconstruction network for the second sub-problem is :

[0028]

[0029] wherein, denotes an image edge sharpness loss term, denotes a corresponding weight coefficient.

[0030] In another aspect, the present application provides a free breathing fast three-dimensional cardiac magnetic resonance imaging system for performing the above method.

[0031] Compared with the prior art, the present application has the following advantages and technical effects:

[0032] 1. The present application jointly solves the respiratory motion field and the motion compensated three-dimensional cardiac magnetic resonance image directly from the under-sampled K-space data by means of implicit neural representation, avoids additional reconstruction and registration of high-fold under-sampled images of each respiratory state, prevents the propagation of motion solving errors, and improves the quality of high-fold under-sampling motion compensation reconstruction;

[0033] 2. The present application combines multiple regularization constraints, improves the signal-to-noise ratio of the reconstructed image, and effectively preserves the details and image edges of the reconstructed image;

[0034] 3. No other training data is needed except for the under-sampled K-space data itself to be reconstructed, which can be generalized to different imaging modalities and under-sampling rates. BRIEF DESCRIPTION OF DRAWINGS

[0035] The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application and are incorporated in and constitute a part of this application. The embodiments of this application and their

[0036] Figure 1 is an implementation flowchart of the free breathing fast three-dimensional cardiac magnetic resonance imaging joint motion estimation and motion compensated image reconstruction method based on implicit neural representation of the embodiment of the present application;

[0037] Figure 2 is a technical principle diagram of the free breathing fast three-dimensional cardiac magnetic resonance imaging joint motion estimation and motion compensated image reconstruction method based on implicit neural representation of the embodiment of the present application;

[0038] Figure 3 FIG. 1 is a schematic diagram of results of an implicit neural representation based free breathing fast three-dimensional cardiac magnetic resonance imaging combined motion estimation and motion compensated image reconstruction method according to an embodiment of the present application. DETAILED DESCRIPTION

[0039] It should be noted that the embodiments and features of the embodiments in the present application can be combined with each other without conflict. The present application will be described in detail below with reference to the accompanying drawings and in conjunction with the embodiments.

[0040] It should be noted that the steps shown in the flowchart of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and although a logical order is shown in the flowchart, in some cases, the steps shown or described herein can be executed in an order different from that shown.

[0041] The present application provides a free breathing fast three-dimensional cardiac magnetic resonance imaging method and system, comprising: acquiring under-sampled K-space (Fourier space) data under each breathing state by using a self-navigation imaging sequence; constructing a loss function including a motion compensation data fidelity term, a motion field smoothness regularization term, an image smoothness regularization term, and an image edge sharpness regularization term constraint as an optimization target; constructing an iterative reconstruction algorithm based on implicit neural representation to minimize the loss function, and jointly solving the respiratory motion field and the motion compensation reconstructed image. The present application directly solves the respiratory motion field and the three-dimensional cardiac magnetic resonance image from the under-sampled K-space data with the help of implicit neural representation, avoids the additional reconstruction and registration of each breathing state image in the traditional method, and combines various regularization term constraints to reduce the imaging time and improve the image signal-to-noise ratio, thereby accelerating the three-dimensional cardiac magnetic resonance imaging data acquisition while ensuring the image reconstruction quality, and without the need for additional training data.

[0042] In order to describe the above technical solutions in detail, the following embodiments are provided for detailed description:

[0043] Embodiment 1:

[0044] Reference Figure 1 and Figure 2 According to the free breathing fast three-dimensional cardiac magnetic resonance imaging method provided by the present application, the method comprises the following steps:

[0045] Step S1: acquiring under-sampled K-space data under each breathing state by using a self-navigation imaging sequence;

[0046] Step S1.1: acquiring free breathing three-dimensional cardiac magnetic resonance imaging data by using a self-navigation sequence;

[0047] Step S1.2: generating a breathing curve according to the self-navigation signal of the acquired data;

[0048] Step S1.3: dividing a plurality of respiratory phase according to the respiratory amplitude of the respiratory curve, each respiratory phase corresponding to a respiratory state from the end of inspiration to the end of expiration;

[0049] Step S1.4: extracting all data collected in each respiratory phase respectively to obtain the undersampled K-space data in the respiratory state.

[0050] Step S2: constructing a loss function including a motion compensation data fidelity term, a motion field smoothness regularization term, an image smoothness regularization term and an image edge sharpness regularization term constraint as an optimization objective.

[0051] Step S3: constructing an iterative reconstruction algorithm based on implicit neural representation to minimize the loss function, jointly solving the respiratory motion field and the motion compensation reconstruction image; using a Half-Quadratic Splitting algorithm to optimize the above loss function, by introducing an auxiliary variable z, the optimization problem is split into two subproblems; using implicit neural representation to solve subproblem 1 (first subproblem), the steps are as follows:

[0052] Step S3.1: constructing a motion solving network and an image reconstruction network, the motion solving network is used to output the motion field matrix corresponding to the spatial coordinates of the voxel points of the motion field, and the image reconstruction network is used to output the image matrix corresponding to the spatial coordinates of the voxel points of the three-dimensional cardiac magnetic resonance image;

[0053] Step S3.2: defining the spatial coordinates of the voxel points of the motion field and the magnetic resonance image in advance;

[0054] Step S3.3: taking the optimization objective of subproblem 1 as the loss function of network training, training the motion solving network and the image reconstruction network based on the defined spatial coordinates;

[0055] Step S3.4: inputting the spatial coordinates of the voxel points of the motion field and the magnetic resonance image into the trained motion solving network and the image reconstruction network respectively to obtain the respiratory motion field estimation result of subproblem 1 and the motion compensation reconstruction image.

[0056] The solution of subproblem 2 (second subproblem) is obtained by sequentially performing image denoising and image sharpening on the motion compensation reconstruction image of subproblem 1. By alternately solving subproblem 1 and subproblem 2, the respiratory motion solving result and the motion compensation reconstruction image of step S3 are obtained.

[0057] The application also provides a free breathing fast three-dimensional cardiac magnetic resonance imaging system, which can be realized by performing the process steps of the free breathing fast three-dimensional cardiac magnetic resonance imaging method, that is, the free breathing fast three-dimensional cardiac magnetic resonance imaging method can be understood by those skilled in the art as a preferred embodiment of the free breathing fast three-dimensional cardiac magnetic resonance imaging system.

[0058] Embodiment 2:

[0059] The application also provides a free breathing fast three-dimensional cardiac magnetic resonance imaging system, which comprises the following modules:

[0060] Module M1: acquiring undersampled K-space data in each breathing state by using a self-navigation imaging sequence;

[0061] Module M1.1: acquiring free breathing three-dimensional cardiac magnetic resonance imaging data by using a self-navigation sequence;

[0062] Module M1.2: generating a breathing curve according to the self-navigation signal of the acquired data;

[0063] Module M1.3: dividing a plurality of breathing phases according to the breathing amplitude of the breathing curve, each breathing phase corresponding to one breathing state from the end of inspiration to the end of expiration;

[0064] Module M1.4: extracting all data acquired in each breathing phase, respectively, to obtain undersampled K-space data in the breathing state.

[0065] Module M2: constructing a loss function comprising a motion compensation data fidelity term, a motion field smoothness regularization term, an image smoothness regularization term, and an image edge sharpness regularization term constraint as an optimization target.

[0066] Module M3: constructing an iterative reconstruction algorithm based on implicit neural representation to minimize the loss function, and jointly solving the respiratory motion field and the motion compensation reconstructed image;

[0067] The above loss function is optimized by using a Half-Quadratic Splitting algorithm, and by introducing an auxiliary variable z, the optimization problem is split into two sub-problems; implicit neural representation is used to solve sub-problem 1, comprising the following modules:

[0068] Module M3.1: constructing a motion solving network and an image reconstruction network, the motion solving network is used to output a motion field matrix corresponding to the spatial coordinates of the voxel points of the motion field, and the image reconstruction network is used to output an image matrix corresponding to the spatial coordinates of the voxel points of the three-dimensional cardiac magnetic resonance image;

[0069] Module M3.2: defining the spatial coordinates of the voxel points of the motion field and the magnetic resonance image in advance;

[0070] Module M3.3: the loss function of network training is the optimization target of sub-problem 1, and the motion solving network and the image reconstruction network are trained based on the defined spatial coordinates;

[0071] Module M3.4: the spatial coordinates of the motion field and the voxel points of the magnetic resonance image are input into the trained motion solving network and image reconstruction network respectively to obtain the respiratory motion field estimation result of sub-problem 1 and the motion compensation reconstruction image.

[0072] The solving result of sub-problem 2 is obtained by sequentially performing image denoising and image sharpening on the motion compensation reconstruction image of sub-problem 1. By alternately solving sub-problem 1 and sub-problem 2, the respiratory motion solving result and the motion compensation reconstruction image of the module M3 are obtained.

[0073] Embodiment 3:

[0074] An embodiment of the application applied to three-dimensional cardiac magnetic resonance coronary imaging will be described in detail below with reference to the accompanying drawings. Figure 1 The flowchart of the embodiment of the application is as follows, Figure 2 The algorithm principle diagram of the application is as follows. Specifically, the embodiment includes the following steps:

[0075] Step S1: acquiring undersampled K-space data in each respiratory state by using a self-navigation imaging sequence:

[0076] Step S1.1: the self-navigation gradient echo coronary imaging sequence acquires free-breathing three-dimensional cardiac magnetic resonance imaging data in the following manner: the sequence uses electrocardiogram gating to apply a T2 preparation pulse after the R wave of the electrocardiogram to enhance the myocardial blood contrast, and the T2 preparation time is 32 ms, after which the imaging data is acquired. The time of applying the T2 preparation pulse and the data acquisition time window are determined according to the cardiac cine image. The sequence uses a pseudo-radial sampling trajectory to acquire K-space data, and in each cardiac cycle, the data of the K-space center point is first acquired as navigation data, which is also used as imaging data, and then other imaging data is acquired from the K-space center to the periphery.

[0077] Step S1.2: extracting the navigation data acquired in each heartbeat, and after one-dimensional inverse Fourier transform, a one-dimensional projection image of the imaging region along the head-to-foot direction is obtained, then a one-dimensional projection image of the end-expiratory heartbeat is selected as a reference image, and a respiratory curve is obtained by correlating all one-dimensional projection images with the reference image;

[0078] Step S1.3: according to the respiratory amplitude of the respiratory curve, four respiratory phase sections are equally divided from bottom to top, corresponding to four respiratory states from the end of inspiration to the end of expiration; wherein, the four respiratory states are equally divided from the end of inspiration (the lowest point of the respiratory curve) to the end of expiration (the highest point of the respiratory curve) according to the respiratory amplitude.

[0079] Step S1.4: dividing the data collected in each heartbeat into the respiratory phase phase where the heartbeat is located, respectively merging all K-space data collected in the four respiratory phase phases to obtain the under-sampled K-space data corresponding to the respiratory state.

[0080] Step S2: constructing a loss function including a motion compensation data fidelity term, a motion field smoothness regularization term, an image smoothness regularization term and an image edge sharpness regularization term constraint as an optimization objective:

[0081] (1) Motion compensation data fidelity term (Ldata): First, the image in each respiratory state is calculated according to the motion compensation reconstruction image x and the solved motion field U in each respiratory state; then, the reconstruction K-space data obtained by coil sensitivity encoding and Fourier transform of the image in each respiratory state is calculated; finally, the square of the 2-norm of the difference between the reconstruction K-space data at the sampling position and the real acquisition K-space data is calculated:

[0082]

[0083] Wherein, represents the K-space data collected in the i-th respiratory state, represents the K-space sampling mask corresponding to the i-th respiratory state, represents the Fourier transform, represents the coil sensitivity, represents the respiratory motion field corresponding to the i-th respiratory state estimated by the motion solving network, represents the motion compensation reconstruction image estimated by the image reconstruction network. (2) Motion field smoothness regularization term (Lsmooth): The smoothness measure of the solved motion field U in each respiratory state is calculated, and the smoothness measure is the square of the 2-norm of the finite difference result:

[0084]

[0085]

[0086] Wherein, represents the finite difference operator.

[0087] (3) Image smoothness regularization term (Lsmooth): The smoothness measure of the motion compensation reconstruction image is calculated, and the smoothness measure is the 1-norm of the finite difference result, that is, the total variation:

[0088]

[0089] (4) Image edge sharpness regularization term (Ledge):​​​ ): To promote the smoothness of the motion-compensated reconstructed image while preserving the image details and high-frequency information such as image edges, the regularization term promotes the edge sharpness of the motion-compensated reconstructed image through image denoising and image sharpening algorithms:

[0090]

[0091] wherein, represents an implicit image edge sharpness loss term.

[0092] The above four losses are weighted and summed by weight coefficients, and then the losses of the four respiratory states obtained are summed to constitute a complete loss function. The weight coefficients are hyperparameters.

[0093] Step S3: Construct an iterative reconstruction algorithm based on implicit neural representation to minimize the loss function, and jointly solve the respiratory motion field and the motion-compensated reconstructed image. Specifically, a Half-Quadratic Splitting algorithm is used to optimize the loss function, and by introducing an auxiliary variable z, the optimization problem is split into two sub-problems:

[0094] The loss function optimized by sub-problem 1 is:

[0095]

[0096] wherein, z represents an auxiliary variable, represents a respiratory phase, i.e., a respiratory state, represents the K-space data acquired in the i-th respiratory state, represents the K-space sampling mask corresponding to the i-th respiratory state, represents a Fourier transform, represents a coil sensitivity, represents the respiratory motion field corresponding to the i-th respiratory state estimated by the motion solving network, represents the motion-compensated reconstructed image estimated by the image reconstruction network, represents a finite difference operator, , , represents a weight coefficient corresponding to different terms.

[0097] The loss function optimized by sub-problem 2 is:

[0098]

[0099] wherein, represents an image edge sharpness loss term, represents a corresponding weight coefficient.​​​

[0100] where, by the principle of Half-Quadratic Splitting algorithm. Subproblem 1 is to solve the optimization problem about the reconstructed image x, and subproblem 2 is to solve the optimization problem about the auxiliary variable z. The reconstructed image x in subproblem 2 is the reconstructed image x obtained after solving subproblem 1, on the basis of which the auxiliary variable z is optimized, and in the next iteration process, the current optimal auxiliary variable z is used to solve subproblem 1.

[0101] Step S3.1: Constructing the joint motion estimation and motion compensation reconstruction model based on implicit neural representation to solve subproblem 1 (first subproblem), as shown in Figure 2 The model includes the following modules:

[0102] (1) Respiratory motion estimation module: The input of this module is 4D coordinates (x, y, z, b) composed of 3D spatial coordinates and 1D respiratory phase coordinates, and the output is the displacement amount of the voxel with spatial coordinates (x, y, z) in the end-expiratory image along the left-right, up-down, and front-back directions when the end-expiratory image is transformed to the bth respiratory phase. First, the input 4D coordinates are hashed to obtain a feature vector, and then the feature vector is input into the motion solving network. The network architecture is a multi-layer perceptron with 3 layers, hidden neurons of 32, 32, 3, and an activation function of Elu except for the last layer, and the output length is 3, corresponding to the displacement along the left-right, up-down, and front-back directions.

[0103] (2) Motion compensation image reconstruction module: The input of this module is 3D spatial coordinates (x, y, z), and the output is the real and imaginary parts of the signal of the reconstructed three-dimensional cardiac magnetic resonance image voxel with spatial coordinates (x, y, z). First, the input 3D coordinates are hashed to obtain a feature vector, and then the feature vector is input into the image reconstruction network. The network architecture is a multi-layer perceptron with 6 hidden layers, hidden neurons of 128, 128, 128, 128, 128, 2, and an activation function of GELU except for the last layer, and the output length is 2, corresponding to the real and imaginary parts of the reconstructed voxel signal.

[0104] Step S3.2: Define the motion field and the spatial coordinates of the magnetic resonance image voxel: The motion field and the magnetic resonance image are equally divided to obtain 3D spatial coordinates; the 4 respiratory states are defined as 1D coordinates (1), (2), (3), and (4); and the 4D coordinates of each respiratory state motion field are composed of 3D spatial coordinates and 1D respiratory state coordinates.

[0105] Step S3.3: training the joint motion estimation and motion compensated reconstruction model with the undersampled K-space data of each respiratory state to be reconstructed: randomly reconstructing the motion field and the motion compensated reconstructed image of one respiratory state in each iteration. First, input the 4D coordinates of the motion field of the respiratory state selected in this iteration into the respiratory motion estimation module to obtain the motion field of the respiratory state; input the 3D spatial coordinates of the magnetic resonance image into the motion compensated image reconstruction module to obtain the real part and the imaginary part of the motion compensated reconstructed image; then, calculate the loss function corresponding to sub-problem 1 in step S3 according to the outputs of the related modules and the undersampled K-space data corresponding to the selected respiratory state, update the learnable parameters of the motion solving network and the image reconstruction network and the hash table by minimizing the loss function until the loss function no longer decreases, and obtain the trained joint motion estimation and motion compensated reconstruction model. In this embodiment, the Adam optimizer is adopted, the initial learning rate of the motion solving network is 0.001, the initial learning rate of the image reconstruction network is 0.01, the learning rate is exponentially decayed, the decay coefficient is 0.1, and a total of 500 iterations are trained.

[0106] Step S3.4: input the 4D coordinates of the voxel points of the motion field of each respiratory state and the 3D spatial coordinates of the voxel points of the magnetic resonance image into the trained motion solving network and the image reconstruction network respectively to obtain the respiratory motion field estimation result of sub-problem 1 and the motion compensated reconstructed image.

[0107] Step S3.5: solving sub-problem 2 (second sub-problem): first performing image denoising on the motion compensated reconstructed image of sub-problem 1, and then performing image sharpening on the denoised image to obtain the solution result of sub-problem 2. In this embodiment, the Non-Local Means algorithm is adopted to realize image denoising, and the Unsharp Masking algorithm is adopted to realize image sharpening.

[0108] Step S3.6: iteratively solving sub-problem 1 and sub-problem 2 according to the above steps to obtain the final motion estimation result and the motion compensated reconstructed image.

[0109] Figure 3 The results of this embodiment are shown in the diagram, in which PSNR is the peak signal-to-noise ratio, and SSIM is the structural similarity index. Among them, the three rows correspond to the reconstruction results of three healthy volunteers; the first column is the image reconstructed by the traditional method NR-SENSE from 3-fold accelerated acquisition data as the reference; the second column is the image reconstructed by the method of the application under 6-fold acceleration; the third and fourth columns are respectively the reconstruction results of different control methods under 6-fold acceleration. Compared with the control methods, the image reconstructed by the method of the application retains more complete details, has no obvious artifacts, and has obvious improvement in the peak signal-to-noise ratio and the structural similarity index relative to the reference image.

[0110] Based on the present application, the respiratory motion field can be estimated from the undersampled K-space data acquired with free breathing and the motion compensated reconstructed image can be obtained. Compared with the prior art, the present application does not need additional registration of the magnetic resonance image containing artifacts, and does not need additional training data, and is expected to improve the clinical usability of three-dimensional cardiac magnetic resonance imaging.

[0111] The above merely provides the preferred embodiments of the present application, but the protection scope of the present application is not limited thereto, and any person skilled in the art can easily think of the changes or replacements within the technical scope disclosed by the present application, which should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. Free-breathing fast three-dimensional cardiac magnetic resonance imaging method, characterized in that, The method comprises: acquiring undersampled K-space data of a heart in different breathing states, wherein the undersampled K-space data is undersampled Fourier space cardiac magnetic resonance data; constructing a motion solving network and an image reconstruction network based on implicit neural representation; the motion solving network adopts a multi-layer perceptron, wherein the motion solving network comprises three hidden layers, the number of hidden neurons of each hidden layer is 32, 32 and 3, the first hidden layer and the second hidden layer adopt an Elu activation function, the length of the output layer is 3, the input is encoded information of a respiratory motion field, wherein the encoded information of the respiratory motion field is obtained by performing hash encoding on spatial coordinates and breathing states of voxel points of the respiratory motion field, and an output result is a respiratory motion field estimation result, wherein the respiratory motion field estimation result is a displacement amount in different directions of each voxel point; the image reconstruction network adopts a multi-layer perceptron, wherein the image reconstruction network comprises six hidden layers, the input is encoded information of a reconstructed image, the encoded information of the reconstructed image is obtained by performing hash encoding on spatial coordinates of voxel points of the reconstructed image, and the output is a real part and an imaginary part of a voxel signal of a voxel point of the reconstructed image; constructing a loss function for joint motion estimation and motion compensated image reconstruction; iteratively updating and reconstructing the motion solving network and the image reconstruction network according to the sampled K-space data in different breathing states through the loss function, to obtain a final respiratory motion field estimation result and a motion compensated reconstructed image of the cardiac magnetic resonance reconstruction; the process of iteratively updating and reconstructing the motion solving network and the image reconstruction network comprises: optimizing the loss function through a semi-quadratic splitting method, and splitting the optimization steps of the loss function on the motion solving network and the image reconstruction network into a first sub-problem and a second sub-problem by introducing an auxiliary variable; for the first sub-problem, obtaining spatial coordinates of voxel points of a reconstructed image and a motion field, reconstructing the spatial coordinates of the voxel points of the reconstructed image and the motion field through the motion solving network and the image reconstruction network to obtain a respiratory motion field estimation result of the first sub-problem and a motion compensated reconstructed image, and optimizing the motion solving network and the image reconstruction network according to the respiratory motion field estimation result of the first sub-problem and the motion compensated reconstructed image; for the second sub-problem, denoising and sharpening the motion compensated reconstructed image of the first sub-problem to obtain a processed image of the second sub-problem, and optimizing the motion solving network and the image reconstruction network according to the processed image of the second sub-problem; iteratively and alternately performing the optimization steps of the first sub-problem and the second sub-problem to obtain a final respiratory motion field estimation result and a motion compensated reconstructed image; wherein the loss function comprises different loss terms of data fidelity, data smoothing and image detail constraint.

2. The method of claim 1, wherein the acquisition process of the undersampled K-space data comprises: ​ acquiring magnetic resonance imaging data of a heart in free breathing by a self-navigation imaging sequence; acquiring, for the magnetic resonance imaging data, data of a center point of K-space in the magnetic resonance imaging data within a collection window of each cardiac cycle as navigation data; generating a breathing curve according to the navigation data, dividing the breathing curve according to a breathing amplitude of the breathing curve to obtain breathing phase, wherein each breathing phase corresponds to a breathing state from end-inspiration to end-expiration; integrating the breathing phase with the magnetic resonance imaging data to obtain under-sampling K-space data in different breathing states.

3. The method of claim 1, wherein the loss function comprises a motion compensated data fidelity term, a motion field smoothness regularization term, an image smoothness regularization term, and an image edge sharpness regularization term; wherein the motion compensated data fidelity term is a loss term between the real under-sampling K-space data and K-space data synthesized by the motion compensated reconstructed image, the motion field smoothness regularization term is a smoothness measurement loss term of the breathing motion field, the image smoothness regularization term is a smoothness measurement loss term of the motion compensated reconstructed image, and the image edge sharpness regularization term is an edge sharpness loss term of the motion compensated reconstructed image.

4. The method of claim 3, wherein the loss function is a weighted sum result of the motion compensated data fidelity term, the motion field smoothness regularization term, the image smoothness regularization term, and the image edge sharpness regularization term.

5. The method of claim 1, wherein the loss function is a weighted sum result of the motion compensated data fidelity term, the motion field smoothness regularization term, the image smoothness regularization term, and the image edge sharpness regularization term.

6. The method of claim 5, wherein the loss function is a weighted sum result of the motion compensated data fidelity term, the motion field smoothness regularization term, the image smoothness regularization term, and the image edge sharpness regularization term.

7. A computer program product comprising computer executable instructions that, when executed by a computer, cause the computer to perform the method of any one of claims 1-6. The loss function for optimizing the motion solving network and the image reconstruction network for the first sub-problem is: wherein z denotes an auxiliary variable, denotes a respiratory phase, i.e. a breathing state, denotes the K-space data acquired for the th breathing state, denotes a K-space sampling mask corresponding to the th breathing state, denotes a Fourier transform, denotes a coil sensitivity, denotes a respiratory motion field corresponding to the th breathing state estimated by the motion solving network, denotes a motion compensated reconstructed image estimated by the image reconstruction network, denotes a finite difference operator, , , denote weight coefficients corresponding to different terms. ​ Loss function for optimizing the motion solving network and the image reconstruction network for the second sub-problem is: wherein, denotes an image edge sharpness loss term, denotes a corresponding weight coefficient.

7. A free-breathing fast three-dimensional cardiac magnetic resonance imaging system, characterized by, ​

Images