Dual spectral limited-angle CT imaging regularized by edge-preserving diffusion and smoothing

Abstract

Limited-angle CT scan is an effective way for nondestructive inspection of planar objects, and various methods have been proposed accordingly. When the scanned object contains high-absorption material, such as metal, existing methods may fail due to the beam hardening of X-rays. In order to overcome this problem, we adopt a dual spectral limited-angle CT scan and propose a corresponding image reconstruction algorithm, which takes the polychromatic property of the X-ray into consideration, makes basis material images free of beam hardening artifacts and metal artifacts, and then helps depress the limited-angle artifacts. Experimental results on both simulated PCB data and real data demonstrate the effectiveness of the proposed algorithm.

Keywords

Planar objects dual spectral limited-angle CT imaging the high-absorption material

1 Introduction

Computed Tomography (CT) is widely used in planar objects’ defect inspection, such as printed circuit boards (PCB), solar panels, slab fossils, etc. Due to the shape restriction of the scanned object, to avoid collision of the scanned object with the X-ray source, we can only acquire projections within a limited-angle range. Even if full-angle projections can be acquired, some may be useless due to the heavy absorption of the X-rays, which also results in limited-angle projections. Image reconstruction from the limited-angle projections is called limited-angle CT (LACT). It is well known that limited-angle projections do not satisfy Tuy’s condition [5]. Images reconstructed with conventional algorithms, such as filtered back projection (FBP) and algebraic iterative reconstruction (ART), usually contain blurring artifacts and streak artifacts along specific directions, which could be more severe if there are metals in the scanned object.

Due to the wide application requirements of LACT, various imaging models and solving algorithms have been proposed. Early methods [29, 32] attempt to incorporate the local or global smoothness prior of projections to recover the complete projections by applying certain “extrapolation” techniques to the limited-angle projections. Recently, many regularization methods have tried to incorporate various priors into the imaging model to eliminate limited-angle artifacts, including geometrical shape, distribution of edges, or densities (gray values) [13, 30] of the image. The commonly used regularization algorithms include total variation (TV) minimization [8, 30], l₀ gradient minimization [17, 27], alternating edge-preserving and smoothing (AEDS) method [12], directional total variation (DTV) algorithm [35], etc., all of which assume that the CT image is piecewise constant and thus has sparse gradients. The regularization reduces the solution space, thereby helping to eliminate image blurring and image aliasing artifacts caused by missing projections. In particular, AEDS and DTV algorithms, which incorporate both the piecewise constant prior to the ideal CT image and the correlation between the limited-angle artifacts and the scanning angular range, are state-of-the-art ones. All the above models, however, have assumed a linear relationship between the reconstructed image and the projections, which generally does not hold for a real CT system. Actually, the X-ray emitted from the source is polychromatic, and the X-ray absorption by the scanned object is energy-dependent. Therefore, images reconstructed with such algorithms using the above linear relationship usually suffer from beam hardening artifacts. The beam hardening artifacts make the reconstructed image not satisfy the piecewise constant assumption, so only suboptimal results can be obtained when the above algorithms are applied to real data. Moreover, if the scanned object contains metal such as copper and iron, the heavy beam hardening artifacts and metal artifacts can blur image edges, in which case the algorithms like AEDS that use edge information would fail to correctly recover blurred edges. In addition, beam hardening artifacts and metal artifacts may also introduce fake edges into the image and affect the final reconstructed results.

It is well known that dual spectral CT (DSCT) imaging can effectively eliminate beam hardening and metal artifacts. By acquiring projections of the scanned object at two different X-ray spectra, DSCT could reconstruct basis material images and synthesize a monochromatic image of the scanned object, which are theoretically free of beam hardening artifacts and metal artifacts [1 , 12].

There are several ways to acquire DSCT projections. The simplest way is to scan the object twice at different tube voltages and tube currents. This method is easy to implement but will increase the scan time and dose. Other typical data acquisition methods for DSCT include the dual-source dual-detector method [16, 23], fast kVp switching method [31, 34] and double-layer detector method [18, 22]. In recent years, the rapidly developing photon-counting detector technology has provided a new way to acquire polychromatic projections [14 , 28].

The key to DSCT imaging is the image reconstruction algorithm. The related study began in the 1970 s, and after over fifty years of exploration and research, scholars have invented a variety of reconstruction methods, which can be roughly classified into three categories: image domain methods [9], projection domain methods [19], and iterative methods [6, 11]. In image domain methods, two basis material images are approximated by linear combinations of the two images reconstructed separately from the two polychromatic DSCT projections. Such basis material images are usually regarded as a first-order approximation of the true images. Generally, beam hardening artifacts from reconstruction methods based on a monochromatic model [1, 11] propagate into basis material images. In projection domain methods, two sets of projections of basis material images are obtained from the two polychromatic projections. Then the basis material images are reconstructed from their projections with conventional image reconstruction methods, such as the FBP or back projection-filtering (BPF). Generally, projection domain methods can obtain better images than image domain methods. But projection domain methods require the DSCT data sets to be geometrically consistent. In other words, two separate flux measurements of the low and high X-ray spectra need to be acquired along each X-ray path [25]. However, this requirement is hard to meet for practical DSCT systems, such as dual-source CT systems. Therefore, several nonlinear optimization-based iterative methods for the DSCT were proposed in [11 , 33]. These methods fit both geometrically consistent and inconsistent projections. The main problem of iterative algorithms is a large amount of computation, but with the rapid development of computing hardware such as GPUs, iterative algorithms are receiving more and more attention, and some CT manufacturers have already applied them to their CT systems.

Recently, scholars have also begun to study dual spectral limited-angle CT (DSLACT) imaging problems. In 2016, Xing et al. [15] proposed a segmental multi-spectral CT algorithm, which divided the angular range of a circular scan trajectory into several arcs to achieve satisfactory limited-angle reconstructed results after adding a prior image. Then, an optimization-based algorithm relying on partial-scanned projections is proposed by Pan et al. [2, 4], which provides new inspiration for DSCT on conventional CT scanners. Later, Zhang et al. [26] tried to eliminate DSLACT artifacts with the guidance of image gradient. Sheng et al. [25] proposed a basis material decomposition method combined with sequence regularization to solve the DSLACT problem, which is mainly for medical images. In 2021 and 2022, Pan’s research group extended the DTV model [35] to DSLACT [3, 36], which achieves promising reconstructed results, especially for very small scanning angular ranges. But as is pointed out in [3], X-ray spectra are not considered in the model, and beam-hardening artifacts could not be eliminated completely.

This paper proposes an image reconstruction model and corresponding reconstruction algorithm called Dual Spectral Alternating Edge-preserving Diffusion and Smoothing (DS-AEDS) for DSLACT, which combines the extended algebraic iterative reconstruction (E-ART) method and the AEDS method. This algorithm could effectively eliminate beam hardening artifacts and metal artifacts in the reconstructed basis material images by making full use of the energy spectrum information, and suppress limited-angle artifacts by using the piecewise constant assumption of the basis material images. The effectiveness of the proposed model and algorithm are verified by numerical experiments on noise-free and noisy data. The beam hardening artifacts, metal artifacts and limited-angle artifacts are well suppressed.

2 DS-AEDS method for DSLACT imaging problem

This section first introduces the LACT imaging problem and AEDS algorithm, and then presents the DSCT imaging model and E-ART algorithm. Finally, the DS-AEDS model and algorithm are proposed for DSLACT imaging.

2.1 LACT imaging problem and AEDS algorithm

For simplicity, we consider the two-dimensional fan-beam LACT system. As illustrated in Fig. 1, Ox₁x₂ as the global coordinate system, the scanned object is located in origin. The fan beam scanning system consists of the X-ray source and the detector rotating around the origin. The symbols D and R represent the distance from the X-ray source to the detector and to the rotation center, respectively, d represents the detector coordinate, and β is the rotation angle of the X-ray source. Let f (x) denote the distribution of the scanned object’s linear attenuation, where x = (x₁, x₂) is the coordinate of a certain point in the scanned object, then fan-beam projections of f (x) can be expressed as:

Fig. 1

Illustration of the limited-angle fan-beam scanning configuration.

$p (d, β) = \int_{- \infty}^{+ \infty} f (R ν (β) + \frac{t}{\sqrt{D^{2} + d^{2}}} (- D ν (β) + ν^{⊥} (β))) dt$ (1) where ν (β) = (cosβ, sinβ), ν^⊥ (β) = (- sin β, cosβ). When $β \in [α, π - α] (0 < α < \frac{π}{2})$ , the formula (1) refers to the limited-angle scanning mode. Mathematically, LACT imaging is to solve f (x) from limited-angle projections.

Problem (1) can be discretized as follows $p_{j} = \sum_{i = 1}^{I} a_{ji} u_{i} + ɛ_{j}, j = 1, 2, \dots, J$ (2) $u_{i} = u_{i_{x} i_{y}}, i = (i_{x} - 1) N_{y} + i_{y}$ where {u_{i
_x
i
_y}|i_x = 1, 2, … , N_x, i_y = 1, 2, … , N_y } is the discrete two-dimensional digital image. The symbols N_x and N_y denote the number of rows and columns of the image, respectively, a_ji is defined as the intersection length between the jth X-ray and the ith pixel, J denotes the total number of rays for the limited-angle scan, and I = N_x × N_y denotes the total number of image pixels. In Equation (2) , we add an extra ɛ_j to denote the errors in p_j, including model error, discretization error, and measurement error. Thus, the LACT imaging problem could be boiled down to solving u from the following linear system: $A u = p + ɛ$ (3) where A = (a_ji) _J×I is the system matrix for the LACT scanning system, p = (p₁, … , p_j, … , p_J) ^T, ɛ = (ɛ₁, … , ɛ_j, … , ɛ_J) ^T.

Due to the strong ill-posedness of LACT imaging problem, severe blurring and streak artifacts along the ending lines of the limited-angular range and tangent to the boundary of some features in the object will be introduced when applying conventional reconstruction algorithms on limited-angle projections Fig. 2 illustrates this phenomenon, where (a) is a simulated planar object and (b) is the image reconstructed with the SART method from limited-angle projections (α = 30°). Obvious blurring could be seen in the reconstructed image, and streak artifacts appear along the directions of 30°and 150°.

Fig. 2

(a) The simulated planar object; (b) Reconstructed image with scanning angular range [30°, 150°].

AEDS is one of the effective algorithms for solving LACT imaging problems. Using edge information of the reconstructed image, AEDS recovers the inner information of the scanned object by edge-preserving diffusion and smoothing operation. The optimization model for AEDS reads, $u^{*} = arg min_{u} {{∥ A u - p ∥}_{2}^{2} + λ_{1} φ (u) + λ_{2} ψ (u)}$ (4) where ${∥ A u - p ∥}_{2}^{2}$ is the data fidelity term, φ (u) denotes l₀ norm of u_x, which is the directional derivative of the image u along the x direction, and ψ (u) denotes l₀ or l₁ norm of u_y, which is the directional derivative of the image u along the y direction. λ₁ and λ₂ are weighting parameters to adjust the regularization strength along x and y directions. The alternating direction algorithm is adopted for solving model (4) in [12].

2.2 DSCT imaging model and E-ART algorithm

Omitting scattered photons, the polychromatic projections could be discretized as: $p_{k, j} = - ln \sum_{m = 1}^{M_{k}} S_{k, m} δ e^{(- φ_{m} A_{j} f - θ_{m} A_{j} g)}, k = 1, 2, j \in J_{k}$ (5) where p_k,j is the polychromatic projection of the kth X-ray spectrum along the jth X-ray path, (S_k,1, … , S_{k,M_k}) is the discrete form of the kth normalized effective X-ray spectrum, and δ is the length of each subinterval. φ_m and θ_m are the sampling values of the mass attenuation coefficients of the two basis material in the kth subinterval. A_j = (a_j,1, a_j,2, … , a_j,I) is the projection vector for the jth ray, a_j,i is the contribution of the ith pixel from f = (f₁, f₂, ⋯ , f_I) ^T and g = (g₁, g₂, ⋯ , g_I) ^T to jth ray; f_i and g_i are the sampling values of f and g on the ith pixel, respectively. I = N_x × N_y is the total number of each basis material image pixels, N_x and N_y are the number of rows and columns of each basis image, respectively. $J_{k}$ indicates the index set of X-ray paths corresponding to the kth spectrum. The DSCT reconstructed problem is converted to solving f and g from Equation (5), which is a nonlinear system.

To solve Equation (5), Zhao et al. proposed the E-ART algorithm, which is an extension of the classic ART algorithm in 2015. It updates the basis material images by projecting the current estimations onto the approximate hyper plane of the nonlinear equation one by one. Let (f⁽ⁿ⁾ ; g⁽ⁿ⁾) be the currently estimated basis material images, the E-ART algorithm calculates the new iteration point (f⁽ⁿ⁺¹⁾ ; g⁽ⁿ⁺¹⁾) according to the polychromatic projection p_k,j as follows: $\begin{matrix} f^{(n + 1)} = f^{(n)} + \frac{Φ_{k, j}^{(n)}}{{(Φ_{k, j}^{(n)})}^{2} + {(Θ_{k, j}^{(n)})}^{2}} \frac{A_{j}^{T} q_{k, j}^{(n)} (p_{k, j} - p_{k, j}^{(n)})}{{∥ A_{j} ∥}_{2}^{2}} \\ g^{(n + 1)} = g^{(n)} + \frac{Θ_{k, j}^{(n)}}{{(Φ_{k, j}^{(n)})}^{2} + {(Θ_{k, j}^{(n)})}^{2}} \frac{A_{j}^{T} q_{k, j}^{(n)} (p_{k, j} - p_{k, j}^{(n)})}{{∥ A_{j} ∥}_{2}^{2}} \end{matrix}$ (6) where $p_{k, j}^{(n)}$ is the calculated polychromatic projection at the point (f⁽ⁿ⁾ ; g⁽ⁿ⁾) with Equation (5), and $\begin{matrix} q_{k, j}^{(n)} = \sum_{m = 1}^{M_{k}} S_{k, m} δ e^{- φ_{m} A_{j} f^{(n)} - θ_{m} A_{j} g^{(n)}} \\ p_{k, j}^{(n)} = - ln \sum_{m = 1}^{M_{k}} S_{k, m} δ e^{- φ_{m} A_{j} f^{(n)} - θ_{m} A_{j} g^{(n)}} \\ Φ_{k, j}^{(n)} = \sum_{m = 1}^{M_{k}} S_{k, m} δ φ_{m} e^{- φ_{m} A_{j} f^{(n)} - θ_{m} A_{j} g^{(n)}} \\ Θ_{k, j}^{(n)} = \sum_{m = 1}^{M_{k}} S_{k, m} δ θ_{m} e^{- φ_{m} A_{j} f^{(n)} - θ_{m} A_{j} g^{(n)}} \end{matrix}$ (7)

2.3 The proposed DS-AEDS model and algorithm

2.3.1 DS-AEDS model

In this section, we extend the AEDS algorithm to solve the DSLACT imaging problem and obtain our DS-AEDS algorithm. The proposed algorithm tries to remove metal and beam hardening artifacts, and simultaneously suppress the limited-angle artifacts. The DS-AEDS optimization model is as follows, $\begin{matrix} (f^{*}, g^{*}) = arg min_{(f, g)} {Γ (f, g)} \\ s . t . f ⊙ g = 0 \\ f ⩾ 0 \\ g ⩾ 0 \end{matrix}$ (8)

where $\begin{matrix} Γ (f, g) = {∥ P_{1} (f, g) - p_{1} ∥}_{2}^{2} + {∥ P_{2} (f, g) - p_{2} ∥}_{2}^{2} + λ_{1} φ_{1} (f) + λ_{2} ψ_{1} (f) \\ + λ_{3} φ_{2} (g) + λ_{4} ψ_{2} (g) \end{matrix}$ $\begin{matrix} φ_{1} (f) = {∥ \nabla_{x} f ∥}_{0} = \sum_{j = 1}^{I} sgn (| f_{j + 1} - f_{j} |) \\ ψ_{1} (f) = {∥ \nabla_{y} f ∥}_{0} = \sum_{j = 1}^{I} sgn (| f_{j + N_{x}} - f_{j} |) \\ φ_{2} (g) = {∥ \nabla_{x} g ∥}_{0} = \sum_{j = 1}^{I} sgn (| g_{j + 1} - g_{j} |) \\ ψ_{2} (g) = {∥ \nabla_{y} g ∥}_{0} = \sum_{j = 1}^{I} sgn (| g_{j + N_{x}} - g_{j} |) \end{matrix}$ (9)

${∥ P_{1} (f, g) - p_{1} ∥}_{2}^{2}$ and ${∥ P_{2} (f, g) - p_{2} ∥}_{2}^{2}$ are data fidelity terms for the acquired high- and low-energy projections p₁ and p₂, respectively. P₁ (f, g) and P₂ (f, g) are the calculated high- and low-energy projections at point (f, g) according to Equation (5). The operators φ_i (·) and ψ_i (·) , i = 1, 2 are the x direction and y direction regularization terms for f and g, respectively, which are of the same meaning as those in the AEDS model. The parameters λ_i (i = 1, 2, 3, 4) are the weighting parameters for regularization terms, controlling the smoothing strength. The f ⊙ g in the constraints denotes the element-wise product of the two basis material images f and g. f ⊙ g = 0 indicates that the scanned object contains only one basis material at each pixel position.

2.3.2 DS-AEDS algorithm

Due to the nonlinearity of DSCT, model (8) is non-convex, only local optimal solutions can be obtained. Here, we adopt a simple alternating minimization technique to solve (8). Let (f⁽ⁿ⁾, g⁽ⁿ⁾) denote the solution after the nth iteration, then (f⁽ⁿ⁺¹⁾, g⁽ⁿ⁺¹⁾) could be obtained by solving the following three subproblems sequentially:

Subproblem 1: $\begin{matrix} (f^{(n + \frac{1}{3})}, g^{(n + \frac{1}{3})}) = arg min_{(f, g)} {\begin{matrix} {∥ f^{(n)} - f ∥}_{2}^{2} + {∥ P_{1}^{(n)} (f, g) - p_{1} ∥}_{2}^{2} + \\ {∥ g^{(n)} - g ∥}_{2}^{2} + {∥ P_{2}^{(n)} (f, g) - p_{2} ∥}_{2}^{2} \end{matrix}} \\ s . t . f ⊙ g = 0 \\ f ⩾ 0 \\ g ⩾ 0 \end{matrix}$ (10)

Subproblem 2: $(f^{(n + \frac{2}{3})}, g^{(n + \frac{2}{3})}) = arg min_{(f, g)} {{∥ f^{(n + \frac{1}{3})} - f ∥}_{2}^{2} + λ_{1} φ_{1} (f) + {∥ g^{(n + \frac{1}{3})} - g ∥}_{2}^{2} + λ_{3} φ_{2} (g)}$ (11)

Subproblem 3: $(f^{(n + 1)}, g^{(n + 1)}) = arg min_{(f, g)} {{∥ f^{(n + \frac{2}{3})} - f ∥}_{2}^{2} + λ_{2} ψ_{1} (f) + {∥ g^{(n + \frac{2}{3})} - g ∥}_{2}^{2} + λ_{4} ψ_{2} (g)}$ (12)

The E-ART algorithm, followed by a thresholding process, is used to solve subproblem 1, which is also an iterative process. In our implementation, we find that one E-ART iteration in each inner iteration is enough to solve our problem. After each E-ART iteration, the negative pixel values are replaced with 0 to satisfy the nonnegative constraint. At last, f ⊙ g = 0 is achieved by setting parts of one basis material image’s pixel values to 0. Suppose that f and g are the basis material images with high and low attenuation, respectively. The specific steps for f ⊙ g = 0 are as follows: firstly, choose a threshold T for f, and the pixels whose values are greater than the selected threshold T are identified as high attenuation material; then, set the corresponding pixels in g to 0. Subproblems 2 and 3, which are quite similar in form, are solved by the alternating direction multipliers method (ADMM) [7]. In addition, since f and g are separable, they could be solved independently and in parallel. But please note that the roles of subproblems 2 and 3 are different, the former recovers image edges and eliminates limited-angle artifacts on the basis material images by edge-preserving diffusion along the x direction, while the latter removes the possible image artifacts caused by the former operations by edge-preserving smoothing along the y direction.

It should be noted that the selection of edge-preserving diffusion operators φ_i (·) (i = 1, 2) and edge-preserving smoothing operators ψ_i (·) (i = 1, 2) affect the reconstructions greatly. For monochromatic imaging, as pointed out in [12], a suitable choice for horizontal edge-preserving diffusion operators φ_i (·) (i = 1, 2) is l₀ norm of image gradient, thanks to its edge-preserving ability. On the other hand, the choice of vertical direction edge-preserving smoothing operators ψ_i (·) (i = 1, 2) depends on the noise level. In our experiment, we adopt the same regularization terms as [12] suggests. Define, the workflow of the DS-AEDS algorithm can be concluded as Algorithm 1.

2.3.4 Parameter discussion

The algorithm parameters of DS-AEDS include λ₁, λ₂, λ₃, λ₄, T, M, N, K and ω. λ₁ and λ₃ control the strength of x direction edge-preserving diffusion for the two different basis materials. And λ₂ and λ₄control the strength of y direction edge-preserving smoothing for the two different basis materials. Forλ_i (i = 1, 2, 3, 4), the larger the value, the stronger the strength. How to determine the threshold T depends on the density of the different materials and the noise level. For example, if the densities of high and low attenuation materials are in [7, 8] and [0.9, 1.2], we recommend that T could be chosen from [2, 6]. The number N is the total iterations, which is determined by the sample structure and noise level. The more complex the sample and the higher the noise level, the more iterations are required. The number M is the iterations to start using the median filter, which is used every K time. M, N and K jointly determine the number of times the median filter is used. ω is the relaxation factor, determined by the noise level of projections. The greater the noise level, the smaller the ω value.

Algorithm 1: DS-AEDS
Input: polychromatic projections p₁, p₂, maximum iteration number λ₁, λ₂, λ₃, λ₄, T, M, N, K
and ω
Initialize: f⁽⁰) = ^′′0, g^(0)′ = ^′′0
Output: f^(N), g^(N), t
for n = 0 : N–1do
Step 1:
$(f^{(n + \frac{1}{3})}, g^{(n + \frac{1}{3})}) = E - ART (f^{(n)}, g^{(n)}, p_{1}, p_{2})$
Let $f^{(n + \frac{1}{3})}$ , guide $g^{(n + \frac{1}{3})}$
$g^{(n + \frac{1}{3})} = {\begin{matrix} 0 & f^{(n + \frac{1}{3})} > T \\ g^{(n + \frac{1}{3})} & e l s e \end{matrix}$
$(f^{(n + \frac{1}{3})}, g^{(n + \frac{1}{3})}) = pos (f^{(n + \frac{1}{3})}, g^{(n + \frac{1}{3})})$
Step 2:
Solving optimization problem:
$(f^{(n + \frac{2}{3})}, g^{(n + \frac{2}{3})}) = arg min_{(f, g)} {{∥ f^{(n + \frac{1}{3})} - f ∥}_{2}^{2} + λ_{1} φ_{1} (f) + {∥ g^{(n + \frac{1}{3})} - g ∥}_{2}^{2} + λ_{3} φ_{2} (g)}$
Step 3:
Solving optimization problem:
$(f^{(n + 1)}, g^{(n + 1)}) = arg min_{(f, g)} {{∥ f^{(n + \frac{2}{3})} - f ∥}_{2}^{2} + λ_{2} ψ_{1} (f) + {∥ g^{(n + \frac{2}{3})} - g ∥}_{2}^{2} + λ_{4} ψ_{2} (g)}$
ifn > M and n % K = =0:
f⁽ⁿ⁺¹⁾ = median (f⁽ⁿ⁺¹⁾)along y direction
g⁽ⁿ⁺¹⁾ = median (g⁽ⁿ⁺¹⁾) along y direction
end if
end for
synthesize f^(N) and g^(N) to a monochromatic image t = s₁f^(N) + s₂g^(N)
where s₁ and s₂ are f and g corresponding mass attenuation coefficients at fixed energy.

3 Experiments

3.1 Simulated experiment

In this section, experiments on a simulated PCB phantom will be carried out to validate the effectiveness of the proposed DS-AEDS algorithm. For simplicity, the numerical experiments are restricted to fan beam CT. The generalization to cone-beam CT is straightforward.

3.1.1 Experiment parameters and conditions

PCB inspection is a typical LACT imaging problem, so a simulated PCB phantom, as shown in Fig. 3, is used in our experiments. The phantom is composed of basis material polyvinyl chloride (PVC) and metal material (copper). The size of the phantom is 12.4mm*2.25 mm. Four different types of the copper route are included in the phantom, whose thicknesses are 0.24 mm, 0.015 mm, 0.03 mm and 0.045 mm, respectively. Due to the existence of high-absorption metal material in the phantom, conventional reconstruction algorithms will introduce severe beam hardening artifacts and metal artifacts in the reconstructed CT image. Furthermore, if an edge information-dependent algorithm is employed, the artifacts may affect the visibility of certain edges in the reconstructed image, making it difficult to exactly recover the inner information of the scanned object. Moreover, beam hardening artifacts and metal artifacts might introduce fake edges to the reconstructed image, affecting the interpretation of correct image.

Fig. 3

The simulated PCB phantom.

Table 1 shows the parameters of the scanning configuration. The resolution of the reconstructed image is 1024*370, with a pixel size of 0.015 mm.

Table 1

The scanning configuration parameters for simulated data

Parameter	Value
Source to object distance	239.046mm
Source to detector distance	697.905mm
Scanning angular range	[30°, 150°]
Scanning angular interval	0 . 25°
Projection views	480
Detector cells	1920
Detector cell size	0.127mm

The polychromatic X-ray spectra are simulated with the software SpectrumGUI, an open-source X-ray spectra simulator (see http://spectrumgui.sourceforge.net/’). We have simulated the X-ray spectra of GE Maxiray 125 X-ray tube with tube voltages 80 and 140 kV, where the latter is filtered with 1 mm copper. The distributions of the X-ray spectra are shown in Fig. 4. The mass attenuation coefficients used in the experiments are retrieved from the National Institute of Standard Technology (NIST) tables of the X-ray mass attenuation coefficient. The length of the sampling interval is 1 keV in the simulation of polychromatic projections.

Fig. 4

The X-ray spectra used in our experiments.

The proposed DS-AEDS is implemented by python and CUDA. All the experiments are performed on the Windows Server 2016 server (NVIDIA Tesla V100 S card with 32 G video memory).

The selection of parameters is critical in each of the above algorithms. The only parameter involved in the SART algorithm is its relaxation parameter which is set to be 1 in our simulated experiments. The parameters of AEDS [12], DTV [3] and Newton+AEDS for noise-free and noisy data are shown in Tables 2, 3 and 4, respectively, where we have used the same notations in the references.

Table 2

Reconstruction parameters used by AEDS algorithm on simulated data

Parameters	λ ₁	λ ₂	K	ω
noise-free	5 × 10^-5	4 × 10^-5	28000	1.0
Poisson noise	5 × 10^-5	4 × 10^-5	28000	0.9

Table 3

Reconstruction parameters used by DTV algorithm on simulated data

Parameters	t _x	t _y	b
noise-free	8000	3125	0.004
Poisson noise	8000	3000	0.004

Table 4

Reconstruction parameters used by the Newton+AEDS algorithm on simulated data

	PVC		Copper
Parameters	λ ₁	λ ₂	λ ₁	λ ₂
noise-free	2 × 10^-4	8 × 10^-5	1 × 10^-3	6 × 10^-5
Poisson noise	2 × 10^-2	1 × 10^-2	1 × 10^-2	1 × 10^-3

As for DS-AEDS, the threshold T for f ⊙ g = 0 is set to 6, according to the basis material density. Our algorithm runs 28000 iterations, in which the median filter is performed once every 250 iterations after 4000 iterations of the DS-AEDS algorithm to accelerate the y direction convergence for the PVC reconstructed result. Because the noise level of the simulated data is relatively low, it is only necessary to apply median filtering on the PVC board to play a role of acceleration. The parameters of our proposed algorithm are listed in Table 5.

Table 5

Reconstruction parameters used by our proposed algorithm on simulated data

Parameters	λ ₁	λ ₂	λ ₃	λ ₄	K	T	M	N	ω
noise-free	2 × 10^-4	8 × 10^-5	4 × 10^-1	5 × 10^-6	28000	6	4000	250	1.0
Poisson noise	2 × 10^-4	7 × 10^-5	4 × 10^-1	5 × 10^-6	28000	6	4000	250	1.0

3.1.2 Experiment results

Due to the polychromatic nature of X-rays as well as the high abortion of the scanned sample, the image reconstructed with SART will be degraded by beam hardening artifacts, as shown in Fig. 5, where Fig. 5(a) is the reconstructed result with SART from 100kVp full-angle polychromatic projections. It can be seen that the gray values at the border parts of the metal are greater than those of the inner part. In addition, the middle part of the thin iron sheet pointed by the red arrow is not visible in the reconstructed result due to beam hardening artifacts. Obviously, there are beam hardening artifacts in the reconstructed result. Figure 5(b) and (c) show the profiles along the red lines in Fig. 5(a), where typical cupping artifacts can be clearly seen.

Fig. 5

(a) The reconstructed result with SART from 100kVp full-angle polychromatic projections; (b) The profile along the top red line in (a); (c) The profile along the bottom red line in (a).

Reconstructed results of the proposed DS-AEDS algorithm and comparative algorithms are presented in this section. Figure 6 and 7 are reconstructed results from noise-free and noisy projections. For noisy data, the noise is the Poisson noise of 10⁷ photons per ray. In both Fig. 6 and 7 , subfigures (a), (b) and (c) are reconstructed results by SART, AEDS [12] and DTV [3] from 100kVp polychromatic projections. Subfigures (d) and (g) are 100 keV monochromatic images reconstructed with Newton+AEDS and the proposed DS-AEDS from 80kVp and 140kVp polychromatic projections. The former is synthesized from the basis material images in subfigures (e) and (f), while the latter is synthesized from the basis material images in subfigures (h) and (i). The Newton+AEDS algorithm is a combination of Newton and AEDS algorithms, which consists of the following two steps: 1) Decompose the limited-angle polychromatic projections directly in the projection domain to obtain line integrals of basis material using the Newton method. 2) Apply the AEDS algorithm on obtained line integrals to reconstruct basis material images. The quantitative measurements including NRMSE, PSNR and SSIM for the reconstructed results are shown in Table 6.

Fig. 6

The reconstructed results from noise-free projections. (a), (b) and (c) are the results of SART, AEDS and DTV from 100kVp polychromatic projections; the display window of (a)-(c) is [– 0.2, 1.8]; (d) is the synthesized 100keV monochromatic image of Newton+AEDS algorithm, whose basis material images are (e) and (f); (g) is DS-AEDS algorithm’s synthesized 100keV monochromatic result, the corresponding basis material images are (h) and (i). For (d) and (g), the display window is [0.015, 0.045]; the display window of (e), (f), (h) and (i) is [– 0.2, 1.8]].

Fig. 7

The reconstructed results from noisy projections, other descriptions and display window levels are the same as those in Fig. 6.

It can be seen from Fig. 6(a)that both beam hardening artifacts and limited-angle artifacts are so severe that the gray value distribution of the metal region in the reconstructed image is not uniform and horizontal image edges are blurred. In addition, due to the effects of the two types of artifacts, the copper routes of 0.015 mm, 0.03 mm and 0.045 mm are completely invisible in the reconstructed result. As shown in Fig. 6(b), AEDS has alleviated the limited-angle artifacts, and several blurred edges are recovered. But, as the polychromatic characteristic of the projections is not considered in AEDS, the model fails to reconstruct all the copper routes, and fake edges are even introduced. The DTV reconstructed result in Fig. 6(c)cannot clearly show the thin copper routes, and the corners indicated by the arrows are not very sharp, due to the effect of beam hardening. Newton+AEDS and the proposed DS-AEDS obtain satisfactory reconstructed results as shown in Fig. 6(d) and (g), where no obvious visual differences exist. But from Table 6, we can see that the quantitative measurements of the results of the proposed algorithm are better than those of Newton+AEDS.

Table 6

The NRMSE, PSNR and SSIM of the reconstructed results are shown in Fig. 6 and 7

	noise-free			Poisson noise
	NRMSE	PSNR	SSIM	NRMSE	PSNR	SSIM
SART	2.7538	11.2621	0.7069	2.7576	11.2500	0.6862
AEDS	3.8101	8.4418	0.6740	3.8849	8.273	0.6896
DTV	3.2973	9.6976	0.6718	3.2826	9.7363	0.6748
Newton+AEDS	0.0393	48.1702	0.9995	0.0746	42.6082	0.9963
DS-AEDS	0.0063	64.060	0.9999	0.0146	56.7528	0.9996

Reconstructed results of the proposed DS-AEDS algorithm and comparative algorithms from noisy data are shown in Fig. 7, from which similar conclusions could be drawn. Similar to noise-free cases, SART, AEDS and DTV still fail to recover the thin copper routes and sharp corners of thick metal. Moreover, due to the high sensitivity to noise, the Newton+AEDS method cannot obtain accurate line integrals of basis material, fake edges are produced by the AEDS method as indicated by the red arrows in Fig. 7(d). But from Fig. 7(g), we can see that our proposed method is robust to noise, limited-angle artifacts and beam hardening artifacts are both well suppressed and all the copper routes are correctly reconstructed.

Figure 8 shows the profiles of the reconstructed results by different algorithms of noise-free and noisy cases. The obvious cupping shape can still be seen in the line profiles of SART, AEDS and DTV’s results, indicating that these methods could not effectively suppress beam hardening artifacts. This is in accordance with our intuition, since they do not take the polychromatic property of the X-ray into consideration. The black dotted lines in Fig. 8 show the results of DS-AEDS. It can be seen that the results of the DS-AEDS match the references very well, even for noisy case, indicating that beam hardening artifacts are effectively eliminated.

Fig. 8

The profiles of the reconstructed results by different algorithms of noise-free and noisy cases. (a) and (b) are the profiles corresponding to noise-free results in Fig. 6(a)(b)(c)(d)(g) along the lines in Fig. 5(a), respectively; (c) and (d) are the profiles corresponding to noisy results in Fig. 7(a)(b)(c)(d)(g) along the lines in Fig. 5(a), respectively.

3.2 Real data experiment

The real data experiments are performed on the industrial CT system of our laboratory, as is shown in Fig. 9(a), which is equipped with a YXLON Y.TU-450-D11 (YXLON GmbH, Hamburg, Germany) X-ray source and an IRay NDT 1616HE detector (IRay Technology, Shanghai, China).

Fig. 9

(a) shows the CT system. (b) is the scanned sample, which is sandwich-shaped, composed of alternately stacked aqua green FR4 plates and thin iron sheets.

3.2.1 Experiment parameters and conditions

The sample is sandwich-shaped, composed of alternately stacked Flame Retardant 4 (FR4) thin plates and thin iron sheets to simulate the circuit board, as shown in Fig. 9(b). The density of FR4, a common circuit board material, is about 1.02g/cm³, and the density of iron is about 7.85g/cm³. The overall length and width of the sample are 80mm and 16.5mm. The thicknesses of FR4 board and iron sheets are shown in Fig. 10. Three sets of scan data are acquired at voltages 250kV, 350kV and 300kV, where the first two are used in Newton+AEDS and DS-AEDS, while the last one is used in SART, DTV and AEDS. The parameters of the scanning configuration are shown in Table 7. The resolution of the reconstructed image is 1024 * 300, with a pixel size of 0.1mm. The parameters for the real data experiments are shown in Tables 8, 9 and 10. The relaxation factor for the SART algorithm is 0.2. For the sake of fairness, AEDS, DTV and DS-AEDS all use 12000 iterations. The number of iterations is smaller than that used in the simulated data because the thinnest iron routes reconstructed in the real data is 2 pixels theoretically (1 pixel in simulated data) and the sample’s physical structure is relatively simple.

Fig. 10

The top view of the scanned sample (mm).

Table 7

The scanning configuration parameters for real data

Parameter	Value
Source to object distance	615mm
Source to detector distance	1205mm
Scanning angular range	[30°, 150°]
Scanning angular interval	0 . 25°
Projection views	480
Detector cells	1024
Detector cell size	0.2mm

Table 8

Reconstruction parameters used by the AEDS algorithm on real data

Parameters	λ ₁	λ ₂	K	ω
	5 × 10^-5	4 × 10^-5	28000	1.0

Table 9

Table9

Parameters	t _x	t _y	b
	2000	100	0.01

Table 10

Reconstruction parameters used by our proposed algorithm on real data

Parameters	λ ₁	λ ₂	λ ₃	λ ₄	K	T	M	N	ω
	2 × 10^-2	1 × 10^-5	5 × 10^-5	1 × 10^-6	12000	1.2	4000	100	0.05

3.2.2 Experiment results

The reconstructed results of real data are shown in Fig. 11 and Fig. 12. Figure 11(a) is the reconstructed result with SART from full-angle projections, Fig. 11(b) and (c)are line profiles along the red lines in Fig. 11(a). From the reconstructed results as well as the line profiles, beam hardening artifacts can be clearly seen. In Fig. 12(a)-(c), are the results of SART, AEDS and DTV, (d) is the 300 keV monochromatic image synthesized from basis material images (e) and (f) that are reconstructed with DS-AEDS. The noise level of the real data is much higher than that of noisy simulated data, so Newton method cannot converge in our experiments. We did not put the Newton method experimental results in Fig. 12.

Fig. 11

(a) The reconstructed result of real data with SART from 300kVp full-angle polychromatic projections; (b) The profile along the top red line in (a); (c) The profile along the bottom red line in (a).

Fig. 12

The reconstructed results from real data. (a) is the result of SART from 300kVp polychromatic projections after 10 iterations; (b) is the result of the AEDS algorithm from 300kVp polychromatic projections; (c) is the result of the DTV algorithm from 300kVp polychromatic projections; (d) is DS-AEDS algorithm’s synthesized 300keV monochromatic result. The display window of (a)-(d) is [0, 0.3]; the corresponding basis material images of (d) are (e) and (f); the display window of (e) is [[0.635, 1.365] and the display window of (f) is [0, 4.0].

From Fig. 12(a), we can see that the limited-angle artifacts and beam hardening artifacts are intertwined in the reconstructed results by SART. For AEDS and DTV, both of them recover the main structures of the sample as shown in Fig. 12(b) and (c). But the beam hardening artifacts indicated by the red arrow cannot be eliminated, and the thin iron sheet indicated by the yellow arrow is not accurately reconstructed. AEDS and DTV with other parameters either fail to remove limited-angle artifacts or introduce fake horizontal edges. From Fig. 12(d), (e) and (f), we can see that our proposed algorithm can effectively remove both hardening artifacts and limited-angle artifacts, and can accurately reconstruct thin iron sheets.

Figure 13 shows the profiles of the reconstruction results of different algorithms along the red lines in Fig. 11(a). From the profiles, a similar conclusion with the simulated experimental results can be drawn, i.e., the proposed method eliminated beam hardening artifacts very well, while other comparison methods failed to suppress beam hardening artifacts as the polychromatic property of X-rays is not considered.

Fig. 13

The profiles of the reconstructed results of real data by different algorithms. (a) and (b) are profiles corresponding to Fig. 12(a)(b)(c)(d). (a)’s profile position is along the top red line in Fig. 11 (a); (b)’s profile position is along the bottom red line in Fig. 11 (a).

4 Conclusion

This paper proposes the DS-AEDS model and algorithm for DSLACT imaging, which simultaneously suppresses limited-angle artifacts and beam hardening artifacts effectively. Preliminary experiments on simulated PCB data verify its validity. It is important to note that although only one PCB phantom’s reconstructed results are illustrated in this paper due to space limitation, multiple different PCB phantoms are tested and positive results are obtained.

In our model, a constraint f ⊙ g = 0 is used, which seems to be too strong for medical imaging. But this constraint usually holds for the industrial nondestructive test, especially for workpieces whose components are generally simple, such as PCB, solar panels, etc. DS-AEDS seems slow since tens of thousands of iterations are needed. We will investigate its acceleration techniques in our future work.

Footnotes

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant 62271330, Grant 61971292, Grant 61871275, Grant 61827809, Grant 61771324; National Key Research and Development Program of China under Grant 2020YFA0712200; and Major Technologies R & D Program of Shenzhen (JSGGZD20220822095600001). The authors would like to thank Professor Peng Zhang and Professor Hongwei Li, for the instructions on constructive suggestions on the modification of the manuscript.

References

Coleman

A.J.

, Sinclair

A beam-hardening correction using dual-energy computed tomography, Physics in Medicine and Biology 30(11) (1985), 1251–1256.

Chen

, Zhang

, Xia

, Sidky

E.Y.

, Pan

X.C.

Algorithm-enabled partial-angular-scan configurations for dual-energy CT, Medical Physics 45(5) (2018), 1857–1870.

Chen

, Zhang

, Xia

, Sidky

E.Y.

, Pan

Dual-energy CT imaging with limited-angular-range data, Physics in Medicine and Biology 66 (2021), 185020.

Chen

, Zhang

, Sidky

E.Y.

, Xia

, Pan

Image reconstruction and scan configurations enabled by optimization based algorithms in multi-spectral CT, Physics in Medicine and Biology 62 (2017), 8763–8793.

Smith

B.D.

Image reconstruction from cone-beam projections: necessary and sufficient conditions and reconstruction methods, IEEE Transactions on Medical Imaging 4(1) (1985), 14–25.

Man, J. Nuyts and P. Dupont, 23 polychromatic algorithm for CT, IEEE Transactions on Medical Imaging 20(10) (2001), 999–1008.

Boyd

, Parikh

, Chu

et al Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers, Foundations & Trends in Machine Learning 3(1) (2010), 1–122.

Sidky

E.Y.

, Pan

Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in Medicine and Biology 53(17) (2008), 4777–4807.

Christ

Exact treatment of the dual-energy method in CT using poly energetic X-ray spectra, Physics in Medicine and Biology 29(12) (1984), 1501–1510.

10.

Park

H.S.

, Hwang

, Seo

J.K.

Metal artifact reduction for polychromatic X-ray CT based on a beam-hardening corrector, IEEE Transactions on Medical Imaging 35(2) (2016), 480–487.

11.

Elbakri

I.A.

, Fessler

J.A.

Statistical image reconstruction for poly energetic X-ray computed tomography, IEEE Transactions on Medical Imaging 21(2) (2002), 89–99.

12.

, Zhao

, Li

, Zhang

An Image Reconstruction Model Regularized by Edge-preserving Diffusion and Smoothing for Limited-angle Computed Tomography, Inverse Problem 35(8) (2019), 085004.

13.

Batenburg

K.J.

, Sijbers DART: A practical reconstruction algorithm for discrete tomography, IEEE Transactions On Image Processing 20(9) (2011), 2542–2553.

14.

Taguchi

, Iwanczyk

J.S.

Vision 20/20: Single photon counting X-ray detectors in medical imaging, Medical Physics 40(10) (2013), 100901.

15.

Shen

, Xing

Multienergy CT acquisition and reconstruction with a stepped tube potential scan, Medical Physics 42 (2014), 282–296.

16.

Petersilka

, Bruder

, Krauss

Technical principles of dual source CT, European Journal of Radiology 68(3) (2008), 362–368.

17.

Storath

, Weinmann

, Frikel

, Unser

Joint image reconstruction and segmentation using the potts model, Inverse Problems 31(2), 2015.

18.

Allec

, Abbaszadeh

, Karim

K.S.

Single-layer and dual layer contrast-enhanced mammography using amorphous selenium flat panel detectors, Physics in Medicine and Biology 56(18) (2011), 5903–5923.

19.

Stenner

, Berkus

, Kachelriess

Empirical dual energy calibration (EDEC) for cone-beam computed tomography, Medical physics 34(9) (2007), 3630–3641.

20.

, Mou

, Tang

Implementation of penalized likelihood statistical reconstruction for polychromatic dual energy CT[A], SPIE Proceedings: Physics of Medical Imaging. San Diego, CA: SPIE (2009), 1–9.

21.

Brooks

R.A.

, Chiro

G.D.

Beam hardening in X-ray reconstructive tomography, Physics in Medicine and Biology 21(3) (1976), 390–398.

22.

Carmi

, Naveh

, Altman

Material separation with dual-layer CT, IEEE Nuclear Science Symposium Conference Record. USA: IEEE (2005), 1876–1878.

23.

Flohr

T.G.

, Mccollough

C.H.

, Bruder

First performance evaluation of a dual-source CT (DSCT) system, European Radiology 16(2) (2006), 256–268.

24.

Barber

W.C.

, Nygarda

, Iwanczyka

J.S.

et al. Characterization of a novel photon counting detector for clinical CT: count rate, energy resolution, and noise performance, Proc SPIE 7258 (2009), 725824–1.

25.

Sheng

W.J.

, Zhao

, Li

M.F.

A sequential regularization based image reconstruction method for limited-angle spectral CT, Physics in Medicine and Biology (2020), 65.

26.

Zhang

, Wang

, Li

, Liang

, Yan

, Hu

Reconstruction method for DECT with one half-scan plus a second limited-angle scan using prior knowledge of complementary support set (Pri-CSS), Physics in Medicine and Biology 2020.

27.

, Zeng

l0 gradient minimization based image reconstruction for limited-angle computed tomography, PloS One 10(7) (2015), e0130793.

28.

Llopart

, Campbell

, Dinapoli

et al. Medipix2: a 64 -k pixel readout chip with 55-μm square elements working in single photon counting mode, IEEE Transactions on Nuclear Science 49(5) (2002), 2279–2283.

29.

Mou

, Tang

, Yu

Comparison on beam hardening correctionof CT based on HL consistency and normal water phantom experiment,In1V, Proc SPIE 6318 (2006), 63181V.

30.

Zhuge

, Palenstijn

W.J.

, Batenburg

K.J.

TVR-DART: A morerobust algorithm for discrete tomography from limited projectiondata with automated gray value estimation, IEEE Transactions OnImage Processing 25(1) (2016), 455–468.

31.

, Langan

D.A.

, Wu

Dual energy CT via fast kVp switching spectrum estimation. SPIE Proceedings: Physics of Medical Imaging, San Diego, CA: SPIE (2009), 1–10.

32.

Xia

, Berger

, Bauer

, Hu

S.Y.

, Aichert

, Maier

An improved extrapolation scheme for truncated CT data using 2D Fourier-based Helgason-Ludwig consistency conditions, International Journal of Biomedical Imaging 2017 (2017), 1867025.

33.

Zhao

, Zhao

, Zhang

An Extended Algebraic Reconstruction Technique (E-ART) for Dual Spectral CT, IEEE Transactions on Medical Imaging 34(3) (2015), 761–768.

34.

Zou

, Silver

M.D.

Analysis of fast kV-switching in dual energy CT using a pre-reconstruction decomposition technique, SPIE Proceedings: Physics of Medical Imaging. San Diego, CA: SPIE, (2008), 1493–1495.

35.

Zhang

, Chen

, Xia

, Sidky

E.Y.

, Pan

Directional-TV algorithm for image 411 reconstruction from limited angular-range data, Medical Image Analysis 2021.

36.

Zhang

, Chen

, Xia

, Sidky

E.Y.

, Pan

Image reconstruction from data over two orthogonal arcs of limited-angular ranges, Medical Physics 49 (2022), 1468–1480.