Fast limited-angle fan-beam ct image reconstruction method based on preconditioning matrix

Abstract

The application relates to a precondition matrix-based fast limited-angle fan-beam CT image reconstruction method and belongs to the technical field of CT image reconstruction. The method comprises the following steps: S1, acquiring projection data; S2, considering the sparsity of an image under tight wavelet frame transformation and image gradient transformation, establishing a double-regularization limited-angle fan-beam CT image reconstruction model to inhibit noise and artifacts and protect a boundary; S3, according to an analytical reconstruction algorithm, designing a precondition matrix, introducing the precondition matrix into the double-regularization limited-angle fan-beam CT image reconstruction model established in the step S2, and solving the model by using an alternating direction multiplier method; and S4, outputting a reconstructed image. The precondition matrix technology is used to integrate the analytical reconstruction algorithm into the regularization iterative reconstruction algorithm, fast convergence of the algorithm is realized, artifacts and noise of the limited-angle fan-beam CT reconstructed image can be effectively inhibited, and therefore the quality of the CT reconstructed image is improved.

Description

Technical Field

[0001] This invention belongs to the field of CT image reconstruction technology, and relates to a fast finite-angle fan-beam CT image reconstruction method based on a precondition matrix. Background Technology

[0002] Limited-angle fan-beam CT, constrained by the size of the scanning scene or the object being examined, can only acquire data within a limited angular range. This results in images reconstructed by traditional analytical reconstruction algorithms containing numerous artifacts, leading to the loss of much important structural information. Methods based on total variation (TV) minimization can suppress artifacts in reconstructed images to some extent, but can also cause oversmoothing. Yu Wei introduced an image gradient L0 minimization method to suppress artifacts and protect boundaries, further improving the quality of the reconstructed image. Wang Chengxiang introduced a compact wavelet transform domain L0 minimization method to avoid the oversmoothing caused by total variation, further suppressing artifacts and noise, thereby improving the quality of the reconstructed image. All of the above methods are based on optimization iterative reconstruction algorithms, requiring significant iteration and reconstruction time to reconstruct a good-quality image, limiting the widespread adoption and promotion of this algorithm in practical commercial CT. Therefore, for limited-angle fan-beam CT, how to quickly reconstruct high-quality images is of great practical significance.

[0003] Patent application CN107978005A discloses "A Finite-Angle CT Image Reconstruction Algorithm Based on Boundary-Preserving Diffusion and Smoothing". This method mainly utilizes the L0 norm of the image gradient to perform boundary-preserving diffusion correction on the horizontal and hammer-shaped aspects of the image. This method can suppress artifacts to a certain extent, but it only introduces a single prior knowledge constraint and requires a large number of iterations and iteration time.

[0004] Patent application CN110717959A discloses a "Method and Apparatus for Reconstructing X-ray Finite-Angle CT Images Based on Curvature Constraints." This method mainly combines image gradient L0 regularization and curvature constraints to suppress artifacts and noise in the reconstructed image. Although this method considers two constraints and can further improve the quality of the reconstructed image, the computational load increases accordingly with the increase of constraints, meaning that a large number of iterations and iteration time are required.

[0005] Patent application CN109697691A discloses a "Finite-Angle Projection Reconstruction Method Based on L0 Norm and Singular Value Threshold Decomposition with Dual Regularization Terms." This method combines image gradient L0 regularization and kernel norm regularization to suppress artifacts and noise in the reconstructed image. Although this method considers two constraints, it can further improve the quality of the reconstructed image. However, with the increase of constraints, a large number of iterations and iteration time are required to ensure a high-quality reconstructed image.

[0006] Patent application CN105590332A discloses "A Fast Algebraic Reconstruction Method for CT Imaging". This method, based on traditional algebraic reconstruction, arbitrarily selects a portion of the hyperplane and further accelerates and adjusts the projection solution vectors obtained by the traditional algebraic reconstruction method to achieve optimality. This method can speed up convergence to some extent, but the iteration time is still relatively long.

[0007] Patent CN105608719A discloses "A Fast CT Image Reconstruction Method Based on Two-Stage Projection Adjustment". This method transforms a non-uniform system of equations into two uniform systems of equations for solution, and utilizes a fast iterative algorithm based on square root error and projection adjustment to accelerate the convergence speed of the algebraic iterative reconstruction algorithm. However, due to the lack of a regularization term, this method results in a large number of artifacts in finite-angle CT reconstructed images.

[0008] This shows that most existing finite-angle fan-beam CT reconstruction methods do not take advantage of the speed of analytical reconstruction algorithms, but simply try to accelerate convergence from the perspective of optimization algorithms. As a result, it is actually difficult to quickly reconstruct high-quality images and thus improve the quality of CT reconstructed images.

[0009] Therefore, a new finite-angle fan-beam CT image reconstruction method is urgently needed to solve the above problems. Summary of the Invention

[0010] In view of this, the purpose of this invention is to provide a fast finite-angle fan-beam CT image reconstruction method based on a precondition matrix. By integrating the analytical reconstruction algorithm into the regularized iterative reconstruction algorithm through the precondition matrix technique, the algorithm can achieve rapid convergence and effectively suppress artifacts and noise in finite-angle fan-beam CT reconstructed images, thereby improving the quality of CT reconstructed images.

[0011] To achieve the above objectives, the present invention provides the following technical solution:

[0012] A fast finite-angle fan-beam CT image reconstruction method based on a preconditioning matrix specifically includes the following steps:

[0013] S1: Obtain projection data;

[0014] S2: Considering the sparsity of images under compact wavelet frame transform and image gradient transform, a double-regularized finite-angle fan-beam CT image reconstruction model is established to suppress noise and artifacts and protect boundaries.

[0015] S3: Based on the analytical reconstruction algorithm, design a precondition matrix and introduce it into the double-regularized finite-angle fan-beam CT image reconstruction model established in step S2. Solve the model using the alternating direction multiplier method.

[0016] S4: Output the reconstructed image.

[0017] Furthermore, in step S2, the established dual-regularized finite-angle fan-beam CT image reconstruction model is as follows:

[0018]

[0019] Where A is the finite-angle CT projection operator, f∈R N×1 The image to be reconstructed is N, which represents the number of reconstructed image sizes. δ ∈R M×1 This is finite-angle CT projection data, where M represents the total number of rays; λ i and Here, is the non-negative regularization parameter, i represents the number of wavelet subbands; W is the compact wavelet frame transform; ||·||0 is a function of the number of non-zero elements in the statistical vector. in The component form is f i′,j′ This represents the (i′, j′)th pixel of the image.

[0020] Furthermore, in step S3, a precondition matrix D is designed and introduced into the double-regularized finite-angle fan-beam CT image reconstruction model established in step S2, that is, model (1) is transformed into the following form:

[0021]

[0022] Where h is an auxiliary variable.

[0023] Furthermore, in step S3, the designed precondition matrix D is:

[0024] D T D=W ham PJ (3)

[0025] Among them, W ham Let P represent the Hamming window function, P represent the filter matrix, and J represent the weighted diagonal matrix.

[0026] Furthermore, in step S3, the alternating direction multiplier method is used to solve the model, specifically including:

[0027] First, model (2) is transformed into an unconstrained optimization problem using the Lagrange augmented function:

[0028]

[0029] Where V represents the Lagrange multiplier, and t > 0 is a parameter introduced by the ADMM algorithm;

[0030] Then, the alternating direction multiplier method is used to solve the problem, and the iterative formula is as follows:

[0031]

[0032] Where k represents the number of iterations;

[0033] To avoid the difficulty of finding the inverse of the system matrix A and the coupling function in the subproblem f of the iterative formula (5), this invention adopts the adjacent alternating linearization method to transform the iterative formula (5) with respect to the subproblems f and h into the adjacent alternating linearization ADMM iterative formula:

[0034]

[0035] in, Let represent the intermediate result of the k-th iteration, and γ,q represent the parameters introduced by the proximity operator;

[0036] Find the optimal solution to the subproblem of iterative formula (6), as shown in the iterative formula below:

[0037]

[0038] Among them, A T This refers to the fan-beam CT backprojection operator and the hard thresholding operator. x represents the threshold. i Let α represent any vector. k Represents the wavelet coefficients after hard thresholding in the k-th iteration;

[0039] In formula (7), V is replaced by a variable. k ←D T V k Then the iterative formula for scaling is:

[0040]

[0041] Formula (8) regarding f k+1 The question is:

[0042]

[0043] The iterative formula for solving equation (9) using the L0 minimization method is as follows:

[0044]

[0045] in, Indicates to The result of the L0 algorithm has the following iterative form: for all image coordinates i′, j′,

[0046]

[0047] Where F represents the Fourier transform, F -1 F(·) represents the inverse Fourier transform. * Represents the complex conjugate of the Fourier transform. Let x and y represent the gradient operators respectively; β represents the control... The similarity parameter, κ (κ>1) represents the parameter controlling the growth rate of β; n represents The number of iterations of the algorithm, The algorithm stops when β is greater than the preset parameter β before iteration. max ,when The algorithm outputs image f after it stops. k+1 .

[0048] The beneficial effects of this invention are as follows: This invention considers the sparsity of images under compact wavelet frame transform and image gradient transform, and suppresses noise and artifacts, as well as protects boundaries, by establishing a double-regularization optimization model. To improve the speed of the reconstruction algorithm, this invention utilizes preconditioning matrix technology, introducing analytical reconstruction algorithms into the iterative reconstruction process. Through this preconditioning matrix technology and double regularization technology, the speed of the reconstruction algorithm is improved while suppressing finite-angle artifacts and noise, significantly improving the quality of CT reconstructed images and promoting the application of iterative reconstruction algorithms in commercial CT.

[0049] Other advantages, objectives, and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination, or may be learned from practice of the invention. The objectives and other advantages of the invention can be realized and obtained through the following description. Attached Figure Description

[0050] To make the objectives, technical solutions, and advantages of the present invention clearer, the preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings, wherein:

[0051] Figure 1 This is a schematic diagram of the geometry of a finite-angle fan-beam CT scanner.

[0052] Figure 2 This is a flowchart of a fast finite-angle fan-beam CT image reconstruction method based on preconditioning matrix technology;

[0053] Figure 3 This is a comparison chart of reconstruction results for a scanning angle range of [0, 140°]. Figure 3 (a) shows the results after 600 iterations without using the preconditioning matrix technique. Figure 3 (b) shows the results of 100 iterations using the precondition matrix technique. Detailed Implementation

[0054] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of the present invention. Unless otherwise specified, the following embodiments and features can be combined with each other.

[0055] The accompanying drawings are for illustrative purposes only and are schematic diagrams, not actual pictures. They should not be construed as limiting the invention. To better illustrate the embodiments of the invention, some parts in the drawings may be omitted, enlarged, or reduced, and do not represent the actual product dimensions. It is understandable to those skilled in the art that some well-known structures and their descriptions may be omitted in the drawings.

[0056] In the accompanying drawings of the embodiments of the present invention, the same or similar reference numerals correspond to the same or similar components. In the description of the present invention, it should be understood that if terms such as "upper," "lower," "left," "right," "front," and "rear" indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, they are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, the terms used to describe positional relationships in the drawings are only for illustrative purposes and should not be construed as limiting the present invention. For those skilled in the art, the specific meaning of the above terms can be understood according to the specific circumstances.

[0057] Figure 1 This is a schematic diagram of the geometry of a finite-angle fan-beam CT scanner, as shown below. Figure 1 As shown, a right-handed Cartesian coordinate system O-x1x2 is established with respect to the ray source S and the rotation center O. (x1,x2) are the coordinates of the point to be reconstructed, T represents the trajectory of the ray source (finite angle Φ), and θ is the angle between the central ray beam of the fan-beam and the x2 axis under the current projection viewpoint, i.e., the projection angle. s is the distance between the ray passing through the point to be reconstructed (x1,x2) at the detector position and the detector center. g δ (θ,s) represents the projection data collected at the projection angle θ and the distance s between the detector position and the detector center.

[0058] Figure 2 A flowchart of a fast finite-angle fan-beam CT image reconstruction method based on preconditioning matrix technology provided by the present invention is shown below. Figure 2 As shown, the method specifically includes the following steps:

[0059] S1: Acquire projection data: Rotate the X-ray source S around the rotation center O along the scanning track T by a finite angle to obtain projection data;

[0060] S2: Establish a finite-angle fan-beam CT image reconstruction model.

[0061] To suppress artifacts and protect boundaries, this invention employs a compact wavelet frame transform to decompose the reconstructed image into low-frequency and high-frequency components, and applies L0 sparse regularization constraints. To achieve a smoother image, L0 regularization constraints are applied to the image gradient. The model established in this invention is as follows:

[0062]

[0063] Where A is the finite-angle CT projection operator, f∈R N×1 The image to be reconstructed is N, which represents the number of reconstructed image sizes. δ ∈R M×1 This is finite-angle CT projection data, where M represents the total number of rays; λ i and Here, is the non-negative regularization parameter, i represents the number of wavelet subbands; W is the compact wavelet frame transform; ||·||0 is a function of the number of non-zero elements in the statistical vector. in The component form is f i′,j′ This represents the (i′, j′)th pixel of the image.

[0064] S3: Finite Angle CT Iterative Reconstruction: Based on the established model (1), variable substitution h and precondition matrix D are introduced to transform model (1) into the following form:

[0065]

[0066] S31: Design the precondition matrix D, which includes the following steps:

[0067] S311: According to Figure 1 The scan geometry and the finite-angle CT analytical reconstruction algorithm are as follows:

[0068]

[0069] Where f(r,φ) is the reconstructed image in polar coordinates, and ρ(s) represents a one-dimensional ramp filter. This indicates a weighted filtering of the projection. To take advantage of the analytical reconstruction algorithm, according to the analytical reconstruction algorithm (3), the precondition matrix D should contain a weighted diagonal matrix J, whose diagonal elements are... In addition, it should include a filter matrix P, which is a one-dimensional ramp filter matrix H with a frequency domain response of |ω|. (P = F)-1 H(ω)F, where F represents the Fourier transform, F -1 (This represents the inverse Fourier transform).

[0070] S312: To suppress noise and smooth the image, a Hamming window function W is introduced. ham (q)=rect(q)(0.54-0.46cos(2πq)),q∈[-0.5,0.5],rect(q) is a rectangular window function.

[0071] S313: The precondition matrix D is designed as: D T D=W ham PJ.

[0072] S32: Finite-angle CT iterative reconstruction: Based on the model (2) established by introducing the precondition matrix D, the alternating direction multiplier method is used to solve model (2); the specific process is as follows:

[0073] First, model (2) is transformed into an unconstrained optimization problem using the Lagrange augmented function:

[0074]

[0075] Then, the alternating direction multiplier method is used to solve the problem, and the iterative format is as follows:

[0076]

[0077] Here, t > 0 is a parameter introduced by the ADMM algorithm.

[0078] To avoid the difficulty of finding the inverse of the system matrix A and the coupling function in the subproblem f of the iterative scheme (5), this invention incorporates the idea of ​​adjacent alternating linearization, transforming the iterative scheme (5) with respect to subproblems f and h into the adjacent alternating linearization ADMM iterative scheme, as follows:

[0079]

[0080] Find the optimal solution to the subproblem of iterative formula (6), which is shown below:

[0081]

[0082] Among them, A T This represents the fan-beam CT backprojection operator.

[0083] (7) Using variable substitution for V k ←D T V k The iterative format for scaling is then:

[0084]

[0085] (8) Regarding f k+1 The question is:

[0086]

[0087] The iteration format for solving equation (9) using the L0 minimization method is as follows:

[0088]

[0089] The iterative form is as follows: for all image coordinates i′, j′,

[0090]

[0091] Where F represents the Fourier transform, F -1 F(·) represents the inverse Fourier transform. * Represents the complex conjugate of the Fourier transform. Let x and y represent the gradient operators respectively; β controls The similarity parameter, κ (κ>1) represents the parameter controlling the growth rate of β; n represents The number of iterations of the algorithm, The algorithm stops when β is greater than the preset parameter β before iteration. max ,when The algorithm outputs image f after it stops. k+1 .

[0092] According to the above method, the finite-angle CT iterative reconstruction in step S32 includes the following 5 sub-steps:

[0093] S321: Hard thresholding under compact wavelet frame transform for artifact and noise suppression. Based on the design of the precondition matrix D in step S31, it can be seen that... It is an iterative analytical reconstruction method, therefore it has a very fast reconstruction speed;

[0094] S322: Inverse transform of compact wavelet frames and nonnegative constraints.

[0095] S323: Perform the following steps on the results of S322: Minimize smoothing to further suppress artifacts and noise, and perform smoothing processing;

[0096] S324: Projection error adjustment;

[0097] S325: Update the dual variable v. When a certain number of iterations N is reached... Tot Stop the iteration if the iteration fails; otherwise, repeat steps S321-S325.

[0098] S4: Output the reconstructed image. When the iterative algorithm in step S32 stops iterating, output the reconstructed image.

[0099] Figure 3 This is a comparison chart of reconstruction results for a scanning angle range of [0, 140°]. Figure 3 (a) shows the results after 600 iterations without using the preconditioning matrix technique. Figure 3 (b) shows the results of 100 iterations using the preconditioning matrix technique. Figure 3 It can be seen that, Figure 3 The precondition matrix technique of this invention can accelerate the reconstruction speed and can reconstruct images faster than methods that do not use this technique.

[0100] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A fast finite-angle fan-beam CT image reconstruction method based on a preconditioning matrix, characterized in that, The method specifically includes the following steps: S1: Obtain projection data; S2: Considering the sparsity of the image under compact wavelet frame transform and image gradient transform, a doubly regularized finite-angle fan-beam CT image reconstruction model is established, expressed as: (1) in, A It is a finite-angle CT projection operator. The image to be reconstructed. N The number of images representing the reconstructed image size. It is finite angle CT projection data. M Indicates the total number of rays; and It is a non-negative regularization parameter. i Indicates the number of wavelet subbands; W It is a compact wavelet frame transform; It is a function of the number of non-zero elements in a vector. ,in The component form is , , The image represents the first 1 pixel; S3: Design the precondition matrix based on the analytical reconstruction algorithm. This is then incorporated into the dual-regularized finite-angle fan-beam CT image reconstruction model established in step S2, that is, model (1) is transformed into the following form: (2) in, h As an auxiliary variable; Designed precondition matrix for: (3) in, This represents the Hamming window function. Represents the filter matrix. Represents a weighted diagonal matrix; The model is solved using the alternating direction multiplier method; S4: Output the reconstructed image.

2. The fast finite-angle fan-beam CT image reconstruction method according to claim 1, characterized in that, In step S3, the alternating direction multiplier method is used to solve the model, specifically including: First, model (2) is transformed into an unconstrained optimization problem using the Lagrange augmented function: (4) in, V Represents the Lagrange multipliers. These are parameters introduced by the ADMM algorithm; Then, the alternating direction multiplier method is used to solve the problem, and the iterative formula is as follows: (5) in, k Indicates the number of iterations; By employing the adjacent alternating linearization method, the iterative formula (5) is transformed with respect to the subproblems. and Transformed into the ADMM iterative formula of adjacent alternating linearization: (6) in, Indicates the first k Intermediate results of the next iteration This represents the parameters introduced by the proximity operator; Find the optimal solution to the subproblem of iterative formula (6), as shown in the iterative formula below: (7) in, Fan-beam CT backprojection operator; hard thresholding operator , Indicates the threshold. Represents any vector, Indicates the first k Wavelet coefficients after hard thresholding in the next iteration; Formula (7) uses variable substitution Then the iterative formula for scaling is: (8) Formula (8) regarding The question is: (9) use L The method of minimizing 0 solves equation (9), and its iterative formula is as follows: (10) in, Indicates to The result of the L0 algorithm, its iterative form is as follows: for all image coordinates , (11) in, Indicates Fourier transform, Indicates the inverse Fourier transform. Represents the complex conjugate of the Fourier transform. , They represent x , y The gradient operator; Indicates control Similarity parameters Indicates control Parameters related to growth rate; express The number of iterations of the algorithm, The algorithm's stopping criteria are Greater than the preset parameters before iteration ,when The algorithm outputs the image after it stops. .

Images