Compressed sensing and Otsu’s method based binary CT image reconstruction technique in non-destructive detection

NUCLEAR ELECTRONICS AND INSTRUMENTATION

Compressed sensing and Otsu’s method based binary CT image reconstruction technique in non-destructive detection

REN Yong，

HE Peng，

WANG Hong-Liang，

CEN Zhong-Jie，

FENG Peng，

WEI Biao

Nuclear Science and Techniques

Vol.26, No.5

Article number 050403

Published in print 20 Oct 2015

Available online 20 Oct 2015

DOI：10.13538/j.1001-8042/nst.26.050403

37301

This paper tries to address the problem of binary CT image reconstruction in non-destructive detection with an algorithm based on compressed sensing (CS) and Otsu’s method, which could reconstruct binary CT image of test object from incomplete detection data. According to binary CT image characteristics, we employ Split-bregman method based on L_1/2 regularization to solve piecewise constant region reconstruction. To improve the reconstructed image quality from incomplete detection data, we utilize a priori knowledge and Otsu’s method as the optimization constraint. In our study, we make numerical simulation to investigate our proposed method, and compare reconstructed results from different reconstruction methods. Finally, the experimental results demonstrate that the proposed method could effectively reduce noise and suppress artifacts, and reconstruct high-quality binary image from incomplete detection data.

Non-destructive detectionComputed tomography (CT)Binary image reconstructionCompressed sensing (CS)

I. INTRODUCTION

X-ray computed tomography (CT), an technique for clinic diagnostics at first, has been developed into industrial CT system for non-destructive detection in aerospace, geology, weapons, metallurgy, etc. [1-5]. Test objects consisting of single composition materials, such as engine, rock specimen, teeth, etc. are used in CT imaging [6, 7], so CT images of these test objects can be considered as binary images that have only two gray values (i.e., black and white). They can be modeled as piecewise constant matrix, and easily sparsified by proper orthogonal transformation in reconstruction [8, 9]. Binary CT image reconstruction is a key technique, for image reconstruction from incomplete projection data by some continuous tomography methods [10, 11]. Also, there are some discrete tomography methods for binary CT image reconstruction from few-view projections [12-15]. In industrial CT systems, scan time is very long due to large size of test object and too many projection views [4, 5]. Test object reconstruction from few-view projections can reduce scan time. In medical CT systems [8], biomedical samples reconstruction from few-view projections is conducive for reducing the radiation dose.

Conventional reconstruction algorithms (i.e., filter back-projection and algebraic reconstruction techniques) cannot reconstruct high-quality CT image from incomplete projection data [1]. Interestingly, compressed sensing (CS) theory can reconstruct high-quality images from less projection data than what is usually considered necessary according to the Nyquist sampling theory [16-18]. The main idea of CS theory is that an image can be reconstructed from a rather limited amount of data as long as an underlying image can be sparsely represented in an appropriate domain and determined from these data [19-21]. Also, some advanced algorithms make use of some priori information in CT image reconstruction to reduce projection data greatly [22, 23]. Thus, it is feasible to reconstruct binary CT image from incomplete projection data using the CS-based reconstruction algorithm and some priori information.

In this paper, we focus on binary CT image reconstruction from incomplete projection data in non-destructive detection, and propose an algorithm based on CS and Otsu’s method. To improve quality of reconstructed image, we employ Split-bregman method with L_1/2 regularization [24], which can produce sparser solution compared with the L₁ regularization method used in CT image reconstruction. Meanwhile, we also use a priori knowledge of gray value information and Otsu’s method [25] as the optimization constraint in reconstruction process, which can segment and extract gray value information of binary image. The remainder of this paper is organized as follows. Our reconstruction algorithm is proposed in Sec. II, images reconstructed using different methods are compared in Sec. III, and in the last section, the image reconstruction process and results are discussed.

II. MATERIALS AND METHODS

In CT imaging, it is of significance to reduce scan time or radiation dose, which, in turn, determine the detection data completeness. Currently, the CS-based reconstruction algorithm is easier to take into account incomplete projection data during the reconstruction process and make a better performance. We propose here an algorithm based on CS theory and Otsu’s method to reconstruct binary CT images from incomplete detection data in non-destructive detection. In the following, we will summarize the proposed algorithm scheme.

It has been widely accepted that a CT image reconstruction can be modeled as a linear system

A f = b,

(1)

where b = (b¹,, bM) RM represents the detection data, f = (f₁,,fM) RN denotes the image object, and A = (aij) is the measurement matrix.

The CT image reconstruction method can be empowered by the CS technique to reduce the necessary datasets and improve the image quality under many less favorable conditions. The CS-based reconstruction method can be expressed as

\min_{f} E (f), subject to A f = b and f > 0,

(2)

where E(f)> is the regularization function. Then, we can use the penalty function method to convert Eq. (2) into an unconstraint optimization problem

\min_{f} E (f) + λ ∥ b - A f ∥_{2}^{2} .

(3)

In CS theory, it is difficult to apply L₀ norm, the most ideal regularization norm, in image reconstruction [26, 27]. Thus, L₀ norm is commonly replaced by L₁ norm in CT image reconstruction [28, 29]. Theoretically, L_1/2 norm is closer to L₀ norm, which can produce sparser solution and reconstruct high-quality CT images [30, 31]. We propose an algorithm based on L_1/2 regularization to solve binary image reconstruction problem [24]: using L_1/2 norm as the regularization norm and gradient transformation as the sparse transformation, Eq. (3) can be redefined as

f = \underset{f}{\arg \min} ∥ ϕ (f) ∥_{1 / 2}^{1 / 2} + β ∥ b - A f ∥_{2}^{2},

(4)

where ϕ(f) is the sparse transformation. Using an intermediate variable d=ϕ (f), Eq. (4) can be converted into

f = \underset{f, d}{\arg \min} {‖ d ‖}_{1 / 2}^{1 / 2} + β {‖ b - A f ‖}_{2}^{2} + ρ {‖ d - ϕ (f) ‖}_{2}^{2} .

(5)

Then we can use Split-bregman method [32, 33] to convert Eq. (5) into two unconstrained optimal problems:

\begin{array}{l} (f^{k + 1}, d^{k + 1}) = & \underset{f, d}{\arg \min} {‖ d ‖}_{1 / 2}^{1 / 2} + β {‖ b - A f ‖}_{2}^{2} \\ + ρ {‖ d - ϕ (f) - b^{k} ‖}_{2}^{2}, \end{array}

(6)

b^{k + 1} = b^{k} + (ϕ (f^{k + 1}) - d^{k + 1}) .

(7)

To solve Eq. (6), we can split it into Eqs. (8) and (9):

f^{k + 1} = \underset{f}{\arg \min} β {‖ b - A f ‖}_{2}^{2} + ρ {‖ d^{k} - ϕ (f) - b^{k} ‖}_{2}^{2},

(8)

d^{k + 1} = \underset{d}{\arg \min} {‖ d ‖}_{1 / 2}^{1 / 2} + ρ {‖ d - ϕ (f^{k + 1}) - b^{k} ‖}_{2}^{2} .

(9)

Here, Split-bregman method based on L_1/2 regularization can reconstruct high-quality binary image from few-view projections [24].

In order to improve the quality of reconstructed binary image from incomplete detection data, we utilize a priori knowledge of the two gray values and a segmentation mechanism (Otsu’s algorithm) in reconstruction process. In Otsu’s algorithm [25, 34], the normalized histogram of reconstructed image is p[m], m∈[m_min, m_max], where [m_min, m_max] is the gray-level range. Setting a threshold T [m_min, m_max] to divide the gray-level range into two classes ([m_min, T] and [T+1, m_max]), the reconstructed binary image can be divide into two gray value regions. The class separability associated with T is defined as

S^{2} (T) = w_{0} {(μ_{0} - μ)}^{2} + w_{1} {(μ_{1} - μ)}^{2} = w_{0} w_{1} (μ_{1} - μ_{0}),

(10)

where w₀(T)=p[m] (m=m_min T), w₁(T)=p[m] (m=T+1 m_max), μ₀(T)=mp[m] (m=m_minT), μ₁(T)=mp[m] (m=T+1m_max), and μ=mp[m] (m=m_minm_max). The purpose of Otsu’s algorithm is to search an optimal threshold ^T’ to maximize S²(T),

T^{'} = \arg \max_{T} [S^{2} (T)],

(11)

and this method can localize the reconstructed image structure information in non-destructive detection. In our algorithm, Otsu’s method calculates the optimal segmentation threshold T, and segment the reconstructed image with two parts. Meanwhile, the priori knowledge (the true gray value) is used to determine the gray value of each part. Finally, we use the segmented image as an intermediate image in next iterative loop.

In the implementation, the whole iteration process of our proposed algorithm can be summarized as follows: [Step 1]

Step 1 Initialize reconstructed image f = 0;

Step 2 Input measured data b, and calculate intermediate image f_ART using ART algorithm:

For i=1,2,⋅,N_angle

f^{k + 1} = f^{k} + λ \frac{a_{i j}}{{‖ A^{i} ‖}^{2}} (b^{i} - A^{i} f^{k}), k = 0,1,...;

(12)

End

Step 3 Positivity constraint for f_ART using Eq. (13):

f = {_{0, f_{i, j} ＜ 0}^{f_{i, j}, f_{i, j} \geq 0};

(13)

Step 4 Update the intermediate image f_SB using the Split-bregman method based on L_1/2 regularization.

For k = 1,2,,N_SB

f^{k + 1} = \underset{f}{\arg \min} β {‖ b - A f ‖}_{2}^{2} + ρ {‖ d^{k} - \nabla f - b^{k} ‖}_{2}^{2};

(14)

d^{k + 1} = \min_{d} {‖ d ‖}_{1 / 2}^{1 / 2} + ρ {‖ d - \nabla f^{k + 1} - b^{k} ‖}_{2}^{2};

(15)

b^{k + 1} = b^{k} + (\nabla f^{k + 1} - d^{k + 1});

(16)

End

Step 5 Update the intermediate image f_Otsu using Otsu’s method, and calculate optimal segmentation threshold ^T’,

T^{'} = \arg \max_{T} {w_{0} (T) w_{1} (T) (μ_{1} (T) - μ_{0} (T))};

(17)

Step 6 Initialize next iteration image f:

f = (1 - δ) f_{SB} + δ f_{Otsu}, 0 < δ < 1;

(18)

Step 7 Go to Step 2 until the stopping criterion is met.

In our algorithm, some key parameters are selected according to experimental analysis. The constraint factor λ in ART method is determined to accomplish initialization reconstruction. Two important parameters of β and ρ in Split-bregman method are determined based on L_1/2 regularization to realize optimization reconstruction. Finally, the scale factor δ is determined to combine Split-bregman method based on L_1/2 regularization and Otsu’s method. For better reconstruction results, reconstruction errors are calculated to obtain optimal parameters of experiments.

III. RESULTS

To demonstrate the feasibility of our proposed method in non-destructive detection, numerical simulations were performed with three binary phantoms, i.e the mandible, turbine blade and limestone phantoms in Fig. 1. The pixel size is 256× 256. The binary mandible phantom was derived from a real mandible CT image, which contains teeth region and bone region; the binary turbine blade phantom having two gray levels (0 and 1) and containing turbine blade and background regions, was derived from a turbine blade CT image; and the binary limestone phantom, derived from a rock specimen CT image, has two gray levels (0 and 1) and contains surrounding rock region, internal porosity region (air) and background region (air).

Fig. 1.

The mandible phantom containing teeth and bone regions (a), the turbine blade phantom consisting of alloy material (b) and the limestone phantom containing surrounding rock and internal porosity (c). Image (a) was from a real mandible CT measurement, and Images (b) and (c) were scanned by industrial CT system.

Here, a typical parallel-beam geometry of the CT system is assumed. To compare and analyze reconstruction results, algebraic reconstruction technique (ART) [35], total variation based algebraic reconstruction technique (ART-TV) [28], Split-bregman method based on L_1/2 regularization (SB-L_1/2) [24], and our method, i.e. Split-bregman method based on L_1/2 regularization and Otsu’s method (SB-Otsu), were used to reconstruct the three phantoms, respectively, in the iteration numbers of 200 for all the reconstruction process. The ART was implemented in Eq. (12). The ART-TV were implemented in two loops, the outer loop implementing ART to reduce data discrepancy, and the inner loop to minimize image’s TV. In the inner loop, the gradient descent method was used:

\begin{array}{l} {‖ \nabla f ‖}_{1} = \sum_{i, j} g_{i, j}, \\ g_{i, j} = \sqrt{{(f_{i, j} - f_{i + 1, j})}^{2} + {(f_{i, j} - f_{i, j + 1})}^{2}}, \end{array}

(19)

where ‖f‖₁ denotes TV of f, f_i,j is the pixel value of the discrete 2D image, and g_i,j is discrete gradient.

f^{(m + 1)} = f^{(m)} － γ ω υ / | υ |,

(20)

where γ is the gradient descent control coefficient, ω= ‖f^(m+1)-f^(m)‖₂ is the scaling coefficient of the gradient descent, $υ = (\partial ∥ \nabla f ∥_{1} / \partial f_{i, j})_{f_{i, j} = f_{i, j} [n, m]}$ is the gradient direction when f_i,j=f_i,j[n,m]. The implementation of Split-bregman method based on L_1/2 regularization was mainly used in Eqs. (14)–(16). Finally, we analyzed key parameters of the reconstructed algorithms, and the values of optimal parameters are λ=0.5, β=1000, ρ=1, δ =0.2 and γ=0.2.

In binary mandible phantom reconstruction, the angular scanning was from 0 to 180 in 30 steps, producing 6 projections. Then we added 0.5% Gaussian noise into projection data, and the reconstructed results using different methods are shown in Fig. 2A. In binary turbine blade phantom reconstruction, the angular scanning was from 0 to 180 in 18 steps, producing 10 projections. We added 0.5% Gaussian noise into projection data and used the 10 projections to reconstruct the turbine blade phantom (Fig. 2B). In binary limestone phantom reconstruction, the angular scanning was from 0 to 180 with 22.5 steps, producing 8 projections. Then we added 0.5% Gaussian noise into projection data, and the reconstructed results using different methods are shown in Fig. 2C.

Fig. 2.

Images of the binary (A) mandible (B) blade and (C)limestone phantoms, reconstructed using (a) ART, (b) ART-TV, (c) SB-L_1/2, and (d) SB-Otsu methods. The display window for the images is [0,1].

Then, the root mean square error (RMSE) in Eq. (21) was used to quantify the reconstructed results,

RMSE = [\sum {(f_{i, j} - {\hat{f}}_{i, j})}^{2} / N_{f}]^{1 / 2},

(21)

where ${\hat{f}}_{i, j}$ is the reconstructed pixel value, f_i,j is the true pixel value, and Nf is the pixel number of the phantom. The RMSE values for the two reconstructed phantoms using different methods are summarized in Table 1, and the iterative process curves for these reconstructed algorithms are shown in Fig. 3.

Fig. 3.

(Color online) RMSE lines of image reconstruction using different algorithms in 200 iterations. (a) mandible phantom reconstructed from 6 projections, (b) turbine blade phantom reconstructed from 10 projections, (c) limestone phantom reconstructed from 8 projections.

RMSE values for the three reconstructed phantoms using different methods

Phantoms	ART	ART-TV	SB-L_1/2	SB-Otsu
Mandible	0.1319	0.0954	0.0753	0.0534
Turbine blade	0.1372	0.1308	0.1287	0.0864
Limestone	0.0728	0.0391	0.0334	0.0275

From Fig. 2, the reconstructed images using ART methods contain a lot of noise and artifacts, while the reconstructed images using ART-TV, SB-L_1/2 and SB-Otsu methods have clearer edges. From Table 1, the SB-Otsu method performed better in binary CT image reconstruction from severe incomplete projection data than other three methods. It is more effective in treating noise and artifacts. From Fig. 3, SB-Otsu algorithm can reconstruct higher quality binary CT images with the same iteration numbers.

IV. DISCUSSIONS AND CONCLUSION

In CS theory, an image can be reconstructed from a rather limited amount of detection data as long as it can be sparsely represented in an appropriate domain and determined from these data. From the results in Sec. III, the CS-based method is effective in treating noise and artifacts for the reconstructed images from incomplete detection data. However, the CS-based algorithm is not omnipotent, the reconstructed images may suffer the loss of some detail information due to severe incomplete detection data. In order to improve reconstructed image quality from severe incomplete detection data, a priori knowledge of gray value information and a segmentation mechanism (Otsu’s method) are introduced into binary CT image reconstruction.

There are still some issues worth further discussion in this paper. First, the amount of projection data used in reconstruction depends on the structure of reconstructed image in our study. More complex the structure of reconstructed image is, more projections reconstructing high-quality image needs. Second, the correlation of projection vectors also determines the quality of reconstructed image. Weaker the correlation of projection vectors is, higher the quality of reconstructed image is. Third, our proposed algorithm combines Split-bregman method based on L_1/2 regularization and Otsu’s method. In reconstruction process, we need to set a proper weighting coefficient δ for Otsu’s algorithm implementation. If the value of weighting coefficient δ is too large, the quality of reconstructed image will be lowered. In our study, we select the value of weighting coefficient δ according to experimental analysis.

Currently, we analyze three binary image phantoms: a mandible phantom, a turbine blade phantom and a limestone phantom. The limestone and turbine blade phantoms are tested by industrial CT, our proposed method can reconstruct these two phantoms from few-view projections to reduce scan time. The mandible phantom is scanned by the medical CT system, our proposed method can reconstruct mandible phantom from few-view projections to reduce radiation dose. This initial methodological study is mainly focused on phantom simulation analysis, the follow-up study will deal with more general settings for industrial and biomedical applications.

In conclusion, we proposed a binary CT image reconstruction algorithm based on compressed sensing and Otsu’s method to reduce scan time or radiation dose in non-destructive detection, and investigated the feasibility and potential of the proposed method. The experimental results demonstrated that our proposed method is very effective to reduce noise and suppress artifacts in binary CT image reconstruction from incomplete projection data. In further work, we will analyze real data reconstruction, and a systematic study is beyond the scope of this initial investigation.

References

[1]

A C Kak and M Slaney. Principles of Computerized Tomographic Imaging. IEEE Press, 1988. DOI: 10.1137/1.9780898719277

[2]

G Wang, H Yu, B D Man.

An outlook on X-ray CT research and development

. Med Phys, 2008, 35: 1051-1064. DOI: 10.1118/1.2836950