Tracking method for region-of-interest of lung CT (Computed Tomography) image

Abstract

The invention provides a method for tracking a region of interest of a lung CT image, and the method comprises the steps: obtaining a lung CT image sequence of a patient at a plurality of time points, carrying out the preprocessing of the CT image sequence, and obtaining a preprocessed CT image sequence; for the preprocessed CT image sequence, extracting lung anatomical structure features in each CT image by adopting a feature extraction algorithm based on an anatomical structure, and constructing feature descriptors; performing shape analysis on the contour of the region of interest, extracting geometric features and topological features of the contour, and establishing a shape model of the region of interest; and based on the shape model, estimating a local deformation field between the regions of interest at different time points by adopting a non-parametric registration method, and calculating a displacement vector of each pixel point to obtain deformation field mapping.

Description

Technical Field

[0001] The present invention relates to the field of information technology, and in particular to a method for tracking a region of interest in a lung CT image. Background Art

[0002] In lung CT image tracking, the patient's respiratory movement and body position changes lead to differences in spatial position and morphology between CT images collected at different time points, which poses a challenge to accurately tracking the region of interest. Traditional registration methods based on grayscale values ​​or texture features are difficult to cope with such complex deformations and are prone to mismatching. Specifically, lung tissue undergoes non-rigid deformation during breathing, and the movement amplitude and direction of different lobes and segments vary. At the same time, the deformation pattern of lesion areas such as tumors is different from that of surrounding normal tissues. This complex local deformation brings difficulties to feature extraction and matching. How to weigh and fuse feature information at different scales, and how to effectively use prior knowledge to constrain and guide the feature matching process are issues that need to be explored in depth. While ensuring tracking accuracy, how to improve the robustness and adaptability of the algorithm so that it can cope with various complex clinical scenarios remains a challenging research topic. Summary of the invention

[0003] The present invention provides a method for tracking a region of interest in a lung CT image, which mainly includes:

[0004] Acquire a lung CT image sequence of a patient at multiple time points, and preprocess the CT image sequence to obtain a preprocessed CT image sequence;

[0005] For the preprocessed CT image sequence, using an anatomical structure-based feature extraction algorithm to extract lung anatomical structure features in each CT image and construct a feature descriptor;

[0006] According to the feature descriptors, an initial registration relationship is established between CT images at different time points, and based on the differences in anatomical structures between different individuals, a population statistics-based method is used to construct an a priori model of the anatomical structure to guide the image registration process;

[0007] In the CT image sequence after image registration, the region of interest is located and segmented, and the contour of the region of interest in each CT image is obtained by using the region growing algorithm or the level set method;

[0008] Performing shape analysis on the contour of the region of interest, extracting geometric features and topological features of the contour, and establishing a shape model of the region of interest;

[0009] Based on the shape model, a non-parametric registration method is used to estimate the local deformation field between the regions of interest at different time points, and a displacement vector of each pixel is calculated to obtain a deformation field map;

[0010] According to the deformation field mapping, deformation compensation is performed on the region of interest at different time points to achieve accurate correspondence between the CT images of the region of interest at different time points. Through deformation compensation, the anatomical structure and the region of interest at different time points are aligned to the same reference coordinate system to complete the tracking task.

[0011] The technical solution provided by the embodiment of the present invention may have the following beneficial effects:

[0012] The invention discloses a tracking method for a region of interest of a lung CT image, and the method is suitable for processing lung CT images acquired at different time points. The image quality is optimized by preprocessing, and the interference of noise and artifacts is eliminated, which lays a solid foundation for subsequent steps; a feature extraction algorithm based on anatomical structure is adopted to accurately capture the characteristics of lung structure, and a stable feature descriptor is constructed, so that good matching performance can be maintained even under complex deformation; in terms of positioning and segmentation of the region of interest (ROI), the application of the regional growing algorithm or the level set method enables the ROI contour to be accurately extracted even under conditions of complex background and low contrast, providing key information for further shape analysis and registration; the establishment of the shape model not only captures the geometric and topological characteristics of the contour, but also provides a morphological benchmark for non-parametric registration, so that the algorithm can more accurately estimate the local deformation field, calculate the pixel displacement, and obtain an accurate deformation field mapping; the deformation compensation method based on the deformation field mapping ensures the accurate correspondence between ROI CT images at different time points, and even in the face of significant changes in the respiratory cycle, the data consistency can be maintained, which greatly improves the accuracy and reliability of tracking. It not only solves the limitations of traditional methods in dealing with complex deformations, but also significantly enhances the robustness and adaptability of the algorithm, enabling it to cope with various clinical scenarios, especially showing great potential in disease dynamic monitoring and pathological development analysis. In general, the technical effects of the present invention include improving the accuracy and efficiency of lung CT image analysis, enabling doctors to more accurately diagnose and monitor the progression of lung diseases, especially in dealing with dynamic changes and pathological development analysis. In addition, through advanced image processing and data analysis technology, highly automated and accurate medical image analysis can be achieved, greatly improving the work efficiency of medical image analysis and the reliability of diagnosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 The present invention is a flow chart of a method for tracking a region of interest in a lung CT image.

[0014] Figure 2It is a schematic diagram of a method for tracking a region of interest in a lung CT image according to the present invention.

[0015] Figure 3 It is another schematic diagram of a method for tracking a region of interest in a lung CT image according to the present invention. DETAILED DESCRIPTION

[0016] The following will describe the technical solutions in the embodiments of the present invention in detail in conjunction with the accompanying drawings in the embodiments of the present invention. The described embodiments are only part of the embodiments of the present invention.

[0017] like Figure 1-3 In this embodiment, a method for tracking a region of interest in a lung CT image may specifically include:

[0018] S101 , acquiring a lung CT image sequence of a patient at multiple time points, and preprocessing the CT image sequence to obtain a preprocessed CT image sequence.

[0019] Acquire a lung CT image sequence of a patient at multiple time points, segment the lung parenchyma region of the CT image, set a threshold range according to the gray value distribution of the CT image, and obtain a segmented lung CT image sequence. Perform format conversion on the segmented CT image sequence, convert the CT image sequence into a DICOM format, and adjust the resolution of the CT image to a preset pixel value. Align the CT images at different time points in the CT image sequence, calculate the similarity between the CT images, and optimize the registration parameters using the gradient descent method to obtain an aligned CT image sequence. Perform noise suppression and contrast enhancement on the aligned CT image sequence, use a Gaussian filter to remove CT image noise and retain edge information, and adjust the gray distribution of the CT image by a histogram equalization method. Perform brightness and contrast standardization on the CT image sequence after adjusting the gray distribution of the CT image. If the brightness value of the CT image is lower than the preset threshold, the brightness of the CT image is enhanced; if the contrast value of the CT image is lower than the preset threshold, the contrast of the CT image is enhanced to obtain a preprocessed CT image sequence.

[0020] Specifically, the lung CT image sequences of patients at multiple time points were obtained, and the threshold segmentation algorithm was used to segment the lung parenchyma area of ​​the CT images at each time point. The threshold range was set according to the image gray value distribution to separate the lung parenchyma area from the surrounding tissue to obtain the segmented lung CT image sequence. The segmented CT image sequence was converted to DICOM format, and the image resolution was adjusted to 512x512 pixels to ensure that all images had consistent formats and resolutions. For the converted CT image sequence, the grayscale-based image registration algorithm was used to align the CT images at different time points, the maximum mutual information criterion method was used to calculate the similarity between the images, and the gradient descent method was used to optimize the registration parameters to obtain the registered CT image sequence. The registered CT image sequence was subjected to noise suppression and contrast enhancement. The Gaussian filter was used to remove image noise and retain edge information. The image gray distribution was adjusted by histogram equalization. The Laplace operator was used to enhance the boundary clarity of the lesion area. Finally, the brightness and contrast were standardized to make the images at different time points comparable, and the preprocessed CT image sequence was obtained. Lung CT image sequences of patients at multiple time points were obtained, and each sequence contained 30-50 cross-sectional images. The threshold segmentation algorithm was used to segment the lung parenchyma region of the CT images, and the grayscale value range was set to -950 to -500 HU as the lung parenchyma threshold. The region growing method was used to expand from the center of the lung until the threshold boundary was reached to obtain the preliminary segmentation result. Morphological operations were applied to fill the small blood vessels and bronchi in the lungs, and a 3x3 structure element was used for closing operations to remove small noise points. The segmented images were converted to a unified DICOM format, and the image resolution was adjusted to 512x512 pixels with a pixel spacing of 1 mm. Grayscale-based registration was performed on the converted image sequences, and the first time point image was selected as the reference image, and the mutual information was calculated as the similarity measure. The Powell optimization algorithm was used to iteratively adjust the affine transformation parameters, including translation, rotation, and scaling, until the mutual information value converged or the maximum number of iterations reached 100 times. A Gaussian filter was applied to the registered images for denoising, with the kernel size set to 5x5 and the standard deviation σ=1.5. The image grayscale range was stretched to 0-4095 by histogram equalization to increase the image contrast. The 3x3 Laplacian operator was used to sharpen the image edges and highlight the contour of the lesion area. Finally, the brightness and contrast were normalized, and the average grayscale value of the image was adjusted to 2048 and the standard deviation was adjusted to 500 to ensure the consistency of the grayscale distribution of the images at different time points, which was convenient for subsequent analysis and comparison.

[0021] S102: for the preprocessed CT image sequence, using an anatomical structure-based feature extraction algorithm to extract lung anatomical structure features in each CT image and construct a feature descriptor.

[0022] The regional growing algorithm is used to segment the lung lobes. The grayscale threshold range of the lung parenchyma area is determined according to the grayscale histogram of the CT image. The seed points of the left and right lungs are selected within the grayscale threshold range. The area is gradually expanded from the seed points to obtain the segmented left and right lung lobes mask images. Based on the segmented lung images, the multi-scale Hessian matrix analysis method is used to enhance the trachea and blood vessels. By comparing the size relationship of the Hessian matrix eigenvalues, it is determined whether the pixel point belongs to a tubular structure. If the eigenvalue meets the preset requirements, it is determined to be a tubular structure and enhanced to obtain the enhanced trachea and blood vessel images. The three-dimensional connected domain analysis method is applied to the enhanced trachea and blood vessel images to extract the three-dimensional structure of the trachea and blood vessels. The volume threshold is set to remove the small volume connected domain. The trachea and blood vessels are distinguished by morphological features to obtain the three-dimensional skeleton of the tracheal tree and blood vessel tree. Based on the three-dimensional skeleton of the tracheal tree and the vascular tree, the geometric and statistical features of each structure are calculated, including the number of branches, branch angles, and diameter change rate of the trachea, the volume fraction, vascular tortuosity, and vascular distribution density of the blood vessels, and the volume, surface area, average density, density standard deviation, shape factor, and curvature of the lung lobes. The calculated features are normalized and combined into feature vectors to construct feature descriptors for each CT image.

[0023] Specifically, for the preprocessed CT image sequence, the region growing algorithm is used to segment the lung lobes. The grayscale threshold range of the lung parenchyma area is determined by calculating the image grayscale histogram. The seed points of the left and right lungs are automatically selected within this range. The region is gradually expanded from the seed points to obtain the segmented left and right lung lobes mask images. The segmented lung images are enhanced for trachea and blood vessels using multi-scale Hessian matrix analysis. The Hessian matrix of the image at different scales is calculated. By comparing the size relationship of the matrix eigenvalues ​​λ1, λ2, and λ3, it is determined whether the pixel belongs to a tubular structure. If |λ1|≤|λ2|≤|λ3|

[0024] If |λ1| approaches 0, it is determined to be a tubular structure and enhanced to obtain enhanced tracheal and vascular images, where the eigenvalues ​​λ1, λ2, and λ3 are the three eigenvalues ​​of the Hessian matrix at a certain point in the 3D image. The enhanced image is subjected to three-dimensional connected domain analysis to extract the three-dimensional structure of the trachea and blood vessels. The volume threshold is set to remove the small volume connected domain. The trachea and blood vessels are distinguished by morphological features. The trachea is a tree-like branching structure, and the blood vessels are distributed in a network-like manner. The three-dimensional skeletons of the tracheal tree and the blood vessel tree are obtained. Based on the extracted anatomical structure, the geometric and statistical features of each structure are calculated. The number of branches, branch angles, and diameter change rate of the trachea are calculated. The volume fraction, vascular tortuosity, and vascular distribution density of the blood vessels are calculated. The volume, surface area, average density, density standard deviation, shape factor, and curvature of the lung lobe are calculated. At the same time, morphological operations and density threshold methods are used to identify lung nodules and interstitial lungs, extract their size, shape, and edge features, and these features are normalized and combined into feature vectors to construct feature descriptors for each CT image. For the preprocessed CT image sequence, the image grayscale histogram was first calculated to determine the grayscale threshold range of the lung parenchyma area from -950HU to -500HU. Within this range, the seed points of the left and right lungs were automatically selected, usually located in the central area of ​​the lung field. The region growing algorithm was used for lung lobe segmentation, and the growth threshold was set to 30HU. The region was gradually expanded until the lung wall boundary was reached. The morphological closing operation was applied to the segmentation results, and the 3x3 structural element was used to fill the small blood vessels and bronchial cavities. Then, a multi-scale Hessian matrix analysis was performed, and three scales σ = 1, 2, and 3mm were selected to calculate the Hessian matrix eigenvalues ​​λ1, λ2, and λ3 of each pixel at each scale. If |λ1|≤|λ2|≤|λ3| and |λ1|<0.1×|λ3|, it was determined to be a tubular structure. For the tubular structure pixels, the response function R = |λ2| / |λ3| was calculated. The smaller the R value, the more obvious the structure. Based on this, the trachea and blood vessels were enhanced. The enhanced images were analyzed by three-dimensional connected domain analysis, and the minimum volume threshold was set to 100 voxels to remove small volume noise points. The trachea and blood vessels were distinguished by calculating the branching pattern and direction of the structure. The trachea showed regular tree-like branches with a trunk diameter of 8-15 mm and a branching angle of 45°-90°; the blood vessels showed irregular network shape with a diameter of 1-5 mm and a more curved direction. The features of each anatomical structure were extracted: the number of tracheal branches, the average branching angle, and the change rate of the tube diameter were calculated; the vascular volume fraction, the average tortuosity, and the length of the blood vessels per unit volume were calculated; the lung lobe volume, surface area, average CT value, standard deviation of CT value, and shape factor (6×π^(1 / 2)×volume / surface area^(3 / 2)) were calculated. At the same time, morphological operations and a density threshold of -600HU were used to identify lung nodules and interstitial tissues, and their maximum diameter, volume, sphericity, and edge sharpness were extracted. Finally, all features were normalized by min-max, and a 128-dimensional feature vector was constructed as the feature descriptor of each CT image.

[0025] S103, establishing an initial registration relationship between CT images at different time points according to the feature descriptors; and constructing a priori model of the anatomical structure to guide the image registration process by using a population statistics-based method based on the differences in anatomical structures between different individuals.

[0026] CT images at two different time points are acquired, and the correspondence between feature points is established using the nearest neighbor search algorithm. The initial matching point pair is obtained by calculating the Euclidean distance of the feature vector. According to the set distance threshold and number of iterations, the random sampling consistency method is used to eliminate abnormal matching points and screen out a stable set of matching points. For the anatomical structure in the group CT image, the covariance matrix is ​​calculated and the eigenvalues ​​and eigenvectors are solved. The first k principal components are selected to construct a statistical shape model and a density distribution model. Based on the initial registration relationship and the prior model, the B-spline free deformation algorithm is used for non-rigid registration; the registration results are quantitatively evaluated, the normalized mutual information and correlation coefficient before and after the registration are calculated, and the average distance error of the corresponding points is calculated by combining the landmark points automatically selected based on the image gradient and anatomical knowledge. If the average distance error is greater than the preset threshold, the non-rigid registration is performed again, the control point density and weight coefficient are adjusted, and the parameters are optimized and re-registered.

[0027] Specifically, based on the feature descriptors, the nearest neighbor search algorithm is used to establish the correspondence between feature points between CT images at different time points. By calculating the Euclidean distance of the feature vector, the feature point pair with the smallest distance is selected as the initial matching point pair to obtain the initial registration relationship. The random sampling consistency algorithm is used to eliminate abnormal matching points, set the distance threshold and the number of iterations, and screen out a stable set of matching points. The principal component analysis of the anatomical structure in the group CT images is performed, the covariance matrix is ​​calculated, the eigenvalues ​​and eigenvectors are solved, the first k principal components are selected, and the statistical shape model and density distribution model are constructed. The shape model is expressed as the average shape plus the linear combination of the principal components, and the density model is expressed as the average density distribution plus the linear combination of the principal components. Based on the initial registration relationship and the prior model, the B-spline free deformation algorithm is used for non-rigid registration, and the control point grid is defined. By minimizing the combined objective function

[0028] E=αEsim(I1,T(I2))+βEshape(T) optimizes the deformation field parameters, where Esim is the image similarity term, which is measured by mutual information, Eshape is the shape constraint term based on the prior model, I1 is the reference image, I2 is the template image, T is the deformation field, T(I2) is the deformed image, α and β are weight coefficients, which respectively determine the importance of the image similarity term and the shape constraint term based on the prior model in the objective function, and they are used to balance the contribution between image similarity and shape constraint. The objective function is iteratively optimized using the gradient descent method to obtain an accurate spatial correspondence. The registration results are quantitatively evaluated, and the normalized mutual information and correlation coefficients before and after the registration are calculated. The average distance error of the corresponding points is calculated by combining the landmark points automatically selected based on the image gradient and anatomical knowledge. If the accuracy does not meet the preset threshold, return to the non-rigid registration step, adjust the control point density and weight coefficient, and perform parameter optimization and re-registration. During the implementation process, we first extracted 128-dimensional feature descriptors for each CT image, and used the kd tree algorithm to perform nearest neighbor search. The search radius was set to 0.1 to find the nearest neighbor of each feature point. Then, we used the random sampling consistency algorithm, set the distance threshold to 2 mm, and the number of iterations to 1000. We randomly selected 3 pairs of points from the initial matching points, calculated the affine transformation matrix, and counted the number of point pairs that met the transformation. We repeated this process to obtain the best matching point set. We performed principal component analysis on the CT images of 100 patients, calculated the covariance matrix, solved the eigenvalues ​​and eigenvectors, and selected the top 10 principal components to explain 95% of the shape variation. We constructed a statistical shape model. where x is a shape vector, is the average shape, P is the principal component matrix, and b is the shape parameter vector. The density distribution model can also be constructed similarly. In non-rigid registration, the cubic B-spline free deformation algorithm is used, a 10x10x10 control point grid is set, and the combined objective function E=0.7NMI(I1, T(I2))+0.3||b||2 is defined, where NMI is the normalized mutual information, ||b||2 is the L2 norm of the shape parameter, and the normalized mutual information measures the degree of information sharing between the two images and normalizes it to a fixed interval, usually [0, 1]. 0.7NMI(I1, T(I2)) measures the similarity between the reference image I1 and the template image T(I2) after the deformation field T transformation. The L-BFGS algorithm is used to optimize the objective function, the maximum number of iterations is set to 100, and the convergence threshold is 1e-6. In the evaluation stage, 50 anatomical landmarks were automatically selected, including the tracheal bifurcation point and the junction of the pulmonary lobe fissure, and the root mean square error of the corresponding distance of the landmarks was calculated, and the threshold was set to 2 mm. If the threshold was exceeded, the control point grid was encrypted to 15x15x15, the weights were adjusted to 0.8 and 0.2, and the registration optimization was re-performed.

[0029] S104 , in the CT image sequence after image registration, locate and segment the region of interest, and use a region growing algorithm or a level set method to obtain the contour of the region of interest in each CT image.

[0030] Obtain a registered CT image sequence, smooth the image using a Gaussian filter according to the CT image sequence, enhance the image contrast through adaptive histogram equalization, and obtain an enhanced CT image sequence. For the enhanced CT image sequence, determine the initial position of the region of interest using the Canny edge detection algorithm based on the image gray value distribution characteristics. According to the initial position, accurately segment the region of interest using a region growing algorithm, and the region growing algorithm includes adding spatial constraints and adaptive thresholds. For the segmented region of interest, optimize the region boundary using morphological operations, extract the precise contour of the region of interest in combination with the Sobel edge detection algorithm, and obtain the final contour of the region of interest in each CT image. If the region growing algorithm cannot obtain the precise contour of the region of interest, the level set method is used as an alternative, and the level set method uses a signed distance function to initialize the level set.

[0031] Specifically, the registered CT image sequence was preprocessed, and the image was smoothed using a Gaussian filter. The kernel size was set to 5x5, and the standard deviation σ=1.5. The image contrast was enhanced by adaptive histogram equalization. The block size was set to 8x8, and the contrast limit parameter was set to 0.02 to obtain the enhanced CT image sequence. In the enhanced CT image sequence, the initial position of the region of interest was automatically located using the image gray value distribution characteristics and the Canny edge detection algorithm. The low threshold of the Canny algorithm was set to 50 and the high threshold was set to 150. The image gradient amplitude and direction were calculated, and the rough outline of the region of interest was preliminarily determined by combining the region growing algorithm. Based on the positioning results, the improved region growing algorithm was used to accurately segment the region of interest. The improvements included adding spatial constraints and adaptive thresholds, selecting multiple seed points from the initial outline, and selecting seed points based on image gradient and gray value. The pixel points with the smallest gradient and gray value within the range of ±10% of the average gray value of the target area were selected. The initial gray value similarity threshold was set to 20, and the spatial connectivity constraint was 26 neighborhoods. The region growing criterion was iteratively updated to obtain the segmented region of interest. For the segmentation results, morphological operations were applied to optimize the region boundaries. 3x3 circular structure elements were used for opening operations to remove small protrusions, and 5x5 circular structure elements were used for closing operations to fill the boundary gaps. The Sobel edge detection algorithm was combined to extract the precise contour of the region of interest, and the final contour of the region of interest in each CT image was obtained. If the region growing algorithm could not obtain satisfactory results, the level set method was used as an alternative. The level set was initialized using the signed distance function. The evolution rate function combined the image gradient and regional statistical information. The number of iterations was set to 200 to obtain the contour of the region of interest. In practical applications, when processing the registered CT image sequence, a 5x5 Gaussian filter was first used for smoothing, and the σ value was set to 1.5, which effectively removed noise while retaining edge information. Then, adaptive histogram equalization was applied to divide the image into 8x8 small blocks, and the contrast limit parameter was set to 0.02, which significantly improved the image contrast. On the enhanced image, the Canny edge detection algorithm was used with a low threshold of 50 and a high threshold of 150 to accurately detect the edge of the region of interest. Combined with the detected edge information, the improved region growing algorithm was used for segmentation. The algorithm selects 10 seed points from the inside of the initial contour, whose gradient values ​​are less than 5 and whose grayscale values ​​are within ±10% of the average grayscale of the target area. The initial similarity threshold is set to 20, and 26 neighborhoods are used as spatial constraints. The preliminary segmentation results are obtained after 50 iterations. Morphological operations are applied to the segmentation results. First, a 3x3 circular structure element is used for opening operations to remove protrusions with an area of ​​less than 10 pixels, and then a 5x5 circular structure element is used for closing operations to fill holes with an area of ​​less than 20 pixels. Finally, the Sobel operator is used to extract the precise contour, and the gradient threshold is set to 100.If the Dice coefficient of the region growing algorithm is lower than 0.85, the method is switched to the level set method. The level set is initialized using the distance transform, the weight of the curvature term in the evolution rate function is set to 0.1, the weight of the region term is 0.6, and the weight of the edge term is 0.3. It is iterated 200 times or until the contour changes less than 0.1%.

[0032] S105 , performing shape analysis on the contour of the region of interest, extracting geometric features and topological features of the contour, and establishing a shape model of the region of interest.

[0033] The geometric features of the contour of the region of interest are obtained, and the geometric features include area, perimeter, major axis length, minor axis length, circularity and eccentricity; the Douglas-Peucker algorithm is used to perform polygon fitting on the contour to obtain the number of vertices and the sum of internal angles of the contour; according to the coordinates of the contour points, a third-order B-spline curve is used for fitting to obtain the curvature distribution of the contour; if the curvature value is greater than a preset value of the average curvature, it is determined to be an inflection point, and if the curvature value is greater than a preset multiple of the average curvature, it is determined to be a cusp; the skeleton of the region of interest is obtained by a medial axis transformation algorithm, and the number of branches and length distribution of the skeleton are calculated; the connectivity of the region of interest is determined according to the Euler number, which is equal to the number of regions minus the number of holes; for the extracted geometric features and topological features, the principal component analysis method is used to reduce the feature dimension; according to the eigenvalues ​​and eigenvectors of the feature covariance matrix, the principal component with a cumulative contribution rate greater than a preset threshold is selected as a shape descriptor; the weight coefficients of the shape descriptor are formed into a vector to obtain a compact shape representation of the region of interest.

[0034] Specifically, the geometric features of the contour of the region of interest are extracted, the area and perimeter are calculated by the coordinates of the contour points, the length of the major axis and the minor axis are obtained by the minimum circumscribed rectangle algorithm, the circularity and eccentricity are calculated, the Douglas-Peucker algorithm is used for polygon fitting, the threshold is set to 0.5% of the contour perimeter, and the number of vertices and the sum of the internal angles of the contour are obtained. The curvature analysis of the contour is performed, and the contour point sequence is fitted with a third-order B-spline curve. The number of control points is set to 1 / 10 of the number of contour points. The curvature value of each point on the curve is calculated, and the mean, variance and maximum value of the curvature are statistically calculated. The inflection point of the contour is extracted by setting the curvature threshold to 2 times the average curvature. Points with curvature greater than 3 times the average curvature are defined as cusps, and the local shape features of the contour are obtained. The topological features of the region of interest are extracted, the skeleton of the region is obtained by the medial axis transformation algorithm, the number of skeleton branches and the length distribution are calculated, the connectivity of the region is calculated by the Euler number, which is equal to the number of regions minus the number of holes, and the thickness distribution of the region is obtained by distance transformation. Based on the extracted geometric and topological features, the shape model of the region of interest was constructed. The principal component analysis method was used to reduce the feature dimension. The eigenvalues ​​and eigenvectors of the feature covariance matrix were calculated. The first K principal components were selected as shape descriptors based on the principle that the cumulative contribution rate was greater than 95%. The weight coefficients of these K principal components were combined into a vector as a compact shape representation of the region of interest. In practical applications, when performing shape analysis on the contour of the region of interest of a pulmonary nodule, the geometric features were first extracted. The area was calculated using the Green formula to obtain 3.14 cm2; the perimeter was calculated using the Euclidean distance between the contour points to obtain 6.28 cm. The minimum circumscribed rectangle algorithm obtained a major axis of 2.5 cm and a minor axis of 1.6 cm. Based on this, the circularity was calculated to be 0.79 and the eccentricity was 0.77. The Douglas-Peucker algorithm was used for polygon fitting, and the threshold was set to 0.0314 cm, which was 0.5% of the perimeter. 12 vertices were obtained, and the sum of the internal angles was 1800°. Then, the curvature analysis was performed, and the contour was fitted using a 3rd order B-spline curve. The number of control points was 63, which is 1 / 10 of the contour points of 630. The average curvature was calculated to be 0.5 cm -1 , curvature variance 0.04cm -2 , maximum curvature 1.2cm -1 . Set the curvature threshold to 1.0cm -1 , which is twice the average curvature, 5 inflection points are extracted; the curvature is greater than 1.5cm -1The points are defined as cusps, and 2 are detected. In the topological feature extraction, the medial axis transformation algorithm obtains a skeleton, which contains 3 main branches, and the longest branch length is 1.8cm. The Euler number is calculated to be 1, including 1 area and 0 holes. The distance transformation obtains an average thickness of 0.8cm and a maximum thickness of 1.2cm. Finally, the shape model is constructed, the principal component analysis is used for dimensionality reduction, and the eigenvalues ​​and eigenvectors of the feature covariance matrix are calculated. The first 6 principal components corresponding to the cumulative contribution rate reaching 95% are selected as shape descriptors. The weight coefficients of these 6 principal components form a vector [-0.35, 0.28, -0.15, 0.10, -0.08, 0.05] as a compact shape representation of the pulmonary nodule.

[0035] S106 , based on the shape model, using a non-parametric registration method to estimate the local deformation field between the regions of interest at different time points, calculating the displacement vector of each pixel point, and obtaining a deformation field map.

[0036] According to the feature points extracted from the shape model, the correspondence between the regions of interest at different time points is established, and the similarity measure is obtained by calculating the Euclidean distance of the feature descriptor. If the similarity measure is less than the preset threshold, it is determined to be a stable matching point pair as the control point for deformation field estimation. The initial deformation field is constructed, and the displacement of each pixel is interpolated and estimated by solving the coefficient matrix of the thin plate spline function to obtain a rough deformation field mapping. The deformation field is optimized using a three-layer multi-resolution strategy, step by step from low resolution to original resolution. At each resolution level, the displacement vector of each pixel is iteratively updated by minimizing the combined objective function of the image mutual information and the sum of squares of the second-order derivatives of the deformation field. The optimized deformation field is post-processed, and the vector median filter is used to remove local abnormal displacements. The topological consistency of the deformation field is maintained by the elastic grid method, and the grid node position is iteratively optimized until convergence or the preset maximum number of iterations is reached. According to the post-processing results, a local deformation field mapping is obtained. The local deformation field mapping is used to characterize the deformation relationship of the region of interest between different time points.

[0037] Specifically, based on the feature points extracted from the shape model, the scale-invariant feature transformation algorithm is used to establish the correspondence between the regions of interest at different time points. The Euclidean distance of the feature descriptor is calculated as the similarity measure, and the threshold is set to 0.7 to screen out stable matching point pairs as control points for deformation field estimation. The thin plate spline interpolation method is used to construct the initial deformation field according to the matched control point pairs, and the coefficient matrix of the thin plate spline function is solved. The displacement of each pixel is interpolated and estimated to obtain a rough deformation field mapping. A three-layer multi-resolution strategy is adopted to gradually optimize the deformation field from 1 / 4, 1 / 2 to the original resolution. At each resolution level, the weight coefficients are set to 0.7 and 0.3 respectively by minimizing the combined objective function of the image mutual information and the sum of squares of the second-order derivatives of the deformation field. The displacement vector of each pixel is iteratively updated using the gradient descent method, and the number of iterations is set to 100. The optimized deformation field was post-processed, and the vector median filter with a 5x5 window was applied to remove local abnormal displacement. The topological consistency of the deformation field was maintained by the elastic grid method. The grid unit size was set to 10x10 pixels, the elastic coefficient was set to 0.1, and the grid node position was iteratively optimized until convergence or the maximum number of iterations reached 50, and finally a fine local deformation field mapping was obtained. In practical applications, when non-parametric registration of CT images of pulmonary nodules at different time points was performed, the scale-invariant feature transformation algorithm was first used to extract feature points, usually about 500 feature points were extracted for each image. The 128-dimensional descriptors of these feature points were calculated, and the Euclidean distance was used as the similarity metric. The matching threshold was set to 0.7, and about 200 pairs of stable matching points were screened out. Subsequently, the initial deformation field was constructed using the thin plate spline interpolation method. The coefficient matrix of the thin plate spline function was solved, and the matrix size was 200x200, corresponding to 200 control points. Each pixel in the 512x512 resolution image was interpolated to obtain the initial deformation field. A three-layer multi-resolution strategy was adopted to optimize the deformation field at 128x128, 256x256, and 512x512 resolutions, respectively. At each resolution level, the weights of the mutual information and deformation field smoothness constraints were set to 0.7 and 0.3, and the gradient descent method was used for optimization, with a learning rate of 0.01 and 100 iterations. For example, at a resolution of 256x256, the displacement updates of 65,536 pixels were calculated at each iteration. After optimization, a vector median filter with a 5x5 window was used to remove abnormal displacements, and the filter window covered 25 adjacent pixels. Finally, the elastic grid method was applied to divide the 512x512 image into 51x51 grid cells, each with a size of 10x10 pixels. The elastic coefficient was set to 0.1, and the grid node positions were iteratively optimized, usually converging within 30-40 iterations. The resulting deformation field accurately describes the local deformation of lung nodules at different time points, with an average displacement error of less than 0.5 pixels.

[0038] S107. According to the deformation field mapping, deformation compensation is performed on the region of interest at different time points to achieve accurate correspondence between the region of interest in the CT images at different time points; through deformation compensation, the anatomical structure and the region of interest at different time points are aligned to the same reference coordinate system to complete the tracking task.

[0039] According to the deformation field mapping, the region of interest at different time points is obtained, and deformation compensation is performed on the region of interest. The CT image at the middle time point with the highest signal-to-noise ratio is determined as the reference coordinate system from the CT images at different time points. The affine transformation matrix is ​​calculated using the least squares method, and the region of interest is globally aligned through the affine transformation matrix; the region of interest is finely aligned using the local deformation field. The mutual information value and correlation coefficient before and after the registration are calculated for the deformation compensation result. If the mutual information value is less than the preset mutual information value threshold or the correlation coefficient is less than the preset correlation coefficient threshold, the gradient descent method is used to locally optimize the deformation compensation result. Shape features, density features and texture features are extracted from the aligned region of interest, and feature correspondences between different time points are established based on the shape features, density features and texture features.

[0040] Specifically, according to the deformation field mapping, the bilinear interpolation algorithm is used to compensate the deformation of the region of interest at different time points. By calculating the four nearest neighbor pixel values ​​of each target pixel in the original image, the weighted average is performed according to the distance weight, and the pixel value in the original image is mapped to the new position to obtain the deformation compensated image. The CT image at the middle time point with the highest signal-to-noise ratio is selected as the reference coordinate system, and the region of interest at other time points is aligned to the reference coordinate system through deformation compensation. The affine transformation matrix is ​​calculated by the least squares method for global alignment, and then the local deformation field is applied for fine alignment. The deformation compensation results are quantitatively evaluated, and the mutual information value and correlation coefficient before and after alignment are calculated. The mutual information value threshold is set to 1.5 and the correlation coefficient threshold is set to 0.95 to judge the alignment quality. If it is lower than the threshold, the gradient descent method is used for local optimization. The learning rate is set to 0.01, and the iterative optimization is performed until the preset threshold is reached or the maximum number of iterations is 100. Based on the aligned ROI, the shape features of the region are extracted, including volume, surface area, maximum diameter, density features, including average CT value, standard deviation and texture features, including energy, contrast and correlation of gray-level co-occurrence matrix, and the feature correspondence between different time points is established. The feature changes at the intermediate time points are estimated by the cubic spline interpolation algorithm to achieve continuous tracking of the ROI. In practical applications, when the CT images of pulmonary nodules at different time points are deformed and tracked, the bilinear interpolation algorithm is first used to resample the pixels according to the deformation field mapping. For example, for an image with a resolution of 512x512, the new value of each pixel is obtained by weighted averaging the adjacent 2x2 pixel blocks in the original image, and the weight is calculated based on the inverse distance. The image at the third time point with the highest signal-to-noise ratio among the five time points is selected as the reference coordinate system, and the images at other time points are globally aligned by calculating the 3x3 affine transformation matrix by the least squares method, and the error threshold is set to 0.5 pixels. The local deformation field is then applied to the sub-pixel level. The deformation compensation results are evaluated, the mutual information value and correlation coefficient are calculated, and the mutual information value threshold is set to 1.5 and the correlation coefficient threshold is set to 0.95. If it is lower than the threshold, the gradient descent method is used for local optimization, the learning rate is 0.01, and the maximum iteration is 100 times. For the aligned region of interest, shape features such as volume, surface area, maximum diameter, density features, including average CT value, standard deviation and texture features are extracted. The gray level co-occurrence matrix is ​​used with a window size of 5x5 to calculate energy, contrast and correlation. For example, for a lung nodule with a diameter of 2 cm, the volume is approximately 4.19 cm 2 The average CT value was -600HU, and the standard deviation was 150HU. The cubic spline interpolation algorithm was used to estimate the feature changes at the four intermediate time points between the five time points. The number of interpolation nodes was set to 20, achieving continuous tracking of nine time points and fully describing the dynamic changes of pulmonary nodules.

[0041] The above description is merely a preferred embodiment of one or more embodiments of the present specification and is not intended to limit one or more embodiments of the present specification. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of one or more embodiments of the present specification shall be included in the scope of protection of one or more embodiments of the present specification.

Claims

1. A method for tracking a region of interest in a lung CT image, characterized in that: The method comprises: Acquire a lung CT image sequence of a patient at multiple time points, and preprocess the CT image sequence to obtain a preprocessed CT image sequence; For the preprocessed CT image sequence, using an anatomical structure-based feature extraction algorithm to extract lung anatomical structure features in each CT image and construct a feature descriptor; According to the feature descriptors, an initial registration relationship is established between CT images at different time points, and based on the differences in anatomical structures between different individuals, a population statistics-based method is used to construct an a priori model of the anatomical structure to guide the image registration process; In the CT image sequence after image registration, the region of interest is located and segmented, and the contour of the region of interest in each CT image is obtained by using the region growing algorithm or the level set method; Performing shape analysis on the contour of the region of interest, extracting geometric features and topological features of the contour, and establishing a shape model of the region of interest; Based on the shape model, a non-parametric registration method is used to estimate the local deformation field between the regions of interest at different time points, and a displacement vector of each pixel is calculated to obtain a deformation field map; According to the deformation field mapping, deformation compensation is performed on the region of interest at different time points to achieve accurate correspondence between the CT images of the region of interest at different time points. Through deformation compensation, the anatomical structure and the region of interest at different time points are aligned to the same reference coordinate system to complete the tracking task.

2. The method according to claim 1, wherein: The step of acquiring a lung CT image sequence of a patient at multiple time points and preprocessing the CT image sequence to obtain a preprocessed CT image sequence includes: Acquire a lung CT image sequence of the patient at multiple time points, perform lung parenchymal region segmentation on the CT image, set a threshold range according to the gray value distribution of the CT image, and obtain a segmented lung CT image sequence; Performing format conversion on the segmented CT image sequence, converting the CT image sequence into a DICOM format, and adjusting the CT image resolution to a preset pixel value; Aligning the CT images at different time points in the CT image sequence, calculating the similarity between the CT images, optimizing the registration parameters using a gradient descent method, and obtaining an aligned CT image sequence; Performing noise suppression and contrast enhancement on the aligned CT image sequence, removing CT image noise using a Gaussian filter and retaining edge information, and adjusting the grayscale distribution of the CT image using a histogram equalization method; Performing brightness and contrast standardization processing on the CT image sequence after adjusting the grayscale distribution of the CT image, and if the brightness value of the CT image is lower than a preset threshold, performing brightness enhancement on the CT image; If the contrast value of the CT image is lower than a preset threshold, the CT image is contrast enhanced to obtain a preprocessed CT image sequence.

3. The method according to claim 1, wherein: The method of extracting lung anatomical structure features in each CT image using an anatomical structure-based feature extraction algorithm for the preprocessed CT image sequence and constructing a feature descriptor includes: The region growing algorithm is used to segment the lung lobes. The grayscale threshold range of the lung parenchyma area is determined according to the grayscale histogram of the CT image. The seed points of the left and right lungs are selected within the grayscale threshold range. The region is gradually expanded from the seed points to obtain the segmented left and right lung lobes mask images. Based on the segmented lung image, the multi-scale Hessian matrix analysis method is used to enhance the trachea and blood vessels. By comparing the size relationship of the Hessian matrix eigenvalues, it is determined whether the pixel point belongs to a tubular structure. If the eigenvalue meets the preset requirements, it is determined to be a tubular structure and enhanced to obtain the enhanced trachea and blood vessel images. Applying a three-dimensional connected domain analysis method to the enhanced trachea and blood vessel images to extract the three-dimensional structures of the trachea and blood vessels, setting a volume threshold to remove small volume connected domains, distinguishing the trachea and blood vessels through morphological features, and obtaining the three-dimensional skeletons of the tracheal tree and the blood vessel tree; Based on the three-dimensional skeletons of the tracheal and vascular trees, the geometric and statistical features of each structure are calculated, including the number of branches, branch angles, and diameter change rate of the trachea, the volume fraction, vascular tortuosity, and vascular distribution density of the blood vessels, and the volume, surface area, average density, density standard deviation, shape factor, and curvature of the lung lobes. The calculated features are normalized and combined into feature vectors to construct feature descriptors for each CT image.

4. The method according to claim 1, wherein: The method of establishing an initial registration relationship between CT images at different time points according to the feature descriptors, and constructing a priori model of the anatomical structure to guide the image registration process by using a population statistics-based method based on the differences in anatomical structures between different individuals, includes: Obtain CT images at two different time points, use the nearest neighbor search algorithm to establish the correspondence between feature points, and obtain the initial matching point pairs by calculating the Euclidean distance of the feature vectors; According to the set distance threshold and number of iterations, the random sampling consistency method is used to remove abnormal matching points and screen out a stable matching point set; For the anatomical structures in the group CT images, the covariance matrix is ​​calculated and the eigenvalues ​​and eigenvectors are solved, the first k principal components are selected, and the statistical shape model and density distribution model are constructed; Based on the initial registration relationship and prior model, the B-spline free deformation algorithm is used for non-rigid registration; Quantitatively evaluate the registration results, calculate the normalized mutual information and correlation coefficient before and after the registration, and calculate the average distance error of the corresponding points by combining the landmark points automatically selected based on the image gradient and anatomical knowledge; If the average distance error is greater than a preset threshold, non-rigid registration is performed again, the control point density and weight coefficient are adjusted, and parameter optimization and re-registration are performed.

5. The method according to claim 1, wherein: The method of locating and segmenting a region of interest in the CT image sequence after image registration, and obtaining the contour of the region of interest in each CT image by using a region growing algorithm or a level set method, comprises: Acquire the registered CT image sequence, use a Gaussian filter to smooth the image according to the CT image sequence, enhance the image contrast by adaptive histogram equalization, and obtain an enhanced CT image sequence; For the enhanced CT image sequence, according to the image gray value distribution characteristics, the Canny edge detection algorithm is used to determine the initial position of the region of interest; According to the initial position, a region growing algorithm is used to accurately segment the region of interest, wherein the region growing algorithm includes adding spatial constraints and adaptive thresholds; For the segmented region of interest, morphological operations are applied to optimize the region boundary, and the precise contour of the region of interest is extracted in combination with the Sobel edge detection algorithm to obtain the final contour of the region of interest in each CT image; If the region growing algorithm cannot obtain the precise outline of the region of interest, the level set method is used as an alternative, and the level set method uses a signed distance function to initialize the level set.

6. The method according to claim 1, wherein: The step of performing shape analysis on the contour of the region of interest, extracting geometric features and topological features of the contour, and establishing a shape model of the region of interest includes: Acquire geometric features of the contour of the region of interest, wherein the geometric features include area, perimeter, major axis length, minor axis length, circularity, and eccentricity; The Douglas-Peucker algorithm is used to perform polygon fitting on the contour to obtain the number of vertices and the sum of the internal angles of the contour. According to the coordinates of the contour points, a third-order B-spline curve is used for fitting to obtain the curvature distribution of the contour; If the curvature value is greater than the preset value of the average curvature, it is determined to be an inflection point; if the curvature value is greater than the preset multiple of the average curvature, it is determined to be a cusp; The skeleton of the region of interest is obtained through the medial axis transformation algorithm, and the number of branches and length distribution of the skeleton are calculated; The connectivity of the region of interest is determined based on the Euler number, which is equal to the number of regions minus the number of holes; For the extracted geometric and topological features, the principal component analysis method is used to reduce the feature dimension; According to the eigenvalues ​​and eigenvectors of the feature covariance matrix, the principal component whose cumulative contribution rate is greater than a preset threshold is selected as the shape descriptor; The weight coefficients of the shape descriptors are combined into a vector to obtain a compact shape representation of the region of interest.

7. The method according to claim 1, wherein: Based on the shape model, a non-parametric registration method is used to estimate the local deformation field between the regions of interest at different time points, and the displacement vector of each pixel is calculated to obtain the deformation field mapping, including: Based on the feature points extracted from the shape model, the correspondence between the regions of interest at different time points is established, and the similarity measure is obtained by calculating the Euclidean distance of the feature descriptors; If the similarity measure is less than a preset threshold, it is determined to be a stable matching point pair and used as a control point for deformation field estimation; Construct the initial deformation field, and interpolate and estimate the displacement of each pixel by solving the coefficient matrix of the thin plate spline function to obtain a rough deformation field mapping; A three-layer multi-resolution strategy is used to optimize the deformation field, step by step from low resolution to original resolution; At each resolution level, the displacement vector of each pixel is iteratively updated by minimizing the combined objective function of the image mutual information and the sum of squares of the second-order derivatives of the deformation field; The optimized deformation field is post-processed and vector median filtering is used to remove local abnormal displacements. The topological consistency of the deformation field is maintained by the elastic mesh method, and the mesh node positions are iteratively optimized until convergence or the preset maximum number of iterations is reached; According to the post-processing results, the local deformation field mapping is obtained; The local deformation field map is used to characterize the deformation relationship of the region of interest between different time points.

8. The method according to claim 1, wherein: According to the deformation field mapping, deformation compensation is performed on the region of interest at different time points to achieve accurate correspondence between the CT images of the region of interest at different time points. Through deformation compensation, the anatomical structure and the region of interest at different time points are aligned to the same reference coordinate system to complete the tracking task, including: Acquire the region of interest at different time points according to the deformation field mapping, and perform deformation compensation on the region of interest; From the CT images at different time points, determine the CT image at the middle time point with the highest signal-to-noise ratio as the reference coordinate system; Calculating an affine transformation matrix using a least squares method, and globally aligning the region of interest using the affine transformation matrix; The local deformation field is used to finely align the region of interest; The mutual information value and correlation coefficient before and after the registration are calculated for the deformation compensation result. If the mutual information value is less than a preset mutual information value threshold or the correlation coefficient is less than a preset correlation coefficient threshold, the gradient descent method is used to locally optimize the deformation compensation result. Shape features, density features and texture features are extracted from the aligned regions of interest, and feature correspondences between different time points are established based on the shape features, density features and texture features.

Images