A spatiotemporal stitching method of dual 4D-CT scans for obtaining whole lung images

Abstract

The present invention discloses a dual 4D-CT scan spatiotemporal stitching method for obtaining a whole lung image, and relates to the technical field of medical image recognition. The method comprises: obtaining a volume package of two 4D-CT scans of the upper lung region and the lower lung region in a time series; performing a spatial correlation analysis on the overlapping area of ​​the two scans to determine the best stitching layer; analyzing the two scans in combination with the time series data according to the determined stitching layer to determine the offset time of the stitching; stitching the images in the volume packages of the two scans according to the determined stitching layer and the offset time. The present invention can accurately locate the stitching position and the time of the scan offset, thereby more stably and reasonably merging the two independent scan data into a new volume package covering the whole lung field of view, providing a relatively more convenient and robust method to obtain complete lung 4D-CT image information.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image recognition, and in particular to a dual 4D-CT scanning spatiotemporal splicing method for obtaining a whole lung image. Background Art

[0002] Currently, to our knowledge, there is no automated technology specifically for stitching upper and lower lung volume scans.

[0003] During natural inhalation and exhalation, the diaphragm and intercostal muscles contract and relax, causing the chest cavity and lung volumes to expand and retract regularly. Therefore, the stitching of 4D-CT lungs involves stitching in two dimensions: space and time. At present, commercialized medical image stitching technologies, such as X-ray-based spine image stitching, usually only involve the spatial merging of multiple sets of cone beam computed tomography scans, because the stitching target is rigid structures such as bones, and these tissues do not change over time. The widely studied multimodal image registration usually refers to comparing images obtained at the same location at different times or using different scanning methods, and aligning them to the same coordinate system. It is not targeted at the lungs, nor does it involve spatial stitching, which is different from the purpose of this technology.

[0004] In the past, 4D-CT stitching was performed under the assumption that the external respiratory monitoring device could obtain two complete overlaps during the default scan without additional correction. Although this respiratory monitoring device can provide patients with a visual respiratory curve during the scan, allowing patients to follow their first scan respiratory curve as closely as possible during the second scan. However, despite this, this method cannot ensure that the respiratory phases of the two acquisitions are completely overlapped. Deviations in the scanning phases caused by various reasons are common, and patients often do not have the opportunity to sample multiple times during the actual detection process, which places higher requirements on the stitching of subsequent images. Manual identification and (manual) image stitching require medical workers to observe the original images at different times and depths. This not only requires the stitcher to have rich image recognition experience and be able to correct the offset time of the two scans, but also a large amount of repeated stitching work will consume huge time and labor costs, and may also result in missing feature changes and reduce stitching accuracy. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a dual 4D-CT scan spatiotemporal stitching method for obtaining a whole lung image, so as to provide a relatively more convenient and robust method to obtain complete lung 4D-CT image information.

[0006] A method for spatiotemporal stitching of dual 4D-CT scans for obtaining a whole lung image, comprising:

[0007] Step S101: obtaining a time series of volume packages of two 4D-CT scans of the upper lung region and the lower lung region;

[0008] Step S102: performing spatial correlation analysis on the overlapping area of ​​the two scans to determine the best stitching layer;

[0009] Step S103: analyzing two scans combined with timing data according to the determined splicing layer to determine the offset time of splicing;

[0010] Step S104: stitching the images in the volume packages scanned twice according to the determined stitching layer and offset time.

[0011] The present invention has the following beneficial effects:

[0012] The present invention targets the volume packages of two scans of the patient's upper lung and lower lung, and can accurately locate the position of splicing and the time of offset through image cross-correlation calculation and lung morphology change curve analysis, so as to more stably and reasonably merge the two independent scan data into a new volume package covering the entire lung field of view. This technology can make up for the splicing problem of segmented scanned images caused by the defects in the hardware conditions of the current scanning equipment to a certain extent. Replacing manual image recognition with image analysis technology saves medical workers a lot of time and energy consumption, and provides reliable basic image information for subsequent segmentation and reconstruction of lung tissue, respiratory strain parameter analysis, pulmonary vascular volume value calculation, diagnosis and monitoring of patient lung diseases, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0014] Figure 1 Schematic diagram of the scan volume and overlapping area of ​​two scans S1 and S2 in the present invention;

[0015] Figure 2 Schematic diagram of the scanning window diameter changing with z1 or z2 in the present invention;

[0016] Figure 3 This is a schematic diagram of the overlapping area of ​​two 4D-CT scans in the present invention;

[0017] Figure 4Schematic diagram of normalized cross-correlation analysis of multi-time spatial scales in the present invention, wherein (a) shows the normalized cross-correlation values ​​at t1=1, t2=1, z2=311:300 and z1=11, (b) shows the normalized cross-correlation values ​​of z2=311:300 and t1=1, z1=11 in the entire time series t, and (c) reflects the total correlation proportion score of z2=311:300 and t1=1, z1=11 in the entire time series t;

[0018] Figure 5 It is a schematic diagram of time scale normalized cross-correlation analysis in the present invention, reflecting the change of the difference in time sequence number between z1=11 at each moment and z2=310 with the maximum correlation in the time series t;

[0019] Figure 6 It is a schematic diagram of lung area deformation curve analysis with time series t in the present invention;

[0020] Figure 7 It is a schematic diagram of the analysis of the sternum motion curve along the time series t in the present invention;

[0021] Figure 8 Schematic diagram of time scale structural similarity analysis in the present invention, wherein (a) corresponds to the similarity value of Vs1 and Vs2 in time series t1, t2, and (b) reflects the change of the difference in time sequence number between Vs1 at each moment and the moment when Vs2 with the maximum similarity value is located in the time series t;

[0022] Fig. 9 Schematic diagram of lung volume deformation curve analysis with time series t in the present invention;

[0023] Fig.10 Schematic diagram of the coronal plane of Y=256 and the sagittal plane of X=256 before and after splicing in the present invention, wherein (a) and (b) are coronal plane images of the upper lung S1 and the lower lung S2 before splicing; (c) is the splicing position with z1=11 and z2=310, and the time offset is Δt * Coronal image of the whole lung after stitching; (d) and (e) are sagittal images of the upper lung S1 and lower lung S2 before stitching; (f) is the stitching position with z1=11 and z2=310, and the time offset is Δt * The sagittal view of the whole lung after stitching;

[0024] Fig.11 This is a flow chart of lung 4D-CT image acquisition and splicing of the present invention;

[0025] Fig.12 This is a flow chart of the spatial stitching and temporal stitching algorithm of the image in the present invention. DETAILED DESCRIPTION

[0026] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0027] It should be clear that the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0028] Chronic obstructive pulmonary disease (COPD) is an obstructive lung disease characterized by persistent airflow limitation. It includes two pathological types, chronic bronchitis and emphysema, which usually worsen progressively and are related to the chronic inflammatory response of the airways and lungs to harmful particles or gases. Although the diagnosis of COPD currently relies mainly on spirometry results and respiratory symptoms, spirometry has limitations, such as lack of regional information, and it is difficult to distinguish it from underlying pathological abnormalities, including emphysema and bronchitis.

[0029] In recent years, quantitative measurement of computed tomography (CT) has gradually become an available means of COPD examination. By performing dual-gas phase registration on the patient's lung CT images during the inspiratory and expiratory phases, the limitations of vital capacity measurement can be compensated to a certain extent. Among them, dynamic ventilation CT scan (DVCT) is an imaging technology specifically designed to evaluate lung function and structure. It has shown important clinical application value in the diagnosis and evaluation of chronic obstructive pulmonary disease, pulmonary fibrosis, and lung transplantation. It combines the high-resolution advantages of multi-slice spiral computed tomography (MDCT) and can carefully observe the tiny structures and lesions of the lungs. Through continuous or intermittent scanning, the dynamic changes of the lungs in different respiratory phases (such as inspiration and exhalation) are captured. Based on this, strain analysis has become a new quantitative measurement method for recording lung deformation during ventilation.

[0030] However, for 4D dynamic ventilation CT scans, due to the limitation of the number of detector rows (256 or 320 rows, single maximum distance = 16 cm), a single scan cannot cover the motion monitoring of the entire lung ventilation process, and dynamic scanning of the whole lung often requires two times to complete. The insufficient scanning width hinders the subsequent reconstruction of the three-dimensional volume of the lung segments and the reconstruction of the pulmonary blood vessels, especially for the quantitative identification of abnormal lung motion in high-risk patients and the observation of potential differences in lung deformation in patients with different degrees of airflow limitation. Relying solely on the external respiratory curve during acquisition to guide the patient without correction often results in the failure to obtain accurate and reasonable whole lung scan images due to the deviation of the respiratory phase. The splicing of dynamic lung scans is a technical problem that needs to be solved urgently.

[0031] Due to the limitation of the number of rows of detectors, the scan area is insufficient and the stitching of scanned images is a pain point of the current 4D-CT technology. Based on this problem, the present invention proposes a spatial and temporal series CT stitching technology for the lungs, that is, a dual 4D-CT scan spatiotemporal stitching method for obtaining a full lung image. Figure 11-12 As shown, including:

[0032] Step S101: obtaining a time series of volume packages of two 4D-CT scans of the upper lung region and the lower lung region;

[0033] This step is to import the scan data. Specifically, two scans S1 (scan of the upper lung area) and S2 (scan of the lower lung area) can be completed to export two 4D-CT volume packages.

[0034] Preferably, in this step, each volume package contains 81 sets of files, each set of files contains 320 images. The 81 sets of files represent the time series of the scans, and the interval between each set of files is 0.1s. The 320 images are slices of different detector scanning positions, and the physical interval between adjacent slices is 0.5mm. Each set is numbered in order from head to toe.

[0035] Preferably, this step includes: importing all images into the software, integrating them into two 4D data sets, defining the horizontal and vertical of the axial plane image as x and y, x, y∈[1,512]. The number of slices is defined as z, z∈[1,320]. The images of the first scan S1 are arranged in sequence from the apex of the lung (z1=320) to the middle of the lung (z1=1), and the images of the second scan S2 are arranged in sequence from the middle of the lung (z2=320) to the base of the lung (z2=1). The time series is defined as t, t∈[1,81].

[0036] As an optional embodiment, before step S101, the following steps are included:

[0037] Step S100: Formulate a 4D-CT segmented scanning plan: obtain a simulated positioning CT scan image, and formulate two 4D-CT scanning plans S1 and S2 in the head-foot direction, S1 is a scan covering the upper lung area, and S2 is a scan covering the lower lung area.

[0038] This step is to obtain 4D-CT images in segments.

[0039] Preferably, the step S100 includes:

[0040] The dynamic volume scanning mode was adopted, the maximum single scan length of the 320-row detector was 16 cm, and the overlap area of ​​the two scans was set to at least 1 cm (set to 1 cm in the subsequent examples), so as to ensure that the two 4D-CT images could cover the entire lung;

[0041] And / or, after setting the 4D-CT scanning parameters, the designed dynamic volume scanning plan is transferred to the workstation, and the imaging device and related parameters are selected according to the plan to perform the first 4D-CT scan S1. In specific implementation, the 4D-CT reconstruction parameters can be: tube current = 60mA; tube voltage = 80kV; rotation time = 0.35s; total scanning time = 8.1s; imaging field of view = 320mm; collimation = 0.5mm; slice thickness = 0.5mm; reconstruction kernel = BODY; reconstruction interval = 0.1 / frame (81 frames in total); iterative reconstruction = DLR-AiCE; slice thickness / interval = 0.5mm / 0.5mm. After the scan is completed, keep the patient's position still, and then move the scanning bed to the position of the scan plan S2, and complete the second 4D-CT scan with the same parameters.

[0042] Step S102: performing spatial correlation analysis on the overlapping area of ​​the two scans to determine the best stitching layer;

[0043] This step is the analysis of the splicing space position:

[0044] First, determine the spatial splicing position of the two 4D-CT scans, and use the cross-correlation algorithm to determine the matching degree of the spatial overlap of the two scans. Since the two scans are acquired successively without the patient changing position on the treatment bed, and the two scan ranges have a 1cm overlap, the theoretical spatial splicing position should be within the overlapping 1cm area.

[0045] Window effect analysis of 4D-CT:

[0046] 4D-CT scanning has a "window effect" phenomenon, that is, at the beginning and end of each scan, the scanning range of the detector fails to fully reach the maximum scale, and the scanning volume presents a cone-shaped distribution, such as Figure 1 shown.

[0047] The scanning window of the axial plane is circular. In the conical area, the window diameter of the layer with the smallest scanning range is 18 pixels. When the window diameter increases to 352 pixels, the complete lung section can usually be displayed, and it continues to increase until the maximum field of view diameter reaches 512 pixels. The relationship between the window diameter and the scanning slice height is as follows: Figure 2 As shown, when z1, z2 = 320 and z1, z2 = 1, the diameter of the window is the smallest (d = 18 pixels), and when z1 or z2 = 306, z1 or z2 = 15, the diameter of the window reaches the maximum (d = 512 pixels).

[0048] Figure 3 The overlapping area of ​​two 4D-CT scans is shown. The dark color represents the largest field of view, the gray color represents the field of view is not the largest but can display the complete lung section, and the light color represents the field of view is not enough to display the complete lung section. If the window range is too small or the field of view is incomplete, the scan information may be incomplete, making it impossible to reconstruct the location. When stitching, you should try to avoid selecting such data. Figure 3 In the figure, for the 1cm overlap area, it includes z1=21:1 in S1 and z2=320:300 in S2. The scanning window diameters of z1=21:15 in the upper lung and z2=306:300 in the lower lung are the maximum. The window diameters of the scanning layers z1=14:10 and z2=311:307 are not the maximum but can display the complete lung section. The scanning layers z1=9:1 and z2=320:312 cannot cover the complete lung section. Therefore, the reference plane for stitching is selected in the range of z1=12:10 and z2=311:309.

[0049] As an optional embodiment, in step S102, the correlation analysis adopts multi-time spatial scale normalized cross-correlation analysis, which is specifically described as follows:

[0050] Normalized cross-correlation analysis of multiple time spatial scales:

[0051] Two-dimensional normalized cross-correlation is an algorithm to measure the similarity between two images. According to the image size, the cross-correlation of the two-dimensional matrix of the image is calculated. The covariance of the two images is divided by the product of the standard deviations of the two signals to normalize the cross-correlation to obtain the correlation coefficient. The calculation formula of the two-dimensional normalized cross-correlation γ(u,v) is:

[0052]

[0053] Among them, f(x,y) is the reference image, t(xu,yv) is the pixel value of the template image at position (u,v) in the reference image, is the mean of the reference image at (u,v), is the mean of the template image.

[0054] Preferably, the multi-time spatial scale normalized cross-correlation analysis is:

[0055] A centered image of a scanning layer in an overlapping area with a certain time sequence number is selected from the images of the first scan (in the aforementioned example, the centered image in the overlapping area whose field of view is not the largest but can display the complete lung cross-section, that is, the centered image z1=11 with z1=12:10) as the stitching reference layer, and all scanning layers in the overlapping area of ​​the second scan within the entire scanning time range that can display the complete lung are used as layers to be analyzed. A normalized cross-correlation analysis is performed on the stitching reference layer and the layer to be analyzed, and the layer to be analyzed with the highest correlation score is used as the best stitching layer corresponding to the stitching reference layer.

[0056] In specific implementation, for example, the layer with t1=1, z1=11 can be selected as the stitching reference section of the upper lung S1, and a mask with a radius of 176 pixels is established with the center of the image (corresponding to the window sizes of z1, z2=311 and z1, z2=10), and normalized cross-correlation analysis is performed with all layers with t2=1:81, z2=311:300 in the lower lung S2, traversing the overlapping area of ​​S2 within the entire scanning time range. Figure 4 (a) shows the normalized cross-correlation values ​​of different z2 and the reference layer when t1=t2=1, and z2=310 has the highest correlation peak corresponding to z1=11. (b) shows the change of normalized cross-correlation values ​​of different z2 and the reference layer in the entire time series. Sort the normalized cross-correlation values ​​at each moment, list the 5 values ​​with the largest correlation and their corresponding layers, and define the correlation score according to the ranking (5 points for the first, 4 points for the second, 3 points for the third, 2 points for the fourth, and 1 point for the fifth), as shown in Figure 4 (c) shows the total score of each layer of z2 in the entire time series. Thus, we verify the accuracy of the scanned spatial position. If the layer with the highest score is z2=310, then z2=310 is the best concatenated layer in S2 corresponding to z1=11 in S1; if the layer with the highest score is z2=z max , then z2=z max The best stitching layer in S2 corresponding to z1=11 in S1. Due to the high precision of the movement of the CT scanning bed, the significance of this method is to verify whether the patient has active large-scale body position movement during the scanning process, which leads to the deviation of the scanning position. In the 4D-CT scanning data of multiple groups of patients we verified, the results of the multi-time spatial scale normalized cross-correlation analysis are consistent with the theoretical scanning position, that is, z2=310 is the best stitching layer for z1=11.

[0057] Figure 4The multi-time spatial scale correlation analysis is shown. (a) shows the normalized cross-correlation value of t1=t2=1, z2=311:300 and z1=11. (b) shows the normalized cross-correlation value of z2=311:300 and t1=1, z1=11 in the entire time series. (c) reflects the total correlation score of each layer of z2=311:300 and t1=1, z1=11 in the entire time series, among which z2=310 is the layer with the highest score.

[0058] Step S103: analyzing two scans combined with timing data according to the determined splicing layer to determine the offset time of splicing;

[0059] This step is the analysis of splicing offset time:

[0060] The lungs of the human body contract and expand repeatedly during the entire respiratory cycle as the heart beats. Since the two scans are collected one after another, even if there is a respiratory curve as a reference during the scan, it is impossible to guarantee that the respiratory phases of the two acquisitions are completely consistent. The time delay of the same respiratory phase during the entire scan must also be considered in the stitching process. After determining the spatial stitching position of the upper and lower lungs S1 and S2, we proposed a method to determine the stitching time offset.

[0061] As an optional embodiment, in step S103, the analysis adopts time-scale normalized cross-correlation analysis, upper and lower lung area deformation curve analysis, sternum motion curve analysis, time-scale structural similarity analysis and / or upper and lower lung volume deformation curve analysis, which are specifically described as follows:

[0062] Time scale normalized cross-correlation analysis:

[0063] Preferably, the time scale normalized cross-correlation analysis is:

[0064] First, in the entire scanning time range, the stitching reference layer z1 and its optimal stitching layer z2 are traversed, and for each time sequence number t 1i Under z1, sort z2 at different times according to the normalized cross-correlation value, and find the time sequence number t of z2 with the largest correlation. 2j , and calculate the time sequence difference with the corresponding z1, which is defined as the delay Δt = t 2j -t 1o ; Then, the interquartile range method is used to exclude outliers for all Δt, and the average value of Δt after excluding outliers is calculated and rounded to get Δt c That is the offset time.

[0065] In specific implementation, for example, the cross-correlation calculation can be performed on the time scale of the entire scanning period, and all z1=11 and z2=310 can be traversed in t1=1:81 and t2=1:81 respectively. 1i Under z1, sort z2 at different times according to the normalized cross-correlation value, and find the time sequence number t of z2 with the largest correlation. 2j , and calculate the time sequence difference with the corresponding z1, which is defined as the delay Δt = t 2j -t 1i .

[0066] Then, the interquartile range method is used to exclude outliers for all Δt, that is, Δt is divided into four equal numerical points from small to large, Q1 is the quartile of the smallest 25%, Q3 is the quartile of the largest 25%, and the interquartile range IQR = Q3-Q1 is calculated, and the lower limit lim of the outlier is set lower =Q1-1.5×IQR, upper limit lim upper =Q3+1.5×IQR, the values ​​in the data that are less than the lower limit or greater than the upper limit are considered as outliers and removed from the data set. Calculate the average value of the filtered Δt and round it off, Δt c The concatenated time delays for upper and lower lung S1 and S2 were derived using the time-normalized cross-correlation algorithm.

[0067] Figure 5 The time scale normalized cross-correlation analysis is shown, reflecting the change in the time series t of the difference between z1=11 and the time z2=310 with the maximum correlation at each moment. The quartile method is used to remove outliers and take the mean. In this case, Δt c =-6.

[0068] Analysis of deformation curves of upper and lower lung areas:

[0069] Preferably, the upper and lower lung area deformation curve analysis is:

[0070] First, for the volume packages of the two scans, coronal images with a depth of 256 are obtained in the entire scanning time range; all coronal images are contrast enhanced and binarized; the images are further eroded and filtered to eliminate the white areas of small blood vessels and trachea in the lung area, fill the small gaps in the image, mark the connected domains of the image, and calculate the area size of the left and right lung areas at different times; the change in the sum of the areas of both lungs is normalized to draw a deformation curve on a time series scale; finally, the horizontal offset of the upper and lower lung curves is calculated using the least squares difference algorithm based on the deformation curve to obtain Δt ad That is the offset time.

[0071] In specific implementation, for example, the volume packages of S1 and S2 can be used to obtain coronal images of Y=256 at t1=1:81 and t2=1:81, respectively. Contrast enhancement and binarization are performed on all coronal images. The image is further eroded and filtered to eliminate scattered white areas such as small blood vessels and trachea in the lung area, fill small gaps in the image, mark the connected domains of the image, and calculate the area size of the left and right lung areas at different times. The change in the sum of the areas of both lungs is normalized to draw the deformation curve on a time series scale. We believe that the deformation characteristics of the lungs are unified within a respiratory cycle, so the horizontal offset of the two curves of the upper and lower lungs corresponds to the time delay of the splicing. Finally, the least squares difference algorithm is used to calculate the horizontal offset Δt of the two curves. ad . Least square difference algorithm formula:

[0072] E(Δt)=Σ t [A1(t)-A2(t+Δt)] 2

[0073] Δt ad =argmin Δt E(At)

[0074] Where A1(t) is the function of the upper lung area over time t, A2(t) is the function of the lower lung area over time t, and Δt ad This is the offset that minimizes E(Δt).

[0075] Figure 6 The lung area deformation curve analysis is shown, reflecting the change curve of the lung area in a patient's breathing cycle. At the peak of the patient's inhalation, the lung expansion corresponds to the largest lung area, and the contraction corresponds to the smallest lung area. The least squares difference algorithm is used to calculate the offset time, Δt in this case ad =-7.

[0076] Sternum motion curve analysis:

[0077] Preferably, the sternum motion curve analysis is:

[0078] First, for the volume packages of the two scans, sagittal images with a width of 256 are obtained in the entire scanning time range; contrast enhancement and binarization are performed on all sagittal images; the images are eroded and filtered to eliminate the white areas of small blood vessels and trachea in the lung area and fill the small gaps in the image; the pixel points of the outermost sternum with a depth of 310 in the sagittal image of the first scan and a depth of 11 in the sagittal image of the second scan are marked, and the changes in the coordinate position in the time series are stored; then the change curve of the horizontal coordinates of the outer side of the sternum in the time series is drawn; finally, the horizontal offset of the motion curve of the upper and lower sternum is calculated using the least square difference algorithm based on the change curve to obtain Δt sd That is the offset time.

[0079] In specific implementation, for example, the sagittal images of X=256 can be obtained for the volume packages of S1 and S2 at t1=1:81 and t2=1:81, respectively. For all sagittal images, contrast enhancement and binarization are performed, that is, segmentation and image processing (erosion, filtering and filling) are performed using the same morphological method as the upper and lower lung area deformation curve analysis. Mark the pixels of the outermost sternum at Y=310 in the sagittal image of S1 and Y=11 in the sagittal image of S2, and store the changes in the coordinate position in the time series. Draw a curve of the change of the horizontal coordinates of the outer side of the sternum on the time series. We believe that the ups and downs of the sternum in a respiratory cycle is a significant expression of the gas phase transformation of the lungs. Due to the bony characteristics, the upper and lower sternum should be completely connected, and the movement trajectory should be coherent and unified. Therefore, the horizontal offset of the movement curve of the upper and lower sternum corresponds to the time offset of the splicing. Finally, the least squares difference algorithm is used to calculate the horizontal offset Δt of the two curves. sd .

[0080] Figure 7 The sternum motion curve analysis is shown, reflecting the horizontal abscissa change curve of the outermost sternum pixel point in a patient's breathing cycle. When the patient is at the peak of inspiration, the chest bulge corresponds to the largest horizontal pixel. The least square difference algorithm is used to calculate the offset time, Δt in this case. sd =-6.

[0081] Time scale structural similarity analysis:

[0082] Structural similarity (SSIM) is another metric used to measure the similarity between two images. It not only involves the difference in pixel intensity of the image, but also considers the structural information such as brightness and contrast of the image, and pays more attention to the changes in the local structure of the image. When used to compare three-dimensional data, we calculate the structural similarity values ​​for multiple layers of images and take the average. The calculation formula for structural similarity SSIM (x, y) is:

[0083]

[0084] Among them, μ x and μ y is the image average, and is the image variance, σ xy is the image covariance, C1 and C2 are two stable constants.

[0085] Preferably, the time scale structural similarity analysis is:

[0086] First, for the volume packages of the two scans, the images in the overlapping area where the field of view is not the largest but can show the complete lung section are selected as the overlapping volume V of the upper and lower lungs to be analyzed. S1 、V S2 ; In the entire scanning time range, traverse all V S1 and V S2 , for each time sequence number t 1i V S1 , according to the structural similarity value, V at different times S2 Sort and find the V with the largest similarity value S2 The time sequence number t 2j , and calculate the corresponding V S1 The time sequence number difference is defined as the offset Δt = t 2j -t 1i ; Then, the interquartile range method is used to exclude outliers for all Δt, and the average value of Δt after excluding outliers is calculated and rounded to get Δt s That is the offset time.

[0087] In a specific implementation, for example, a structural similarity calculation can be performed on the time scale of the entire scanning period, and for the volume packages of S1 and S2, z1=10, 11, 12 and z2=309, 310, 311 are selected as the overlapping volumes V of the upper and lower lungs to be analyzed, respectively. S1 、V S2 . Traverse all V in t1 = 1:81, t2 = 1:81 S1 With V S2 For each time series t 1i V S1 , according to the structural similarity value, V at different times S2 Sort and find V with the largest similarity value S2 The time sequence number t 2j , and calculate the corresponding V S1 The time sequence number difference is defined as the delay Δt = t 2j -t 1i .

[0088] Then, the interquartile range method is used to exclude outliers for all Δt, that is, the values ​​in the data that are less than the lower limit or greater than the upper limit are removed from the data set. The average value of the filtered Δt is calculated and rounded to the nearest integer, Δt s Time offsets of the splicing using a temporal structural similarity algorithm for upper and lower lung S1 and S2.

[0089] Figure 8 The time scale structural similarity analysis is shown, where (a) corresponds to V at t1 and t2. S1 With V S2 The similarity value of (b) reflects the V at each moment S1 , with which V has the largest similarity value S2 The difference between the two time series at the moment of time changes in the time series t. Use the quartile method to remove outliers and take the mean. In this case, Δt s =-6.

[0090] Analysis of upper and lower lung volume deformation curves:

[0091] Preferably, the upper and lower lung volume deformation curve analysis is:

[0092] First, for the volume packages of the two scans, the entire scanning time range is traversed respectively, and the starting coronal layer is selected. The depth Y of the starting coronal layer i ∈[160,180], with a certain depth interval, a total of at least 15 coronal layers are analyzed to cover the complete lung volume; all coronal images are contrast enhanced and binarized, and the images are further eroded and filtered to eliminate the white areas of small blood vessels and trachea in the lung area, fill the small gaps in the image, mark the connected domains of the image, and calculate the area sizes of the left and right lung areas at different times; the sum of the areas of at least 15 coronal layers of both lungs at the same time is defined as the volume eigenvalue at that time, and the volume eigenvalues ​​of all time series are calculated; the changes in the volumes of both lungs are normalized to draw the deformation curve on the time series scale; finally, the least square difference algorithm is used to calculate the horizontal offset of the two curves based on the deformation curve to obtain Δt vd That is the offset time.

[0093] In specific implementation, for example, lung volume analysis characterized by the sum of multiple coronal lung areas can be used to obtain more accurate results than single-layer area analysis. For the volume packages of S1 and S2, traverse in t1 = 1:81 and t2 = 1:81 respectively, and select the starting coronal layer Y i , Y i∈[160,180], with a certain depth (10 in this case) as an interval, a total of at least 15 (19 in this case) coronal layers are analyzed to cover the complete lung volume. For all coronal images, contrast enhancement and binarization are performed, and segmentation and image processing are performed using the same morphological method as the upper and lower lung area deformation curve analysis to mark the connected domains. The sum of the areas of both lungs in the 19 coronal layers at the same moment is defined as the volume eigenvalue at that moment, and the volume eigenvalues ​​of all time series are calculated. The changes in the volumes of both lungs are normalized to draw the deformation curves on the time series scale. Finally, the least squares difference algorithm is used to calculate the horizontal offset Δt of the two curves. vd .

[0094] Fig. 9 The lung volume deformation curve analysis is shown, reflecting the change curve of the lung volume in a patient's breathing cycle. At the peak of inspiration, the lung expansion corresponds to the largest lung volume, and the exhalation contraction corresponds to the smallest lung volume. The offset time calculated using the least squares difference algorithm, Δt in this case vd =-3.

[0095] As another optional embodiment, step S103 includes:

[0096] Step S1031: According to the determined splicing layer, the two scans are combined with the time series data for analysis, that is, the time scale normalized cross-correlation analysis, the upper and lower lung area deformation curve analysis and the sternum motion curve analysis are used to obtain the offset time Δt c , Δt ad , Δt sd ;

[0097] Step S1032: The offset time Δt c , Δt ad , Δt sd Perform the first offset time Δt * determination.

[0098] This step is to calculate the offset time Δt obtained by the three methods. c , Δt ad , Δt sd In the specific implementation, if the absolute value of the difference between any result and another result is ≤5, a horizontal line is drawn under the two offset times, indicating that the two results are close; in this way, if the results of the three algorithms are all crossed out, the average method is used to obtain the corrected offset time Δt * , Δt * =(Δt c +Δt ad +Δt sd ) / 3; if Δt cIf it is not crossed out, the structural similarity algorithm is used for further verification; if Δt ad If no line is drawn, the upper and lower lung volume deformation curve analysis is used for further verification; if Δt sd If no line is drawn, the result is discarded and the offset time Δt * =(Δt c +Δt ad ) / 2; if the differences between the three are all greater than 5, it is considered that the multi-segment 4D-CT scan is not conducive to splicing, and manual verification should be performed to check whether the original scan image and the respiratory monitoring curve are synchronized. This step is the first offset time Δt * Determination, the cases using structural similarity algorithm and lung volume deformation curve analysis need to be followed by a second offset time Δt * determination.

[0099] As another optional embodiment, the step S1032 includes:

[0100] Step S1033: According to the determined splicing layer, the two scans are analyzed in combination with the time series data, that is, the time scale structural similarity analysis and / or the upper and lower lung volume deformation curve analysis are used to obtain a newly obtained offset time result Δt s or Δt vd ;

[0101] Step S1034: respectively calculate the corresponding Δt s , Δt ad , Δt sd , or Δt c , Δt vd , Δt sd Perform the second offset time Δt * determination.

[0102] This step is the second offset time Δt * In specific implementation, the newly obtained offset time result Δt s (or Δt vd ), respectively for the corresponding Δt s , Δt ad , Δt sd (or Δt c , Δt vd , Δt sd ) is used for judgment. If the absolute value of the difference between any result and the other result is ≤5, a horizontal line is drawn under the two offset times to indicate that the two results are close. In this way, if the results of the three algorithms are all crossed out, the average method is used to obtain the corrected offset time Δt * , Δt * =(Δt s +Δtad +Δt sd ) / 3(or Δt * =(Δt c +Δt vd +Δt sd ) / 3); if Δt s If no line is drawn, the result is discarded and the offset time Δt * =(Δt ad +Δt sd ) / 2; if Δt vd If no line is drawn, the result is discarded and the offset time Δt * =(Δt c +Δt sd ) / 2.

[0103] Step S104: stitching the images in the volume packages scanned twice according to the determined stitching layer and offset time.

[0104] This step is data integration and image stitching. Preferably, step S104 includes:

[0105] The 4D-CT scan data of the two scans are integrated according to the corresponding time sequence numbers. The image data with the time sequence number M in S1 is selected, and the image data with the time sequence number M+Δt in S2 is selected. * Image data, and satisfy 1≤M≤81 and 1≤M+Δt * ≤81; for each group in S1, select the stitching reference layer as the stitching section, and for each group in S2, select the best stitching layer as the stitching section; compare the image data from the lung apex scanning layer to the stitching reference layer in the S1 time series numbered M with the corresponding S2 time series numbered M+Δt * The image data from the best stitching layer to the lung bottom scanning layer are integrated and stored in a new folder. The folders are numbered in chronological order, and the layers in each new folder are numbered in order from head to toe.

[0106] In specific implementation, for example, the 4D-CT scan data of the upper lung S1 and the lower lung S2 can be integrated with the corresponding time sequence numbers, and the image data with the time sequence number M in S1 is selected, and the image data with the time sequence number M+Δt in S2 is selected. * Image data, and satisfy 1≤M≤81 and 1≤M+Δt * ≤81. For each group in S1, select z1=11 as the splicing section, and for each group in S2, select z2=310 as the splicing section; compare the image data of z1=320:11 in the S1 time series numbered M with the corresponding S2 time series numbered M+Δt *The image data of z2=310∶1 in the image are integrated and stored in a new folder. The folders are numbered in chronological order, and the layers in each new folder are numbered in order from head to toe. So far, each folder at a certain moment will contain 620 images, which contains all the information of the dynamic 4D-CT scan of the whole lung of the patient in one respiratory cycle.

[0107] Fig.10 The Y=256 coronal plane and X=256 sagittal plane before and after splicing are shown. (a) and (b) are coronal plane images of the upper lung S1 and the lower lung S2 before splicing; (c) is the splicing position with z1=11 and z2=310, and the time offset is Δt * (d) and (e) are sagittal images of the upper lung S1 and lower lung S2 before splicing; (f) is the splicing position z1 = 11 and z2 = 310, and the time offset is Δt * Sagittal view of the whole lung after splicing.

[0108] Fig.10 The two scans S1 and S2 are shown. The coronal plane (Y=256) and sagittal plane (X=256) of the upper and lower lungs without time offset before stitching. After the spatial stitching position is confirmed and the first offset time is determined, z1=11 and z2=310 are used as the stitching section, and Δt c , Δt ad , Δt sd The average value of Δt * As the offset time, the stitched whole lung coronal and sagittal images were integrated.

[0109] The process of segmented acquisition and splicing of 4D-CT lung images is as follows: Fig.11 As shown, the algorithm flow chart of spatial splicing and temporal splicing is as follows Fig.12 shown.

[0110] In summary, the existing 4D-CT splicing for the lungs, in addition to locating the spatial splicing position of multiple scans, also needs to determine the time offset of lung deformation. The previous 4D-CT splicing depends on the breathing curve provided by the external respiratory monitor during scanning, and controls the breathing rhythm of the two scans by the patient's subjective consciousness, and assumes that the two curves completely overlap. Although the monitor can provide real-time respiratory phase display, it is difficult for patients to make two completely identical breaths. Without advance training and repeated experiments, the phase offset of the two breaths is ubiquitous. To eliminate this offset, medical personnel can only manually distinguish from a large amount of image data, which consumes a lot of energy and time. The present invention aims to alleviate the pain points of this medical technology to a certain extent. Through the provided method, including the cross-correlation calculation of the image and the analysis of the lung deformation curve, the spatial position and time offset of the splicing can be determined more quickly and accurately, providing a more stable and automated 4D-CT splicing process. The output volume package containing the whole lung scan provides reliable image information for subsequent three-dimensional reconstruction of the lungs and blood vessels, strain analysis, etc.

[0111] For the volume packages of the two scans of the patient's upper and lower lungs, the normalized cross-correlation calculation of the images, the analysis of the lung deformation curve, and the analysis of the sternal motion curve can accurately locate the stitching position and the offset time, so as to more stably and reasonably merge the two independent scan data into a new volume package covering the entire lung field of view. This technology can to a certain extent make up for the problem of stitching segmented scan images caused by the defects in the hardware conditions of the current scanning equipment. The stitching process also enhances the reliability of the respiratory monitoring curve during scanning. Replacing manual image recognition with image analysis technology saves a lot of time and energy consumption for medical workers, and provides reliable basic image information for subsequent segmentation and reconstruction of lung tissue, analysis of respiratory strain parameters, calculation of pulmonary vascular volume values, diagnosis and monitoring of patient lung diseases, etc.

[0112] The present invention is equipped with a complete image cross-correlation algorithm (multi-time spatial scale normalized cross-correlation analysis, time scale normalized cross-correlation analysis), which can complete the splicing of two 4D-CT scans of the lungs in spatial position and time axis based on the maximum correlation calculation of the axial plane; and, equipped with a complete image processing method for lung deformation curve analysis (image contrast enhancement, binarization, filtering and corrosion, connected domain marking, coronal lung area calculation), which can complete the splicing of two 4D-CT scans of the lungs in time axis based on the offset of the bilateral lung area or volume change curve of the coronal plane.

[0113] The present invention aims to compensate to some extent for the technical problem that a complete 4D-CT image of the lung cannot be obtained through a single scan due to the current CT equipment. Its potential applications are:

[0114] (1) In the radiology departments of major hospitals, chest CT scans are an important prerequisite for accurately determining lung deformation and diagnosing lung diseases, including tumors, pulmonary fibrosis, chronic obstructive pulmonary disease, lung tissue segmentation, and pulmonary vascular reconstruction;

[0115] (2) Based on the correlation calculation of the upper and lower lungs during splicing and the overlap of the deformation curves, feedback can be given on the scanning protocol (such as the scanning overlap area), optimizing the supporting product program or helping to develop a better scanning plan to obtain scanning data with closer respiratory phases for the patient.

[0116] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed by the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims

1. A dual 4D-CT scan spatiotemporal stitching method for obtaining a whole lung image, characterized in that: include: Step S101: obtaining a time series of volume packages of two 4D-CT scans of the upper lung region and the lower lung region; Step S102: performing spatial correlation analysis on the overlapping area of ​​the two scans to determine the best stitching layer; Step S103: analyzing two scans combined with timing data according to the determined splicing layer to determine the offset time of splicing; Step S104: stitching the images in the volume packages scanned twice according to the determined stitching layer and offset time; Wherein, before step S101, the process includes: Step S100: Formulate a 4D-CT segmented scanning plan: obtain a simulated positioning CT scan image, and formulate two 4D-CT scanning plans S1 and S2 in the head-foot direction, S1 is a scan covering the upper lung area, and S2 is a scan covering the lower lung area; In step S102, the correlation analysis adopts multi-time spatial scale normalized cross-correlation analysis, where: The multi-time spatial scale normalized cross-correlation analysis is as follows: a central image of a scanning layer of a certain time sequence number and an overlapping area is selected from the first scanned image as a stitching reference layer, all images of the scanning layer of the overlapping area of ​​the second scan within the entire scanning time range that can show the complete lungs are used as the layers to be analyzed, a normalized cross-correlation analysis is performed on the stitching reference layer and the layers to be analyzed, and the layer to be analyzed with the highest correlation ratio score is used as the best stitching layer corresponding to the stitching reference layer; In step S103, the analysis adopts time-scale normalized cross-correlation analysis, upper and lower lung area deformation curve analysis, sternum motion curve analysis, time-scale structural similarity analysis and / or upper and lower lung volume deformation curve analysis, wherein: The time scale normalized cross-correlation analysis is as follows: first, traverse the splicing reference layer in the entire scanning time range z 1 and its best stitching layer z 2. For each time sequence number Next z 1. According to the normalized cross-correlation value, z 2 Sort and find the most relevant z 2 The time sequence number , and calculate the corresponding z The time difference of 1 is defined as the delay ; then for all The interquartile range method is used to exclude outliers. Calculate the average and round to the nearest integer. That is, the offset time; The deformation curve analysis of the upper and lower lung areas is as follows: first, for the volume packages of the two scans, coronal images with a depth of 256 are obtained in the entire scanning time range; all coronal images are contrast enhanced and binarized; the images are further eroded and filtered to eliminate the white areas of small blood vessels and trachea in the lung area, fill the small gaps in the image, mark the connected domains of the image, and calculate the area sizes of the left and right lung areas at different times; the changes in the sum of the areas of both lungs are normalized to draw the deformation curve on the time series scale; finally, the horizontal offset of the upper and lower lung curves is calculated using the least squares difference algorithm based on the deformation curve to obtain That is, the offset time; The sternum motion curve analysis is as follows: first, for the volume packages of the two scans, a sagittal image with a width of 256 is obtained in the entire scanning time range; for all sagittal images, contrast enhancement and binarization are performed; the image is eroded and filtered to eliminate the white areas of small blood vessels and trachea in the lung area and fill the small gaps in the image; the pixel points of the outermost sternum with a depth of 310 in the sagittal image of the first scan and a depth of 11 in the sagittal image of the second scan are marked, and the changes in the coordinate position in the time series are stored; then a change curve of the horizontal coordinates of the outer side of the sternum in the time series is drawn; finally, the horizontal offset of the motion curve of the upper and lower sternum is calculated based on the change curve using the least square difference algorithm to obtain That is, the offset time; The time-scale structural similarity analysis is as follows: first, for the volume packages of the two scans, the images of the overlapping areas that can cover the complete lung section are selected as the overlapping volumes of the upper and lower lungs to be analyzed. Vs 1. Vs 2; Throughout the entire scanning time range, traverse all Vs 1 and Vs 2. For each time sequence number Next Vs 1. According to the structural similarity value, Vs 2 Sort and find the one with the largest similarity value Vs 2 The time sequence number , and calculate the corresponding Vs The time sequence number difference of 1 is defined as the offset ; then for all The interquartile range method is used to exclude outliers. Calculate the average and round to the nearest integer. That is, the offset time; The upper and lower lung volume deformation curve analysis is as follows: first, for the volume packages of the two scans, respectively, traverse the entire scanning time range, select the starting coronal surface layer, and the depth of the starting coronal surface layer , at a certain depth interval, a total of at least 15 coronal layers are analyzed to cover the complete lung volume; all coronal images are contrast enhanced and binarized, and the images are further eroded and filtered to eliminate the white areas of small blood vessels and trachea in the lung area, fill the small gaps in the image, mark the connected domains of the image, and calculate the area sizes of the left and right lung areas at different times; the sum of the areas of both lungs in at least 15 coronal layers at the same time is defined as the volume eigenvalue at that time, and the volume eigenvalues ​​of all time series are calculated; the changes in the volumes of both lungs are normalized to draw the deformation curve on the time series scale; finally, the horizontal offset of the two curves is calculated based on the deformation curve using the least squares difference algorithm to obtain That is the offset time.

2. The method according to claim 1, characterized in that The step S100 includes: The dynamic volume scanning mode is adopted, the maximum single scanning length of the 320-row detector is 16 cm, and the overlapping area of ​​the two scans is set to 1 cm, so as to ensure that the two 4D-CT images can cover the entire lung; And / or, after setting the 4D-CT scanning parameters, the designed dynamic volume scanning plan is transferred to the workstation, and the first 4D-CT scan is performed after the imaging device and related parameters are selected according to the plan. After the scan is completed, the patient's position is kept still, and then the operating table is manually translated along the head and feet direction to the position of the scanning plan S2, and the second 4D-CT scan is completed with the same parameters.

3. The method according to claim 1, characterized in that In step S101, each volume package contains 81 sets of files, each set of files contains 320 images; the 81 sets of files represent the time series of scanning, and the interval between each set of files is 0.1s; the 320 images are slices of different detector scanning positions, the physical interval between adjacent slices is 0.5mm, and each set is numbered in order from head to toe; And / or, the step S101 includes: All images were imported into the software and integrated into two 4D data sets, with the horizontal and vertical axes of the axial images defined as x and y, x,y [1,512]; Define the number of slices as z , z [1,320], the first scan image is from the apex of the lung, z 1=320 to the middle of the lung, i.e. z1=1. The second scan is from the middle of the lung, i.e. z 2=320 to the lung base, i.e. z2=1; the time series is defined as t , t [1,81].

4. The method according to claim 3, characterized in that The step S103 includes: Step S1031: According to the determined splicing layer, the two scans are combined with the time series data for analysis, that is, the time scale normalized cross-correlation analysis, the upper and lower lung area deformation curve analysis and the sternum motion curve analysis are used to obtain the offset time ; Step S1032: The offset time First offset time Determine if the absolute value of the difference between any result and another result , then draw a horizontal line under the two offset times, indicating that the two results are close; in this way, if the results of the three algorithms are crossed out, the average method is used to obtain the corrected offset time , ;like If it is not crossed out, it will be further checked using the structural similarity algorithm; if If no line is drawn, the upper and lower lung volume deformation curve analysis is used for further verification; If it is not crossed out, the result is discarded and the time is offset If the differences between the three are all greater than 5, it is considered that the multi-segment 4D-CT scan is not conducive to stitching, and manual verification should be performed to check whether the original scan image and the respiratory monitoring curve are synchronized.

5. The method according to claim 4, characterized in that The step S1032 includes: Step S1033: According to the determined splicing layer, the two scans are combined with the time series data for analysis, that is, the time scale structural similarity analysis and / or the upper and lower lung volume deformation curve analysis are used to obtain a newly obtained offset time result. or ; Step S1034: Respectively ,or Perform the second offset time determination.

6. The method according to claim 5, characterized in that The step S1034 includes: If the absolute value of the difference between any result and another result , then draw a horizontal line under the two offset times, indicating that the two results are close; in this way, if the results of the three algorithms are crossed out, the average method is used to obtain the corrected offset time , ,or ;like If it is not crossed out, the result is discarded and the time is offset ;like If it is not crossed out, the result is discarded and the time is offset .

7. The method according to claim 6, characterized in that The step S104 includes: The 4D-CT scan data of the two scans are integrated according to the corresponding time sequence numbers. The image data with the time sequence number M in S1 is selected, and the image data with the time sequence number M+ in S2 is selected. image data, and satisfy and For each group in S1, the stitching reference layer is selected as the stitching section, and for each group in S2, the best stitching layer is selected as the stitching section; the image data from the lung apex scanning layer to the stitching reference layer in the S1 time series numbered M is compared with the image data from the S2 time series numbered M+ The image data from the best stitching layer to the lung bottom scanning layer are integrated and stored in a new folder. The folders are numbered in chronological order, and the layers in each new folder are numbered in order from head to toe.

Images