Whole-cycle lung function magnetic resonance imaging method

Abstract

The invention discloses a complete-cycle lung function magnetic resonance imaging method which comprises the following steps: in a free breathing state of a subject, acquiring lung coronal magnetic resonance images at different time points at the same anatomical position by using a gradient echo sequence, and selecting a lung aorta region R1 and a region R2 containing lung parenchyma and liver tissue; after the average signal time sequences of the R1 and the R2 are respectively subjected to high-pass filtering to remove the influence of respiration and low-pass filtering to remove the influence of heart blood flow, fitting is carried out in a cosine function form to obtain the phase of the signal, and the average signal time sequences in the R1 and the R2 are respectively reordered as RS1 and RS2 in a heart cycle and a respiration cycle according to the sizes of phase values. According to the method, the lung image is subjected to segmentation and registration and then is reordered according to RS1 and RS2, respiration and heart cycle parameters of each local area of the lung image are calculated, lung ventilation, lung perfusion, ventilation peak reaching time and perfusion peak reaching time images are generated, the local function of the lung can be non-invasively evaluated, and the method has very high clinical value.

Description

Technical Field

[0001] The present invention relates to the technical field of magnetic resonance imaging, and particularly to a method for whole-cycle pulmonary functional magnetic resonance imaging under free breathing. Background Art

[0002] Briefly introduce the value of local functional imaging of pulmonary ventilation and perfusion. The most traditional assessment of pulmonary ventilation function is vital capacity testing, which can obtain pulmonary ventilation function indexes such as forced vital capacity (FVC), forced expiratory volume in one second (FEV1), etc. These indexes play important roles in the diagnosis and treatment detection of pulmonary diseases. However, vital capacity measurement can only give global indexes and cannot observe local pulmonary ventilation function. In fact, the ventilation function of the lungs may be asymmetric between the left and right, with large differences between the upper and lower parts, etc. The local pulmonary ventilation function can be measured by imaging methods. Currently, the internationally recognized gold standard is the nuclear medicine method, such as scintigraphy and single photon emission computed tomography (SPECT). However, these methods require patients to inhale radioactive drugs, which is not beneficial to the health of patients, cannot be used for regular monitoring after treatment, and have low temporal and spatial resolutions. Magnetic resonance imaging (MRI) has no radiation and high resolution. One idea is to combine gas inhalation for MRI scanning to obtain pulmonary ventilation imaging, such as the relatively advanced hyperpolarized gas imaging of Xe-129. However, it has problems such as expensive gas preparation equipment, long preparation time, and complex operation, and it is difficult to be popularized in clinics in the short term. It is also possible to measure pulmonary function by inhaling oxygen to change the MRI contrast of the lungs, which still requires additional devices.

[0003] The perfusion detection of pulmonary blood vessels has great value in the clinical processes of pulmonary diseases such as the diagnosis of chronic thromboembolic pulmonary hypertension and the prediction of lung cancer metastasis. The local perfusion assessment of the lungs depends on imaging methods. Currently, the first choice is SPECT / CT imaging, but as mentioned before, it has radiation. In addition, contrast agents can be intravenously injected into patients, and then enhanced CT or MRI scans are performed. The contrast agent infiltrates into tissues along with blood vessels and changes the local relaxation time, thereby changing the image signal value. However, contrast agents have slight toxicity, so they are not applicable to groups such as pregnant women and children and cannot be used frequently. Some non-contrast-agent perfusion imaging methods attempt to solve this problem, such as magnetic resonance arterial spin labeling (ASL) method, but under the existing magnetic resonance hardware conditions, the signal-to-noise ratio of this method is not high, and its feasibility for pulmonary imaging with relatively low signal-to-noise ratio is particularly low.

[0004] Therefore, non-invasive assessment of local lung function, including local ventilation and perfusion functions, by imaging methods has high clinical value. Summary of the Invention

[0005] The object of the present invention is to propose a method for magnetic resonance functional imaging of the lungs, which is a non-invasive lung function evaluation method based on magnetic resonance imaging, and has the advantage of being able to dynamically and radiation-free synchronously quantify the ventilation and perfusion functions of the lungs.

[0006] The specific technical solution for achieving the object of the present invention is as follows:

[0007] A method for magnetic resonance imaging of the lungs over a full cycle, comprising the following specific steps:

[0008] Step S10, while the subject is breathing freely, continuously acquire a series of dynamic coronal magnetic resonance images of the lungs at a time resolution of 3 times per second or higher.

[0009] Step S20, determine two reference regions of interest for lung function image processing, including the visible lung parenchyma region on all images and the reference tissues for determining the cardiac cycle, i.e., the pulmonary aorta region R1, and the reference tissues for determining the respiratory cycle, i.e., the region R2 containing the lung parenchyma and liver tissue; respectively average the signal values within the reference regions of interest R1 and R2 to obtain the time series of average signal values for region R1 and region R2.

[0010] Step S30, for the series of dynamic coronal magnetic resonance images of the lungs acquired, segment the lung tissue region, then register the segmented lung tissue in the time domain, and then extract the dynamic signal change information of all local blocks of the registered lung tissue. Each local block can obtain a time series of signal values.

[0011] Step S40, for the time series of average signal values of region R1 obtained based on step S20, use the periodicity of the signal to rearrange its order within one cycle as RS1, rearrange the segmented and registered dynamic lung tissue images obtained based on step S30 into one cardiac cycle according to the order RS1, and perform data fitting based on the time series of signal values of each local block to obtain cardiac cycle parameters, including signal fluctuation amplitude, reference signal amplitude, and regional phase.

[0012] Step S50, for the time series of average signal values of region R2 obtained based on step S20, use the periodicity of the signal to rearrange its order within one cycle as RS2, rearrange the segmented and registered dynamic lung tissue images obtained based on step S30 into one respiratory cycle according to the order RS2, and perform data fitting based on the time series of signal values of each local block to obtain respiratory cycle parameters, including signal fluctuation amplitude, reference signal amplitude, and local phase.

[0013] Step S60, generate lung ventilation images, lung perfusion images, ventilation peak time images, and perfusion peak time images according to the results of steps S20 - S50.

[0014] Furthermore, the segmented lung tissue region specifically includes: implementing semi-automatic lung tissue segmentation using the idea of Non-Local Means, and then training the model using the nnUNet framework to achieve automatic lung tissue segmentation.

[0015] Furthermore, the time-domain registration of the segmented lung tissue specifically includes: selecting an image at the middle position of the respiratory cycle of the lungs as the reference image, and registering all images to this reference image; detecting the feature points of the images by combining the use of SURF (speeded-up robust features), BRIEF (binary robust independent elementary feature), ORB (oriented FAST and rotated BRIEF), (binary robust invariant scalable keypoints, BRISK), KAZE, and AKAZE (accelerated KAZE) methods, and registering the images using the matching feature points, with the registration method using affine transformation.

[0016] Furthermore, in step S40, the signal is rearranged in a cycle to RS1 using the periodicity of the signal. The segmented and registered dynamic lung tissue images obtained based on step S30 are rearranged in a cardiac cycle according to the order RS1, and the cardiac cycle parameters are obtained by fitting the data according to the signal value time series of each local block. Specifically: performing high-pass filtering on the signal value time series of each local block obtained in step S30 and the average signal value time series of region R1 obtained based on step S20 to remove the influence of respiratory motion; then fitting the average signal value time series of region R1 after high-pass filtering in a cardiac cycle in the form of a cosine function, and reordering the mean signals within region R1 collected at different time points according to the fitted phase difference, that is, the phase at each moment in the cardiac cycle, to obtain a new signal order RS1; reordering the signal value time series of each local block of the lung tissue after high-pass filtering in a cardiac cycle according to the new signal order RS1, so as to achieve a data sample time resolution higher than the acquisition resolution and obtain the cardiac cycle parameters.

[0017] Further, in step S50, the signal is rearranged in one period according to its periodicity to obtain RS2, and the segmented and registered dynamic lung tissue images obtained based on step S30 are rearranged in one respiratory cycle according to the order RS2. The respiratory cycle parameters are obtained by fitting the time series of signal values of each local block. Specifically: the time series of signal values of each local block of the lung tissue obtained in step S30 and the time series of average signal values of region R2 obtained based on step S20 are subjected to low-pass filtering to remove the influence of the heartbeat, and then the low-pass filtered signals are fitted in one respiratory cycle in the form of a cosine function. According to the fitted phase difference, that is, based on the phase at each moment in the respiratory cycle, the mean signals in region R2 collected at different time points are reordered to obtain a new signal order RS2; the time series of signal values of each local block of the lung tissue after low-pass filtering are reordered in one respiratory cycle according to the new signal order RS2, so as to achieve a data sample time resolution higher than the acquisition resolution and obtain the respiratory cycle parameters.

[0018] Advantages of the present invention: Due to its non-invasive and non-radiative properties, as well as the characteristic of being able to quantify functional parameters, the present invention can replace CT and nuclear medicine scans, avoiding the cumulative radiation risk for children, pregnant women, and long-term follow-up patients. It synchronously evaluates ventilation and perfusion functions and can more accurately distinguish diseases. Brief Description of the Drawings

[0019] Figure 1 is the overall implementation flowchart provided by an example of the present invention;

[0020] Figure 2 is the flowchart of post-processing the cardiac cycle based on the collected magnetic resonance images provided by an example of the present invention;

[0021] Figure 3 is the flowchart of post-processing the respiratory cycle based on the collected magnetic resonance images provided by an example of the present invention;

[0022] Figure 4 is the re-ordered standard cardiac cycle curve provided by an example of the present invention;

[0023] Figure 5 is the fitted standard respiratory cycle curve provided by an example of the present invention;

[0024] Figure 6 is the functional magnetic resonance image of the lungs of a healthy volunteer output by an example of the present invention. Detailed Embodiments

[0025] To better elaborate the purpose, technical features, and application potential of the present invention, the following examples in combination with the accompanying drawings and implementation steps are used to further illustrate the present invention. The specific embodiments described below are only used to better explain the present invention and are not used to limit the protection scope of the present invention.

[0026] Referring to Figure 1 , the pulmonary functional magnetic resonance imaging method involved in the present invention includes the following steps:

[0027] Step S10, dynamically acquire a series of pulmonary magnetic resonance images in a free-breathing state; specifically, in the state where the subject breathes freely, continuously acquire a series of dynamic pulmonary magnetic resonance imaging at a time resolution of 3 times per second or higher.

[0028] In this example, the sequence used for magnetic resonance imaging is a two-dimensional gradient echo sequence. By setting a short repetition time (TR ~ 3 ms), a high time resolution is achieved; 200 images are acquired for each layer, taking 1 minute; by setting a short echo time (TE ~ 1 ms), signal attenuation of the lung parenchymal tissue is avoided, and a relatively high signal-to-noise ratio of the lung tissue is achieved; the field of view FOV = 500 * 500 mm 2 Includes the entire thorax; a slice thickness of 15 mm is selected to achieve a relatively high signal-to-noise ratio.

[0029] In another example of the present invention, an implementation method for three-dimensional magnetic resonance acquisition in step S10 is demonstrated; the rationality for this is that, when the hardware conditions of the MRI scanner permit, three-dimensional acquisition has a better signal-to-noise ratio; the magnetic resonance sequence uses a three-dimensional gradient echo sequence, sets a small flip angle (FA ~ 3.5°), and a short repetition time (TR ~ 1.9 ms); partial Fourier acquisition mode (6 / 8) is enabled for acceleration; each three-dimensional slab contains 3 layers, each layer with a thickness of 6 mm; the echo time is set to the shortest (TE ~ 0.9 ms) to achieve a relatively high signal-to-noise ratio of the lung tissue; the field of view FOV = 500 * 500 mm 2 Includes the entire thorax; the final acquisition time for each slab is about 1 minute, and images at 200 time points are acquired.

[0030] Step S20, draw the regions of interest of the pulmonary reference tissue; specifically, first determine the regions of the lung tissue on all the images acquired in S10; subsequently, on the magnetic resonance pulmonary image at a certain moment, outline the boundaries of the region R1 of the reference tissue for determining the cardiac cycle and the region R2 of the reference tissue for determining the respiratory cycle, and use the boundaries for all magnetic resonance pulmonary images at all moments.

[0031] In this example, all lung parenchymal tissues were determined as the regions of interest (ROIs) through an automatic segmentation algorithm. The average signal value of the ROIs in the image at each moment was calculated automatically. The magnetic resonance image at the moment when the signal value was at the median among all images was selected to draw the regions R1 and R2 of the reference tissue. Subsequently, the central pulmonary artery was selected for R1, and the position of the diaphragm was selected for R2. The upper 1 / 3 was lung tissue, and the lower 2 / 3 was liver tissue. Their boundaries were determined by manual delineation.

[0032] The specific objects and segmentation methods of the ROIs can be adjusted according to the specific clinical or research scenarios, and are not limited to the semi-automatic and automatic segmentation of the lung parenchymal range or the manual segmentation of the reference tissue region proposed in this example.

[0033] As an example, for lung images containing tumors, since it is considered that tumors do not have the possibility of ventilation from a clinical perspective, the tumors can be excluded by reverse delineation within the ROI of the lung tissue.

[0034] Step S30: Segment the lung tissues from a series of collected images; The data used for the training of the automatic segmentation model was sourced from 25 semi-automatically segmented cases, with a total of 5000 layers of lung images. Each data set contained 200 dynamically acquired coronal magnetic resonance images of the lungs. The 200 images reflected the characteristics of the up-and-down periodic movement of the lungs during the breathing process. To improve the segmentation efficiency, a semi-automatic segmentation method was adopted during manual segmentation. Manual segmentation was only performed on the positions with the largest and smallest volumes within the lung movement cycle. The regions in the middle positions were completed using a 3D interpolation algorithm based on the idea of Non-Local Means, that is, the pixel blocks corresponding to the current pixel were found in the maximum and minimum position layers of the lungs. By comparing the similarity of the corresponding pixel blocks in the maximum and minimum position layers with the pixel block of the current layer pixel, it was determined whether the current layer pixel block belonged to the lung tissue, thus realizing the semi-automatic segmentation of the lung tissue. The nnUNet framework was used to train the model for the semi-automatically segmented lung tissue images. Among them, the parameter epoch was set to 250, and 5-fold training was adopted to finally obtain a 2D segmentation model to achieve the automatic segmentation of the lung tissue.

[0035] Perform time-domain registration on a series of segmented lung tissue images; Specifically, an image of the lung in the middle position after segmentation in step S20 was selected as the reference image, and all the lung tissue images after segmentation in step S20 were registered to this reference image. The dynamic signal change information of all local blocks of the registered lung tissue was extracted, and each local block could obtain a time series of signal values.

[0036] In this example, the feature points of the image are detected by combining SURF (speeded-up robust features) (SURF), BRIEF (binary robust independent elementarv feature), ORB (oriented FAST and rotatedBRIEF), (binary robust invariant scalable key points, BRISK), KAZE, AKAZE (accelerated KAZE) methods, and the images are registered using the matched feature points. The registration method uses affine transformation.

[0037] Step S40, rearrange the signal value time series of each local block of the lung obtained based on step S30, obtain the signal value change curve of each local block of the lung in a cardiac cycle, and obtain the cardiac cycle parameters by calculation, including signal fluctuation amplitude, reference signal amplitude, and local phase. Specifically, first convert the original signal in the time domain (the time series of the spatial average signal value of region R1 and the time series of the signal value of each local block of the lung obtained in S30) into a frequency domain signal, apply high-pass filtering to eliminate the low-frequency fluctuation of the signal caused by breathing, and then convert the filtered frequency domain signal back to the time domain signal; fit the time domain signal of region R1 of the reference tissue to a series of periodic cardiac cycle curves, and record its frequency ω cardic ,; then, based on the cosine function curve as the basis for fitting the cardiac cycle, the R1 time domain signal is fitted to obtain the phase of the signal at each moment in the cardiac cycle, and the mean signals in the region R1 collected at different time points are reordered to obtain a new signal order RS1.

[0038] For each local block position (x, y) of the lung tissue, the new signal intensity time series S can be obtained by reordering based on the new signal sequence RS1. cardic (x, y, t); using the interpolation algorithm, reconstruct a standard cardiac cycle to obtain the signal fluctuation amplitude, baseline signal strength, and phase relative to the reference tissue R1. The fitting parameters included.

[0039] In this example, step S40 first removes the time point data of the first 20 images that have not reached a steady state to ensure data stability. The time domain signals of the remaining 180 data points are converted into frequency domain signals using fast Fourier transform (FFT), high-pass filtered with a cutoff frequency of 0.8 Hz, and then converted back into new time domain signals using inverse fast Fourier transform (IFFT); the pulmonary aorta, i.e., the spatial average signal value time series S of region R1, is removed. R1The lowest 5% and the highest 5% of the signals (x, y, t) are used to avoid the influence of outliers and improve the stability of the algorithm. In this example, the cosine function curve is selected as the basis for fitting the cardiac cycle, based on the signal S R1 (x, y, t), and the difference between the maximum value max(S R1 (x, y, t)) and the minimum value min(S R1 (x, y, t)) at all times is used to obtain the signal fluctuation amplitude A R1 (x, y). The mean value of the maximum and minimum values is used to obtain the reference signal amplitude I R1 (x, y). Subsequently, based on the dynamic change of the time series of the spatially averaged signal value in the region R1 over time, the R1 signal at each moment is classified into the rising period or the falling period for phase fitting. For the signal in the rising period, the following formula is followed: For the signal in the falling period, the following formula is followed: The phase of the signal S R1 (x, y, t) is calculated and obtained Each S R1 (x, y, t) is reordered according to the phase value to form a cardiac cycle, and the sequence order is RS1. Based on this new signal order, for each local block position (x, y) of the lung tissue, a new signal intensity time series S cardic (x, y, t) can be obtained. Using the Nadaraya-Watson non-parametric regression algorithm based on the Gaussian kernel to perform interpolation at 30 equally spaced time nodes, the fitting parameters including the signal fluctuation amplitude A cardic (x, y), the reference signal intensity I cardic (x, y), and the phase are obtained, as well as the dynamic lung perfusion images at 30 time nodes. Attached Figure 2 is a flowchart for processing the heartbeat cycle based on the collected magnetic resonance images provided by an example of the present invention;

[0040] Attached Figure 4 is the signal intensity curve of a cardiac cycle after reordering the time series of the spatially averaged signal value of the pulmonary aorta provided by the example.

[0041] Step S50: Rearrange the time series of signal values for each local block of the lungs obtained based on Step S30 to obtain the signal value change curve of each local block of the lungs in one respiratory cycle, and calculate the respiratory cycle parameters, including the signal fluctuation amplitude, the reference signal amplitude, and the local phase at each moment. Specifically, first convert the original signals in the time domain (the time series of spatially averaged signal values in region R2 and the time series of signal values for each local block of the lungs obtained in S30) into frequency domain signals, apply low-pass filtering to eliminate the high-frequency signal fluctuations caused by blood inflow, and then convert the filtered frequency domain signals back into time domain signals; fit the time domain signal of region R2 of the reference tissue to a series of periodic respiratory cycle curves and record its frequency ω resp ; Subsequently, based on the cosine function curve as the basis for respiratory cycle fitting, fit the R2 time domain signal to obtain the phase of the signal at each moment within the respiratory cycle, and reorder the mean signals within the R2 region collected at different time points to obtain a new signal order RS2.

[0042] Optionally, before calculating the signal fluctuation amplitude in Step S50, outliers with excessive inhalation and exhalation in the lungs in all images can be removed by methods such as variance test to improve the algorithm stability.

[0043] Optionally, at the end of Step S50, interpolation can be used, and Nadaraya-Watson non-parametric regression based on Gaussian kernel is used for interpolation, and interpolation is performed at 60 equally spaced time points within one respiratory cycle;

[0044] In this example, in Step S50, the time point data of the first 20 images that have not reached a steady state are first removed to ensure data stability. The time domain signal of the remaining 180 data points is converted into a frequency domain signal using the fast Fourier transform FFT, low-pass filtering with a cut-off frequency of 0.6 Hz is used, and then it is converted back into a new time domain signal S resp (x, y, t); Remove the signals with the lowest 5% of the lung inhalation expansion amount (S resp (x, y, t)) and the highest 5% of the lung exhalation contraction amount (S resp (x, y, t)) from the time series of spatially averaged signal values in region R2 to avoid the influence of outliers and improve the algorithm stability. In this example, a simple cosine function curve is selected as the basis for lung respiratory cycle fitting. Based on the change of the time series of spatially averaged signal values S R2 (x, y, t) in region R2, the difference between the maximum value max(S R2 (x, y, t)) and the minimum value min(S R2 (x, y, t)) at all moments is used to obtain the signal fluctuation amplitude A R2 (x, y), and the mean value of the maximum and minimum values is used to obtain the reference signal amplitude I R2(x, y); Based on the dynamic change of the signal over time after spatial averaging in the R2 region, classify the lung images at each moment as the exhalation or inhalation phase for phase fitting; for the images in the inhalation phase, follow the formula: For the images in the exhalation phase, follow the formula: Calculate the phase For each s R2 (x, y, t) is reordered according to the phase value to form a cardiac cycle, and the sequence order is RS2. Based on this new signal order, for each local block position (x, y) of the lung tissue, a new signal intensity time series S resp (x, y, t) can be obtained; use the Nadaraya-Watson non-parametric regression algorithm based on a Gaussian kernel to perform interpolation at 60 equally spaced time nodes, and obtain fitting parameters including the signal fluctuation amplitude A resp (x, y), the reference signal intensity I resp (x, y), the phase and so on, as well as the dynamic lung respiration images at 60 time nodes. Attached Figure 3 is a flowchart for respiratory cycle processing based on the collected magnetic resonance images provided by an example of the present invention;

[0045] Attached Figure 5 is the signal intensity curve of one respiratory cycle after the mean signal in the region R2 provided by the example is rearranged in the RS2 order.

[0046] Step S60, according to the results of the previous steps, draw parameter maps such as lung ventilation images, lung perfusion images, and time-to-peak images. Specifically, use the results of step S30 to segment and register the lung tissue region, and extract the signal value time series of each local block of the registered lung tissue; then, according to the cardiac cycle parameters obtained in step S40, calculate the perfusion parameter map and the perfusion time-to-peak map; according to the respiratory cycle parameters obtained in step S50, calculate the ventilation parameter map and the ventilation time-to-peak map;

[0047] Optionally, in step S60, based on the perfusion time-to-peak map and the ventilation time-to-peak map, and based on any reasonable threshold, determine whether the perfusion and ventilation time-to-peak of each lung tissue match, and calculate the perfusion-ventilation matching situation map.

[0048] Optionally, in step S60, based on the dynamic change of the ventilation volume, using the ventilation fraction change rate value and the ventilation fraction value as two dimensions, draw a ventilation cycle loop diagram.

[0049] In this example, the perfusion-weighted signal value satisfies the cosine function: Based on the change of QWI(x, y, t), record the time when each position reaches the maximum value, and output it as the perfusion time-to-peak map Q-TTP(x, y); by Calculate the pulmonary perfusion image. The ventilation function signal value satisfies the cosine function: The quantitative ventilation V(x, y) of the lung tissue at a certain local block position (x, y) is: V(x, y) = ∝(S ex (x, y) - S in (x, y)) = ∝A resp (x, y); where S ex (x, y) is the signal intensity at the end of expiration, that is, the highest position of the respiratory cycle signal, and S in (x, y) is the signal intensity at the end of inspiration, that is, the lowest position of the respiratory cycle signal; calculate the local ventilation map of the lung tissue at a certain position (x, y): Calculate the ventilation fraction of the lung tissue at a certain position (x, y) at a certain moment t Based on the change of VF(x, y, t), record the proportion of the time interval from the end of inspiration to the end of expiration in the respiratory cycle, and output it as the ventilation time-to-peak map V-TTP(x, y); for the lung tissue with Q-TTP(x, y) greater than 200 ms, it is determined as perfusion lag; for the lung tissue with V-TTP greater than 60% of the respiratory cycle, it is determined as ventilation lag; based on the above determination, classify the lung tissue into three labels: ventilation-perfusion matching, perfusion lag, and ventilation lag, and thus output the perfusion-ventilation matching situation map V / Q-Match(x, y); taking the ventilation fraction change rate value and the ventilation fraction value as two dimensions, draw the ventilation cycle loop diagram; Figure 6 Figure is the functional magnetic resonance image (QWI, VR, Q-TTP, V-TTP) of the lungs of healthy volunteers output by the example of the present invention.

Claims

1. A full-cycle pulmonary function magnetic resonance imaging method, characterized in that The steps of the method include: Step S10, continuously acquire a series of dynamic coronal lung magnetic resonance images at a time resolution of 3 times per second or higher while the subject is breathing freely; Step S20, determine two reference regions of interest for lung function image processing, including the visible lung parenchyma region on all images, and the reference tissues for determining the cardiac cycle, i.e., the pulmonary aorta region R1, and the reference tissues for determining the respiratory cycle, i.e., the region R2 containing lung parenchyma and liver tissue; respectively average the signal values within the reference regions of interest R1 and R2 to obtain the average signal value time series of region R1 and region R2; Step S30, segment the lung tissue region from the acquired series of dynamic coronal lung magnetic resonance images, then register the segmented lung tissue in the time domain, and then extract the dynamic signal change information of all local blocks of the registered lung tissue. Each local block can obtain a signal value time series; Step S40, for the average signal value time series of region R1 obtained based on Step S20, use the periodicity of the signal to rearrange its order within one cycle as RS1, rearrange the segmented and registered dynamic lung tissue images obtained based on Step S30 into one cardiac cycle according to the order RS1, and perform data fitting according to the signal value time series of each local block to obtain cardiac cycle parameters, including signal fluctuation amplitude, reference signal amplitude, and regional phase; Step S50, for the average signal value time series of region R2 obtained based on Step S20, use the periodicity of the signal to rearrange its order within one cycle as RS2, rearrange the segmented and registered dynamic lung tissue images obtained based on Step S30 into one respiratory cycle according to the order RS2, and perform data fitting according to the signal value time series of each local block to obtain respiratory cycle parameters, including signal fluctuation amplitude, reference signal amplitude, and local phase; Step S60, generate lung ventilation images, lung perfusion images, ventilation peak time, and perfusion peak time images according to the results of Steps S20 - S50.

2. The full-cycle pulmonary function magnetic resonance imaging method according to claim 1, wherein The segmentation of the lung tissue region specifically includes: implementing semi-automatic lung tissue segmentation using the idea of Non - Local Means, and then training the model using the nnUNet framework to achieve automatic segmentation of the lung tissue.

3. The full-cycle pulmonary function magnetic resonance imaging method according to claim 1, wherein The registration of the segmented lung tissue in the time domain specifically includes: selecting an image at the middle position of the respiratory cycle of the lung as the reference image, and registering all images to this reference image; detecting the feature points of the images by using methods such as SURF, BRIEF, ORB, BRISK, KAZE, and AKAZE, and registering the images using the matching feature points. The registration method uses affine transformation.

4. The full-cycle pulmonary function magnetic resonance imaging method according to claim 1, characterized in that, As described in step S40, the signals are rearranged in a cycle with the order RS1 by using the periodicity of the signals. The segmented and registered dynamic lung tissue images obtained based on step S30 are rearranged in a cardiac cycle according to the order RS1, and cardiac cycle parameters are obtained by performing data fitting based on the time series of the signal values of each local block. Specifically: the time series of the signal values of each local block obtained in step S30 and the time series of the average signal values of region R1 obtained based on step S20 are subjected to high-pass filtering to remove the influence of respiratory motion; then the time series of the average signal values of region R1 after high-pass filtering are fitted in a cardiac cycle in the form of a cosine function, and according to the fitted phase difference, that is, the phase at each moment in the cardiac cycle, the average signals in region R1 collected at different time points are reordered to obtain a new signal order RS1; the time series of the signal values of each local block of the lung tissue after high-pass filtering are reordered in a cardiac cycle according to the new signal order RS1, so as to achieve a time resolution of data samples higher than the acquisition resolution and obtain cardiac cycle parameters.

5. The full-cycle pulmonary function magnetic resonance imaging method according to claim 1, characterized in that, As described in step S50, the signals are rearranged in a cycle with the order RS2 by using the periodicity of the signals. The segmented and registered dynamic lung tissue images obtained based on step S30 are rearranged in a respiratory cycle according to the order RS2, and respiratory cycle parameters are obtained by performing data fitting based on the time series of the signal values of each local block. Specifically: the time series of the signal values of each local block of the lung tissue obtained in step S30 and the time series of the average signal values of region R2 obtained based on step S20 are subjected to low-pass filtering to remove the influence of heartbeat, and then the signals after low-pass filtering are fitted in a respiratory cycle in the form of a cosine function, and according to the fitted phase difference, that is, based on the phase at each moment in the respiratory cycle, the average signals in region R2 collected at different time points are reordered to obtain a new signal order RS2; the time series of the signal values of each local block of the lung tissue after low-pass filtering are reordered in a respiratory cycle according to the new signal order RS2, so as to achieve a time resolution of data samples higher than the acquisition resolution and obtain respiratory cycle parameters.

Images