A limited-angle CT image reconstruction algorithm based on directional rebound total variation

Abstract

The application discloses a limited-angle CT image reconstruction algorithm based on directional rebound total variation, and belongs to the technical field of computer tomography imaging. In view of the problems of poor directional modeling adaptability, excessive regularization smoothing, and high calculation cost of traditional ADMM solution of existing DTV and STV reconstruction algorithms, the application provides a limited-angle CT image reconstruction method based on directional rebound total variation regularization, fuses a mixed regularization framework constructed by directional gradient prior and weakly convex rebound optimization, matches a DCA convex difference decomposition and a FL-ADMM complete linearization fast solution architecture, suppresses strip-like artifacts under the condition of extremely limited-angle sampling, completely retains fine textures and weak boundaries, greatly reduces iteration operation amount, and improves reconstruction speed.

Description

Technical Field

[0001] This invention belongs to the field of computed tomography imaging technology, specifically involving a high-precision reconstruction algorithm for CT tomographic images under limited angle projection data. It is applicable to scenarios such as linear trajectory CT, breast CT, short exposure low-dose clinical CT, and limited scanning industrial defect detection CT, and can be integrated into CT post-processing workstations, industrial imaging software, and embedded reconstruction processing terminals. Background Technology

[0002] Computed tomography (CT) is a mainstream non-destructive imaging technique widely used in industrial defect detection and clinical medical diagnosis. However, in real-world scanning scenarios, limitations such as the size of the object being measured, the restricted scanning trajectory of linear CT scans, the need to shorten exposure time to suppress motion artifacts, and the limited acquisition range in breast scans restrict the acquisition of projection data to a finite angle of less than 180°. Traditional analytical reconstruction algorithms, such as filtered back projection (FBP), suffer from insufficient imaging information due to missing projection data. This results in reconstructed images with numerous stripe artifacts and severely blurred edges, failing to meet the accuracy requirements for diagnosis and detection.

[0003] With the widespread application of compressed sensing theory in CT reconstruction, iterative regularization algorithms based on optimization have become the mainstream solution for finite-angle reconstruction problems. Iterative regularization algorithms based on total variation (TV) utilize the sparse prior of image gradients to suppress strip artifacts and preserve edge structures even with incomplete projection data, thus gaining widespread use. However, existing improved TV-type algorithms have significant drawbacks: 1.1. Directed Total Variation (DTV) algorithm applies gradient sparsity constraints along the horizontal and vertical directions respectively to fully exploit the prior of the edge direction of medical images and suppress directional artifacts. However, it only uses ℓ1 norm regularization, which inevitably leads to problems such as over-smoothing, loss of texture details, and blurring of weak boundaries. 2.2. Springback Total Variation (STV) model introduces a weak convex springback penalty term and combines the advantages of ℓ1 norm stability and ℓ2 norm detail preservation to alleviate the defects of ℓ1 regularization over-smoothing. However, it adopts an isotropic gradient design and has no adaptive perception capability for multi-directional anatomical edges. In finite angle scenarios with severe angle loss, it has insufficient modeling capability for complex tissue structures and limited reconstruction accuracy.

[0004] Meanwhile, the traditional Alternating Direction Multiplier Method (ADMM) requires inverting large system matrices when solving iterative optimization models, resulting in high computational overhead and long reconstruction time, which is difficult to meet the needs of real-time clinical imaging and rapid industrial detection. In summary, existing technologies cannot simultaneously achieve multi-directional edge adaptive modeling, weak convexity regularization detail preservation, and fast solution with low computational cost, and both the effect and efficiency of finite-angle CT reconstruction have obvious shortcomings. Summary of the Invention

[0005] To address the problems of poor adaptability of directional modeling, excessive smoothing of regularization, and high computational cost of traditional ADMM in existing DTV and STV reconstruction algorithms, this invention provides a finite-angle CT image reconstruction algorithm based on directional rebound total variation. It integrates directional gradient prior and weak convex rebound optimization to construct a hybrid regularization framework, and combines DCA convexity decomposition and FL-ADMM fully linearized fast solution architecture. Under extremely limited angle sampling conditions, it suppresses strip artifacts, fully preserves fine textures and weak boundaries, and significantly reduces the amount of iterative computation, thereby improving reconstruction speed.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: A finite-angle CT image reconstruction algorithm based on directional rebound total variation includes the following steps: Step 1: Obtain CT projection data g under finite angle scanning conditions and establish the discrete projection linear imaging equation: g = Au ; in, g For discrete projection data, A For the system matrix, u The vector of the CT tomographic image to be reconstructed; Step 2: Construct a directional rebound total variation regularized objective optimization model. This model is based on the directional total variation gradient framework. The weak convex rebound correction term of the rebound total variation is embedded into the gradient constraints in the horizontal and vertical directions, respectively, and data fidelity constraints are set. Step 3: Solve the target optimization model to obtain the reconstructed CT tomographic image.

[0007] Furthermore, the directional rebound total variation regularized objective optimization model described in step 2 is expressed as: , ; In the formula: D x This is the discrete gradient transformation matrix in the horizontal direction of the image; D y This is the discrete gradient transformation matrix in the vertical direction of the image; α 1 represents the horizontal gradient regularization weight. α 2 represents the vertical gradient regularization weight; α x , α y It is a bidirectional weakly convex balance coefficient used for regulation. Norm sparsity constraints and Weights for norm detail preservation; εLet || be the tolerance threshold for the projection residual. Au - g ||2≤ ε for The norm sphere data fidelity constraint limits the range of residuals between the reconstructed projection and the measured projection.

[0008] Furthermore, in step 3, the convexity algorithm is used to decompose the target optimization model, and the non-convex terms in the objective function are linearly approximated, transforming the original non-convex minimization problem into a convex sub-optimization problem with multiple rounds of outer layer iterations.

[0009] Furthermore, in step 3, the convexity difference algorithm is used to linearly approximate the second-order ℓ2 term in the objective function.

[0010] Furthermore, for each round of convexity outer layer iteration generating convex sub-optimization problem, the fully linearized alternating direction multiplier method is used for solution, including: introducing auxiliary variable splitting constraints, constructing an augmented Lagrangian function, and sequentially completing the closed-loop iterative update of each variable; performing a first-order Taylor linear expansion on the image variable sub-problem to avoid the system matrix inversion operation.

[0011] Furthermore, the auxiliary variables include: gradient auxiliary variables. , and residual auxiliary variables .

[0012] Furthermore, the variables include image variables. u Gradient auxiliary variables, residual auxiliary variables, and Lagrange multipliers.

[0013] Furthermore, the maximum number of outer and inner iterations is set, and the normalized target error, normalized data error, and normalized total variation error are calculated in real time. The iteration is terminated when all three indicators are lower than the preset convergence threshold or the maximum number of iterations is reached, and the reconstructed image is output.

[0014] Furthermore, the system matrix is ​​a two-dimensional Radon transform system matrix constructed based on the pixel-driven projection method.

[0015] Compared with the prior art, the present invention has the following advantages: This invention pioneers a hybrid regularization framework for DSTV-oriented weak convexity, combining the advantages of DTV bidirectional gradient orientation modeling with the weak convexity and bounce penalty characteristics of STV, while simultaneously solving the problems of DTV algorithms. It overcomes the dual shortcomings of over-smoothing regularization and the inability of STV isotropic gradients to adapt to multi-directional edges; in extremely limited angle sampling scenarios of 19°~50°, it significantly suppresses strip artifacts, fully preserves weak boundaries and subtle anatomical textures, and the reconstruction PSNR and SSIM metrics are significantly better than the existing mainstream DTV and STV algorithms.

[0016] The system employs a dual-layer solution architecture of DCA and FL-ADMM. DCA decomposes the non-convex optimization model to ensure stable convergence of the iteration. FL-ADMM fully linearizes the quadratic terms, eliminating the need for matrix inversion in large systems, significantly reducing the computational overhead of iteration, and greatly improving the reconstruction speed to meet the needs of real-time clinical imaging and rapid industrial detection.

[0017] The algorithm is robust and has been validated through multiple experiments using Shepp-Logan phantoms, Forbild fine phantoms, and real clinical brain CT scans. It can achieve high-precision reconstruction in noise-free, limited-angle projection scenarios, and its parameters are adjustable to adapt to different imaging devices and different scanning angle scenarios.

[0018] It is highly versatile and can be directly compiled and integrated into medical CT post-processing workstations, linear track industrial CT non-destructive testing software, and breast low-dose CT imaging systems, adapting to various practical imaging conditions such as limited angle scanning, short exposure, and low dose. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 A comparison image of the Shepp-Logan phantom from all angles of reverse crime reconstruction; Figure 2 The convergence curves of NOE, NDE, and NTVE for the Shepp-Logan motif iteration; Figure 3 A comparison image of the Forbild phantom from all angles of reverse crime reconstruction; Figure 4 Visual comparison images of DTV, STV, and DSTV reconstructed from different finite angles (19°, 20°, 30°, 40°, 50°) of the Shepp-Logan phantom; Figure 5 Quantitative evaluation curves (NOE, PSNR, SSIM) were reconstructed for the Shepp-Logan phantom at different finite angles. Figure 6 Visual comparison images of DTV, STV, and DSTV reconstructed from the Forbild phantom at different finite angles (50°, 55°, 60°, 65°, 70°); Figure 7 Quantitative evaluation curves (NOE, PRSE, SSIM) were reconstructed for the Forbild phantom at different finite angles. Figure 8 Comparison of DTV, STV, and DSTV reconstruction effects under different limited angles (50°, 60°, 70°, 80°, 90°) in clinical brain CT scans; Figure 9 Quantitative evaluation curves for reconstruction under different limited angles in clinical brain CT scans; Figure 10 This is a magnified comparison of details from a clinical brain CT scan at a limited angle of 50°. Detailed Implementation

[0021] To gain a deeper understanding of this invention, we will provide a comprehensive and detailed description. However, this invention has various implementations and is not limited to the specific examples listed herein. These examples are presented to enhance a full understanding of the disclosure of this invention.

[0022] We will first introduce the applicable scenarios and usage methods of the DCA algorithm, soft thresholding algorithm, projection algorithm, and FL-ADMM algorithm.

[0023] The minimization problem is formulated as a convexity optimization problem: (1); The Convex Function Difference (DCA) algorithm handles the non-convex structure of the objective function through sequential convex approximation. By... By linearizing the solution, an iterative solution scheme is derived: (2).

[0024] Soft thresholding algorithms are widely used in compressed sensing, image reconstruction, and other fields. Soft thresholding (also known as shrinkage operators) is one of the most fundamental core operators for solving incomplete data optimization problems. It primarily reduces small-amplitude gradient coefficients to zero while preserving large-amplitude edge-related gradients, thus applying gradient sparsity to CT image regularization. One-dimensional shrinkage algorithm: (3), x It is a vector of size N. It is the Σ1 norm of a vector. It is the Σ2 norm of a vector.

[0025] Assumption This represents the first vector. i There are elements. The solution to this optimization problem is: , (4), S It is a one-dimensional contraction operator. It is a standard symbolic function; its value of 1 represents a positive number, its value of 0 represents zero, and its value of -1 represents a negative number. Operators Choose the maximum value between its two parameters.

[0026] Projection on convex set and The situation. The Projection of Convex Sets (POCS) operator is a fundamental component of convex optimization, implementing convex constraints by mapping a point to its nearest feasible point. Formally, given a convex set C, a vector... The projection onto C is defined as the solution to the following constrained least squares problem: , (5) The projection operator obtained from this is denoted as This directly yields the optimal point: (6).

[0027] Projected to On a sphere, if a convex set has a radius of... of Ball, denoted as Then the projection has a simple closed-form solution: (7). If the input vector Already inside the ball If , then the projection will keep it unchanged. If it is located outside the sphere, the projection will shrink it to the upper edge of the sphere's boundary. A point indicating direction.

[0028] The FL-ADMM algorithm is a general, simple, efficient, and accurate method for solving convex optimization models for image reconstruction, regardless of whether they are unconstrained or constrained. The FL-ADMM algorithm is suitable for solving optimization problems of the following forms: , (8), of which x and y It is a multidimensional vector, where A and B are matrices representing linear transformations. and It is a convex function, but not necessarily a smooth function.

[0029] The corresponding augmented Lagrangian function is: The iterative optimization framework is as follows: ; The ADMM algorithm divides the original problem into three simple subproblems. The formulas for these subproblems (10a) and (10b) are similar, so we will only discuss the linearization of formula (10a). If Based on Taylor expansion, any quadratic function can be linearized at a given point. Therefore, the linearized form of formula (10a) is: ;in, , The closed-form solution can be expressed as: (11).

[0030] This invention provides a finite-angle CT image reconstruction algorithm based on directional rebound total variation: It employs a constrained total variation (TV) model, with the regularization term based on a bounce-penalty model, and incorporates the directionality of the directional total variation (DTV) algorithm. This is a combination of... and Based on the norm regularization method, the direction of the model is constrained to promote sparse solutions and maintain the stability of the solutions, thereby enabling the effective recovery of unknown signals under incomplete observation conditions.

[0031] The DSTV model is described in detail below: , (12); among which, u Represents the image to be recovered. , These are TV weight parameters that control the horizontal and vertical directions respectively. , These are the control norms in the horizontal and vertical directions, respectively. sum norm The positive parameter of the weight, To measure the goodness of fit between the observed data and the reconstructed projection, ε This represents the data tolerance describing the level of projection noise. (Matrix) and The sizes are N×N, representing the distances along the path. x shaft and y Discrete gradient transformation of the axes. The mathematical definition is as follows: , ; DSTV Algorithm: The optimization problem corresponding to formula (12) can be reformulated as follows: (15) here, express The indicator function of the ball is used to reconstruct the data. Au With measurement data g The differences between them are limited to a specified tolerance range. This tolerance is defined as follows: (16);

[0032] According to the general framework of DCA In formula (15), respectively level The second item of direction and the second term in the vertical direction Replace with line type item and The iterative scheme of DCA is expressed as: (17) The FL-ADMM algorithm is used to solve the derived DCA subproblem by introducing auxiliary variables. , , ; To avoid confusion, the external iteration variable... u k Denote as constant u a It remains fixed during the internal iteration process. Equation (17) can thus be expressed as: , (18) The corresponding augmented Lagrangian function is: (19) Here, β 1, β 2, β 3 > 0 is the penalty parameter. λ 1, λ 2, λ 3 represents the Lagrange multipliers used for enforcing constraints. These multipliers are updated in each iteration. The FL-ADMM framework is represented as:

[0033] Update subproblem: According to formula (6), by re-pointing u p By performing a Taylor expansion, subproblem (20.a) can be linearized into a strongly convex quadratic optimization problem, making it easier to solve. Therefore, we can obtain:

[0034] Assuming the gradient of the objective function in formula (19) is 0, it can be directly derived that... u Iterative update formula:

[0035] here, , , .

[0036] Update subproblem , ; The subproblem in formula (20.b) is solved in the same way as that in formula (20.c). In this embodiment, we take the solution of the subproblem in formula (20.b) as an example, which can be simplified as follows:

[0037] make , Formula (23) can be expressed as:

[0038] The solution can be obtained using the one-bit soft threshold method. :

[0039] here, S It is a one-dimensional contraction operator. It is a standard symbolic function, with a value of 1 representing a positive number, 0 representing 0, and -1 representing a negative number. It is an operator that selects the maximum value.

[0040] Similarly, let , ,get : .

[0041] Update subproblem z: The subproblem in equation (20.d) can be transformed into the following form: ,

[0042] That is, to Projected onto a radius of of On the norm sphere.

[0043]

[0044] in, P This represents the projection operator, which means Already located If it is within the norm sphere, then it is projected onto the boundary of the sphere by scaling.

[0045] The method of the present invention will now be tested and verified.

[0046] Example 1: Finite-Angle Reconstruction Experiment of Shepp-Logan Phantom In this embodiment, the Shepp-Logan digital phantom is used as the reconstruction object, and the phantom resolution is set to 256×256 pixels. The acquisition of full-angle projection data follows a standardized process: a detector with a length of 256 is used to uniformly collect projection information at 256 angles within the range of 0~π radians (0~180°), and the system matrix A is defined through a precise pixel-driven modeling method.

[0047] Limited-angle experiments were conducted within limited angle ranges of 19°, 20°, 30°, 40°, and 50°, and the DSTV algorithm was compared with the STV and DTV algorithms. The angle sampling interval for each reconstruction was approximately 0.7° (180° / 256).

[0048] Taking the reconstruction of a limited angular range of 20° as an example, the imaging conditions are as follows: image size is 256×256, each pixel size is 1×1; projection size is 256×29, that is, there are 29 projections, each projection has 256 measurements; detector partition size is 1. Parameter settings are: α 1 = 1.5, α 2 = 0.2, α x =1.5, α y =0.2, β 1 = β 2 = 10, β 3 = 0.01. The reconstruction uses noise-free projection data acquired within a limited angular range, ε = 0. Maximum number of external iterations. MaxOit =3, maximum number of internal iterations MaxIt =5000.

[0049] Experimental results show that, within a smaller angular range, the reconstruction accuracy of the DSTV algorithm is higher than that of the STV and DTV algorithms (e.g., ...). Figure 4 (As shown). Within the angle range θ = 19°, the DSTV algorithm produces only minor artifacts, which do not affect the observation of the Shepp-Logan phantom, while DTV and STV both exhibit significant artifacts. Within the angle range θ = 20°, DSTV successfully reconstructs an image visually indistinguishable from the original, while DTV and STV only achieve a relatively complete reconstruction of the Shepp-Logan image structure when the angle range increases to 40°.

[0050] Quantitative comparisons of normalized target error (NOE), peak signal-to-noise ratio (PSNR), and structural similarity index (SSIM) further validate the above conclusions (e.g. Figure 5 As shown): In the angular range θ = 20°, only the NOE value of DSTV is above 10. -3The magnitude is significant, and DTV and STV only fall below this critical value when the angle range θ = 40°. As the angle range increases, the evaluation metrics of all algorithms improve, but the NOE, PSNR, and SSIM values ​​of the DSTV algorithm remain the best at the same angle.

[0051] In the anti-crime verification experiment, Figure 1 (a) and Figure 1 (b) The intuitive comparison shows that the reconstructed image has a very high degree of structural consistency with the original phantom; Figure 1 (c) shows that the profile curves of the centerline of the original image and the reconstructed image almost perfectly overlap. Figure 2 The iterative convergence curves show that the three core metrics, NOE, NDE, and NTVE, rapidly decreased to 10 after approximately 4000 iterations. -4 It is below the order of magnitude and maintains a stable convergence state in subsequent iterations.

[0052] Example 2: Finite-angle reconstruction experiment of Forbild fine phantom This embodiment uses the Forbild fine phantom as the reconstruction object. Because this phantom contains a series of fine structures on the right side, it better meets the needs of reproducing minute details in actual clinical imaging. The phantom resolution is set to 256×256 pixels.

[0053] The Forbild phantom was reconstructed using DSTV, STV, and DTV models at different angle ranges (e.g. Figure 6 (As shown). Even within the angular range θ = 50°, DSTV retains fine structural details and generates visually reliable reconstruction results. In contrast, DTV and STV exhibit significant artifacts and distortion. Quantitative comparisons of NOE, PRSE, and SSIM metrics all show that the DSTV algorithm achieves the best results (e.g., Figure 7 (As shown).

[0054] In the anti-crime verification experiment, Figure 3 (a) and Figure 3 (b) The intuitive comparison shows that the reconstructed image has a very high structural consistency with the original phantom. Even when the local area is magnified for observation, no obvious loss of detail or artifact interference is found. Figure 3 (c) shows that the fine structure on the right side of the phantom is completely preserved in the reconstructed image, and the overlap of the cross-sectional curves is close to 100%. Figure 4 The iterative convergence curves show that the three metrics, NOE, NTVE, and NDE, decrease rapidly with the number of iterations, all dropping to 10 after approximately 4000 iterations. -5 The following is an example of convergence, which is achieved subsequently.

[0055] Example 3: Limited-angle reconstruction experiment of clinical brain CT images This embodiment uses real clinical brain CT image data for experiments to evaluate the effectiveness and robustness of the algorithm under real-world conditions.

[0056] Reconstruction of clinical brain CT images using DSTV, DTV, and STV models at different angular ranges (e.g.) Figure 8 (As shown). The maximum number of iterations is set to 20,000, and the parameters are set as follows: α 1 = 1.5, α 2 = 1.5, α x =1.5, α y =0.4, β 1 = β 2 = 1, β 3 = 0.5.

[0057] The results show that DSTV outperforms other algorithms (such as...) across all evaluation metrics. Figure 9 As shown). Within the angular range θ = 50°, a detailed comparison of the fine structural region is performed (e.g., Figure 10 As shown in the image, DSTV preserves edges more clearly and reproduces subtle features more faithfully. In contrast, STV and DTV exhibit significant blurring when handling fine structures and boundaries.

[0058] Contents not described in detail in this specification are prior art known to those skilled in the art. Although illustrative specific embodiments of the invention have been described above to facilitate understanding by those skilled in the art, it should be understood that the invention is not limited to the scope of the specific embodiments. Various modifications are readily apparent to those skilled in the art as long as they fall within the spirit and scope of the invention as defined and determined by the appended claims, and all inventions utilizing the concept of this invention are protected.

Claims

1. A finite-angle CT image reconstruction algorithm based on directional rebound total variation, characterized in that, Includes the following steps: Step 1: Obtain CT projection data under finite angle scanning conditions and establish the discrete projection linear imaging equation: g = Au; Where g is discrete projection data, A is the system matrix, and u is the CT tomographic image vector to be reconstructed; Step 2: Construct a directional rebound total variation regularized objective optimization model. This model is based on the directional total variation gradient framework. The weak convex rebound correction term of the rebound total variation is embedded into the gradient constraints in the horizontal and vertical directions, respectively, and data fidelity constraints are set. Step 3: Solve the target optimization model to obtain the reconstructed CT tomographic image.

2. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 1, characterized in that, The directional rebound total variation regularized objective optimization model described in step 2 is expressed as follows: , In the formula: This is the discrete gradient transformation matrix in the horizontal direction of the image; This is the discrete gradient transformation matrix in the vertical direction of the image; For horizontal gradient regularization weights, These are the vertical gradient regularization weights; , It is a bidirectional weakly convex balance coefficient used for regulation. Norm sparsity constraints and Weights for norm detail preservation; This is the tolerance threshold for the projected residual. for The norm sphere data fidelity constraint limits the range of residuals between the reconstructed projection and the measured projection.

3. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 1, characterized in that, In step 3, the convexity algorithm is used to decompose the target optimization model, and the non-convex terms in the objective function are linearly approximated, transforming the original non-convex minimization problem into a convex sub-optimization problem with multiple rounds of outer layer iterations.

4. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 3, characterized in that, In step 3, the convexity algorithm is used to process the second-order convexity of the objective function. The terms are approximated linearly.

5. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 3, characterized in that, For each round of outer convexity iteration generating the convex sub-optimization problem, the fully linearized alternating direction multiplier method is used for solution, including: introducing auxiliary variable splitting constraints, constructing an augmented Lagrangian function, and sequentially completing the closed-loop iterative update of each variable; performing a first-order Taylor linear expansion on the image variable sub-problem to avoid the system matrix inversion operation.

6. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 5, characterized in that, The auxiliary variables include: gradient auxiliary variables. , and residual auxiliary variables .

7. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 5, characterized in that, The variables include image variables, gradient auxiliary variables, residual auxiliary variables, and Lagrange multipliers.

8. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 3, characterized in that, Set the maximum number of outer and inner iterations, calculate the normalized target error, normalized data error, and normalized total variation error in real time, and terminate the iteration when all three indicators are below the preset convergence threshold or the maximum number of iterations is reached, and output the reconstructed image.

9. The finite-angle CT image reconstruction algorithm based on directional rebound total variation according to claim 1, characterized in that, The system matrix is ​​a two-dimensional Radon transform system matrix constructed based on the pixel-driven projection method.

Images