A 4D-CT-based method for visualizing lung perfusion

Abstract

This invention discloses a 4D-CT-based method for visualizing lung perfusion, belonging to the field of intelligent medical technology. The method includes: acquiring 4D-CT scan data of a patient's lungs; acquiring the patient's cardiac waveform based on the scan data; acquiring the patient's lung waveform based on the scan data; performing correlation analysis on the cardiac and lung waveforms to obtain a correlation coefficient; and outputting a time-series color image of lung perfusion based on the correlation coefficient. This invention can directly utilize conventional breath-hold 4D-CT volumetric data, intuitively presenting the propagation characteristics of lung perfusion without the need for additional hardware or examination procedures. This method does not require the injection of contrast agents, completely avoiding the risks associated with contrast agents, and has better safety and repeatability. Furthermore, the overall calculation process is highly automated and has a short detection time, improving examination and evaluation efficiency to a certain extent.

Description

Technical Field

[0001] This invention relates to the field of intelligent medical technology, and in particular to a method for visualizing lung perfusion based on 4D-CT. Background Technology

[0002] Pulmonary perfusion describes the size and distribution of blood flow within the lungs. When we inhale, oxygen from the air is stored in tiny alveoli in the lungs, while simultaneously, blood pumped by the heart flows through the capillaries surrounding these alveoli, facilitating efficient gas exchange. Matching ventilation and perfusion is essential for proper oxygenation. Pulmonary perfusion assessment is crucial for the diagnosis and functional evaluation of lung diseases such as pulmonary embolism, chronic obstructive pulmonary disease (COPD), and pulmonary hypertension. Currently, the most commonly used clinical methods for lung perfusion imaging include computed tomography pulmonary angiography (CTPA) and radionuclide lung ventilation / perfusion imaging (V / Q lung imaging). Both rely on exogenous contrast agents or radioactive tracers and have the following limitations: CTPA requires the injection of iodine-containing contrast agents, which may cause allergic reactions and kidney damage, and is not suitable for patients with renal insufficiency or pregnancy; V / Q lung imaging is complex to operate, has low spatial resolution, requires exposure to radioactive materials, and has limited sensitivity to subsegmental embolism.

[0003] In recent years, with the improvement of equipment hardware and analysis software, perfusion analysis methods based on dynamic X-ray examination have emerged. This technology assesses perfusion by analyzing the temporal changes of lung pixels during the cardiac cycle. Although it does not require contrast agents, it has fundamental limitations: 1. Insufficient image quality and spatial information: It is a two-dimensional projection imaging, with overlapping anterior and posterior anatomical structures, resulting in low spatial resolution and making it difficult to accurately locate small lesions or distinguish blood vessels from adjacent tissues.

[0004] 2. Limited quantitative accuracy: Quantification based on pixel changes in planar projection is greatly affected by tissue overlap, patient position, and respiratory movements, resulting in insufficient stability and accuracy of quantitative results, and is prone to errors in complex cases.

[0005] 3. Separation of functional and anatomical information: The perfusion function images obtained lack a corresponding high-resolution three-dimensional anatomical background and need to be interpreted in comparison with other images (such as CT), making it impossible to achieve a true integrated "structure-function" assessment.

[0006] In summary, existing technologies either rely on exogenous contrast agents or are limited by low-dimensional, static imaging modalities, making it impossible to achieve high-resolution, three-dimensional, dynamic, and anatomically and functionally integrated lung perfusion assessment under non-invasive conditions. Therefore, a novel clinical solution is urgently needed. Summary of the Invention

[0007] The purpose of this invention is to provide a lung perfusion visualization method based on 4D-CT, so as to at least partially solve the above-mentioned technical problems.

[0008] To solve the above-mentioned technical problems, the present invention provides the following technical solution: A 4D-CT-based method for visualizing lung perfusion includes: Step S100: Acquire the patient's lung 4D-CT scan data; Step S200: Obtain the patient's cardiac waveform based on the scan data; Step S300: Obtain the patient's lung waveform based on the scan data; Step S400: Perform correlation analysis on the heart waveform and lung waveform to obtain the correlation coefficient; Step S500: Output the time series lung perfusion color image based on the correlation coefficient.

[0009] Based on the above technical solution, the present invention has at least the following beneficial effects: This invention presents a 4D-CT-based lung perfusion visualization method that visualizes lung perfusion in the form of a color dynamic perfusion map. It can directly utilize conventional breath-hold 4D-CT volumetric data, without requiring additional hardware or examination procedures. Based on cardiac reference waveform extraction, heart rate adaptive bandpass filtering, and pixel grayscale cross-correlation analysis, it automatically obtains the intensity distribution and phase delay of perfusion signals related to cardiac activity within the thoracic region of interest (ROI), generating a time-progressing color overlay image to intuitively present the propagation characteristics of lung perfusion. Compared to existing perfusion imaging methods that rely on radioactive tracers or contrast agents, this method does not require contrast agent injection, completely avoiding contrast agent-related risks and offering better safety and repeatability. Simultaneously, the overall computational process is highly automated and has a short detection time, improving examination and assessment efficiency to a certain extent. It is suitable for applications requiring rapid assessment or large-scale screening, such as emergency departments, providing reliable image information for screening diseases such as pulmonary embolism or chronic thromboembolic pulmonary hypertension. Attached Figure Description

[0010] The accompanying drawings in this application are intended to supplement the textual description in the specification with graphics, and to further explain the technical solution of this application. They do not constitute an undue limitation on this application.

[0011] Figure 1This is a flowchart illustrating the lung perfusion visualization method based on 4D-CT of the present invention. Figure 2 This is a schematic diagram illustrating the principle of the lung perfusion visualization method based on 4D-CT of the present invention; Figure 3 This is a schematic diagram illustrating the algorithm principles for the heart and lung waveforms in this invention. Figure 4 This is a schematic diagram of obtaining the heart waveform in this invention, wherein (a) shows the selection of the heart mask, (b) shows the waveform signal of the corresponding initial heart waveform after high-pass filtering and delinearization, and (c) shows the heart rate estimation result and the normalized heart waveform after band-pass filtering. Figure 5 This is a schematic diagram of lung waveform acquisition in this invention, wherein (a) shows the generation of the thoracic mask, with the arrow pointing to any point in the lung area, and (b) shows the lung waveform at that point after being filtered and normalized by the same adaptive passband range pb as the heart waveform; Figure 6 This is a graph showing the changing trends of the heart waveform at the same moment and the lung waveform at any point in the present invention. Figure 7 This is a dynamic lung perfusion color image from the present invention, wherein (a) shows a second background image, and (b) shows the diastolic phase (offset frame number). FL (c) Perfusion images at (=0), with maximum correlation intensity in the cardiac region; (c) shows the systolic phase (offset frame number). FL The perfusion image at =4) shows the maximum correlation intensity shifted to the pulmonary vessels. Detailed Implementation

[0012] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0013] In the description of this invention, it should be understood that the terms "center," "lateral," "longitudinal," "front," "rear," "left," "right," "upper," "lower," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limiting the scope of protection of this invention.

[0014] In recent years, four-dimensional computed tomography (4D-CT) has been increasingly widely used in clinical practice. Its core feature is the continuous imaging of a three-dimensional CT volume containing moving structures over a period of time, thereby obtaining a dynamic volumetric dataset. With the development of large-field-of-view CT technology, 4D-CT is no longer limited to radiotherapy planning but has also demonstrated significant value in clinical diagnosis. For example, strain analysis based on 4D dynamic ventilation CT can record lung tissue deformation during ventilation, thus enabling the monitoring of pathophysiological changes in COPD.

[0015] This invention proposes a method for imaging and visualizing pulmonary blood flow perfusion based on 4D-CT scan images. The physiological principle behind this method lies in the changes in X-ray permeability during the cardiac cycle. Each heartbeat causes periodic changes in blood flow within the heart and pulmonary artery, which can be further manifested as fluctuations in the grayscale values ​​of CT images. Specifically, during end-diastole, the heart has a larger blood volume, making it less permeable to X-rays, while the pulmonary artery has a smaller blood volume, making it more permeable to X-rays, resulting in a lower grayscale value in the image. During end-systole, the heart's blood flow decreases, making it more permeable to X-rays, while the pulmonary artery's blood flow increases, decreasing X-ray permeability, corresponding to an increase in the grayscale value in the image. Continuous 4D-CT sequences obtained under breath-hold conditions contain pulmonary blood flow information caused by cardiac pumping. By extracting the changes in grayscale values ​​over time, a correlation between cardiac and pulmonary waveforms can be established and visualized as a color dynamic perfusion map. The advantage of this technology is that it can directly utilize existing breath-hold 4D-CT data without requiring additional examinations. Meanwhile, this method is non-invasive, requiring no contrast agents or radioactive nuclides, making it suitable for a wide range of people with virtually no contraindications. Furthermore, the entire scan takes only about 8 seconds, significantly improving examination efficiency, especially beneficial for rapid assessment in emergency situations.

[0016] This invention provides a lung perfusion visualization method based on 4D-CT, such as... Figure 1-3 As shown, it includes: Step S100: Acquire the patient's lung 4D-CT scan data; Step S200: Obtain the patient's cardiac waveform based on the scan data; Step S300: Obtain the patient's lung waveform based on the scan data; Step S400: Perform correlation analysis on the heart waveform and lung waveform to obtain the correlation coefficient; Step S500: Output the time series lung perfusion color image based on the correlation coefficient.

[0017] Each step of the present invention will be described in detail below: Step S100: Acquire the patient's lung 4D-CT scan data This step involves acquiring 4D-CT volumetric data. The specific steps are as follows: (1) 4D-CT scan plan: Dynamic volumetric scanning mode was adopted, with the patient's heart as the center, and 320-slice detectors covering a maximum scan length of 16cm. Total scan time = 8.1s; number of scan frames = 10fps; collimation = 0.5mm; reconstruction nucleus = FC15 (mediastinum); iterative reconstruction = AIDR3D (mildsetting); slice thickness / interval = 0.5mm / 0.5mm. During the scan, the patient needs to hold their breath to maintain the lung volume at its maximum. That is to say, the scan data is preferably the scan data under the patient's breath-holding state.

[0018] (2) Importing Scan Data: After scanning, export the 4D-CT volume package. Each volume package contains 81 sets of files, and each set contains 320 images. The 81 sets of files represent the time series of the scan, with an interval of 0.1s between each set. The 320 images are slices from different detector scanning positions, with a physical interval of 0.5mm between adjacent slices, numbered in head-to-toe order along the axial plane. Import all images into the software and integrate them into a 4D dataset. The horizontal and vertical axes of the axial plane image are defined as X and Y, x , y ∈[1,512]. The number of layers in the slice is defined as Z. z ∈[1,320], according to the apex of the lung ( z =320) to the base of the lungs ( z Arranged in the order of 1. Define the time scale as T. t ∈[1,81].

[0019] Step S200: Obtain the patient's cardiac waveform based on the scan data. As an optional embodiment, step S200 includes: Step S201: Generate a cardiac mask based on the scan data; Furthermore, step S201 includes: Step S2011: Based on the scan data, obtain a three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart; This step defines a three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart. Specifically, step S2011 may include: Step A1: On the axial image of the intermediate slice layer of the scan data at any given time, measure the distance from the anterior edge to the posterior edge of the heart along the anteroposterior direction, and define this distance as the thickness of the heart in the Y-axis direction; During breath-holding (or within one heartbeat cycle), the change in thoracic cavity volume with time T can be considered negligible. Since the scan is centered on the heart, the specific implementation can be... t =1, z =160 layers were used to confirm the thickness of the heart in the Y-axis direction. h c1 The reference cross-section. On the axial image of this layer, the distance from the anterior edge to the posterior edge of the heart is measured along the anteroposterior direction, and this distance is defined as the thickness of the heart in the Y-axis direction. .

[0020] Step A2: Determine the sampling range based on the thickness of the heart in the Y-axis direction, and average the scan data along the Y-axis within this range; In specific implementation, To determine the sampling range, the scanned data, i.e., the volumetric data, is averaged along the Y-axis within this range, thus obtaining a new three-dimensional array. This array contains the average grayscale information of the cardiac region in the coronal plane over time.

[0021] Step A3: Based on the ratio of slice thickness to axial plane pixel spacing, the averaged three-dimensional array is scaled proportionally in the Z direction to obtain the three-dimensional array of the thoracic coronal plane containing the heart thickness.

[0022] In practice, to ensure that the pixel spacing of the image has a consistent physical meaning in different directions, the Z-axis can be scaled proportionally based on the ratio of the slice thickness to the pixel spacing of the axial plane, resulting in an adjusted array. The adjusted array is the three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart.

[0023] Step S2012: Average the three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart along the time dimension to obtain a first time mean image; In practice, the three-dimensional array obtained in the previous step will be used as a reference. The first time-mean image is obtained by averaging along the time dimension. .

[0024] Step S2013: Use the first time mean image as the first background image, and select the heart ROI on it; In specific implementation, As the primary background image, manually draw the heart's ROI (Region of Interest) on it. The ROI should cover as much of the visible area of ​​the heart as possible.

[0025] Step S2014: Generate a binary mask for the heart based on the heart ROI.

[0026] In practice, a binary mask for the heart is generated based on the heart's ROI. The inner ROI mask value is 1, and the outer ROI mask value is 0.

[0027] Thus, through the above steps S2011-S2014, the cardiac mask can be generated effectively.

[0028] Step S202: Draw the initial heart waveform using the heart mask; Furthermore, step S202 includes: The cardiac binary mask is applied to the three-dimensional array of the coronal plane of the thoracic cavity, which includes the thickness of the heart, for each time... t The intensity of all pixels within the heart region (ROI) is averaged to obtain the global average grayscale waveform of the heart region.

[0029] In specific implementation, in a three-dimensional array Cardiac binary masking is used above That is, when When =1, the selected pixel is located within the ROI. When =0, the pixel is outside the ROI. For each time... The global average grayscale waveform of the heart region is obtained by averaging the intensity of all pixels within the heart region ROI. : ; in, Represents each time frame Total number of pixels within the ROI. This is the global average waveform. As a representative signal of overall changes in the cardiac region, it is used as input for subsequent heart rate estimation and the construction of filtered reference signals.

[0030] Step S203: Perform high-pass filtering on the initial cardiac waveform to obtain the drift-free signal; In practice, to reduce low-frequency baseline changes caused by breath-holding residual motion and image intensity drift, the global average grayscale waveform of the heart region is adjusted. Perform the first high-pass processing. Specifically, set the drift time window. And based on the sampling rate =10Hz, convert it to sliding window length The local baseline is estimated using a moving average and subtracted from the original waveform to obtain the drift-free signal. ; This step effectively reduces low-frequency variations below 0.5Hz, thereby preserving higher-frequency components related to heartbeats.

[0031] Step S204: Perform delinearization processing on the drift-free signal to obtain a delinearized signal; In practice, for signals that have already had slow drift removed This may still retain a global trend of slow rise or fall from the beginning to the end of the sequence, which could affect the stability of subsequent frequency domain heart rate estimation. Therefore, the waveform... Further delinearization processing is performed. The overall trend line is fitted using the least squares method and then subtracted from the original waveform. ; in, and These represent the slope and intercept of the fitted trend, respectively. The signal after delinearization. This allows the waveform to fluctuate stably around the zero mean, thereby improving the robustness of subsequent heart rate peak detection and bandpass filtering.

[0032] Furthermore, step S204 may also include: The delinearized signal is subjected to spectral analysis to estimate the heart rate.

[0033] In practice, the cardiac waveform after drift removal and delinearity trend removal processing is... Spectral analysis is performed to estimate the heart rate. First, the signal length is... Zero-padded to a value not less than The length of the smallest power of 2 (when T When the length is 81, =128), and calculate L The point is subjected to a Fast Fourier Transform (FFT) to obtain the frequency domain representation. Further squaring yields the power spectrum. ; Constructing the frequency axis And search for the frequency point corresponding to the maximum power. Let the index corresponding to this maximum power be . Then the main peak frequency (i.e., heart rate, Hz) is It can also be further converted into heart rate. .

[0034] Step S205: Perform bandpass filtering on the delinearized signal.

[0035] This step involves adaptive bandpass filtering and cardiac reference waveform construction. Specifically, step S205 may include: Step B1: Construct an adaptive passband range centered on the heart rate, design a third-order Butterworth bandpass filter, and perform zero-phase filtering on the delinearized signal to obtain the bandpass-processed cardiac signal; In practice, to further extract the periodic components related to heartbeat, heart rate is used as the starting point. Construct an adaptive passband range around the center Its lower and upper limits are defined as follows: ; in, Based on this, a third-order Butterworth bandpass filter was designed and used to process the de-stressed cardiac waveform. Zero-phase filtering is performed to avoid phase delay, resulting in a bandpass-processed cardiac signal. .

[0036] Step B2: Standardize the cardiac signal and take its negative sign to obtain the cardiac waveform.

[0037] In specific implementation, Standardize the waveform and take a negative sign to obtain the cardiac reference waveform (i.e., the cardiac waveform mentioned above): ; in, and They are respectively The mean and standard deviation. The final result is... Used for subsequent correlation analysis with lung waveforms.

[0038] Thus, through the above steps S201-S205, the patient's cardiac waveform can be obtained relatively well.

[0039] Figure 4 The example in the image shows (a) the selection of the cardiac mask and (b) the corresponding initial cardiac waveform. (Represented by dashed lines) Waveform signal after high-pass filtering and delinearization. (represented by solid lines), and (c) the results of heart rate estimation and the normalized cardiac reference waveform after bandpass filtering. .

[0040] Step S300: Based on the scan data, acquire the patient's lung waveform. As an optional embodiment, step S300 includes: Step S301: Generate a thoracic mask based on the scan data; Furthermore, step S301 includes: Step S3011: Based on the scan data, obtain a three-dimensional array of the coronal plane of the thoracic cavity including lung thickness; This step defines a three-dimensional array of the coronal plane of the thoracic cavity containing lung thickness. Specifically, step S3011 may include: Step C1: On the axial image of the intermediate slice layer of the scan data at any given time, measure the distance of the maximum anteroposterior diameter of the lung field along the anteroposterior direction, and define this distance as the thickness of the lung in the Y-axis direction; During breath-holding, the thoracic cavity volume can be considered to change with time. T The changes are negligible. Since the scan is centered on the heart, the specific implementation can be chosen... t =1, z A layer of 160 was used to confirm the thickness of the lung in the Y-axis direction. The reference section. On the axial image of this layer, the maximum anteroposterior diameter of the lung field is measured along the anteroposterior direction, and this distance is denoted as... .

[0041] Step C2: Determine the sampling range based on the thickness of the lung in the Y-axis direction, and average the scan data along the Y-axis within this range; In specific implementation, To determine the sampling range, the scanned data, i.e., the volumetric data, is averaged along the Y-axis within this range, thus obtaining a new three-dimensional array. This array contains the average grayscale information of the lung regions in the coronal plane over time.

[0042] Step C3: Based on the ratio of slice thickness to axial plane pixel spacing, the averaged three-dimensional array is scaled proportionally in the Z direction to obtain the three-dimensional array of the thoracic coronal plane containing lung thickness.

[0043] In practice, to ensure that the pixel spacing of the image has a consistent physical meaning in different directions, the Z-axis can be scaled proportionally based on the ratio of the slice thickness to the pixel spacing of the axial plane, resulting in an adjusted array. The adjusted array is the three-dimensional array of the coronal plane of the thoracic cavity containing lung thickness.

[0044] Step S3012: Average the three-dimensional array of the coronal plane of the thoracic cavity containing lung thickness along the time dimension to obtain a second time-averaged image; In practice, the three-dimensional array obtained in the previous step will be used as a reference. A second time-mean image is obtained by averaging along the time dimension. .

[0045] Step S3013: Use the second time-averaged image as the second background image, and select the thoracic ROI on it; In specific implementation, As a second background image, the thoracic ROI is obtained on a two-dimensional plane, which covers the bilateral lung fields and mediastinal region (including the heart and great vessels).

[0046] Step S3014: Generate a binary mask for the thoracic cavity based on the ROI of the thoracic cavity.

[0047] In practice, the thoracic ROI is used in the second time-averaged image. Automatic segmentation is performed. Specifically, a mask for the torso region is first obtained through threshold segmentation and morphological operations (including hole filling and closing operations). Subsequently in The internal air region is segmented based on intensity threshold to obtain the double lung mask. And preserve the main connected regions of both lungs. To obtain the external boundary of the thoracic cavity including the mediastinum, Calculate the convex hull of the object and its union Intersecting the lines yields the final binary mask for the thoracic cavity. The automatic segmentation results undergo visual quality control, with minor manual corrections made when necessary.

[0048] Thus, through the above steps S3011-S3014, the thoracic mask can be generated effectively.

[0049] Step S302: Draw the initial lung waveform using the pleural mask; Furthermore, step S302 includes: A local time series is constructed for each effective pixel location within the thoracic cavity ROI, i.e., an adjustable time series is set centered on each effective pixel. The local query window is used to average the pixel intensity within the window at each time step, resulting in a local grayscale waveform corresponding to that location. This allows extraction of a waveform of length [value missing] from each location within the thoracic ROI. The local time series forms a spatially distributed set of lung signals.

[0050] In practice, after obtaining the binary mask of the thoracic cavity... Subsequently, to characterize the dynamic changes at different locations within the thoracic cavity region, the effective pixel location within the ROI can be defined. Construct a local time series. Specifically, for each valid pixel... Set an adjustable one at the center. Local query window At each time step, the average pixel intensity within the window is taken to obtain the local grayscale waveform corresponding to that location: ; in, Indicator window The coordinates of any pixel within the range, This represents the number of valid pixels within the window. If the local window intersects with the ROI boundary, then only those satisfying the condition are considered. The pixels are averaged to reduce the influence of tissues outside the ROI on the local waveform. Through this processing, a waveform of length [value missing] can be extracted from each location within the thoracic ROI. Local time series This forms a spatially distributed set of lung signals, which is used for the dynamic characterization of subsequent lung perfusion.

[0051] Step S303: Perform bandpass filtering on the initial lung waveform.

[0052] This step involves bandpass filtering and lung reference waveform construction. Specifically, step S303 may include: Step D1: Construct an adaptive passband range centered on the heart rate, design a third-order Butterworth bandpass filter, and apply zero-phase filtering to each local grayscale waveform to obtain a set of local waveforms after bandpass filtering. In practice, to extract cardiac-related dynamic components within the thoracic region of interest (ROI), the already obtained local lung signal set can be analyzed. Perform bandpass filtering. The passband range can be an adaptive passband obtained from heart rate estimation. A third-order Butterworth bandpass filter was designed based on this passband, and zero-phase filtering was applied to each local waveform to obtain the set of local waveforms after bandpass filtering. .

[0053] Step D2: Standardize the set of local waveforms after bandpass filtering and take a negative sign to obtain the lung waveform.

[0054] In practice, the local waveform set after bandpass filtering is standardized, and a negative sign is used to unify the waveform direction in order to eliminate the influence of amplitude differences at different spatial locations. ; in, and They are respectively The mean and standard deviation are calculated. The standardized local signals are then spatially converged to obtain a reference waveform encompassing any location in both lungs. It can be used for subsequent correlation analysis.

[0055] Thus, by going through the above steps S301-S303, the patient's lung waveform can be obtained relatively well.

[0056] Figure 5The examples in the image show (a) the chest mask (arrows indicate any point in the lung region) and (b) the pixel location indicated by the arrow passing through the same adaptive passband range as the heart waveform. Filtered and normalized lung reference waveform .

[0057] Figure 6 The example in the text shows the heart signal at the same moment and the lung signal at any point, with the two showing opposite trends.

[0058] Step S400: Perform correlation analysis on the heart waveform and lung waveform to obtain the correlation coefficient. As an optional embodiment, step S400 includes: Step S401: Reverse the sign of the heart waveform to align it with the direction of the lung waveform; Continuous images acquired during breath-holding contain periodic changes in pulmonary blood volume caused by cardiac pumping, manifested as subtle temporal grayscale fluctuations in the lungs and mediastinum. Generally, during cardiac systole, the blood volume within the heart decreases, X-ray attenuation within the heart decreases, and the average grayscale of the cardiac region of interest (ROI) decreases; simultaneously, pulmonary artery blood flow increases, leading to enhanced local attenuation within pulmonary vessels and an increase in average grayscale value; the opposite trend is observed during diastole. The local grayscale waveform within the lungs should have a certain temporal correlation with the cardiac reference waveform, and there may be a phase delay. To quantitatively characterize the temporal correlation between the cardiac signal and the local signals within the lungs, normalized cross-correlation analysis was used to analyze the correlation between the cardiac reference waveform and waveforms at various locations within the thoracic ROI. Considering that the cardiac ROI and pulmonary signals may exhibit opposite directions of change during systole and diastole, the cardiac reference waveform was sign-flipped to align with the direction of the pulmonary waveform.

[0059] ; This corresponds to the cardiac reference waveform constructed above.

[0060] Step S402: For a lung waveform at any location within the thoracic ROI, calculate its time delay compared to the sign-flipped heart waveform. The correlation coefficient is given below, where the correlation coefficient is based on heart rate. Heart cycle Set offset time Corresponding to the number of offset frames , The sampling rate; in the offset frame number FL Cross-correlation calculations were performed on each effective pixel location within the thoracic ROI to obtain a correlation coefficient sequence.

[0061] In practice, a local reference waveform at any location within the thoracic ROI is used. Calculate its relationship with Correlation coefficient: ; in, and This represents the mean of the corresponding sequence. Given the periodicity of cardiac signals, the cross-correlation function may exhibit approximate peaks at multiple cycles. To obtain a stable and unique delay estimate, the delay search range is limited to one cardiac cycle. Based on heart rate... Heart cycle Set offset time Corresponding to the number of offset frames At the offset frame number Cross-correlation calculations were performed on each effective pixel location within the thoracic cavity ROI to obtain a correlation coefficient sequence. . It represents the coupling strength between local lung signals and cardiac reference signals, and this index can be used to characterize the temporal propagation characteristics of dynamic perfusion in the lungs.

[0062] Step S500: Output the time series lung perfusion color image based on the correlation coefficient. As an optional embodiment, step S500 includes: Step S501: Construct a color map based on the correlation coefficient sequence, wherein the second time mean image is used as the grayscale background, and the result is restricted to the thoracic cavity binary mask, while only the positively correlated components are retained for perfusion display; In this step, to visualize lung perfusion signals, the correlation coefficient strength is obtained based on cross-correlation analysis. Construct a color map. Using the second time-mean image. As a grayscale background, the results are restricted to a binary mask of the thoracic cavity. Internally, since the aforementioned steps have already performed sign flipping on the cardiac reference waveform to ensure that it changes in the same direction as the signal in the lungs, only the positively correlated components are retained for perfusion visualization.

[0063] Step S502: At the same time, FL When the correlation coefficients of pixels at different locations are the same, they are normalized. The normalized correlation coefficients are then used as perfusion weights and mapped onto the red intensity. This is then superimposed on the second background image to obtain a color overlay image of lung perfusion. In practice, it will be implemented at the same time (i.e. When they are the same, they correspond to pixels at different positions. Normalization is performed, that is: ; The pixel position with the highest correlation after normalization Other pixel positions will proportionally fall within [0,1]. As an infusion weight, it is mapped to the red intensity and superimposed on the second background image. This process yields a color overlay image of lung perfusion. The redder the color, the stronger its correlation with the cardiac reference signal.

[0064] Step S503: Calculate the offset frame number using the same method. FL The perfusion weights at various locations within the perfusion area are determined, and time-series perfusion images are generated.

[0065] In practice, the same method is used to further calculate the number of offset frames. FL The perfusion weights at various locations within the lungs are calculated, and time-series perfusion images are generated. The variations in the red areas and the distribution of the vascular tree-like pattern represent the propagation process of cardiac-related perfusion in the pulmonary vessels and lung parenchyma.

[0066] Figure 7 The example in the image shows color images of lung perfusion at two time points, where (a) shows the second background image and (b) represents the lung perfusion during diastole. FL The perfusion image at =0 shows that the pixels with the highest correlation intensity within the entire thoracic ROI are all located in the cardiac region itself, and are colored red in the image; (c) indicates the perfusion image during cardiac systole, corresponding to FL The perfusion image at =4 shows a large influx of blood into the pulmonary artery. The pixels with the highest correlation intensity move from the heart to the pulmonary vessels, which is represented by the flow in the red area.

[0067] In summary, this invention visualizes lung perfusion using color dynamic perfusion imaging. Its advantages include: This method can directly utilize conventional breath-hold 4D-CT volumetric data, automatically obtaining the intensity distribution and phase delay of perfusion signals related to cardiac activity within the thoracic region of interest (ROI) based on cardiac reference waveform extraction, heart rate adaptive bandpass filtering, and pixel grayscale cross-correlation analysis, without requiring additional hardware or examination procedures. It then generates a time-lapsed color overlay image, intuitively presenting the propagation characteristics of lung perfusion. Compared to existing perfusion imaging methods that rely on radioactive tracers or contrast agents, this method eliminates the need for contrast agent injection, completely avoiding contrast agent-related risks and offering better safety and repeatability. Furthermore, the overall computational process is highly automated and has a short detection time, improving examination and assessment efficiency to some extent. It is suitable for applications requiring rapid assessment or large-scale screening, such as emergency departments, providing reliable image information for screening diseases such as pulmonary embolism or chronic thromboembolic pulmonary hypertension.

[0068] The above embodiments are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A lung perfusion visualization method based on 4D-CT, characterized in that, include: Step S100: Acquire the patient's lung 4D-CT scan data; Step S200: Obtain the patient's cardiac waveform based on the scan data; Step S300: Obtain the patient's lung waveform based on the scan data; Step S400: Perform correlation analysis on the heart waveform and lung waveform to obtain the correlation coefficient; Step S500: Output the time series lung perfusion color image based on the correlation coefficient.

2. The method according to claim 1, characterized in that, In step S100, the scan data is the scan data taken while the patient is holding their breath.

3. The method according to claim 1, characterized in that, Step S200 includes: Step S201: Generate a cardiac mask based on the scan data; Step S202: Draw the initial heart waveform using the heart mask; Step S203: Perform high-pass filtering on the initial cardiac waveform to obtain the drift-free signal; Step S204: Perform delinearization processing on the drift-free signal to obtain a delinearized signal; Step S205: Perform bandpass filtering on the delinearized signal.

4. The method according to claim 3, characterized in that, Step S201 includes: Step S2011: Based on the scan data, obtain a three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart; Step S2012: Average the three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart along the time dimension to obtain a first time mean image; Step S2013: Use the first time mean image as the first background image, and select the heart ROI on it; Step S2014: Generate a binary mask for the heart based on the heart ROI.

5. The method according to claim 4, characterized in that, Step S2011 includes: On the axial image of the intermediate slice layer of the scan data at any given time, the distance from the anterior edge to the posterior edge of the heart is measured along the anteroposterior direction, and this distance is defined as the thickness of the heart in the Y-axis direction; The sampling range is determined based on the thickness of the heart in the Y-axis direction, and the scan data is averaged along the Y-axis within this range; Based on the ratio of slice thickness to axial plane pixel spacing, the averaged three-dimensional array is scaled proportionally in the Z direction to obtain the three-dimensional array of the thoracic coronal plane containing the heart thickness. And / or, step S202 includes: The heart binary mask is applied to the three-dimensional array of the coronal plane of the thoracic cavity containing the thickness of the heart. For each time t, the intensity of all pixels within the heart ROI is averaged to obtain the global average grayscale waveform of the heart region. And / or, step S204 includes: Spectral analysis is performed on the delinearized signal to estimate the heart rate; And / or, step S205 includes: An adaptive passband range was constructed with the heart rate as the center, a third-order Butterworth bandpass filter was designed, and the signal after delinearization was subjected to zero-phase filtering to obtain the bandpass processed cardiac signal. The cardiac signal is standardized and its sign is negative to obtain the cardiac waveform.

6. The method according to claim 1, characterized in that, Step S300 includes: Step S301: Generate a thoracic mask based on the scan data; Step S302: Draw the initial lung waveform using the pleural mask; Step S303: Perform bandpass filtering on the initial lung waveform.

7. The method according to claim 6, characterized in that, Step S301 includes: Step S3011: Based on the scan data, obtain a three-dimensional array of the coronal plane of the thoracic cavity including lung thickness; Step S3012: Average the three-dimensional array of the coronal plane of the thoracic cavity containing lung thickness along the time dimension to obtain a second time-averaged image; Step S3013: Use the second time-averaged image as the second background image, and select the thoracic ROI on it; Step S3014: Generate a binary mask for the thoracic cavity based on the ROI of the thoracic cavity.

8. The method according to claim 7, characterized in that, Step S3011 includes: On the axial image of the intermediate slice layer of the scan data at any given time, the distance between the anteroposterior diameter of the lung field is measured along the anteroposterior direction, and this distance is defined as the thickness of the lung in the Y-axis direction. The sampling range is determined based on the thickness of the lung in the Y-axis direction, and the scan data is averaged along the Y-axis within this range; Based on the ratio of slice thickness to axial plane pixel spacing, the averaged three-dimensional array is scaled proportionally in the Z direction to obtain the three-dimensional array of the thoracic coronal plane containing lung thickness. And / or, step S302 includes: A local time series is constructed for each effective pixel location within the thoracic cavity ROI, i.e., an adjustable time series is set centered on each effective pixel. A local query window is used to average the pixel intensity within the window at each time step, resulting in a local grayscale waveform corresponding to that location. This allows extraction of a waveform of length [value missing] from each location within the thoracic region of interest (ROI). The local time series, thus forming a spatially distributed set of lung signals; And / or, step S303 includes: An adaptive passband range was constructed with the heart rate as the center, a third-order Butterworth bandpass filter was designed, and zero-phase filtering was applied to each local grayscale waveform to obtain a set of local waveforms after bandpass filtering. The set of local waveforms after bandpass filtering is standardized and negative to obtain the lung waveform.

9. The method according to any one of claims 1-8, characterized in that, Step S400 includes: Step S401: Reverse the sign of the heart waveform to align it with the direction of the lung waveform; Step S402: For a lung waveform at any location within the thoracic ROI, calculate its time delay compared to the sign-flipped heart waveform. The correlation coefficient is given below, where the correlation coefficient is based on heart rate. Heart cycle Set offset time Corresponding to the number of offset frames , The sampling rate; in the offset frame number FL Cross-correlation calculations were performed on each effective pixel location within the thoracic ROI to obtain a correlation coefficient sequence.

10. The method according to claim 9, characterized in that, Step S500 includes: Step S501: Construct a color map based on the correlation coefficient sequence, wherein the second time mean image is used as the grayscale background, and the result is restricted to the thoracic cavity binary mask, while only the positively correlated components are retained for perfusion display; Step S502: At the same time, FL When the correlation coefficients of pixels at different locations are the same, they are normalized. The normalized correlation coefficients are then used as perfusion weights and mapped onto the red intensity. This is then superimposed on the second background image to obtain a color overlay image of lung perfusion. Step S503: Calculate the offset frame number using the same method. FL The perfusion weights at various locations within the perfusion area are determined, and time-series perfusion images are generated.

Images