Low-dose CT image reconstruction method guided by CT images of group patients

Abstract

The invention provides a low-dose CT image reconstruction method guided by CT images of group patients, and belongs to the technical field of medical image processing. According to the method, for each target region-of-interest structure of a reconstructed image, a similar example sample is selected from a conventional dose CT image of a group patient, then a priori subspace adaptive to the region-of-interest structure is learned by adopting a PCA method, and adaptive sparse modeling is performed on the priori subspace by adopting a weighted Laplacian model. And on the basis of the prior model, constructing a low-dose CT image reconstruction target function, and carrying out iterative optimization on the target function by adopting an alternating direction multiplier method (ADMM). According to the method, a priori subspace adaptive to each target area is learned from a high-quality example sample, and adaptive sparse modeling is performed on the subspace, so that the purpose of accurately capturing a local sparse structure is achieved, and high-quality reconstruction of the low-dose CT image with good texture retention capability is further realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of medical image processing, and more specifically, relates to a low-dose CT image reconstruction method guided by CT images of group patients (PI-MBIR for short). Background Art

[0002] Computed Tomography (CT) has become an important and versatile imaging modality, valued for its high resolution in screening and diagnostic applications. However, the X-rays used in CT scans produce ionizing radiation, which has caused great health concerns. In response to these risks, significant advances have been made in low-dose CT (LDCT) scanning technology, which effectively reduces radiation exposure to patients. Low-dose CT generally requires reducing the tube current and / or voltage of the X-ray source, but this reduction often results in increased image noise and artifacts, which may hinder accurate diagnosis and pose challenges in the clinical setting. How to minimize radiation exposure while maintaining high-quality CT images remains a core area of ​​medical imaging research.

[0003] Several noise suppression strategies have been proposed to address the problem of noise artifacts in LDCT. These methods include sinogram-domain filtering, image-domain denoising, and model-based iterative reconstruction (MBIR). The sinogram-domain filtering method directly smoothes the raw data and then reconstructs the CT image using the traditional filtered back projection (FBP) algorithm. Although this method is computationally efficient, the resulting CT images often show excessive smoothing of subtle tissue structures. In contrast, image-domain post-processing techniques focus on reducing noise and artifacts in reconstructed CT images, but they do not fully consider the physical properties of the noise source. MBIR combines the advantages of sinogram methods and image-domain methods by integrating prior knowledge of raw data statistics and CT image features into the optimization process. It optimizes the energy function by considering the physical system model, raw data statistics, and prior models of the object to be reconstructed. This optimization is achieved through an iterative algorithm and represents a significant advance in noise artifact mitigation in LDCT. The present invention focuses on the application of MBIR in LDCT imaging.

[0004] The development of accurate model-based iterative reconstruction in low-dose CT imaging relies on the use of reasonable prior models that are able to effectively represent the inherent structural details of CT images. Several previous models have been introduced to improve MBIR. For example, Liu et al. introduced an adaptive weighted total variation constraint method based on anisotropic edge features between adjacent image voxels, which provided a more nuanced approach for image reconstruction. Niu et al. combined penalized weighted least squares with generalized total variation techniques to improve LDCT image quality. Zhang et al. developed an LDCT reconstruction method based on adaptive non-local means (NLM) regularization, which focused on local image features and combined spatial adaptation to achieve more efficient reconstruction. Liu et al. proposed the use of Hessian Schatten penalty for CBCT reconstruction.

[0005] Extensive research has highlighted the significant progress achieved in LDCT imaging by introducing new MBIR methods. However, it is important to recognize that despite intensive studies on the statistical properties of raw data and the application of many innovative a priori models, the imaging results are still suboptimal and have not yet reached the stringent standards required for clinical applications. In this context, a major challenge is "photon starvation," which results in a significant loss of original information about specific human tissues in the raw data. As reported in previous studies, conventional MBIR methods mainly rely on these incomplete data to build imaging models. However, this endeavor is particularly challenging due to the inherent limitations of raw data under photon starvation conditions.

[0006] Whether for individual patients or from different patient populations, CT images often have similar radiological features such as texture and structure. In low-dose CT iterative reconstruction imaging, image quality enhancement can not only rely on inherent data, but also use knowledge mining to extract high-quality prior information from normal dose (NDCT) sample images. This additional prior information makes up for the lack of inherent data and potentially improves the quality of LDCT images. Recent research in this field has begun to explore such methods. A notable strategy is to use the patient's own previous NDCT images of the same or adjacent slices as exemplar images. This method uses image registration to extract prior information such as pixel intensity from the sample image using the self-similarity of the image. This information is then used to refine the regularization of the corresponding pixels in the LDCT image, thereby improving its quality. Although these methods show promise in improving LDCT imaging, they rely on the patient's own previous images of the same or adjacent slices, which limits their applicability in clinical practice.

[0007] Another research approach uses NDCT images from a patient population as exemplar images for MBIR in LDCT imaging. This approach uses feature learning techniques to extract exemplar information to capture the structural features of these exemplar images and lay the foundation for MBIR in LDCT imaging. For example, Mou et al. implemented dictionary learning in LDCT reconstruction. They focused on the sparse constraints of dictionary coefficients and used global or adaptive dictionaries to demonstrate the effectiveness of the method and pioneered its application in CT reconstruction. Chen et al. proposed a patch-based dictionary learning method specifically for improving the quality of abdominal tumor images in LDCT. Chen et al. developed a stationary nonlocal concentrated sparse representation (SNCSR) model that combines stationary PCA sub-dictionaries, nonlocal concentrated sparse representations, and relative total variation. Yu et al. proposed a CT reconstruction method using weighted adaptive nonlocal dictionaries. Zheng et al. introduced a CT reconstruction method called PWLS-ULTRA, which is regularized by a union of learned transformations and pre-trained using CT images or image patches of the dataset. Zhang et al. learned tissue-specific Markov random field (MRF) coefficients from previous NDCT scans and used this data as prior knowledge for LDCT reconstruction. Gao et al. studied the feasibility of extracting tissue-specific MRF coefficients from the NDCT database to reconstruct LDCT images. These strategies allow the use of NDCT images of a patient population as exemplar images, thereby expanding their scope of application. However, it usually involves training all exemplar images to obtain a limited set of prior features that may not adequately represent individual variability. This may challenge the adaptive representation of features in specific regions of interest in the target image, potentially limiting the ability to accurately represent texture details. Summary of the invention

[0008] 1. Problem that the invention aims to solve The present invention aims to overcome the deficiencies of the prior art and introduces a low-dose CT image reconstruction method guided by group patient CT images (PI-MBIR). The method addresses the challenges associated with the inherent limitations of raw data in a low-dose scanning environment by utilizing prior information extracted from CT images of different groups of patients. The PI-MBIR method selects similar example samples from the conventional dose CT images of the group of patients for each target region of interest structure of the reconstructed image, and then uses the PCA method to learn the prior subspace adapted to the region of interest structure, reflecting its unique structural and texture characteristics. This targeted approach ensures that the reconstruction process is accurately targeted to the specific settings of each region of interest, thereby improving the preservation of texture details. The method also proposes an adaptive sparse modeling strategy, in particular a weighted Laplace distribution model of the prior subspace. The proposed prior model is fine-tuned to match the signal-to-noise ratio of each transformed band. This adjustment allows stronger constraints to be imposed on bands with low signal-to-noise ratios and lighter constraints to bands with high signal-to-noise ratios, thereby allowing a more accurate characterization of the prior structure of the CT image. Based on the prior model, the objective function of low-dose CT image reconstruction is constructed, and the alternating direction multiplier method framework is used to iteratively optimize the objective function.

[0009] 2. Technical solution To achieve the above objectives, the present invention develops an innovative group patient CT image-guided MBIR method for LDCT imaging, which mainly includes the following steps: Step 1: Construction of a sample image dataset of group patients. First, conventional-dose CT images from different groups of patients are collected. Then, two-dimensional image blocks of size 5×5 (referred to as sample image blocks) are automatically extracted from the conventional-dose CT images using a sliding window method, and the image blocks are merged to construct a sample image dataset.

[0010] Step 2: Screening of similar example samples. This step selects example samples for each target image block of the reconstructed image for subsequent prior subspace learning. Specifically, the Euclidean distance is used to measure the similarity between the target image block and each image block in the example image dataset, and the N image blocks with the highest similarity to the target image block are selected as example samples.

[0011] Step 3: Prior subspace transformation basis learning. For each target image block of the reconstructed image, the principal component transformation basis is learned using the example samples selected in step 2 using the principal component analysis (PCA) method. The prior subspace that is suitable for the target image block to be reconstructed can be constructed using the principal component transformation basis.

[0012] Specifically, for the i target image patch, given its example sample , where each sample is a column vector of length M, N is the number of example samples, and the principal component transformation basis is learned by applying the principal component analysis (PCA) method through the following steps: a) By subtracting the mean The example samples are centralized, where is the mean of the samples. In the present invention, in order to better align with the target image block structural features, we use a weighted average strategy to calculate the sample mean: (1) in, is the weight set according to the similarity of the image blocks: (2) In the above formula is the Euclidean distance between the target image block and its kth example sample, h is a constant parameter set empirically, and Z is a normalization factor.

[0013] b) Calculate the covariance matrix of the example sample after centering: , and decompose the covariance matrix into: (3) , in Contains the eigenvalues, sorted in descending order. It is composed of the corresponding eigenvectors. Construct the required principal component transformation basis.

[0014] Step 4: Adaptive sparse modeling of prior subspace. In low-dose CT imaging tasks, the noise in CT images exhibits non-stationary characteristics. Different image blocks often have different statistical characteristics. At the same time, the transformation subspaces (referred to as transformation band subspaces) constructed by different principal component transformation bases will also show significant differences. The present invention proposes a new weighted Laplace model to perform adaptive sparse modeling of the prior subspace. The model can perform adaptive sparse modeling of transformation coefficients: (4) in, is the basis of the j-th transform band subspace, denoted by . is a small constant added for numerical stability. represents the expectation of the j-th transformed band subspace of the coefficients, is the weight parameter, estimated using the example sample as follows: (5) in for The jth value of .

[0015] In this study, we propose an adaptive setting This method takes into account the energy characteristics of both signal and noise in the PCA transform subspace during the setting process, that is, the signal energy is mainly concentrated in the first few bands of the PCA subspace, while the noise energy is distributed more evenly. In order to assign appropriate weights, we adopted the following strategy: (6) in is the variance of the noise in the local image patch, is a small positive constant. The function represents the statistical expectation, represents the signal energy of the target image block in the jth transformed band in the PCA subspace. Here, by carefully adjusting the weights To match the signal-to-noise ratio (SNR) of each transformed band subspace. Therefore, applying stronger constraints to bands with low SNR and weaker constraints to bands with high SNR helps suppress noise while maintaining fine tissue texture.

[0016] It is worth noting that although the weight By signal The energy of is determined, but the image to be reconstructed itself is not directly usable. Since the transform subspace is derived from clean example samples, the decomposed eigenvalues ​​mainly capture the signal variability of each transform band. According to the PCA principle, we estimate for: (7) in is the j-th eigenvalue of .

[0017] Step 5: Based on the prior subspace adaptive sparse modeling, the low-dose CT image reconstruction objective function is constructed, and the alternating direction multiplier method (ADMM) is used to iteratively optimize the objective function.

[0018] 3. Beneficial effects Compared with the prior art, the technical solution provided by the present invention has the following beneficial effects: (1) This paper proposes a patient-based CT image-guided low-dose CT image reconstruction method (PI-MBIR). PI-MBIR effectively addresses the challenges associated with the inherent limitations of raw data in a low-dose scanning environment by leveraging exemplar information extracted from high-quality images of different patient groups.

[0019] (2) In the present invention, the PI-MBIR method selects similar example samples from the conventional dose CT images of a group of patients for each target ROI structure of the reconstructed image, and then uses the PCA method to learn the prior subspace adapted to the ROI structure, reflecting its unique structural and texture characteristics. This targeted approach ensures that the reconstruction process is precisely targeted to the specific settings of each ROI, thereby improving the preservation of texture details.

[0020] (3) The present invention also proposes an adaptive sparse modeling strategy, specifically a weighted Laplace distribution model of the prior subspace. The proposed prior model is fine-tuned to match the signal-to-noise ratio of each transformed band. This adjustment allows stronger constraints to be imposed on bands with low signal-to-noise ratios and lighter constraints to bands with high signal-to-noise ratios, thereby allowing a more accurate characterization of the prior structure of the CT image.

[0021] (4) In the efficient CT reconstruction algorithm designed by the present invention, an optimization method suitable for the reconstruction algorithm is developed using the alternating direction method of multipliers (ADMM) framework. The algorithm uses an alternating method to solve the joint optimization problem. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 The specific implementation steps of the PI-MBIR method of the present invention are summarized; DETAILED DESCRIPTION

[0023] The present invention is a low-dose CT image reconstruction method guided by CT images of a group of patients. Figure 1 , describe in detail the specific implementation process of the present invention.

[0024] Iterative reconstruction of low-dose CT focuses on reconstructing noisy line-integrated measurements Estimated expected output This task is inherently challenging due to the presence of noise in these measurements. To address this problem, a Maximum A Posteriori (MAP) estimation framework is used to estimate the Reconstruct the desired CT image : (8) in, is a weighting parameter that controls the balance between likelihood and prior probability. Assuming that each line integral measurement of each projection partition is almost independent and can be modeled as a non-stationary Gaussian distribution, the above equation can be reformulated as: (9) In the formula, is the system matrix, where each element Reflects the contribution of voxel j to projection ray i. is a diagonal matrix where each element represents the inverse of the variance of the corresponding measurement in y.

[0025] The second term in equation (9) is the penalty model, which incorporates prior knowledge about the image. Developing a model that effectively captures extensive prior structural knowledge of the target CT image is crucial to improving the performance of the reconstruction framework. In the present invention, the local subspace is constructed as prior knowledge by adaptively learning the principal component transformation basis of the noise-free example images of each ROI in the target image.

[0026] Taking the logarithm of equation (4), the second term of equation (9) is transformed into: (10) In this expression, is considered as an L1 sparse regularization constraint, and Controls the strength of this constraint.

[0027] By integrating equations (9) and (10), the CT maximum a posteriori probability reconstruction problem is formulated as a penalized weighted least squares (PWLS) optimization problem: (11) To solve this optimization function (11), we use the alternating direction method of multipliers (ADMM). We first introduce the auxiliary variable in equation (11) ,make Therefore, the scaled augmented Lagrangian dual form of equation (11) becomes: (12) In this formula, represents the Lagrange multiplier, and γ is a positive penalty scalar. The optimization strategy consists of an alternating approach where we iteratively optimize the variables and This is done by minimizing equation (12) while keeping other variables constant. At the end of each iteration, the Lagrange multiplier The update is performed by gradient ascent. The detailed steps for optimizing Eq. (12) are outlined in the following subsections and Figure 1 Summarized in .

[0028] 1) Sub-questions: Sub-problem: Given and , and the related sub-problems are expressed as: (13) In the formula is the second term of the above formula, and the superscript t is the number of iterations. The scale is large, and it is not feasible to directly solve problem (13). Therefore, in the present invention, the ordered subset accelerated relaxed linearized augmented Lagrangian method (OS-LALM) is used to solve (13), as shown in Figure 1 shown.

[0029] Specifically, The gradient in (13) is: (14) yes The diagonal matrix of the Hessian is defined as: (15) in is the identity matrix. yes The diagonal matrix of is expressed as: (16) in represents an all-ones vector of appropriate dimension. Figure 1 Where Q is the number of iterations, is the number of ordered subsets, (value range is [1,2)) is the relaxation parameter, , , Corresponding to , , The subvector of the mth subset of . As each iteration decreases: (17) 2) Sub-problem: Given and ,about The sub-problem is expressed as: (18) pass Defined as , formula (18) becomes: (19) Formula (19) has a complex structure and cannot be solved directly. For simplicity, we use represents each image block, where .because is from So the strategy is to first Make an estimate and then aggregate using To reconstruct The challenge is to Medium estimate To this end, let Indicates from Extracted from Corresponding image blocks at the same position. According to the law of large numbers, at each iteration, there is a very large probability (limited to 1) to obtain the following equation: (20) I, M and S represent The size of , the image block size, and the total number of pixels in the image block. Since I and S are both large enough for the law of large numbers (a standard 512×512 clinical CT image is usually >260,000), equation (19) can be transformed into: (twenty one) in is a constant parameter. and Respectively and The transformation coefficients of . Using the orthogonality of principal component analysis, we can derive: (twenty two) Therefore, formula (21) can be restated as: (twenty three) This problem can be solved by using the closed soft threshold method: (twenty four) in , , (25) but Calculated as: (26) After obtaining estimates of all overlapping image patches, relocate the image patches to their original locations, average the overlapping pixels, and reconstruct the complete image. : (27) 3) Subproblem: In the ADMM framework, the Lagrange multiplier Update as follows (28)

[0030] The present invention is a low-dose CT image reconstruction method guided by CT images of a group of patients. The method learns a priori subspace adapted to each target area from high-quality example samples. At the same time, for different transformed band subspaces, a weighted Laplace model is proposed to perform adaptive sparse modeling on the priori subspace. The model allows stronger constraints to be imposed on bands with low signal-to-noise ratios and lighter constraints to bands with high signal-to-noise ratios, thereby allowing a more accurate characterization of the priori structure of the CT image.

Claims

1. The present invention discloses a novel group patient CT image guided low-dose CT image reconstruction method (PI-MBIR for short), characterized by the following steps: Step 1: Construction of a sample image dataset of group patients. First, conventional dose CT images from different groups of patients are collected. Then, two-dimensional image blocks of size 5×5 (referred to as sample image blocks) are automatically extracted from the conventional dose CT images using a sliding window method, and the image blocks are merged to construct a sample image dataset; Step 2: Screening of similar example samples. This step selects example samples for each target image block of the reconstructed image for subsequent prior subspace learning. Specifically, the Euclidean distance is used to measure the similarity between the target image block and each image block in the example image dataset, and the N image blocks with the highest similarity to the target image block are selected as example samples; Step 3: Prior subspace transformation basis learning. For each target image block of the reconstructed image, use the example samples selected in step 2 and use the principal component analysis (PCA) method to learn the principal component transformation basis. The principal component transformation basis can be used to construct a prior subspace that is suitable for the target image block to be reconstructed; Step 4: Adaptive sparse modeling of prior subspace. In low-dose CT imaging tasks, the noise in CT images exhibits non-stationary characteristics. Different image blocks often have different statistical characteristics. At the same time, the transformation subspaces (referred to as transformation band subspaces) constructed by different principal component transformation bases will also show significant differences. The present invention proposes a new weighted Laplace model to perform adaptive sparse modeling of the prior subspace; Step 5: Based on the prior subspace adaptive sparse modeling, the low-dose CT image reconstruction objective function is constructed, and the alternating direction multiplier method (ADMM) is used to iteratively optimize the objective function.

2. The method for low-dose CT image reconstruction guided by CT images of a group of patients according to claim 1, characterized in that: In step 2: Step 2.1: By representing each two-dimensional example image patch as a column vector, the i-th example image patch in the dataset is , where M represents the total number of pixels in the image block. It is assumed that the noise in the local neighborhood of the low-dose CT image is additive, independent and identically distributed. Under this assumption, the noisy version of the target image block to be reconstructed can be expressed as: ； in The image patch variable representing the noise, ,in Represents the noise of each pixel in the image block. The main characteristics of noise are that its mean value is zero and its variance is ; To calculate the noisy target image patch And each example image patch in the dataset The Euclidean distance between , the present invention uses the following equation: ； In order to select example samples similar to the target image block from the example image block dataset, the present invention selects image blocks that meet the following conditions as candidate example samples: ； Where T is a given threshold value, and the present invention takes 4×10 3 ; Step 2.2: When processing the selected candidate example samples, two situations are encountered. The first situation occurs when the number of candidate samples N is too small, which may not provide a robust estimate of the local prior subspace and may introduce inaccuracies in the reconstructed CT image. The second situation occurs when N is too large, which significantly increases the computational burden. A large number of experiments have shown that for an image block vector of size M, selecting at least M most similar candidate samples can ensure a more robust estimate of the prior subspace. At the same time, it is observed that the reconstruction performance of more than 3M candidate samples is slightly improved. Therefore, based on these findings, the following criteria are proposed to select the number of training samples in the present invention: (1) If , all candidate samples are used as exemplar samples; (2) If , using 3M noise-similar image patches from low-dose CT images; (3) If , 3M candidate samples are selected as example samples.

3. The method for low-dose CT image reconstruction guided by CT images of a group of patients according to claim 1, characterized in that: In step 3: Step 3.1: For the i-th target image patch, given its example sample , where each sample is a column vector of length M, N is the number of example samples, and the principal component transformation basis is learned by applying the principal component analysis (PCA) method through the following steps: First, by subtracting the mean Centralize the example samples. is the mean of the samples. In the present invention, in order to better align with the target image block structural features, we adopt a weighted average strategy to calculate the sample mean: ； in, is the weight set according to the similarity of the image blocks: ； In the above formula is the Euclidean distance between the target image block and its kth example sample, h is a constant parameter set empirically, and Z is a normalization factor; Step 3.2: Calculate the covariance matrix of the centered example sample: , and decompose the covariance matrix into: ,in Contains the eigenvalues, sorted in descending order. It is composed of the corresponding eigenvectors. Construct the required principal component transformation basis.

4. The method for low-dose CT image reconstruction guided by CT images of a group of patients according to claim 1, characterized in that: In step 4: Step 4.1: The present invention introduces a new weighted Laplace model for modeling for each different transform band subspace. The model can perform adaptive sparse modeling on the transform coefficients, as follows: ； in, is the basis of the j-th transform band subspace, denoted by . is a small constant added for numerical stability. represents the expectation of the j-th transformed band subspace of the coefficients, is the weight parameter, estimated using the example sample as follows: ； in for The jth value of ; Step 4.2: In this study, we propose an adaptive setting This method takes into account the energy characteristics of both signal and noise in the PCA transform subspace during the setting process, that is, the signal energy is mainly concentrated in the first few bands of the PCA subspace, while the noise energy is distributed more evenly. In order to assign appropriate weights, we adopted the following strategy: ； in is the variance of the noise in the local image patch, is a small positive constant. The function represents the statistical expectation, represents the signal energy of the target image block in the jth transformed band in the PCA subspace. Here, by carefully adjusting the weights To match the signal-to-noise ratio (SNR) of each transformed band subspace. Therefore, applying stronger constraints to bands with low SNR and weaker constraints to bands with high SNR helps suppress noise while maintaining good tissue texture; Step 4.3: It is worth noting that although the weight By signal The energy is determined, but the image to be reconstructed itself is not directly usable. Since the transform subspace is derived from clean example samples, the decomposed eigenvalues ​​mainly capture the signal variability of each transform band. According to the PCA principle, we estimate for: ； in for The j-th eigenvalue of .

5. The method for low-dose CT image reconstruction guided by CT images of a group of patients according to claim 1, characterized in that: In step 5: Step 5.1: Based on the prior subspace adaptive sparse modeling, construct the low-dose CT image reconstruction objective function: ； In the formula is the line integral measurement, is the reconstructed CT image, is the system matrix, where each element Reflects the contribution of voxel j to projection ray i. is a diagonal matrix where each element represents the inverse of the variance of the corresponding measurement in y; Step 5.2: Use the alternating direction multiplier method (ADMM) to iteratively optimize the objective function. First, introduce the auxiliary variable ,make . Therefore, the scaled augmented Lagrangian dual form of the total function becomes: ； In this formula, represents the Lagrange multiplier, and γ is a positive penalty scalar. The optimization strategy used is an alternating method that iteratively optimizes the variables and This is done by minimizing equation (10) while keeping other variables constant. At the end of each iteration, the Lagrange multiplier Update via gradient ascent; Step 5.2: Sub-problem: Given and ,and The relevant sub-problems are expressed as: ； In the formula is the second term of the above formula, and the superscript t is the number of iterations. Regulation The modulus is large, and it is not feasible to directly solve problem (11). Therefore, in our invention, the ordered subset accelerated relaxed linearized augmented Lagrangian method (OS-LALM) is used to solve (11); Step 5.3: Sub-problem: Given and ,about The sub-problem is expressed as: ； Step 5.4: Subproblem: In the ADMM framework, the Lagrange multiplier Updated as follows: 。

Images