A method for tracking a region of interest of a lung CT image

Abstract

The application provides a tracking method for a lung CT image region of interest, comprising: acquiring a lung CT image sequence of a patient at multiple time points, pre-processing the CT image sequence to obtain a pre-processed CT image sequence; for the pre-processed CT image sequence, extracting lung anatomical structure features in each CT image by using an anatomical structure-based feature extraction algorithm, and constructing a feature descriptor; performing shape analysis on a region of interest contour, extracting geometric features and topological features of the contour, and establishing a shape model of the region of interest; based on the shape model, estimating a local deformation field between regions of interest at different time points by using a non-parametric registration method, calculating a displacement vector of each pixel point, and obtaining a deformation field mapping.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of information technology, and particularly relates to a tracking method for a region of interest in a lung CT image. BACKGROUND

[0002] In the tracking of lung CT images, the respiratory motion and body position change of a patient result in spatial position and shape differences between CT images collected at different time points, which brings challenges to the accurate tracking of a region of interest. Traditional registration methods based on gray value or texture features are difficult to cope with such complex deformation and are prone to false matching. Specifically, the lung tissue undergoes non-rigid deformation during the breathing process, and the motion amplitude and direction of different lung lobes and segments are different. At the same time, the deformation mode of a lesion region such as a tumor is different from that of the surrounding normal tissue. Such complex local deformation brings difficulties to feature extraction and matching. How to balance and fuse feature information at different scales, and how to effectively use prior knowledge to constrain and guide the feature matching process, are problems that need to be further explored. How to improve the robustness and adaptability of the algorithm while ensuring the tracking accuracy, so that it can cope with various complex clinical scenarios, is still a challenging research topic. SUMMARY

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

[0004] Obtaining a lung CT image sequence of a patient at multiple time points, pre-processing the CT image sequence to obtain a pre-processed CT image sequence;

[0005] For the pre-processed CT image sequence, a feature extraction algorithm based on anatomical structure is used to extract the lung anatomical structure features in each CT image, and a feature descriptor is constructed;

[0006] According to the feature descriptor, an initial registration relationship is established between CT images at different time points, and based on the difference of anatomical structure between different individuals, a prior model of anatomical structure is constructed by using a group statistical method to guide the image registration process;

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

[0008] The shape of the region of interest contour is analyzed, the geometric features and topological features of the contour are extracted, and a shape model of the region of interest is established;

[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, calculate the displacement vector of each pixel point, and obtain the deformation field mapping;

[0010] According to the deformation field mapping, the deformation compensation of the regions of interest at different time points is performed, the accurate correspondence between the regions of interest at different time points is realized, the anatomical structures and the regions of interest at different time points are aligned to the same reference coordinate system through deformation compensation, and the tracking task is completed.

[0011] The technical scheme provided by the embodiment of the present application can include the following beneficial effects:

[0012] The present application discloses a tracking method for a region of interest of a lung CT image, which is suitable for processing lung CT images obtained at different time points. Through pre-processing to optimize image quality, the interference of noise and artifacts is eliminated, laying a solid foundation for subsequent steps; a feature extraction algorithm based on anatomical structure is used to accurately capture the lung structure features and construct stable feature descriptors, which can maintain good matching performance even under complex deformation; in the positioning and segmentation of the region of interest (ROI), the application of region growing algorithm or level set method enables accurate extraction of the ROI contour even under complex background and low contrast conditions, providing key information for further shape analysis and registration; the establishment of the shape model not only captures the geometric and topological features of the contour, but also provides a morphological reference for non-parametric registration, enabling the algorithm to more accurately estimate the local deformation field, calculate the pixel point displacement, and obtain accurate deformation field mapping; the deformation compensation method based on the deformation field mapping ensures the accurate correspondence of the ROI between CT images at different time points, maintains data consistency even in the face of significant changes within the respiratory cycle, and greatly improves the accuracy and reliability of tracking. Not only does it solve the limitations of traditional methods in dealing with complex deformation, but also significantly enhances the robustness and adaptability of the algorithm, enabling it to handle various clinical scenarios, especially in dynamic monitoring and pathological development analysis. In general, the technical effects of the present application 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 handling dynamic changes and pathological development analysis. In addition, through advanced image processing and data analysis techniques, 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 DRAWINGS

[0013] Fig. 1 A flowchart of a tracking method for a region of interest of a lung CT image according to the present application.

[0014] Fig. 2A schematic diagram of a lung CT image region of interest tracking method according to the present application.

[0015] Fig. 3 Another schematic diagram of a lung CT image region of interest tracking method according to the present application. DETAILED DESCRIPTION

[0016] The technical solutions in the embodiments of the present application will be described clearly and in detail below with reference to the drawings in the embodiments of the present application. The described embodiments are only some of the embodiments of the present application.

[0017] As Figs. 1-3 , the lung CT image region of interest tracking method according to the present application can specifically include the following steps.

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

[0019] A lung CT image sequence of a patient at multiple time points is acquired, lung parenchyma region segmentation is performed on the CT image, a threshold range is set according to the gray value distribution of the CT image, and a segmented lung CT image sequence is obtained. The segmented CT image sequence is converted in format, the CT image sequence is uniformly converted into a DICOM format, the resolution of the CT image is adjusted to a preset pixel value. The CT images at different time points in the CT image sequence are aligned, the similarity between the CT images is calculated, the registration parameters are optimized by using a gradient descent method, and an aligned CT image sequence is obtained. The aligned CT image sequence is subjected to noise suppression and contrast enhancement, a Gaussian filter is used to remove CT image noise and retain edge information, and the CT image gray value distribution is adjusted by a histogram equalization method. The CT image sequence after adjusting the CT image gray value distribution is subjected to brightness and contrast standardization processing, if the brightness value of the CT image is lower than a preset threshold, the brightness of the CT image is enhanced, if the contrast value of the CT image is lower than a preset threshold, the contrast of the CT image is enhanced, and a preprocessed CT image sequence is obtained.

[0020] Specifically, lung CT image sequences were acquired from patients at multiple time points. A threshold segmentation algorithm was used to segment the lung parenchyma region of each CT image at each time point. A threshold range was set based on the image grayscale distribution to separate the lung parenchyma from surrounding tissue, resulting in a segmented lung CT image sequence. The segmented CT image sequence was then converted to DICOM format, with images of varying formats uniformly converted. The image resolution was adjusted to 512x512 pixels to ensure consistent format and resolution across all images. A grayscale-based image registration algorithm was then applied to the converted CT image sequence to align the CT images at different time points. Image similarity was calculated using the maximum mutual information criterion, and registration parameters were optimized using gradient descent to obtain the registered CT image sequence. Noise suppression and contrast enhancement were performed on the registered CT image sequence. A Gaussian filter was used to remove image noise while preserving edge information. Histogram equalization was used to adjust the image grayscale distribution. The Laplacian operator was applied to enhance the boundary definition of the lesion region. Finally, brightness and contrast were normalized to ensure comparability between images at different time points, resulting in the preprocessed CT image sequence. Lung CT image sequences were acquired at multiple time points, each containing 30-50 cross-sectional images. A threshold segmentation algorithm was used to segment the lung parenchyma, with a grayscale threshold ranging from -950 to -500 HU. A region growing algorithm was used to expand the region from the center of the lung until the threshold boundary was reached, yielding a preliminary segmentation result. Morphological operations were applied to fill in small pulmonary vessels and bronchi, and a 3x3 structuring element was used for closing to remove small noise. The segmented images were converted to a unified DICOM format and resized to 512x512 pixels with a pixel pitch of 1 mm. Grayscale-based registration was performed on the converted image sequences, with the first time point image selected as the reference image. Mutual information (MI) was calculated as the similarity metric. Affine transformation parameters, including translation, rotation, and scaling, were iteratively adjusted using the Powell optimization algorithm until the MI converged or the maximum number of iterations reached 100. A Gaussian filter was applied to the registered images for denoising, with a kernel size of 5x5 and a standard deviation of σ = 1.5. Histogram equalization was used to stretch the image grayscale range to 0–4095 to increase image contrast. A 3x3 Laplacian operator was used to sharpen image edges and highlight the contours of the lesion area. Finally, brightness and contrast were normalized, with the mean grayscale value adjusted to 2048 and the standard deviation to 500 to ensure consistent grayscale distribution across images at different time points for ease of subsequent analysis and comparison.

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

[0022] The lung lobes are segmented by using a region growing algorithm, the gray scale threshold range of the lung parenchyma region is determined according to the gray scale histogram of the CT image, the seed points of the left and right lungs are selected in the gray scale threshold range, and the region is gradually expanded from the seed points to obtain the mask image of the segmented left and right lung lobes. Based on the segmented lung image, the trachea and blood vessels are enhanced by using a multi-scale Hessian matrix analysis method, whether a pixel point belongs to a tubular structure is judged by comparing the size relationship of the eigenvalues of the Hessian matrix, if the eigenvalues meet the preset requirements, the tubular structure is determined and is subjected to enhancement processing to obtain the enhanced trachea and blood vessel image. The three-dimensional structure of the trachea and blood vessels is extracted by applying a three-dimensional connected domain analysis method to the enhanced trachea and blood vessel image, the small volume connected domains are removed by setting a volume threshold, the trachea and blood vessels are distinguished by morphological characteristics to obtain the three-dimensional skeletons of the trachea tree and the blood vessel tree. Based on the three-dimensional skeletons of the trachea tree and the blood vessel tree, the geometric characteristics and statistical characteristics of each structure are calculated, including calculating the branch number, branch angle and tube diameter change rate of the trachea, calculating the volume fraction, blood vessel tortuosity and blood vessel distribution density of the blood vessels, and calculating the volume, surface area, average density, density standard deviation, shape factor and curvature of the lung lobes. The calculated characteristics are normalized and combined into a feature vector to construct the feature descriptor of each CT image.

[0023] Specifically, for the preprocessed CT image sequence, the lung lobes are segmented by using a region growing algorithm, the gray scale threshold range of the lung parenchyma region is determined by calculating the image gray scale histogram, the seed points of the left and right lungs are automatically selected in the range, and the region is gradually expanded from the seed points to obtain the mask image of the segmented left and right lung lobes. The trachea and blood vessels are enhanced by using a multi-scale Hessian matrix analysis on the segmented lung image, the Hessian matrix of the image under different scales is calculated, whether a pixel point belongs to a tubular structure is judged by comparing the size relationship of the eigenvalues λ1, λ2, λ3 of the matrix, if |λ1|≤|λ2|≤|λ3|

[0024] and if |λ1| approaches 0, it is determined as a tubular structure and enhanced to obtain the enhanced trachea and blood vessel images, wherein the eigenvalues λ1, λ2, λ3 are three eigenvalues of the Hessian matrix at a point of the 3D image. The enhanced images are subjected to three-dimensional connected domain analysis, and the three-dimensional structures of the trachea and blood vessels are extracted. A volume threshold is set to remove small volume connected domains. The trachea and blood vessels are distinguished by morphological features. The trachea is a tree-like branch structure, and the blood vessels are a network-like distribution. The three-dimensional skeletons of the trachea tree and blood vessel tree are obtained. Based on the extracted anatomical structures, the geometric and statistical features of each structure are calculated. The branch number, branch angle, and tube diameter change rate of the trachea are calculated. The volume fraction, vessel tortuosity, and vessel distribution density of the blood vessels are calculated. The volume, surface area, average density, density standard deviation, shape factor, and curvature of the lung lobes are calculated. Meanwhile, the lung nodules and lung interstitium are identified by morphological operation and density threshold method, and the size, shape, and edge features thereof are extracted. After normalization processing, the features are combined into a feature vector to construct a feature descriptor of each CT image.

[0025] S103, establishing an initial registration relationship between the CT images at different time points according to the feature descriptors; based on the difference of the anatomical structures between different individuals, a prior model of the anatomical structures is constructed by using a population statistical method to guide the image registration process.

[0026] CT images at two different time points are acquired, a nearest neighbor search algorithm is used to establish a corresponding relationship of feature points, and an initial matching point pair is obtained by calculating the Euclidean distance of the feature vectors. According to the set distance threshold and the number of iterations, the random sample consensus method is used to eliminate abnormal matching points, and a stable matching point set is selected. For the anatomical structures in the population CT images, the covariance matrix is calculated and the eigenvalues and eigenvectors are solved, the first k principal components are selected, and a statistical shape model and a density distribution model are constructed. Based on the initial registration relationship and the prior model, a B-spline free-form deformation algorithm is used for non-rigid registration; the registration result is quantitatively evaluated, the normalized mutual information and the correlation coefficient before and after registration are calculated, and the average distance error of the corresponding points is calculated combined with the landmark points automatically selected based on the image gradient and anatomical knowledge. If the average distance error is greater than a preset threshold, the non-rigid registration is performed again, the control point density and the weight coefficient are adjusted, and parameter optimization and re-registration are performed.

[0027] Specifically, according to the feature descriptors, a nearest neighbor search algorithm is used to establish a corresponding relationship of feature points between CT images at different time points, and the Euclidean distance of the feature vectors is calculated to select the feature point pair with the smallest distance as the initial matching point pair, thereby obtaining an initial registration relationship. The random sample consensus algorithm is used to eliminate abnormal matching points, the distance threshold and the number of iterations are set, and a stable matching point set is selected. Principal component analysis is performed on the anatomical structures in the population CT images, the covariance matrix is calculated, the eigenvalues and eigenvectors are solved, the first k principal components are selected, and a statistical shape model and a density distribution model are constructed. The shape model is represented as the average shape plus the linear combination of the principal components, and the density model is represented as the average density distribution plus the linear combination of the principal components. Based on the initial registration relationship and the prior model, a B-spline free-form deformation algorithm is used for non-rigid registration, a control point grid is defined, and a combined objective function

[0028] E = aEsim(I1, T(I2)) + bEshape(T) where Esim is the image similarity term, which is measured by mutual information, Eshape is the shape constraint term based on prior model, I1 is the reference image, I2 is the template image, T is the deformation field, T(I2) is the deformed image, a and b are the weight coefficients, which determine the importance of the image similarity term and the shape constraint term in the objective function, respectively, and are used to balance the contributions of image similarity and shape constraint. Gradient descent method is used to iteratively optimize the objective function to obtain the accurate spatial correspondence. The registration results are quantitatively evaluated, and the normalized mutual information and correlation coefficient before and after registration are calculated, and the average distance error of the corresponding points is calculated combined with the automatically selected landmark points based on 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 the weight coefficient, and perform parameter optimization and re-registration. In the implementation process, first, extract 128-dimensional feature descriptors from each CT image, use the k-d tree algorithm for nearest neighbor search, set the search radius to 0.1, and find the nearest neighbor of each feature point. Then use the random sample consensus algorithm, set the distance threshold to 2mm, and the iteration number to 1000, randomly select 3 pairs of points from the initial matching points, calculate the affine transformation matrix, count the number of point pairs that meet the transformation, and repeat the process to obtain the best matching point set. Perform principal component analysis on the CT images of 100 patients, calculate the covariance matrix, solve the eigenvalues and eigenvectors, select the first 10 principal components, and explain 95% of the shape variation. Construct a statistical shape model where x is the shape vector, is the mean shape, P is the principal component matrix, and b is the shape parameter vector. A density distribution model can also be constructed similarly. In non-rigid registration, a cubic B-spline free-form deformation algorithm is used, a control point grid of 10x10x10 is set, and a 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, the normalized mutual information measures the degree of information sharing between two images and is normalized to a fixed interval, usually [0, 1]. The 0.7NMI(I1, T(I2)) part measures the similarity between the reference image I1 and the template image T(I2) transformed by the deformation field T. The L-BFGS algorithm is used to optimize the objective function, with a maximum iteration number of 100 and a convergence threshold of 1e-6. In the evaluation stage, 50 anatomical landmarks are automatically selected, including the tracheal bifurcation point, the lung lobe fissure junction, etc., and the root mean square error of the corresponding distance of the landmarks is calculated, with a threshold of 2mm. If the threshold is exceeded, the control point grid is encrypted to 15x15x15, the weight is adjusted to 0.8 and 0.2, and the registration optimization is performed again.

[0029] S104, in the image registration CT image sequence, locate and segment the region of interest, using region growing algorithm or level set method, get the region of interest contour in each CT image.

[0030] The registered CT image sequence is obtained, the image is smoothed by using a Gaussian filter according to the CT image sequence, the image contrast is enhanced by adaptive histogram equalization, and the enhanced CT image sequence is obtained. According to the enhanced CT image sequence, the initial position of the region of interest is determined by using the Canny edge detection algorithm according to the image gray value distribution characteristics. According to the initial position, the region growing algorithm is used to accurately segment the region of interest, and the region growing algorithm includes adding spatial constraint and adaptive threshold. For the segmented region of interest, the morphological operation is applied to optimize the region boundary, and the accurate contour of the region of interest is extracted by combining the Sobel edge detection algorithm, and the final contour of the region of interest in each CT image is obtained. If the region growing algorithm cannot obtain the accurate 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] In detail, the pre-processing of the registered CT image series was performed, and a Gaussian filter was used to smooth the images with a kernel size of 5x5 and a standard deviation of 1.5. Adaptive histogram equalization was used to enhance the contrast of the images, and the images were divided into 8x8 blocks with a contrast limit of 0.02. The initial position of the region of interest was located automatically using the image gray value distribution characteristics and the Canny edge detection algorithm, with a low threshold of 50 and a high threshold of 150. The rough outline of the region of interest was determined by calculating the image gradient amplitude and direction and combining the region growing algorithm. Based on the positioning results, the improved region growing algorithm was used to accurately segment the region of interest, including adding spatial constraints and adaptive threshold. Multiple seed points were selected from the initial outline, and the seed points were selected based on the image gradient and gray value. The pixel points with the smallest gradient and gray value within the target region average gray value ±10% range were selected, and the initial gray value similarity threshold was set to 20. The spatial connectivity constraint was 26 neighbors, and the region growing criteria were iteratively updated to obtain the segmented region of interest. Morphological operations were applied to optimize the region boundary, and a 3x3 circular structure element was used for opening operation to remove small protrusions, and a 5x5 circular structure element was used for closing operation to fill the boundary gaps. The accurate outline of the region of interest was extracted using the Sobel edge detection algorithm, and the final outline of the region of interest in each CT image was obtained. If the region growing algorithm cannot obtain satisfactory results, the level set method is used as an alternative. The signed distance function is used to initialize the level set, and the evolution speed function is combined with the image gradient and region statistical information. The iteration number is set to 200 to obtain the outline of the region of interest. In practical applications, the registered CT image series is processed first using a 5x5 Gaussian filter for smoothing with a sigma value of 1.5 to effectively remove noise while preserving edge information. Then adaptive histogram equalization is applied, and the image is divided into 8x8 blocks with a contrast limit of 0.02 to significantly improve the image contrast. On the enhanced image, the Canny edge detection algorithm is used with a low threshold of 50 and a high threshold of 150 to accurately detect the edges of the region of interest. Combined with the detected edge information, the improved region growing algorithm is used for segmentation. The algorithm selects 10 seed points from the initial outline, and the gradient value of these points is less than 5 and the gray value is within the target region average gray value ±10% range. The initial similarity threshold is set to 20, and the 26-neighborhood is used as a spatial constraint. After 50 iterations, the preliminary segmentation result is obtained. Morphological operations are applied to the segmentation result, and a 3x3 circular structure element is used for opening operation to remove protrusions with an area less than 10 pixels, and a 5x5 circular structure element is used for closing operation to fill holes with an area less than 20 pixels. Finally, the Sobel operator is used to extract the accurate outline with a gradient threshold of 100.If the Dice coefficient of the region growing algorithm is lower than 0.85, switch to the level set method. The level set initialization uses distance transform, the curvature term weight in the evolution speed function is set to 0.1, the region term weight is 0.6, the edge term weight is 0.3, and iteration is performed 200 times or until the contour changes less than 0.1%.

[0032] S105, shape analysis is performed on the region of interest contour, geometric features and topological features of the contour are extracted, and a shape model of the region of interest is established.

[0033] Geometric features of the region of interest contour are obtained, the geometric features include area, perimeter, major axis length, minor axis length, circularity and eccentricity; a Douglas-Peucker algorithm is used for polygon fitting of the contour to obtain the number of vertices and the internal angle sum of the contour; a 3rd order B-spline curve is used for fitting according to the contour point coordinates 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 as an inflection point, and if the curvature value is greater than a preset multiple of the average curvature, it is determined as a sharp point; a skeleton of the region of interest is obtained through a medial axis transformation algorithm, the number of branches and the length distribution of the skeleton are calculated; the connectivity of the region of interest is judged according to the Euler number, the Euler number is equal to the number of regions minus the number of holes; principal component analysis method is used to reduce the feature dimension for the extracted geometric features and topological features; according to the eigenvalues and eigenvectors of the feature covariance matrix, principal components with a cumulative contribution rate greater than a preset threshold are selected as shape descriptors; weight coefficients of the shape descriptors are combined into a vector to obtain a compact shape representation of the region of interest.

[0034] Specifically, geometric feature extraction is performed on the contour of the region of interest, the area and perimeter are calculated through the contour point coordinates, the length of the long axis and the short axis are obtained by using the minimum circumscribed rectangle algorithm, the circularity and the eccentricity are calculated, the Douglas-Peucker algorithm is used for polygon fitting, the threshold is set to 0.5% of the contour perimeter, the number of contour vertices and the internal angle sum are obtained. Curvature analysis is performed on the contour, a 3-order B-spline curve is used to fit the contour point sequence, 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, the mean, variance and maximum of the curvature are counted, the inflection points of the contour are extracted by setting the curvature threshold to 2 times the average curvature, the points with curvature greater than 3 times the average curvature are defined as sharp points, and the local shape features of the contour are obtained. Topological features of the region of interest are extracted, the skeleton of the region is obtained by using the medial axis transformation algorithm, the number of skeleton branches and the length distribution are calculated, the connectivity of the region is calculated by using the Euler number, and the Euler number is equal to the number of regions minus the number of holes. The thickness distribution of the region is obtained by distance transformation. Based on the extracted geometric and topological features, a shape model of the region of interest is constructed, principal component analysis method is used to reduce the feature dimension, eigenvalues and eigenvectors of the feature covariance matrix are calculated, the first K principal components are selected as shape descriptors according to the principle that the cumulative contribution rate is greater than 95%, and the weight coefficients of the K principal components are combined into a vector, which is used as the compact shape representation of the region of interest. In practical application, when the shape of the contour of the region of interest of the lung nodule is analyzed, geometric features are first extracted. The area is calculated by Green's formula, and the area is 3.14 cm2; the perimeter is calculated by accumulating the Euclidean distance between the contour points, and the perimeter is 6.28 cm. The minimum circumscribed rectangle algorithm obtains the long axis of 2.5 cm and the short axis of 1.6 cm, and the circularity is 0.79 and the eccentricity is 0.77. Douglas-Peucker algorithm is used for polygon fitting, and the threshold is set to 0.0314 cm, which is 0.5% of the perimeter, and 12 vertices are obtained, and the internal angle sum is 1800°. Then curvature analysis is performed, a 3-order B-spline curve is used to fit the contour, the number of control points is 63, which is 1 / 10 of the number of contour points. The average curvature is calculated to be 0.5 cm -1 , the curvature variance is 0.04 cm -2 , and the maximum curvature is 1.2 cm -1 . The curvature threshold is set to 1.0 cm -1 , which is 2 times the average curvature, and 5 inflection points are extracted; the curvature greater than 1.5 cm -1The points defined as cusps are detected as 2. In the topology feature extraction, the skeleton is obtained by the axis transform algorithm, containing 3 main branches, and the longest branch length is 1.8 cm. The Euler number is calculated as 1, in which 1 region and 0 hole. The average thickness is 0.8 cm and the maximum thickness is 1.2 cm by the distance transform. Finally, the shape model is constructed, the dimension is reduced by the principal component analysis, the eigenvalues and eigenvectors of the characteristic covariance matrix are calculated, and the first 6 principal components corresponding to the cumulative contribution rate of 95% are selected as the shape descriptors. The weight coefficients of the 6 principal components form a vector [-0.35, 0.28, -0.15, 0.10, -0.08, 0.05], which is used as the compact shape representation of the lung nodule.

[0035] In S106, 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, calculate the displacement vector of each pixel point, and obtain the deformation field mapping.

[0036] According to the feature points extracted from the shape model, the corresponding relationship 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 as a stable matching point pair, which is used as the control point for deformation field estimation. An initial deformation field is constructed, and the displacement of each pixel point is estimated 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, from low resolution to original resolution step by step. At each resolution level, the displacement vector of each pixel point is iteratively updated by minimizing the combined objective function of the image mutual information and the square sum of the second derivative of the deformation field. The optimized deformation field is post-processed, and the local abnormal displacement is removed by vector median filtering. 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 iteration number is reached. According to the post-processing result, the local deformation field mapping is obtained. The local deformation field mapping is used to represent the deformation relationship between the regions of interest at different time points.

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

[0038] S107, according to the deformation field mapping, performing deformation compensation on the region of interest at different time points to realize accurate correspondence of the region of interest between CT images at different time points; through the deformation compensation, aligning the anatomical structure and the region of interest at different time points 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 the region of interest is subjected to deformation compensation. The CT image at the intermediate time point with the highest signal-to-noise ratio is determined from the CT images at different time points as a reference coordinate system. An affine transformation matrix is calculated by using a least square method, and the region of interest is globally aligned by using the affine transformation matrix; and the region of interest is finely aligned by using a local deformation field. The mutual information value and the correlation coefficient before and after 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, then 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 a feature correspondence relationship between different time points is established according to the shape features, the density features and the texture features.

[0040] In particular, according to the deformation field mapping, the deformation compensation is performed on the regions of interest at different time points by using the bilinear interpolation algorithm, the four nearest neighbor pixel values of each target pixel point in the original image are calculated, the pixel values in the original image are mapped to the new positions according to the distance weight for weighted average, and the deformation compensated image is obtained. The CT image at the middle time point with the highest signal-to-noise ratio is selected as the reference coordinate system, the regions of interest at other time points are aligned to the reference coordinate system through deformation compensation, the affine transformation matrix is calculated by using the least square method for global alignment, and then the local deformation field is applied for fine alignment. The deformation compensation result is quantitatively evaluated, the mutual information value and the correlation coefficient before and after registration are calculated, the mutual information value threshold is set to 1.5, the correlation coefficient threshold is set to 0.95, the alignment quality is judged, 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 or the maximum iteration number 100 is reached. Based on the aligned region of interest, 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 the gray level co-occurrence matrix, the corresponding relationship between the features at different time points is established, the feature changes at the four intermediate time points are estimated by using the cubic spline interpolation algorithm, and the continuous tracking of the region of interest is realized. In practical application, when the deformation compensation and tracking are performed on the CT images of the lung nodules at different time points, first, the pixel resampling is performed by using the bilinear interpolation algorithm 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 average of the adjacent 2x2 pixel blocks in the original image, and the weight is calculated based on the inverse distance. The CT 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, the affine transformation matrix of 3x3 is calculated by using the least square method for global alignment of the images at other time points, and the error threshold is set to 0.5 pixels. Then, the local deformation field is applied for accurate sub-pixel level. The deformation compensation result is evaluated, the mutual information value and the correlation coefficient are calculated, 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 number is 100. For the aligned region of interest, the shape features such as volume, surface area and maximum diameter are extracted, the density features include average CT value, standard deviation and texture features, and the gray level co-occurrence matrix is used, the window size is 5x5, and the energy, contrast and correlation are calculated. For example, for a lung nodule with a diameter of 2 cm, the volume is about 4.19 cm 2 , the average CT value is -600 HU, and the standard deviation is 150 HU. The feature changes at the four intermediate time points between the five time points are estimated by using the cubic spline interpolation algorithm, the interpolation node number is set to 20, the continuous tracking of the nine time points is realized, and the dynamic change process of the lung nodule is completely described.

[0041] The above description is only the preferred embodiment of one or more embodiments of the specification, and is not used to limit one or more embodiments of the specification. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of one or more embodiments of the specification should be included in the protection range of one or more embodiments of the specification.

Claims

1. A method for tracking a region of interest in a lung CT image, characterized by, The method comprises: acquiring a lung CT image sequence of a patient at multiple time points, pre-processing the CT image sequence to obtain a pre-processed CT image sequence; for the pre-processed CT image sequence, a feature extraction algorithm based on an anatomical structure is used to extract the lung anatomical structure features in each CT image to construct a feature descriptor; according to the feature descriptor, an initial registration relationship between CT images at different time points is established, and based on the difference of anatomical structures between different individuals, a prior model of the anatomical structure is constructed by using a population statistical method to guide the image registration process; in the CT image sequence after image registration, the region of interest is located and segmented, and a region growing algorithm or a level set method is used to obtain the region of interest contour in each CT image; the shape of the region of interest contour is analyzed, the geometric features and topological features of the contour are extracted, and a shape model of the region of interest is established; 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, the displacement vector of each pixel point is calculated, and a deformation field mapping is obtained; according to the deformation field mapping, the regions of interest at different time points are deformed and compensated, the accurate correspondence between the regions of interest at different time points is realized, the anatomical structures and the regions of interest at different time points are aligned to the same reference coordinate system through deformation compensation, and the tracking task is completed.

2. The method of claim 1, wherein, The method comprises: acquiring a lung CT image sequence of a patient at multiple time points, pre-processing the CT image sequence to obtain a pre-processed CT image sequence; for the pre-processed CT image sequence, a feature extraction algorithm based on an anatomical structure is used to extract the lung anatomical structure features in each CT image to construct a feature descriptor; according to the feature descriptor, an initial registration relationship between CT images at different time points is established, and based on the difference of anatomical structures between different individuals, a prior model of the anatomical structure is constructed by using a population statistical method to guide the image registration process; in the CT image sequence after image registration, the region of interest is located and segmented, and a region growing algorithm or a level set method is used to obtain the region of interest contour in each CT image; the shape of the region of interest contour is analyzed, the geometric features and topological features of the contour are extracted, and a shape model of the region of interest is established; 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, the displacement vector of each pixel point is calculated, and a deformation field mapping is obtained; according to the deformation field mapping, the regions of interest at different time points are deformed and compensated, the accurate correspondence between the regions of interest at different time points is realized, the anatomical structures and the regions of interest at different time points are aligned to the same reference coordinate system through deformation compensation, and the tracking task is completed.

3. The method of claim 1, wherein, The method comprises the following steps: for the preprocessed CT image sequence, anatomical structure feature extraction algorithm is adopted to extract lung anatomical structure features in each CT image, and a feature descriptor is constructed, including: a region growing algorithm is adopted to segment lung lobes, a lung parenchyma region gray threshold range is determined according to a CT image gray histogram, seed points of left and right lungs are selected in the gray threshold range, and a region is gradually expanded from the seed points to obtain a segmented left and right lung lobe mask image; based on the segmented lung image, a multi-scale Hessian matrix analysis method is used for trachea and blood vessel enhancement, whether a pixel point belongs to a tubular structure is judged by comparing the size relationship of Hessian matrix eigenvalues, if the eigenvalues meet preset requirements, the tubular structure is determined and enhancement processing is performed, and an enhanced trachea and blood vessel image is obtained; a three-dimensional connected domain analysis method is applied to the enhanced trachea and blood vessel image, three-dimensional structures of the trachea and blood vessels are extracted, a volume threshold is set to remove small volume connected domains, the trachea and blood vessels are distinguished through morphological features, and a three-dimensional skeleton of a trachea tree and a blood vessel tree is obtained; based on the three-dimensional skeleton of the trachea tree and the blood vessel tree, geometric features and statistical features of each structure are calculated, including calculating branch number, branch angle and tube diameter change rate for the trachea, calculating volume fraction, blood vessel tortuosity and blood vessel distribution density for the blood vessels, and calculating volume, surface area, average density, density standard deviation, shape factor and curvature for the lung lobes; the calculated features are normalized and combined into a feature vector to construct a feature descriptor of each CT image.

4. The method of claim 1, wherein, According to the feature descriptor, an initial registration relationship is established between CT images at different time points, based on the difference of anatomical structures between different individuals, a prior model of anatomical structures is constructed by using a population statistics-based method to guide the image registration process, including: acquiring CT images at two different time points, using a nearest neighbor search algorithm to establish a feature point correspondence relationship, and calculating an initial matching point pair by calculating the Euclidean distance of the feature vector; according to the set distance threshold and the number of iterations, using a random sample consensus method to remove abnormal matching points and selecting a stable matching point set; for the anatomical structures in the population CT image, a covariance matrix is calculated and eigenvalues and eigenvectors are solved, the first k principal components are selected, a statistical shape model and a density distribution model are constructed; based on the initial registration relationship and the prior model, a B-spline free form deformation algorithm is used for non-rigid registration; the registration result is quantitatively evaluated, the normalized mutual information and the correlation coefficient before and after registration are calculated, the corresponding point average distance error is calculated by combining the landmark points selected automatically based on image gradients 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 the weight coefficient are adjusted, parameter optimization and re-registration are performed.

5. The method of claim 1, wherein, The CT image sequence after image registration is positioned and segmented, and a region growing algorithm or a level set method is used to obtain the contour of the region of interest in each CT image, including: obtaining the CT image sequence after registration, smoothing the image using a Gaussian filter based on the CT image sequence, enhancing the contrast of the image through adaptive histogram equalization to obtain an enhanced CT image sequence; for the enhanced CT image sequence, the initial position of the region of interest is determined using a Canny edge detection algorithm based on the image gray value distribution characteristics; based on the initial position, a region growing algorithm is used to accurately segment the region of interest, the region growing algorithm includes adding spatial constraints and adaptive threshold; for the segmented region of interest, morphological operations are applied to optimize the region boundary, and the accurate contour of the region of interest is extracted in combination with a 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 accurate contour of the region of interest, a 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 of claim 1, wherein, The shape of the region of interest contour is analyzed, the geometric and topological features of the contour are extracted, and a shape model of the region of interest is established, including: obtaining the geometric features of the region of interest contour, the geometric features including area, perimeter, major axis length, minor axis length, circularity and eccentricity; the Douglas-Peucker algorithm is used for polygon fitting of the contour to obtain the number of contour vertices and the internal angle sum; a 3rd order B-spline curve is used for fitting based on the contour point coordinates to obtain the curvature distribution of the contour; if the curvature value is greater than the average curvature by a preset value, it is determined as a inflection point, and if the curvature value is greater than the average curvature by a preset multiple, it is determined as a sharp point; the skeleton of the region of interest is obtained through the center axis transformation algorithm, and the branch number and length distribution of the skeleton are calculated; the connectivity of the region of interest is judged according to the Euler number, and the Euler number is equal to the number of regions minus the number of holes; for the extracted geometric and topological features, principal component analysis is used to reduce the feature dimension; according to the eigenvalues and eigenvectors of the feature covariance matrix, the principal components with a cumulative contribution rate greater than a preset threshold are selected as shape descriptors; 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 of claim 1, wherein, The non-parametric registration method is used to estimate the local deformation field between the regions of interest at different time points based on the shape model, and the displacement vector of each pixel point is calculated to obtain the deformation field mapping, including: establishing the corresponding relationship between the regions of interest at different time points according to the feature points extracted from the shape model, and obtaining the similarity measure by calculating the Euclidean distance of the feature descriptors; if the similarity measure is less than a preset threshold, it is determined as a stable matching point pair as the control point for deformation field estimation; an initial deformation field is constructed, and the displacement of each pixel point 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 by using a three-layer multi-resolution strategy, which is gradually performed from low resolution to original resolution; at each resolution level, the displacement vector of each pixel point is iteratively updated by minimizing the combined objective function of the image mutual information and the square sum of the second derivative of the deformation field; the optimized deformation field is post-processed, and the local abnormal displacement is removed by using the vector median filter; the topological consistency of the deformation field is maintained by using the elastic grid method, and the grid node position is iteratively optimized until convergence or the preset maximum iteration number is reached; the local deformation field mapping is obtained according to the post-processing result; and the local deformation field mapping is used to represent the deformation relationship between the regions of interest at different time points.

8. The method of claim 1, wherein, According to the deformation field mapping, the deformation compensation of the regions of interest at different time points is performed to realize the accurate correspondence of the regions of interest between the CT images at different time points, and the anatomical structures and the regions of interest at different time points are aligned to the same reference coordinate system through deformation compensation to complete the tracking task, including: obtaining the regions of interest at different time points according to the deformation field mapping, and performing deformation compensation on the regions of interest; determining the CT image at the intermediate time point with the highest signal-to-noise ratio as the reference coordinate system from the CT images at different time points; calculating the affine transformation matrix by using the least squares method, and globally aligning the regions of interest by using the affine transformation matrix; finely aligning the regions of interest by using the local deformation field; calculating the mutual information value and the correlation coefficient before and after registration 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, then the gradient descent method is used to locally optimize the deformation compensation result; extracting shape features, density features and texture features from the aligned regions of interest, and establishing the feature correspondence relationship between different time points according to the shape features, the density features and the texture features.

Images