X-ray computer layered imaging reconstruction method and device

Abstract

The invention relates to an X-ray computer layered imaging reconstruction method and device, image reconstruction is carried out based on an X-ray computer layered imaging reconstruction model, and gradient sparse terms in the x, y and z directions of an image and low-rank constraint terms of the image are arranged in the image reconstruction model. The method can effectively improve the image quality of layered imaging reconstruction of the X-ray computer.

Description

Technical Field

[0001] The present invention relates to a method and apparatus for X-ray computed laminographic reconstruction. Background Art

[0002] Computed Tomography (CT) is a relatively mature non-destructive testing method and has played an important role in fields such as industry and medicine. As shown in Figure 1 (a), in a typical cone-beam CT system, an object is placed between an X-ray source and a flat detector. The X-ray beam generated from the X-ray source is attenuated by the object and collected by the detector. However, for planar objects with a large aspect ratio such as printed circuit boards (PCBs), wings, and solar panels, CT has some limitations. First, CT requires a 360° scan of the sample to be detected. However, for certain specific detection scenarios, such as limited detection space or the sample to be detected being difficult to rotate, a full-angle scan cannot be achieved. In this case, some of the collected projection data is missing, resulting in a low imaging quality. Second, even if a full-angle scan can be performed, the energy attenuation of the X-ray beam along the length direction of the detected plate-like component is obvious and it is even difficult to penetrate, resulting in serious image artifacts.

[0003] Computed Laminography (CL) provides a feasible method to solve the above problems. With the emergence of flat-panel detectors, cone-beam CL reconstruction has become a current research hotspot. Typical cone-beam CL scanning geometries are divided into the following four categories: planar type, swing type, rotation type, and more complex geometries. As shown in Figure 1 (b), the rotation type CL system shown detects a plate-like object, and the angle θ between the rotation axis and the central X-ray is less than 90°. This scanning geometry is not limited by the scale of the planar object, making the length of most X-ray beam penetration paths close to the height of the sample and reducing the working X-ray beam energy. Tiny structures in the object are more likely to be detected in the projection image. In addition, the imaging field of view is extended, allowing the scanning of large-sized objects.

[0004] Although CL scanning solves the problem that traditional CT is not suitable for planar object detection, its imaging characteristics still make the obtained projection data incomplete. Specifically, due to the lack of X-ray beams penetrating the long side of the object, we cannot obtain projection information that directly reflects the model hierarchy. According to the "visible boundary and invisible boundary" theory, if the edge of the scanned object is tangent to the straight line in the projection dataset, the edge can be easily reconstructed from these projections. Otherwise, the edge should be difficult to reconstruct. This theory indicates that applying traditional reconstruction methods to CL data is likely to cause interlayer aliasing and cone angle artifacts in the reconstructed image, specifically manifested as blurring and structural degradation along the vertical direction of the reconstructed image. Therefore, high-quality CL image reconstruction has become an urgent technical problem in the field of nondestructive testing.

[0005] Modern image reconstruction algorithms are mainly divided into analytical reconstruction algorithms and iterative reconstruction algorithms. The analytical method is based on the Radon transform and attempts to restore complete projection data through image processing techniques, such as the Filtered BackProjection (FBP) algorithm and the Feldkamp-Davis-Kress (FDK) algorithm. The analytical method is simple to implement and has a fast reconstruction speed. However, it performs poorly under incomplete projection data. In this case, the iterative method, which mainly aims to solve linear equations, shows significant advantages. Classical iterative methods include the algebraic reconstruction technique (ART), the simultaneous algebraic reconstruction technique (SART), etc. They do not use any prior constraints and are only applicable to occasions with good projection quality. With the development of the Compressed Sensing (CS) theory, recent methods attempt to incorporate various priors into the iterative reconstruction framework to improve the reconstruction quality. It is worth noting that CS-based imaging methods utilize the sparsity of the image to be reconstructed in various transform domains (such as wavelet, gradient, Fourier transform domains, etc.). For example, Sajid et al. designed a spherical sinusoidal geometry and demonstrated the effectiveness of the Total Variation (TV) minimization method for sparse-view CL. Liu et al. used the Truncated Adaptive-Weight Total Variation (TAwTV) to improve the overall CL image. Lu et al. proposed an anisotropic constraint on the sparsity of the image gradient along three orthogonal directions to reduce the interlayer aliasing and blurring observed in the reconstruction of existing CL algorithms. However, the existing methods mentioned above have the drawback that they only focus on the local structure of the image and largely ignore the useful information from distant voxels, making it difficult to effectively solve the problems of image cone angle artifacts and structural degradation. Summary of the Invention

[0006] The object of the present invention is to provide an X-ray computed tomography reconstruction method and apparatus, which can effectively improve the image quality of X-ray computed tomography reconstruction.

[0007] Based on the same inventive concept, the present invention has two independent technical solutions:

[0008] 1. An X-ray computed tomography reconstruction method, which performs image reconstruction based on an X-ray computed tomography reconstruction model. It is characterized in that, in the image reconstruction model, there are gradient sparse terms in the x, y, and z directions of the image, and a low-rank constraint term of the image is provided.

[0009] Furthermore, the gradient sparse terms are set based on a gradient sparse model, and the gradient sparse model is as follows.

[0010] The total variation (TV) minimization problem of the image is expressed as:

[0011]

[0012] Wherein, is the degraded observed value, is the target image to be restored, and TV(X) represents the total variation of the reconstructed image, and its definition is:

[0013] Where is the gradient operator, and ‖·||1 is the l1 norm. respectively represent the gradient operators along the x, y, and z directions of the image.

[0014] Furthermore, the low-rank constraint term is set based on a tensor low-rank model, and the tensor low-rank model is as follows.

[0015] The tensor nuclear norm (TNN) minimization problem of the image is expressed as:

[0016]

[0017] Wherein, and are third-order real-valued tensors, and ||X|| * represents the tensor nuclear norm (TNN) of X, and λ is the regularization parameter;

[0018] Let X = U·S·V * be the tensor singular value decomposition (t-SVD) of X, then ||X|| * is defined as:

[0019]

[0020] where \(r = \text{rank}\) t \(\text{tubal - rank}(X)\) represents the tubal rank of \(X\), and is defined as:

[0021] \(\text{rank}\) t (X)=\(\#\{i|S(i,i,1)\neq0\}\).

[0022] where \(\#\{\cdot\}\) is a counting operator.

[0023] Furthermore, image reconstruction is performed based on the following X - ray computed tomography reconstruction model

[0024]

[0025] In the formula represents the system matrix, and its element \(a\) ij is calculated from the intersection length of the \(i\) - th X - ray beam passing through the \(j\) - th voxel;

[0026] represents the three - dimensional image vector to be reconstructed; \(i = 1,2,\cdots,N_1\), \(j = 1,2,\cdots,N_2\), \(k = 1,2,\cdots,N_3\);

[0027] represents the projection data vector;

[0028] and are the gradient operators in the \(x\), \(y\), and \(z\) directions respectively, and are defined as:

[0029]

[0030] \(\|\cdot\|_1\) represents the \(l_1\) norm, and is defined as:

[0031]

[0032] \(\|f\|\) * represents the tensor nuclear norm. Let \(f = U\cdot S\cdot V\) * be the tensor singular value decomposition of \(f\), then \(\|f\|\) * is defined as:

[0033]

[0034] where \(r = \text{rank}\) t \(\text{tubal - rank}(f)\) represents the tubal rank of \(f\), and is defined as:

[0035] \(\text{rank}\) t (X)=\(\#\{i|S(i,i,1)\neq0\}\)

[0036] where \(\#\{\cdot\}\) is a counting operator;

[0037] is a data fidelity term, and are the gradient sparsity terms in the x, y, and z directions respectively, ||f|| * is the low-rank constraint term;

[0038] λ1, λ2, λ3, and λ4 are regularization factors used to control the weighting of each term.

[0039] Furthermore, the reconstruction model is solved based on the Chambolle-Pock algorithm.

[0040] Furthermore, the reconstruction model is solved based on the Chambolle-Pock algorithm. The specific method is as follows:

[0041] Step 1: Obtain the projection data g, and set the parameters λ1, λ2, λ3, λ4, and N iter numerical value;

[0042] Step 2: Calculate L = ||K||2, τ = 1 / L, σ = 1 / L;

[0043] Set γ = 1, n = 0, u0 = 0, y0 = 0, p0 = 0, q0 = 0, r0 = 0;

[0044] Step 3: Calculate

[0045] Step 4: Calculate

[0046] Step 5: Calculate

[0047] Step 6: Calculate

[0048] Step 7: Calculate

[0049] Step 8: Calculate

[0050] Step 9: Calculate n = n + 1; if n is less than N iter , then return to Step 3;

[0051] Step 10: Output the reconstructed image x n ;

[0052] where N iter is the total number of iterations of the algorithm;

[0053] L = ||K||2 is a constant, which is set to the l2 norm of matrix K, i.e., the largest singular value of matrix K;

[0054] $n$ represents the number of algorithm iterations, where $n = 0, 1, \ldots, N$ iter ;

[0055] $I$ represents an image with all voxel values being 1;

[0056] The superscript $T$ represents transpose;

[0057] represents the solution to the tensor nuclear norm minimization problem, obtained through the tensor singular value thresholding operator;

[0058] In the formula, $x = f$, $y = Af$,

[0059] Furthermore, in the PCB image reconstruction, for noiseless projection data, the parameters are set as $\lambda_1=\lambda_2 = 0.008$, $\lambda_3 = 0.0005$, $\lambda_4 = 0.001$; for noisy projection data, the parameters are set as $\lambda_1=\lambda_2 = 0.01$, $\lambda_3 = 0.0005$, $\lambda_4 = 0.001$. In the workpiece image reconstruction, for noiseless projection data, the parameters are set as $\lambda_1=\lambda_2 = 0.02$, $\lambda_3 = 0.001$, $\lambda_4 = 0.002$; for noisy projection data, the parameters are set as $\lambda_1=\lambda_2 = 0.025$, $\lambda_3 = 0.001$, $\lambda_4 = 0.002$.

[0060] Furthermore, for the PCB image reconstruction, the total number of algorithm iterations $N$ iter is set to 15000.

[0061] Furthermore, for the workpiece image reconstruction, the total number of algorithm iterations $N$ iter is set to 5000.

[0062] 2. An X-ray computed tomography reconstruction device for performing the method according to claims 1 - 9 above.

[0063] The beneficial effects of the present invention are:

[0064] The present invention performs image reconstruction based on an X-ray computed tomography reconstruction model, characterized in that in the image reconstruction model, there are gradient sparse terms in the $x$, $y$, and $z$ directions of the image, and a low-rank constraint term of the image. The gradient sparse terms are set based on the gradient sparse model, and the gradient sparse model is as follows,

[0065] The total variation TV minimization problem of the image is expressed as:

[0066]

[0067] where, is the degraded observed value, Let \(X\) be the target image to be restored, and \(TV(X)\) represents the total variation of the reconstructed image, which is defined as:

[0068] where \(\nabla\) is the gradient operator, and \(\|\cdot\|_1\) is the \(l_1\) norm, \(\nabla_x\), \(\nabla_y\), and \(\nabla_z\) respectively represent the gradient operators along the \(x\), \(y\), and \(z\) directions of the image.

[0069] Set the low-rank constraint term based on the tensor low-rank model, and the tensor low-rank model is as follows.

[0070] The problem of minimizing the tensor nuclear norm \(TNN\) of the image tensor is expressed as:

[0071]

[0072] where and are third-order real-valued tensors, \(\|X\| * represents the tensor nuclear norm \(TNN\) of \(X\), and \(\lambda\) is the regularization parameter;

[0073] Let \(X = U\cdot S\cdot V * be the tensor singular value decomposition \(t - SVD\) of \(X\), then \(\|X\| * is defined as:

[0074]

[0075] where \(r = rank t (X)\) represents the tubal rank of \(X\), which is defined as:

[0076] rank t (X)=\#{i|S(i,i,1)\neq0}\).

[0077] where \(\#{·}\) is the counting operator.

[0078] The image reconstruction of the present invention is carried out based on the following X-ray computed tomography reconstruction model.

[0079]

[0080] In the formula, \(A\) represents the system matrix, and its element \(a ij is calculated from the intersection length of the \(i\)-th X-ray beam passing through the \(j\)-th voxel;

[0081] \(x\) represents the three-dimensional image vector to be reconstructed; \(i = 1,2,\cdots,N_1\), \(j = 1,2,\cdots,N_2\), \(k = 1,2,\cdots,N_3\);

[0082] \(b\) represents the projection data vector.

[0083] and are the gradient operators in the x, y, and z directions respectively, and their definitions are as follows:

[0084]

[0085] ||·||1 represents the l1 norm, and its definition is:

[0086]

[0087] ||f|| * represents the tensor nuclear norm. Let f = U·S·V * be the tensor singular value decomposition of f, then ||f|| * is defined as:

[0088]

[0089] where r = rank t (f) represents the tubal rank of f, and its definition is:

[0090] rank t (X) = #{i|S(i,i,1)≠0}

[0091] where #{·} is the counting operator;

[0092] is the data fidelity term, and are the gradient sparsity terms in the x, y, and z directions respectively, and ||f|| * is the low-rank constraint term;

[0093] λ1, λ2, λ3, and λ4 are regularization factors, which are used to control the weighting of each term respectively.

[0094] The present invention proposes a computerized laminography (CL) reconstruction method (model) based on anisotropic gradient sparsity and low-rank regularization (AGSLR). The local transform sparsity is described by the one-dimensional total variation (TV) along three orthogonal directions, and the global spatial correlation is characterized by the tensor nuclear norm (TNN). The global spatial correlation and local smoothness of the object are fully considered, and the anisotropic TV and tensor nuclear norm regularization are integrated into a unified framework to complement each other, thus better characterizing the inherent structure of the image and effectively solving the problems of image cone angle artifacts, structural degradation, and noise suppression, and effectively improving the image quality of X-ray computerized laminography reconstruction. In the present invention, one-dimensional TV constraints with different intensities are added along the x, y, and z directions of the reconstructed image to characterize the sparsity in the gradient domain of the CL image, and the low-rank correlation descriptor tensor nuclear norm (TNN) is used to constrain the non-local similarity of the image, not only eliminating image artifacts and noise through the TV regularizer but also inheriting the advantages of the TNN in image feature recovery.

[0095] The present invention solves the reconstruction model based on the Chambolle-Pock algorithm, and the specific method is as follows.

[0096] Step 1: Obtain the projection data g, and set the parameters λ1, λ2, λ3, λ4, and N iter numerical values;

[0097] Step 2: Calculate L = ||K||2, τ = 1 / L, σ = 1 / L;

[0098] Set γ = 1, n = 0, u0 = 0, y0 = 0, p0 = 0, q0 = 0, r0 = 0;

[0099] Step 3: Calculate

[0100] Step 4: Calculate

[0101] Step 5: Calculate

[0102] Step 6: Calculate

[0103] Step 7: Calculate

[0104] Step 8: Calculate

[0105] Step 9: Calculate n = n + 1; if n is less than N iter , then return to Step 3;

[0106] Step 10: Output the reconstructed image x n ;

[0107] Wherein, N iter is the total number of iterations of the algorithm;

[0108] L = ||K||2 is a constant, which is set to the l2 norm of matrix K, that is, the largest singular value of matrix K;

[0109] n represents the iteration number of the algorithm, n = 0, 1,..., N iter ;

[0110] I represents an image with all voxel values being 1;

[0111] The superscript T represents transpose;

[0112] represents the solution of the tensor nuclear norm minimization problem, which is obtained through the tensor singular value thresholding operator;

[0113] Wherein, x = f, y = Af,

[0114] The present invention solves the reconstruction model through the above calculation method to obtain a reconstructed image. Experiments show that the calculation method proposed by the present invention can effectively restore the image edge and suppress aliasing artifacts, and this method is more powerful than other existing calculation methods in terms of accuracy.

[0115] When the present invention solves the reconstruction model through the above calculation method,

[0116] In PCB image reconstruction, for noiseless projection data, the parameters are set as λ1 = λ2 = 0.008, λ3 = 0.0005, λ4 = 0.001; for noisy projection data, the parameters are set as λ1 = λ2 = 0.01, λ3 = 0.0005, λ4 = 0.001. In workpiece image reconstruction, for noiseless projection data, the parameters are set as λ1 = λ2 = 0.02, λ3 = 0.001, λ4 = 0.002; for noisy projection data, the parameters are set as λ1 = λ2 = 0.025, λ3 = 0.001, λ4 = 0.002.

[0117] For PCB image reconstruction, the total number of iterations N of the algorithm iter is set to 15000; for workpiece image reconstruction, the total number of iterations N of the algorithm iter is set to 5000. Through the setting of the above parameters, the present invention further ensures the quality of the X-ray computed tomography reconstruction image. BRIEF DESCRIPTION OF THE DRAWINGS

[0118] Figure 1 is the working principle diagram of an existing computed tomography (CT) system and a computed laminography (CL) system;

[0119] Figure 2 is the PCB model and workpiece model images;

[0120] Figure 3 is the reconstructed image of the PCB model;

[0121] Figure 4 is the one-dimensional contour map of the horizontal section at z = 16 of the reconstructed PCB model;

[0122] Figure 5 is the reconstructed image of the workpiece model;

[0123] Figure 6 is the absolute difference image between the reconstructed result of the workpiece model and the real image;

[0124] Figure 7 is the ROI enlarged view of the reconstructed result of the workpiece model without noise;

[0125] Figure 8 is the ROI enlarged view of the reconstructed result of the workpiece model with added noise;

[0126] Figure 9 is the RMSE and PCC curve graph of the AGSLR-CP algorithm of the present invention. Detailed implementation manners

[0127] The present invention will be described in detail below in conjunction with the various implementation manners shown in the drawings. However, it should be noted that these implementation manners do not limit the present invention, and any equivalent transformation or substitution in function, method, or structure made by those of ordinary skill in the art according to these implementation manners shall fall within the protection scope of the present invention.

[0128] Example 1:

[0129] X-ray computed tomography reconstruction method

[0130] The X-ray computed tomography (CL) reconstruction problem can be represented by the following discrete-to-discrete linear system:

[0131] Af = g, (1)

[0132] where, represents the three-dimensional image vector to be reconstructed. We also adopt this form of f, i = 1, 2, …, N1, j = 1, 2, …, N2, k = 1, 2, …, N3; the conversion between the two forms can be written as n = (k - 1) × N1 × N2 + (j - 1) × N1 + i. is the projection data vector. is the system matrix, and its element a ijCalculated from the intersection length of the \(i\)-th X-ray beam passing through the \(j\)-th voxel. Image reconstruction is to find the reconstructed image \(f\) by solving Equation (1). However, due to the limitation of computer memory, the solution of Equation (1) usually cannot be directly obtained through the inverse matrix. In the prior art, generally, by adding prior constraints on the quantity \(f\) to be solved, it is transformed into the following objective function:

[0133]

[0134] Among them, the first term is the data fidelity term that constrains the reconstructed image and the projection data, and the second term is the regularization term defined based on prior knowledge. \(\lambda\) is the regularization parameter, and a suitable value needs to be set to balance the proportion of the fidelity term and the regularization term. In order to reconstruct an image without interlayer aliasing and obvious structural degradation from incomplete CL data, \(\varPhi(\cdot)\) is the key to CL reconstruction.

[0135] The present invention proposes an X-ray computed tomography reconstruction method, which performs image reconstruction based on an X-ray computed tomography reconstruction model. In the image reconstruction model, gradient sparse terms in the \(x\), \(y\), and \(z\) directions of the image are provided, and a low-rank constraint term of the image is provided.

[0136] (1) Setting the gradient sparse term based on the gradient sparse model

[0137] The gradient sparse model is as follows,

[0138] The total variation TV minimization problem of the image is expressed as:

[0139]

[0140] Among them, is the degraded observed value, is the target image to be restored, and \(TV(X)\) represents the total variation of the reconstructed image, and its definition is:

[0141]

[0142] Among them is the gradient operator, and \(\|\cdot\|_1\) is the \(l_1\) norm, respectively represent the gradient operators along the \(x\), \(y\), and \(z\) directions of the image.

[0143] (2) Setting the low-rank constraint term based on the tensor low-rank model

[0144] The tensor low-rank model is as follows:

[0145] The tensor nuclear norm TNN minimization problem of the image is expressed as:

[0146]

[0147] Among them, and is a third-order real-valued tensor, ||X|| * denotes the tensor nuclear norm TNN of X, and λ is the regularization parameter;

[0148] Let X = U·S·V * be the tensor singular value decomposition t-SVD of X, then ||X|| * is defined as:

[0149]

[0150] where r = rank t (X) represents the tubal rank of X, which is defined as:

[0151] rank t (X) = #{i|S(i,i,1)≠0}. (7)

[0152] where #{·} is the counting operator. TV is isotropic, and minimizing TV tends to equally penalize all image gradients in different directions.

[0153] where the closed-form solution D λ (Y) can be obtained by the existing tensor singular value thresholding (t-SVT) operator. Algorithm 1 below gives the process of t-SVT, where fft(·,[],3) and ifft(·,[],3) denote the discrete Fourier transform (DFT) and inverse discrete Fourier transform (IDFT) along the third dimension of the tensor, respectively. denotes the integer greater than or equal to t. Γ λ (S) is defined as Γ λ (S)(i,i) = max(S(i,i) - λ, 0), i = 1,..., min(N1,N2). conj(·) represents the complex conjugate of a matrix.

[0154]

[0155] (III) The X-ray computed tomography reconstruction model of the present invention

[0156] The ideal image to be reconstructed is usually constant in regions or has low-level variations, which makes the local structures in the image sparse or compressible under certain transformations. Theoretically, it is easily proven to be piecewise constant in the x, y, and z directions. Additionally, the imaging characteristics of CL result in different structural recovery capabilities along different directions. Therefore, in the present invention, one-dimensional TV terms with different intensities are added along three orthogonal directions of the reconstructed image to constrain the gradient sparse features of the reconstructed image in each direction in an anisotropic form. Since the horizontal edges in the usually reconstructed CL image can be relatively accurately reconstructed while the vertical edges are prone to distortion, the constraint intensity should mainly focus on the x and y directions, and the constraint intensity in the z direction is the weakest.

[0157] In addition to having gradient sparsity, the reconstructed CL image often has a large number of similar structures, such as pads, vias, and wires in a PCB, and spars, struts, and rivets in a wing. The data matrix formed by the stacking of these similar blocks usually has the property of low rank, while noise does not have the low-rank property. Therefore, imposing a low-rank constraint on the image matrix can achieve the effect of noise reduction while restoring the image structure and details. Therefore, regarding the CL image as a third-order tensor, TNN can be used as a low-rank constraint term to describe the global spatial correlation characteristics.

[0158] Based on the above analysis, the reconstructed CL image usually exhibits obvious regional constancy characteristics and non-local self-similarity. Therefore, the present invention proposes an X-ray computed laminography reconstruction model (AGSLR) based on anisotropic gradient sparsity and low-rank regularization, as follows:

[0159]

[0160] In the formula, represents the system matrix, and its element a ij is calculated from the intersection length of the i-th X-ray beam passing through the j-th voxel;

[0161] represents the three-dimensional image vector to be reconstructed; i = 1, 2, …, N1, j = 1, 2, …, N2, k = 1, 2, …, N3;

[0162] represents the projection data vector;

[0163] and are the gradient operators in the x, y, and z directions respectively, and their definitions are:

[0164]

[0165] ||·||1 represents the l1 norm, and its definition is:

[0166]

[0167] ||f|| * represents the tensor nuclear norm. Let f = U·S·V * be the tensor singular value decomposition of f, then ||f|| * is defined as:

[0168]

[0169] where r = rank t (f) represents the tubal rank of f, and its definition is:

[0170] rank t (X) = #{i|S(i,i,1)≠0}

[0171] where #{·} is the counting operator;

[0172] is the data fidelity term, and are the gradient sparsity terms in the x, y, and z directions respectively, ||f|| * is the low-rank constraint term;

[0173] λ1, λ2, λ3, and λ4 are regularization factors, which are used to control the weight of each term respectively.

[0174] (IV) Solving the reconstruction model based on the Chambolle-Pock algorithm

[0175] To solve the minimization problem involved in the X-ray computer tomography reconstruction model (AGSLR) of the present invention, the present invention uses the CP (Chambolle-Pock) algorithm to solve it.

[0176] 1. CP algorithm framework

[0177] The general form of the original minimization problem applicable to the CP algorithm is as follows:

[0178]

[0179] And the following dual maximization problem:

[0180]

[0181] where x and y are two finite-dimensional real vectors in spaces X and Y, K is a linear transformation from X to Y, F and G are convex functions, and F* and G* are their convex conjugates. The superscript T represents the transpose of a matrix.

[0182] For a convex function H(z) with a variable z ∈ Z, its convex conjugate can be obtained by the following operation:

[0183]

[0184] Among them, <·,·> Z represents the inner product in the vector space Z.

[0185] The framework of the CP algorithm is shown in Algorithm 2 below. Among them, L is the l2 norm of matrix K, that is, the largest singular value of matrix K. τ and σ are non - negative parameters, both set to 1 / L. N iter is the total number of iterations of the algorithm. θ ∈ [0,1], set to 1 here. When θ = 1 and satisfies , the algorithm can achieve convergence with a speed of O(1 / n). prox σ [F * and prox τ [G] are two proximal mappings.

[0186]

[0187] For the convex function H(z), its proximal mapping is calculated as follows:

[0188]

[0189] The key to deriving an instance of the CP algorithm for a specific optimization model is to calculate the convex conjugate function F * and the two proximal mapping operators prox σ [F * and prox τ [G].

[0190] 2. Derivation of the AGSLR - CP algorithm instance of the present invention

[0191] To derive an instance of the CP algorithm, we constructed the following associations for equations (8) and (9): x = f, y = Af,

[0192] F(y,p,q) = F1(y)+F2(p)+F3(q)+F4(r), (14)

[0193]

[0194] G(x) = λ4||x|| * , (16)

[0195] The terms F1(y), F2(p), F3(q) and F4(r) in equation (14) are all convex functions, so the function F(y,p,q) is convex. G(x) is also a convex function.

[0196] According to formula (11), the convex conjugates of F1(y), F2(p), F3(q) and F4(r) are obtained as:

[0197]

[0198]

[0199] In formulas (5.23) to (5.25), the definition of the indicator function δ Box(a) (·) is as follows:

[0200]

[0201] where, ||·|| ∞ norm is the maximum value of the absolute values of the components in the vector, and Box(a) consists of vectors whose all components are not greater than a.

[0202] Derive F1 * (y), F4 * (r) and the proximal mappings of G(x) are as follows:

[0203]

[0204] The voxels in the I in formulas (23) to (25) are all 1 in the image. Substitute formulas (22) to (26) into Algorithm 1 to obtain Algorithm 3, which is the CP algorithm instance of the AGSLR model of the present invention.

[0205] (III) The present invention solves the reconstruction model based on the Chambolle-Pock algorithm

[0206] The specific method (Algorithm 3) is as follows

[0207] Step 1: Obtain the projection data g, and set the parameters λ1, λ2, λ3, and λ4, N iter numerical value;

[0208] Step 2: Calculate L = ||K||2, τ = 1 / L, σ = 1 / L;

[0209] Set γ = 1, n = 0, u0 = 0, y0 = 0, p0 = 0, q0 = 0, r0 = 0;

[0210] Step 3: Calculate

[0211] Step 4: Calculate

[0212] Step 5: Calculate

[0213] Step 6: Calculate

[0214] Step 7: Calculate

[0215] Step 8: Calculate

[0216] Step 9: Calculate n = n + 1; if n is less than N iter , then return to Step 3;

[0217] Step 10: Output the reconstructed image x n ;

[0218] wherein, N iter is the total number of algorithm iterations;

[0219] L = ||K||2 is a constant, which is set to the l2 norm of matrix K, that is, the largest singular value of matrix K;

[0220] n represents the algorithm iteration number, n = 0, 1,..., N iter ;

[0221] I represents an image with all voxel values being 1;

[0222] The superscript T represents transpose;

[0223] represents the solution of the tensor nuclear norm minimization problem, which is obtained through the tensor singular value thresholding operator;

[0224] wherein, x = f, y = Af,

[0225] In PCB image reconstruction, for noiseless projection data, the parameters are set as λ1 = λ2 = 0.008, λ3 = 0.0005, λ4 = 0.001; for noisy projection data, the parameters are set as λ1 = λ2 = 0.01, λ3 = 0.0005, λ4 = 0.001. In workpiece image reconstruction, for noiseless projection data, the parameters are set as λ1 = λ2 = 0.02, λ3 = 0.001, λ4 = 0.002; for noisy projection data, the parameters are set as λ1 = λ2 = 0.025, λ3 = 0.001, λ4 = 0.002.

[0226] For PCB image reconstruction, the total number of algorithm iterations N iter is set to 15000.

[0227] For workpiece image reconstruction, the total number of algorithm iterations N iter is set to 5000.

[0228]

[0229] The beneficial effects of the present invention will be further described below in conjunction with experiments.

[0230] The AGSLR-CP algorithm was evaluated using PCB models and workpiece models in the experiment, and was compared with SART, ASD-POCS, POCS-AwDaRTV, and POCS-AAR algorithms. To quantitatively analyze the reconstruction accuracy of each method, image quality evaluation metrics of the experimental results were calculated and listed, including Root Mean Squared Error (RMSE), Pearson Correlation Coefficient (PCC), Mean Structural Similarity (MSSIM), and Universal Quality Index (UQI). All experiments were conducted on a personal computer with a 2.9 GHz Intel Core i7-10700 CPU processor and an NVIDIA GTX 1050Ti graphics card, and were coded in MATLAB R2020a.

[0231] (I) Data acquisition

[0232] We used PCB models and workpiece models for the experiment. The resolution of the PCB model is 256×256×75, and the voxel size is 0.5×0.5×0.5 mm 3 . The PCB model consists of 10 evenly stacked wiring layers, each wiring layer has a thickness of 3 voxels, and the interval between wiring layers is 2 voxels. The workpiece model is a flat cylinder composed of various geometries with different attenuation coefficients. The resolution of the workpiece model is 350×350×50, and the voxel size is 0.5×0.5×0.5 mm 3 . Figure 2 (a) and 2(b) show the horizontal section (upper left), coronal section (upper right), and sagittal section (lower left) of the PCB model and the workpiece model, where the yellow crosshair indicates the position of the other two sections. Two Regions of Interest (ROIs) were selected in the workpiece model, and the enlarged views are shown in Figure 2 (c) and 2(d). The structure of the CL scanning configuration is shown in Figure 1 (b). The scanning geometric parameters are shown in Table 1.

[0233] Table 1 Scanning geometric configuration

[0234]

[0235] In an actual CL imaging system, the noise mainly comes from quantum noise (Poisson noise model) and electronic noise (Gaussian noise model). Therefore, in order to verify the robustness of the proposed method to noise, Poisson noise and Gaussian noise are added to the noiseless data. In the PCB experiment, the incident intensity of Poisson noise is 6×10 4 , and the average value of Gaussian noise is 0, and the standard deviation is 0.5% of the maximum projection value. In the workpiece experiment, Poisson noise with an incident intensity of 2×10 5 is added to the noiseless data, as well as Gaussian noise with an average value of 0 and a standard deviation of 0.5% of the maximum projection value.

[0236] (II) PCB model experiment

[0237] In the experiment, the initial value of each method is set to f = 0. We tested a series of parameters of these methods within a reasonable range and obtained the parameters with the best image reconstruction quality for different situations. For noiseless projections, the parameter settings of the proposed method are λ1 = λ2 = 0.008, λ3 = 0.0005, λ4 = 0.001. For noisy projections, λ1 = λ2 = 0.01, λ3 = 0.0005, λ4 = 0.001. The maximum number of iterations N iter of the SART, ASD-POCS, POCS-AwDaRTV, and POCS-AAR algorithms are all 400, and the N iter of AGSLR-CP is 15000 to ensure that all algorithms obtain their stable solutions.

[0238] Figure 3 shows the reconstruction results of each method for the PCB model. For the reconstruction result of SART, the horizontal section can only show the general outline of the wiring layer but the details are not clear, and the structures in the coronal section and sagittal section are severely aliased. And with the addition of noise, the image quality decreases significantly. In the result of ASD-POCS, although the addition of the TV regularization term suppresses the cone angle artifacts and noise to a certain extent, the structures of other layers can still be faintly seen in the horizontal section, and the interlayer aliasing and blurring are not alleviated. In the reconstruction result of POCS-AwDaRTV, the edges and shapes in the horizontal section are basically restored, and each wiring layer can be easily distinguished, but the gaps between adjacent wiring layers are not effectively reconstructed. For the result of POCS-AAR, the edges in the horizontal section are clearly reconstructed, and the adhesion artifacts between the wiring layers are greatly suppressed, but the overall image is darker. In the result of the AGSLR-CP of the present invention, the structure is closer to the real model, and satisfactory results are obtained even with the addition of noise. This shows that the algorithm of the present invention has strong denoising ability. To clearly show the superiority of the method of the present invention, Figure 4Shows the one-dimensional profiles of the 128th row and 128th column of the reconstructed image horizontal slice. It can be seen that the resulting profiles obtained by the AGSLR-CP algorithm of the present invention fit better with the original model, indicating that this algorithm performs well in preserving the image structure and details. To quantitatively compare the reconstruction accuracies of various methods, Table 2 statistically analyzes the evaluation metrics of the above reconstructed images, where the optimal values are highlighted in bold. We can find that the method of the present invention has advantages in all evaluation metric values, indicating that the method of the present invention can restore the image sequence more effectively than other methods. Figure 3 Is the reconstruction result of the PCB model. From left to right, each column is the reconstruction result of the PCB model obtained by the SART, ASD-POCS, POCS-AwDaRTV, POCS-AAR, and AGSLR-CP methods respectively. The first row and the second row are the reconstruction results of noiseless and noisy projection data respectively. The displayed gray level is [0, 1]. Figure 4 Is the one-dimensional profile of the horizontal section at z = 16 of the reconstructed PCB model. Among them, (a) the 128th column profile in the noiseless case; (b) the 128th row profile in the noiseless case; (c) the 128th column profile in the noisy case; (d) the 128th row profile in the noisy case.

[0239] Table 2 Comparison of RMSE, PCC, MSSIM, and UDI of the reconstructed PCB model by different algorithms

[0240]

[0241] (III) Workpiece model experiment

[0242] In the experiment, for noiseless data, the parameter settings of the proposed method are λ1 = λ2 = 0.02, λ3 = 0.001, λ4 = 0.002. For noisy projection, λ1 = λ2 = 0.025, λ3 = 0.001, λ4 = 0.002. The parameters of other algorithms are determined according to the error and trial techniques. The maximum number of iterations N of the SART, ASD-POCS, POCS-AwDaRTV, and POCS-AAR algorithms iter Are all 150, and the N of AGSLR-CP iter Is 5000.

[0243] The reconstruction results of each method are as Figure 5 Shown, and the absolute residual image between it and the original image is as Figure 6As shown, there are a large number of aliasing artifacts and noises in the reconstruction results of SART, and the visual effect is poor. The image quality of ASD-POCS is significantly improved, and the edges of the horizontal section are effectively restored. However, even in the case of no noise, the edges in the vertical direction are still blurred. POCS-AwDaRTV, POCS-AAR, and AGSLR-CP all achieve good reconstruction results. There are no obvious noises and artifacts in the images, and the interlayer aliasing is significantly suppressed. The image quality of these three methods is very close, and it is difficult to intuitively distinguish their differences. By observing the absolute difference images, it can be found that the reconstructed image of AGSLR-CP is closest to the original model, indicating that this algorithm performs excellently in artifact elimination and structure preservation. Figure 7 and Figure 8 The enlarged views of the ROIs of the reconstruction results of each method under the conditions of no noise and noisy are compared respectively. It can be seen that the edges of some details in the reconstruction results of POCS-AwDaRTV are blurred. Although POCS-AAR can reconstruct the edges of the image structure relatively clearly, its reconstruction accuracy still needs to be improved. The edges of the details in the reconstructed image of the proposed AGSLR-CP are clear and accurate, demonstrating its excellent ability in restoring fine structures. In addition, even in the face of noise interference, AGSLR-CP can still maintain high accuracy and robustness. By comparing the quantitative evaluation indexes of the reconstruction results of each algorithm listed in Table 3, the above analysis and conclusions can be further verified. Figure 5 are the reconstruction results of the workpiece model. From left to right, each column is the reconstruction result of the PCB model obtained by SART, ASD-POCS, POCS-AwDaRTV, POCS-AAR, and AGSLR-CP methods. The first row and the second row are the reconstruction results of noiseless and noisy projection data respectively. The display gray level is [0,1]. Figure 6 are the absolute difference images between the reconstruction results of the workpiece model and the real images. From left to right are the absolute difference images of SART, ASD-POCS, POCS-AwDaRTV, POCS-AAR, and AGSLR-CP methods respectively. The first row and the second row are the absolute difference images of noiseless and noisy projection data respectively. The gray level window is [0,0.3]. Figure 7 are the enlarged views of the reconstruction of ROI1 (the first row) and ROI2 (the second row) in the workpiece model by SART, ASD-POCS, POCS-AAR, POCS-AwDaRTV, and AGSLR-CP algorithms (from left to right) with noiseless data. The display window is [0,1]. The display gray level is [0,1]. Figure 8It is an enlarged view of the noise data reconstruction of ROI1 (the first row) and ROI2 (the second row) in the workpiece model by the SART, ASD-POCS, POCS-AAR, POCS-AwDaRTV, and AGSLR-CP algorithms (from left to right). The display window is [0,1].

[0244] Table 3 Comparison of RMSE, PCC, MSSIM, and UDI of different algorithms for reconstructing the workpiece model

[0245]

[0246] (IV) Conclusion

[0247] Taking the above-mentioned noise-free PCB model and workpiece model data as an example, as Figure 9 shown, the RMSE and PCC curves during the iterative process of the AGSLR-CP algorithm of the present invention are plotted. Figure 9 (a) is the RMSE and PCC curve graph corresponding to the noise-free PCB model; Figure 9 (b) is the RMSE and PCC curve graph of the noise-free workpiece model. When the algorithm parameters are reasonably selected, the CP algorithm is proven to be convergent mathematically and has the convergence of a first-order method. From Figure 9 it can be seen that during the iterative process, the RMSE shows a monotonically decreasing trend on a global scale, while the PCC value monotonically increases, and finally both converge to a stable position, indicating that the intermediate result image gradually approaches the reference image, and the algorithm can effectively optimize the objective function to obtain a satisfactory solution.

[0248] In order to eliminate the cone angle artifacts and structural degradation caused by cone-beam CL reconstruction, the present invention proposes a reconstruction model based on anisotropic gradient sparsity and low-rank regularization (AGSLR). This model fully considers the global spatial correlation and local smoothness of the target, integrates anisotropic TV and tensor nuclear norm regularization into a unified framework, and complements each other, thus better characterizing the inherent structure of the image. Among them, TV uses the local gradient information of neighboring pixels, and the tensor nuclear norm describes the overall properties of the image. During the iterative process, the AGSLR-CP algorithm of the present invention is used to solve the optimization problem. Numerical experiments show that the method of the present invention can effectively restore image edges and suppress aliasing artifacts, and this method is more powerful than existing methods in terms of accuracy.

[0249] Example 2:

[0250] X-ray computed tomography reconstruction device

[0251] The X-ray computed tomography reconstruction device is used to perform the above-mentioned X-ray computed tomography reconstruction method.

[0252] The series of detailed descriptions listed above are only specific descriptions of the feasible embodiments of the present invention, and they are not intended to limit the protection scope of the present invention. Any equivalent embodiments or changes made without departing from the technical spirit of the present invention should be included within the protection scope of the present invention.

[0253] For those skilled in the art, it is obvious that the present invention is not limited to the details of the above-mentioned exemplary embodiments, and the present invention can be implemented in other specific forms without departing from the spirit or basic characteristics of the present invention. Therefore, from any point of view, the embodiments should be regarded as exemplary and non-limiting. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, it is intended to embrace all changes falling within the meaning and scope of the equivalent elements of the claims in the present invention.

Claims

1. An X-ray computed tomography reconstruction method, which performs image reconstruction based on an X-ray computed tomography reconstruction model, characterized in that, In the image reconstruction model, there are gradient sparse terms in the x, y, and z directions of the image, and a low-rank constraint term of the image is provided.

2. The X-ray computed tomography reconstruction method according to claim 1, wherein: The gradient sparse terms are set based on the gradient sparse model, and the gradient sparse model is as follows. The total variation TV minimization problem of the image is expressed as: Among them, is a degenerate observation value, is the target image to be restored, and TV(X) represents the total variation of the reconstructed image, which is defined as: where is the gradient operator, and ‖·‖1 is the l1 norm, representing the gradient operators along the x, y, and z directions of the image, respectively.

3. The X-ray computed tomography reconstruction method according to claim 2, characterized in that: The low-rank constraint term is set based on the tensor low-rank model, and the tensor low-rank model is as follows. The tensor nuclear norm TNN minimization problem of the image is expressed as: Among them, and are third-order real-valued tensors, ||X|| * represents the tensor nuclear norm TNN of X, and λ is the regularization parameter; Let \(X = U\cdot S\cdot V\) * be the tensor singular value decomposition t-SVD of \(X\), then \(\|X\|\) * is defined as: where r = rank t (X) represents the tubal rank of X, which is defined as: rank t (X) = #{i | S(i, i, 1) ≠ 0} Where, #{·} is the counting operator.

4. The X-ray computed tomography reconstruction method according to claim 3, wherein: Image reconstruction is performed based on the following X-ray computed tomography reconstruction model. In the formula, represents the system matrix, and its element a ij is calculated from the intersection length of the i-th X-ray beam passing through the j-th voxel; Represents the three-dimensional image vector to be reconstructed; i = 1, 2, …, N1, j = 1, 2, …, N2, k = 1, 2, …, N3; represent the projection data vector; and are the gradient operators in the x, y, and z directions respectively, and their definitions are as follows: ||·||1 represents the l1 norm, and its definition is: ||f|| * denotes the tensor nuclear norm. Let f = U·S·V * be the tensor singular value decomposition of f, then ||f|| * is defined as: where r = rank t (f) represents the tubal rank of f, which is defined as: rank t (X) = #{i|S(i, i, 1) ≠ 0} Where, #{·} is the counting operator; is a data fidelity term, and are gradient sparsity terms in the x, y, and z directions respectively, and ||f|| * is a low-rank constraint term; λ1, λ2, λ3, and λ4 are regularization factors, which are respectively used to control the action weights of each term.

5. The X-ray computed tomography reconstruction method according to claim 4, characterized in that: The reconstruction model is solved based on the Chambolle-Pock algorithm.

6. The X-ray computed tomography reconstruction method according to claim 5, characterized in that: The reconstruction model is solved based on the Chambolle-Pock algorithm, and the specific method is as follows. Step 1: Obtain projection data g, and set parameters λ1, λ2, λ3, λ4, and N iter numerical values; Step 2: Calculate L = ||K||2, τ = 1 / L, σ = 1 / L; Set γ = 1, n = 0, u0 = 0, y0 = 0, p0 = 0, q0 = 0, r0 = 0; Step 3: Calculate Step 4: Calculate Step 5: Calculate Step 6: Calculate Step 7: Calculate Step 8: Calculate Step 9: Calculate n = n + 1 ； If n is less than N iter , then return to Step 3; Step 10: Output the reconstructed image x n ; Where N iter is the total number of algorithm iterations; L = ||K||2 is a constant, which is set to the l2 norm of matrix K, that is, the largest singular value of matrix K. n represents the number of algorithm iterations, n = 0, 1, ..., N iter ; I represents an image with all voxel values being 1. The superscript T represents the transpose. denotes the solution to the tensor nuclear norm minimization problem, obtained by the tensor singular value thresholding operator; 7. The X-ray computed tomography reconstruction method according to claim 6, wherein: In PCB image reconstruction, for noiseless projection data, the parameters are set as λ1 = λ2 = 0.008 , λ3 = 0.0005 , λ4 = 0.001 ； For noisy projection data, the parameters are set as λ1 = λ2 = 0.01 , λ3 = 0.0005 , λ4 = 0.001 。 In workpiece image reconstruction, for noiseless projection data, the parameters are set as λ1 = λ2 = 0.02, λ3 = 0.001 , λ4 = 0.002 ； For noisy projection data, the parameters are set as λ1 = λ2 = 0.025, λ3 = 0.001, λ4 = 0.

002.

8. The X-ray computed tomography reconstruction method according to claim 7, characterized in that: For PCB image reconstruction, the total number of algorithm iterations N iter is set to 15000.

9. The X-ray computed tomography reconstruction method according to claim 7, characterized in that: For workpiece image reconstruction, the total number of algorithm iterations N iter is set to 5000.

10. An X-ray computerized tomography reconstruction device, characterized in that, It is used to execute the method described in any one of claims 1-9.

Images