research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

A modified discrete algebraic reconstruction technique for multiple grey image reconstruction for limited angle range tomography

CROSSMARK_Color_square_no_text.svg

aNational Synchrotion Radiation Laboratory, University of Science and Technology of China, 3#419, No.42 Hezuohua South Road, Hefei, Anhui 230029, People's Republic of China
*Correspondence e-mail: ychtian@ustc.edu.cn

Edited by P. A. Pianetta, SLAC National Accelerator Laboratory, USA (Received 22 July 2015; accepted 11 January 2016; online 20 February 2016)

The `missing wedge', which is due to a restricted rotation range, is a major challenge for quantitative analysis of an object using tomography. With prior knowledge of the grey levels, the discrete algebraic reconstruction technique (DART) is able to reconstruct objects accurately with projections in a limited angle range. However, the quality of the reconstructions declines as the number of grey levels increases. In this paper, a modified DART (MDART) was proposed, in which each independent region of homogeneous material was chosen as a research object, instead of the grey values. The grey values of each discrete region were estimated according to the solution of the linear projection equations. The iterative process of boundary pixels updating and correcting the grey values of each region was executed alternately. Simulation experiments of binary phantoms as well as multiple grey phantoms show that MDART is capable of achieving high-quality reconstructions with projections in a limited angle range. The interesting advancement of MDART is that neither prior knowledge of the grey values nor the number of grey levels is necessary.

1. Introduction

Tomography techniques have been widely adopted for obtaining three-dimensional images of physical objects non-invasively. Fourier-based methods and algebraic reconstruction methods are two popular types of principal algorithms for image reconstruction. Fourier-based methods, which are directly related to the Radon transform, are still commonly used because of their small computational burden (Mueller et al., 1999[Mueller, K., Yagel, R. & Wheller, J. J. (1999). IEEE Trans. Med. Imaging, 18, 519-537.]). Algebraic reconstruction methods, which formulate the reconstruction problem as a linear system of equations (Zeng, 2010[Zeng, G. L. (2010). Medical Image Reconstruction: A Conceptual Tutorial. Berlin: Springer.]; Herman, 2009[Herman, G. T. (2009). Fundamentals of Computerized Tomography: Image Reconstruction from Projections. Springer.]), can achieve images with better quality than Fourier-based methods, especially when the number of projections is limited or the angles are in a limited range.

A limited range of tilt angle occurs frequently in electron tomography (Milne & Subramaniam, 2009[Milne, J. L. & Subramaniam, S. (2009). Nat. Rev. Microbiol. 7, 666-675.]), soft and hard X-ray nano-tomography (Duke et al., 2013[Duke, E. M., Razi, M., Weston, A., Guttmann, P., Werner, S., Henzler, K., Schneider, G., Tooze, S. A. & Collinson, L. M. (2013). Ultramicroscopy, 143, 77-87.]; Tian et al., 2008[Tian, Y. C., Li, W. J., Chen, J., Liu, L. H., Liu, G., Tkachuk, A., Tian, J. P., Xiong, Y., Gelb, J., Hsu, G. & Yun, W. B. (2008). Rev. Sci. Instrum. 79, 103708.]). A flat sample holder occludes the light when the tilt angle of the tomography reaches a high angle (usually ±75°). The limited spacing for specimen holders between the pole pieces of the objective lens may also prevent the rotation of the sample holders. These would lead to a `missing wedge' in the collected data. In this case it is very difficult to obtain a quantitative measure of the object based on the reconstruction. In other words, the boundaries and the grey values of the reconstructions do not match well with the physical objects.

Discrete tomography has evolved into a powerful, robust and flexible method for reconstructing images with a small number of projections when the physical objects contain only a few discrete grey levels (Herman & Kuba, 2007[Herman, G. T. & Kuba, A. (2007). Advances in Discrete Tomography and its Applications. Berlin: Springer.]). Numerous discrete tomography reconstruction algorithms have been proposed in the last few decades. Schüle et al. presented a primal-dual subgradient algorithm for reconstructing binary images from projections with few-view or a limited range of angles (Weber et al., 2004[Weber, S., Schule, T., Hornegger, J. & Schnorr, C. (2004). Lect. Notes Comput. Sci. 3322, 38-51.]; Schüle et al., 2005[Schüle, T., Schnörr, C., Weber, S. & Hornegger, J. (2005). Discrete Appl. Math. 151, 229-243.]). There are many evolutionary algorithms for discrete tomography (Batenburg, 2005[Batenburg, K. J. (2005). Discrete Appl. Math. 151, 36-54.]; Weber et al., 2006[Weber, S., Nagy, A., Schulle, T., Schnorr, C. & Kuba, A. (2006). Lect. Notes Comput. Sci. 4245, 146-156.]). However, most of these algorithms are limited in achieving a high-quality reconstruction of binary images.

Recently, Batenburg et al. proposed the discrete algebraic reconstruction algorithm (DART), which can address binary images and images that contain three or more grey levels with prior knowledge of the grey levels of each material (Batenburg & Sijbers, 2007[Batenburg, K. & Sijbers, J. (2007). IEEE International Conference on Image Processing (ICIP 2007), pp. IV-133-IV-136. IEEE.], 2011[Batenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542-2553.]; Batenburg et al., 2010[Batenburg, K. J., Sijbers, J., Poulsen, H. F. & Knudsen, E. (2010). J. Appl. Cryst. 43, 1464-1473.]). When samples consist of a few discrete materials, it is possible to obtain accurate reconstructions even if the number of projections is limited. Furthermore, thresholding is already applied during the reconstruction, thus DART can achieve segmented images directly. DART has been steadily developed into a powerful algorithm for many experimental datasets, such as reconstructing catalyst particles for bamboo-like carbon nanotubes from transmission electron microscopy (Bals et al., 2007[Bals, S., Batenburg, K. J., Verbeeck, J., Sijbers, J. & Van Tendeloo, G. (2007). Nano Lett. 7, 3669-3674.]), quantitative analysis of zeolite through discrete electron tomography (Bals et al., 2009[Bals, S., Batenburg, K. J., Liang, D. D., Lebedev, O., Van Tendeloo, G., Aerts, A., Martens, J. A. & Kirschhock, C. E. A. (2009). J. Am. Chem. Soc. 131, 4769-4773.]) and reconstruction of a diamond from a series of micro-computed tomography projections (Batenburg & Sijbers, 2011[Batenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542-2553.]).

A common assumption of DART is that all of the grey levels of the reconstruction are known a priori; therefore, it is important to know the appropriate assumption of the grey levels and the thresholds for obtaining high-quality reconstructions. Although it is easy to know the grey levels in simulation experiments, obtaining a priori knowledge in practice is difficult. To reduce the requirements of prior knowledge, a semi-automatic discrete grey level selection (DGLS) method was proposed to estimate the grey levels of an image (Batenburg et al., 2011[Batenburg, K. J., van Aarle, W. & Sijbers, J. (2011). Pattern Recognit. Lett. 32, 1395-1405.]). A homogeneous area in the reconstruction is selected manually to calculate the grey value of each grey level. Afterwards, a projection distance minimization DART (PDM-DART) method was proposed to estimate the grey level parameters adaptively; the method required that only the number of grey levels should be known in advance (van Aarle et al., 2012[Aarle, W. van, Batenburg, K. J. & Sijbers, J. (2012). IEEE Trans. Image Process. 21, 4608-4621.]). Some algorithms keep a specific component of the reconstruction fixed, and the rest of the components changed during the process. For example, the BgART algorithm could detect the background and keep the background value fixed for every iteration (Messaoudi et al., 2013[Messaoudi, C., Aschman, N., Cunha, M., Oikawa, T., Sorzano, C. O. & Marco, S. (2013). Microsc. Microanal. 19, 1669-1677.]). This algorithm facilitates segmentation of objects and denoising of the reconstruction.

The series of algorithms of DART are effective when the number of grey levels is small (Batenburg & Sijbers, 2011[Batenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542-2553.]). In some actual applications, such as for integrated circuit chips, marine sediment and shale, which have multiple phases and various compositions (Kanitpanyacharoen et al., 2013[Kanitpanyacharoen, W., Parkinson, D. Y., De Carlo, F., Marone, F., Stampanoni, M., Mokso, R., MacDowell, A. & Wenk, H.-R. (2013). J. Synchrotron Rad. 20, 172-180.]; Uramoto et al., 2014[Uramoto, G.-I., Morono, Y., Uematsu, K. & Inagaki, F. (2014). Limnol. Oceanogr. Methods, 12, 469-483.]), prior knowledge of grey levels, even simply the number of grey levels, is not easy to acquire. As the number of grey levels increases, the reconstruction quality decreases (van Aarle et al., 2012[Aarle, W. van, Batenburg, K. J. & Sijbers, J. (2012). IEEE Trans. Image Process. 21, 4608-4621.]). Hence, it is necessary to develop a method of reconstructing multiple grey images when the range of the tilt angle is limited and without prior knowledge of the grey levels.

In this paper a modified DART (MDART) algorithm is proposed to reconstruct multiple grey images. In our method the available projection data were reconstructed using an algebraic reconstruction method. Then, the reconstruction images were divided into discrete regions depending on threshold values that were selected manually. The segmentation images were used as an input of the proposed method. Each independent region was a research object. The grey values of each discrete region were estimated by solving the algebraic reconstruction linear equations. Thus, the number of unknown variables in the equations was decreased by several orders of magnitude. Then, updating boundary pixels and correcting the grey values of each region were iteratively alternated until the stop criterion or the maximum iteration number was met. This method has been tested on both binary and multiple grey images while the angle of projection was limited.

2. Notation and concepts

Let [{\bf{x}}] = [({{x_j}}) \in {{\bb{R}}^{n}}] represent a two-dimensional unknown image, which denotes the object function with n pixels. A finite set of m measured projection data [{\bf{p}}] = [({{p_i}}) \in {{\bb{R}}^m}] is defined by

[\textstyle\sum\limits_{j\,=\,1}^{n} {w}_{ij}\,{x}_{j} = {p}_{i}, \qquad i=1,2,\ldots,m,\eqno(1)]

where wij can be interpreted as the contribution of the jth pixel xj to the ith detector value pi. [{\bf{W}}] = [({{w_{ij}}}) \in {{\bb{R}}^{m \times n}}] is the system matrix. Therefore, the tomography reconstruction algorithm tries to find the unknown image [{\bf{x}}] by solving the linear system of equations

[{\bf{Wx}}={\bf{p}}.\eqno(2)]

In practice, numerous reconstruction algorithms have been proposed to solve equation (2)[link] (Kak & Slaney, 2001[Kak, A. C. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics.]; Zeng, 2010[Zeng, G. L. (2010). Medical Image Reconstruction: A Conceptual Tutorial. Berlin: Springer.]). For the experiments in this paper, the simultaneous algebraic reconstruction technique (SART) (Andersen & Kak, 1984[Andersen, A. & Kak, A. (1984). Ultrason. Imaging, 6, 81-94.]) TVM (SART-TVM) will be used to reconstruct the images to obtain a preliminary result, and SART will also play an important role in MDART. Therefore, SART and TVM are briefly reviewed. SART is a linear algorithm that finds a solution [{\bf{x}}] such that the difference between the projections from the original image and the reconstruction image [||{\bf{Wx}}-{\bf{p}}||] is minimal. Almost all iteration algorithms start with an initial guess [{\bf{x}}] = [{\bf{x}}^{(0)}], and a new reconstruction [{\bf{x}}^{(k)}] is iteratively calculated from [{\bf{x}}^{(k-1)}]. SART carries out updating operations using all projections at an angle each time as follows,

[{x}_{j}^{(k)} = {x}_{j}^{(k-1)} + \lambda\,{{1}\over{{\sum_{i\,=\,1}^{h}}{w}_{ij}}} \,\,\sum_{i\,=\,1}^{h} {{ {w}_{ij}\left({p}_{i}-\textstyle\sum_{t\,=\,1}^{n} {w}_{it\,}{x}_{t}^{(k-1)}\right) }\over{ \textstyle\sum_{t\,=\,1}^{n} {w}_{it} }}, \eqno(3)]

for j = [1\ldots n], where [\lambda] is a relaxation parameter and h is the number of total detectors at each angle.

TVM-based algorithms find the approximate solution as the following optimization problem (Yang et al., 2015[Yang, X., Hofmann, R., Dapp, R., van de Kamp, T., dos Santos Rolo, T., Xiao, X., Moosmann, J., Kashef, J. & Stotzka, R. (2015). Opt. Express, 23, 5368-5387.]),

[x^*= {\rm{arg}}\mathop {\min }\limits_x \left[||Wx-p||_2^2+\gamma{\rm{TV}}(x)\right],\eqno(4)]

where [\gamma] is a regularization parameter indicating the importance of the two terms in the equation and [{\rm{TV}}(\ldots)] represents the TVM norm, which is the l1-norm of the gradient image as defined by Sidky et al. (2006[Sidky, E. Y., Kao, C.-M. & Pan, X. (2006). J. X-ray Sci. Technol. 14, 119-139.]):

[\eqalignno{{\rm{TV}}(x)= {}& \sum\limits_{i,j,k} \Big\{ {{\big[{x\left({i + 1,j,k}\right) - x\left({i,j,k} \right)} \big]}^2} \cr&+ {{\big[{x\left({i,j + 1,k}\right) - x\left({i,j,k} \right)} \big]}^2} \cr&+ {{\big[{x\left({i,j,k + 1} \right) - x\left({i,j,k} \right)} \big]}^2} \Big\}^{1/2}.&(5)}]

3. The modified discrete algebraic reconstruction technique

The main limitations in traditional discrete tomography methods are the requirements of prior knowledge of the greyscales and the limitation in the number of discrete greyscales. In order to overcome these limitations the MDART algorithm was proposed for reconstruction of incomplete tomographic data of multiple discrete greyscale phantoms without any prior knowledge. A flow chart of the algorithm is shown in Fig. 1[link]. In this approach, we first perform tomographic reconstruction using a traditional algorithm (SART-TVM) to provide an initial guess for further calculation. The initial slice is then segmented based on the greyscales of the image resulting in a number of isolated regions. Since the tomographic data were incomplete, the boundary and the greyscales in all of the regions might not be accurate and need to be improved through iterations. These regions were treated as research objects and their original greyscales disregarded. We then perform a forward Radon transform of the slice with unknown greyscales in all of the regions. The result of the forward Radon transform needs to match the experimentally measured projection data, which enables us to recover the greyscales in the regions by solving the linear equations [{\bf{Wx}}_s] = [{\bf{p}}] using the LSRQ method (Paige & Saunders, 1982[Paige, C. C. & Saunders, M. A. (1982). ACM Trans. Math. Softw. 8, 43-71.]; Barrett et al., 1994[Barrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C. & Van der Vorst, H. (1994). Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Philadelphia: Siam.]). Due to the uncertainty of the region boundaries, the greyscales on the boundary pixels are further updated by evaluating the difference between the experimentally measured projection data and the Radon transform of the slice with the boundary pixels set to zero. The boundary pixels were smoothed to regularize the boundary reconstruction. The updated intensity in each boundary pixel is used to identify the neighbouring region that has the closest greyscale. This boundary pixel is subsequently merged into the chosen region, effectively shifting the region boundaries. The above-described procedure is performed iteratively until the stop criterion is met, indicating successful reconstruction of the slice.

[Figure 1]
Figure 1
Flow chart of the MDART algorithm.

There are some main differences with the previous methods. The procedure for identifying regions is added to the algorithm. Each independent region is chosen as a research object, instead of the greyscale. The greyscale of each region was calculated by solving the linear equations. Then the boundary pixel is merged into one of its neighbouring regions that has the closest greyscale.

A simple example is presented to demonstrate the process and details of the algorithm. Fig. 2[link] indicates the various steps of the algorithm. Fig. 2(a)[link] is a binary image with two grey levels (only 0 and 1). Suppose that there were only 121 projections that were equally distributed between 0° and 120°. Fig. 2(b)[link] shows the initial guess through SART-TVM reconstruction after 100 iterations. It is obvious that the grey values and the shape of the reconstruction do not match well with the original phantom.

[Figure 2]
Figure 2
Sequence of operations of the various steps of the MDART algorithm.

Then, the initial reconstruction was segmented into discrete regions depending on the grey values. Many pieces of software (such as Amira, ImageJ) could be used to segment the reconstruction with friendly interactive interfaces. To make segmentation work more efficiently, a group of threshold values was chosen to ensure that homogeneous adjacent pixels were divided into the same region as much as possible. For a group of three-dimensional data, a few slices that contain the most compositions are adopted to choose the threshold values. The threshold values should be appropriate for all of the slices. There is no strict requirement for the selection of the threshold, but it would affect the convergence rate of the algorithm. The threshold values [\tau] = [({\tau}_{1}\ldots{\tau}_{l})] were selected depending on the manual selection. Then, the initial reconstruction was segmented according to the threshold function,

[s\left({x}_{i}\right) = \left\{\matrix{ 1\hfill & \quad ({x}_{i}\le{\tau}_{1})\hfill \cr 2\hfill & \quad ({{\tau}_{1}\,\lt\,{x}}_{i}\le{\tau}_{2})\hfill \cr \,\vdots\hfill & \quad \quad\quad\vdots\hfill \cr l\hfill & \quad ({{\tau}_{l-1}\,\lt\,{x}}_{i}\le{\tau}_{l})\hfill \cr l+1\hfill & \quad \left({x}_{i}\,\gt\,{\tau}_{l}\right)\hfill }\right..\eqno(6)]

A group of threshold values [\tau] = (225, 200, 165, 70, 20, 10)/255 is chosen depending on the features of Fig. 2(b)[link]. The segmentation result is shown in Fig. 2(c)[link].

The procedure for identifying regions was performed afterwards. If the difference between two adjacent regions was less than the threshold value , the two regions were combined. 143 regions are identified and presented in Fig. 2(d)[link].

In this algorithm, each independent region was chosen as a research object, instead of the grey values. Thus, there were 143 unknown variables in the new matrix [{\bf{x}}_s]. The values of each independent region were calculated by solving the linear equations [{\bf{Wx}}_{s}] = [{\bf{p}}] using the LSQR method (Paige & Saunders, 1982[Paige, C. C. & Saunders, M. A. (1982). ACM Trans. Math. Softw. 8, 43-71.]; Barrett et al., 1994[Barrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C. & Van der Vorst, H. (1994). Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Philadelphia: Siam.]). The image endowed with new values is shown in Fig. 2(e)[link].

The boundary region N of the image was the set of all pixels that had at least one different pixel in their eight-connected neighbourhoods. The boundary of the image with the new values is shown in Fig. 2(f)[link]. Subsequently, the values of the pixels were updated by the SART iterations, while the other pixels were fixed. Thus, the number of variables in the linear equations was obviously reduced compared with the conventional SART. Fig. 2(g)[link] shows the result of the boundary reconstruction after ten SART iterations. To regularize the boundary reconstruction, smoothing must be applied to the boundary pixels. Let b denote the average of eight-connected neighbourhoods of a pixel. The boundary pixels were smoothed with b as follows,

[x_i^{\,(t)}= 0.7x_i^{\,(t)} + 0.3b. \eqno(7)]

After that, the pixels in N were merged into one of their neighbouring regions that had the closest greyscale.

Subsequently, iterations of the algorithm were performed until the stop criterion was met. In this study, the sum of the squared difference between two successive iteration results was used as a measurement of the degree of convergence. The iterations were stopped once this parameter was smaller than 500 or the maximum iteration number was met.

The final reconstruction of the phantom is shown in Fig. 2(h)[link]. It was almost the same as the original phantom image, except for a few pixels with different values. It reveals that MDART is effective while the angle range of the projections is limited.

4. Experiments and discussion

In this section, a series of experiments are presented to demonstrate the reconstruction performance of the proposed method and compared with some commonly used methods. For all algorithms, the six phantom images shown in Fig. 3[link] were used. Phantom 1 and phantom 2 are binary images (Batenburg & Sijbers, 2011[Batenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542-2553.]). Phantom 3 and phantom 6 are multiple grey phantoms that contain 15 and 14 grey levels, respectively. Their grey-level distribution is shown in Figs. 3(d) and 3(h)[link], respectively. Phantom 4 is the well known Shepp–Logan phantom, which contains six grey levels. Phantom 5 is composed of 49 homogeneous discs. Every circle has a different value, and a total of 50 grey levels are contained in the phantom. The size of all of the phantoms is 256 × 256 pixels.

[Figure 3]
Figure 3
Simulated phantoms used for the simulation experiments: (a), (b) are binary phantoms; (c) is a multiple grey phantom, and its grey values are shown in (d); (e) Shepp–Logan phantom; (f) multiple grey circle phantom; (g) multiple grey Shepp–Logan phantom, and (h) shows the grey values.

Five reconstruction algorithms were compared in this section. We provide a brief review of them as follows:

(i) FBP. Filtered back projection (FBP) is a Fourier-based method that is related to the inverse Radon transform. In the experiments, we used a Ram-Lak filter and linear interpolation in the projection domain as the MATLAB defaults.

(ii) SART. The SART algorithm, which has been described in §2[link], was run with 500 iterations in the experiments. A positivity constraint was used after each iteration to obtain a better result.

(iii) SART-TVM. This algorithm had also been described in §2[link], and there are many studies that present details of the algorithm (Sidky et al., 2006[Sidky, E. Y., Kao, C.-M. & Pan, X. (2006). J. X-ray Sci. Technol. 14, 119-139.]; Yu & Wang, 2009[Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791-2805.]; Liang et al., 2013[Liang, Z., Guan, Y., Liu, G., Bian, R., Zhang, X., Xiong, Y. & Tian, Y. (2013). Proc. SPIE, 8851, 885113.]). In this work, we refer to the original article (Sidky et al., 2006[Sidky, E. Y., Kao, C.-M. & Pan, X. (2006). J. X-ray Sci. Technol. 14, 119-139.]) for details except using SART to replace ART. We used a = 0.2 (gradient descent parameter), Ngrad = 20 (the total number of gradient descent steps) and 500 iterations for the analysis.

(iv) DART. The DART algorithm program is included in the ASTRA toolbox (Palenstijn et al., 2013[Palenstijn, W. J., Batenburg, K. J. & Sijbers, J. (2013). The ASTRA Tomography Toolbox, in Proceedings of the 13th International Conference on Computational and Mathematical Methods in Science and Engineering (CMMSE 2013), 24-27 June 2013, Cabo de Gata, Almeria, Spain.]) and is available online. 500 iterations were performed to obtain convergent reconstructions. The grey levels and threshold values were applied as prior knowledge, but they were not adopted in the four other algorithms. The 500 SIRT algebraic reconstruction method has been used for the initial reconstruction and 15 in the successive ones. The fix probability is 0.99, and equation (7)[link] is used to smooth the boundary reconstruction.

(v) MDART. Details of MDART were described in §3[link]. The threshold value of the regions that need to be identified was 0.01 in the experiments.

The total number of different pixels between the reconstruction and original image was used as the performance metric. Because prior knowledge of the grey levels was not used in all of the algorithms, except for DART, the segmentation process was not performed after the iteration step. Therefore, a certain degree of difference in the pixel values between the reconstruction and original image was acceptable. We define the computing method for the total number of different pixels as follows:

[K= \left|\left\{i\,\Big|\left[\left|{x}_{i}^{\,r}-{x}_{i}\right| \,\,\gt\,\, \max(\omega\Delta\rho,0.003)\right]\right\}\right|,\eqno(8)]

where [\Delta\rho] is the minimum interval between the two grey levels in the original images; [\omega] is an indicator between 0 and 1 ([\omega] = 0.03 was used in the experiment); and [ \omega\Delta\rho] is the threshold value of the errors. However, ordinary image formats are usually 256 grey levels. When the difference of a pixel between the reconstruction and original image is less than a grey level it can be ignored. Thus, when [\omega\Delta\rho] < 0.003, the threshold is set to 0.003.

If prior knowledge of the grey levels [{\boldrho}] = [({\rho}_{1},{\rho}_{2},\ldots,{\rho}_{l})] is known in advance, the reconstructions can be segmented as follows:

[s\left({x}_{i}\right) = \left\{ \matrix{ {\rho}_{1},\hfill & \hfill {x}_{i}\,\le\, {{({\rho}_{1}+{\rho}_{2})}/{2}} \cr & \vdots & \cr {\rho}_{i},\hfill & \quad\hfill {{({{\rho}_{i-1}+{\rho}_{i})}/{2}} \,\,\lt\,\, x}_{i}\,\le\, {{({\rho}_{i}+{\rho}_{i+1})}/{2}} \cr & \vdots \cr {\rho}_{l},\hfill & \hfill {x}_{i} \,\,\gt\,\,{{( {\rho}_{l-1}+{\rho}_{l} )}/{ 2}} } \right. \eqno(9)]

The relative error of the value of the grey levels was defined as follows:

[{{\left|{\rho }_{i}^{\,r}-{\rho }_{i}\right|}\over{\max\left(\{\rho _{i}\}\right)}}\times 100%, \qquad i=1\ldots l.\eqno (10)]

4.1. Binary phantoms reconstruction

In this section, we present the reconstruction results of binary phantoms based on the projections with a limited angular range. As mentioned above, for the DART algorithm, prior knowledge of the grey levels is necessary, whereas, for other techniques, it is not required. Fig. 4[link] shows the reconstructions of phantom 1 with an angle range of 130° and phantom 2 with an angle range of 140° by FBP, SART, TVM, DART and MDART. The angular step is 1°. Obviously, when the angular range of the projections is limited, a missing wedge has a significant influence on the reconstructions using FBP and SART. The pixel error distribution caused by using five algorithms as a function of the angular range is shown in Fig. 5[link], and the total number of all the phantoms is 65536. The results of TVM, DART and MDART are more accurate as the angle range increases. Combining Figs. 4[link] and 5[link], the reconstruction of MDART is better than TVM for the same projections. If prior knowledge of the grey level is applied, DART can obtain very accurate reconstructions of the binary phantoms with a small angle range. It can be concluded that, for binary images of a limited angle range, MDART can achieve high-quality reconstructions without prior knowledge of the grey levels.

[Figure 4]
Figure 4
Reconstructions of binary images using FBP (column 1), SART (column 2), TVM (column 3), DART (column 4) and MDART (column 5). Phantom 1 has an angular range of α = 130° and phantom 2 has an angular range of α = 140°.
[Figure 5]
Figure 5
Pixel error K as a function of the angle range of the projections for FBP, SART, TVM, DART and MDART.

4.2. Multiple grey phantoms reconstruction

For most discrete tomography reconstruction algorithms, it is a challenge to reconstruct a multiple-level grey image. Conventional continuous reconstruction methods cannot achieve high-quality results when the angle range is limited because of the `missing wedge' of data. Fig. 6[link] shows the reconstructed results of phantoms 3, 4, 5 and 6 using the five methods. Many artifacts and deformations of morphology could be found in the reconstructions of FBP and SART. The results of TVM showed a very undesirable staircase effect as the number of grey levels increases. DART also suffered from the missing wedge influence. The MDART technique showed good behaviour for both morphology and grey value. The pixel error K as a function of the angle range of the projections for FBP, SART, TVM, DART and MDART is shown in Fig. 7[link], where erroneous pixels number are calculated using equation (8)[link]. The angular step in all cases is 1° between the projections. It clearly showed that, even without using prior knowledge of the grey levels, MDART performs better than FBP, SART, TVM and DART when processing multiple grey phantoms in a limited angle range.

[Figure 6]
Figure 6
Reconstructions of multiple grey phantoms using FBP (column 1), SART (column 2), TVM (column 3), DART (column 4) and MDART (column 5). Phantom 3 has an angular range of α = 130°, phantom 4 has an angular range of α = 147°, phantom 5 has an angular range of α = 135° and phantom 6 has an angular range of α = 145°.
[Figure 7]
Figure 7
Pixel error K as a function of the angle range of the projections for FBP, SART, TVM, DART and MDART.

Fig. 8[link] shows the number of error pixels of the reconstructions that have been segmented by equation (9)[link] as a function of the angle range of the projections. When the results obtained by FBP, SART and TVM have been segmented, all of the reconstructed results obviously improved. However, it was obvious that the number of error pixels of MDART was still much less than for the four other methods. Note that segmentation was not applied to the reconstructions of DART because they had already been segmented in the algorithm.

[Figure 8]
Figure 8
Pixel error K as a function of the angle range of the projections for FBP, SART, TVM, DART and MDART. The results have been segmented with prior knowledge of the grey levels.

The limitation of the number of grey levels was also investigated. Phantoms, resembling phantom 5, were tested with an increasing number of discs and grey levels. For each phantom, varying angle ranges of projections were simulated. The pixel error K of the reconstructed images as a function of the angle range is shown in Fig. 9[link]. As the number of grey levels increased, the quality of the reconstructions decreased.

[Figure 9]
Figure 9
Pixel error K of the reconstructed images as a function of the angle range of the projections for seven phantoms with different numbers of grey levels; the total number of pixels in the phantoms is 65536.

The relative errors of the grey value of phantoms 3, 4, 5 and 6 using MDART were calculated using equation (10)[link] in a different angle range, as shown in Fig. 10[link]. For the four multiple grey phantoms, the relative errors decreased to a low level with the angle range increased. Most of the relative errors are less than 1% when, for the angle range of phantom 3, α > 135°; for phantom 4, α > 137°; for phantom 5, α > 120°; and for phantom 6, α > 150°. The precision is good enough to distinguish different compositions. It shows that, for most of the phantoms, an angle range of 150° can usually result in a precise reconstruction. MDART is suitable for multiple grey phantom reconstruction in a limited angle range.

[Figure 10]
Figure 10
Relative error of the grey value of the MDART reconstructions of phantoms 3, 4, 5 and 6 in different angle range.

4.3. Convergence

To understand the impact of the initial value of the segmentation on the algorithm, a multiple grey image, phantom 3, with an angle range of 130° is tested.

To segment the reconstruction of the SART-TV result into discrete regions, a group of threshold values was chosen. In Table 1[link] there are ten groups of threshold values with different values. The threshold values of the first eight groups are chosen depending on the features of the reconstruction. The ninth group, which is roughly segmented, is a case where the third hole in the figure is not identified, as shown in Fig. 11(d)[link]. The tenth group starts with a hole in the middle that spreads outward. The eleventh group is segmented without threshold values, but it is manually marked with some obvious feature regions, as shown in Fig. 11(f)[link]. Some of the segmentation results are shown in Figs. 11(b)–11(e)[link], and the others are shown in Figs. S1 and S2 of the supporting information . Using different segmentation results as an input to MDART leads to a different number of required iterations and different final error, as shown in Table 1[link]. The values of the error pixels are slightly different. The number of error pixels as a function of the iteration number is shown in Fig. 11(g)[link]. The threshold values are related to the convergence rate of the algorithm. For three-dimensional data, a few typical slices are used to select the threshold values, and then the threshold values are applied to the entire three-dimensional data. Overall, the method of selecting the threshold values using software (such as Amira, ImageJ, etc.) is convenient and robust.

Table 1
The algorithm progress of ten groups of threshold values; manual segmentation is used to segment the SART-TV reconstruction

  Segment number
  1 2 3 4 5 6 7 8 9 10 11
Threshold values (/255) 240 235 230 220 225 230 195 190 178 240  
225 220 215 195 180 200 135 50 130 230  
205 195 185 155 150 160 70 10 90 215  
185 160 155 125 120 120 40   55 175  
160 140 130 85 85 80 10   10 125  
135 120 100 60 60 35       65  
115 90 65 35 30 9       45  
90 60 35 10 10         25  
65 35 10             10  
30 10                  
10                    
Pixel error 201 223 603 149 152 642 596 144 597 155 446
Iterations 40 73 38 99 59 41 26 6 110 344 220
[Figure 11]
Figure 11
Segment results of the SART-TV reconstruction of phantom 3 and the convergence rate of different segment results. (a) The reconstruction of SART-TV with an angular range of α = 130°. (b) Segmented carefully with 11 threshold values. (c) Segmented with three threshold values. (d) One of the holes in the middle is not identified. (d) One of the holes in the middle spreads outward. (f) According to the reconstruction, a few obvious features of the regions are marked out manually. (g) The convergence rates of different segment results as an initial input into the MDART.

In the past decade, many important tomography reconstruction methods have been proposed. Equally sloped tomography with oversampling reconstruction (Miao et al., 2005[Miao, J., Förster, F. & Levi, O. (2005). Phys. Rev. B, 72, 052103.]) overcomes the limitation of interpolation error of the points from a polar grid to a Cartesian grid in Fourier space. This method has been applied to many tomography fields. Liu et al. proposed an iterative algorithm to reconstruct the refractive index information from data collected by an analyzer-based imaging setup (Liu et al., 2007[Liu, Y. J., Zhu, P. P., Chen, B., Wang, J. Y., Yuan, Q. X., Huang, W. X., Shu, H., Li, E. R., Liu, X. S., Zhang, K., Ming, H. & Wu, Z. Y. (2007). Phys. Med. Biol. 52, L5-13.]). It is effective and accurate compared with the standard analytical method. For discrete tomography, several competing algorithms also have been presented. DART has proved to be a successful algorithm for objects consisting of only a few different compositions with prior knowledge of the grey levels (Batenburg & Sijbers, 2011[Batenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542-2553.]). PDM-DART has been applied in tomography that requires prior knowledge only of the number of grey levels (van Aarle et al., 2012[Aarle, W. van, Batenburg, K. J. & Sijbers, J. (2012). IEEE Trans. Image Process. 21, 4608-4621.]). A partially discrete algebraic reconstruction technique (PDART) could reconstruct and segment the densest nanoparticles and allow the rest of the reconstruction to vary freely (Roelandts et al., 2011[Roelandts, T., Batenburg, K. J. & Sijbers, J. (2011). Proceedings of the 11th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology & Nuclear Medicine (Fully 3D), Potsdam, Germany.]). The BgART algorithm detects the background and keeps the background fixed at a uniform value for every iteration (Messaoudi et al., 2013[Messaoudi, C., Aschman, N., Cunha, M., Oikawa, T., Sorzano, C. O. & Marco, S. (2013). Microsc. Microanal. 19, 1669-1677.]). This facilitates denoising the reconstruction and segmentation of objects without prior knowledge of the materials. The MDART algorithm treats each region as a research object, and the value of each region is calculated by solving the linear system of equations. It reconstructs images with multiple compositions without prior knowledge of the grey values.

5. Conclusions

In this paper, a modified DART (MDART) algorithm for reconstructing an image with multiple grey and discrete regions is presented. MDART focuses on reconstructing images with a limited angle range and without prior knowledge of the grey values. By using each independent region instead of pixels as a research object, the number of unknowns in the algebraic reconstruction linear equations is decreased by several orders of magnitude.

The results show that MDART is a powerful algorithm for compensating for the `missing wedge' effects that present a major challenge for quantitative analysis of an object when using tomography methods. MDART can achieve a higher-quality reconstruction without prior knowledge of the grey levels for both binary and multiple grey phantoms. The number of erroneous pixels in the multiple grey reconstructions using MDART is smaller than the numbers obtained by other reconstructed methods. Although the number of grey levels increases to 50, there is no obvious reduction in the accuracy. Additionally, the grey level values calculated using MDART are precise enough when compared with the original images. Thus, MDART can be used as an alternative method to quantitatively analyze objects with multiple compositions, such as integrated circuit chips, the cathode and anode of a solid oxide fuel cell, marine sediment, shale, nanoparticles, etc. In future work, we intend to optimize the algorithm to make it suitable for different types of experimental data with discrete compositions obtained from soft X-ray microscopy, hard X-ray microscopy, electron microscopy, etc.

Supporting information


Acknowledgements

This work was supported by grants from the Major State Basic Research Development Program of China (973 Program) (No. 2012CB825804), the PhD Programs Foundation of the Ministry of Education of China (No. 2012340213002) and the National Science Foundation of China (No. 11275204, No. 11475175, No. 11405175).

References

First citationAarle, W. van, Batenburg, K. J. & Sijbers, J. (2012). IEEE Trans. Image Process. 21, 4608–4621.  PubMed Google Scholar
First citationAndersen, A. & Kak, A. (1984). Ultrason. Imaging, 6, 81–94.  CrossRef CAS PubMed Google Scholar
First citationBals, S., Batenburg, K. J., Liang, D. D., Lebedev, O., Van Tendeloo, G., Aerts, A., Martens, J. A. & Kirschhock, C. E. A. (2009). J. Am. Chem. Soc. 131, 4769–4773.  CrossRef PubMed CAS Google Scholar
First citationBals, S., Batenburg, K. J., Verbeeck, J., Sijbers, J. & Van Tendeloo, G. (2007). Nano Lett. 7, 3669–3674.  Web of Science CrossRef CAS Google Scholar
First citationBarrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C. & Van der Vorst, H. (1994). Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Philadelphia: Siam.  Google Scholar
First citationBatenburg, K. J. (2005). Discrete Appl. Math. 151, 36–54.  CrossRef Google Scholar
First citationBatenburg, K. & Sijbers, J. (2007). IEEE International Conference on Image Processing (ICIP 2007), pp. IV-133–IV-136. IEEE.  Google Scholar
First citationBatenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542–2553.  Web of Science CrossRef PubMed Google Scholar
First citationBatenburg, K. J., Sijbers, J., Poulsen, H. F. & Knudsen, E. (2010). J. Appl. Cryst. 43, 1464–1473.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationBatenburg, K. J., van Aarle, W. & Sijbers, J. (2011). Pattern Recognit. Lett. 32, 1395–1405.  CrossRef Google Scholar
First citationDuke, E. M., Razi, M., Weston, A., Guttmann, P., Werner, S., Henzler, K., Schneider, G., Tooze, S. A. & Collinson, L. M. (2013). Ultramicroscopy, 143, 77–87.  CrossRef PubMed Google Scholar
First citationHerman, G. T. (2009). Fundamentals of Computerized Tomography: Image Reconstruction from Projections. Springer.  Google Scholar
First citationHerman, G. T. & Kuba, A. (2007). Advances in Discrete Tomography and its Applications. Berlin: Springer.  Google Scholar
First citationKak, A. C. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics.  Google Scholar
First citationKanitpanyacharoen, W., Parkinson, D. Y., De Carlo, F., Marone, F., Stampanoni, M., Mokso, R., MacDowell, A. & Wenk, H.-R. (2013). J. Synchrotron Rad. 20, 172–180.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationLiang, Z., Guan, Y., Liu, G., Bian, R., Zhang, X., Xiong, Y. & Tian, Y. (2013). Proc. SPIE, 8851, 885113.  Google Scholar
First citationLiu, Y. J., Zhu, P. P., Chen, B., Wang, J. Y., Yuan, Q. X., Huang, W. X., Shu, H., Li, E. R., Liu, X. S., Zhang, K., Ming, H. & Wu, Z. Y. (2007). Phys. Med. Biol. 52, L5–13.  CrossRef PubMed CAS Google Scholar
First citationMessaoudi, C., Aschman, N., Cunha, M., Oikawa, T., Sorzano, C. O. & Marco, S. (2013). Microsc. Microanal. 19, 1669–1677.  CrossRef CAS PubMed Google Scholar
First citationMiao, J., Förster, F. & Levi, O. (2005). Phys. Rev. B, 72, 052103.  Web of Science CrossRef Google Scholar
First citationMilne, J. L. & Subramaniam, S. (2009). Nat. Rev. Microbiol. 7, 666–675.  CrossRef PubMed CAS Google Scholar
First citationMueller, K., Yagel, R. & Wheller, J. J. (1999). IEEE Trans. Med. Imaging, 18, 519–537.  CrossRef PubMed CAS Google Scholar
First citationPaige, C. C. & Saunders, M. A. (1982). ACM Trans. Math. Softw. 8, 43–71.  CrossRef Google Scholar
First citationPalenstijn, W. J., Batenburg, K. J. & Sijbers, J. (2013). The ASTRA Tomography Toolbox, in Proceedings of the 13th International Conference on Computational and Mathematical Methods in Science and Engineering (CMMSE 2013), 24–27 June 2013, Cabo de Gata, Almeria, Spain.  Google Scholar
First citationRoelandts, T., Batenburg, K. J. & Sijbers, J. (2011). Proceedings of the 11th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology & Nuclear Medicine (Fully 3D), Potsdam, Germany.  Google Scholar
First citationSchüle, T., Schnörr, C., Weber, S. & Hornegger, J. (2005). Discrete Appl. Math. 151, 229–243.  Google Scholar
First citationSidky, E. Y., Kao, C.-M. & Pan, X. (2006). J. X-ray Sci. Technol. 14, 119–139.  Google Scholar
First citationTian, Y. C., Li, W. J., Chen, J., Liu, L. H., Liu, G., Tkachuk, A., Tian, J. P., Xiong, Y., Gelb, J., Hsu, G. & Yun, W. B. (2008). Rev. Sci. Instrum. 79, 103708.  Web of Science CrossRef PubMed Google Scholar
First citationUramoto, G.-I., Morono, Y., Uematsu, K. & Inagaki, F. (2014). Limnol. Oceanogr. Methods, 12, 469–483.  CrossRef Google Scholar
First citationWeber, S., Nagy, A., Schulle, T., Schnorr, C. & Kuba, A. (2006). Lect. Notes Comput. Sci. 4245, 146–156.  CrossRef Google Scholar
First citationWeber, S., Schule, T., Hornegger, J. & Schnorr, C. (2004). Lect. Notes Comput. Sci. 3322, 38–51.  CrossRef Google Scholar
First citationYang, X., Hofmann, R., Dapp, R., van de Kamp, T., dos Santos Rolo, T., Xiao, X., Moosmann, J., Kashef, J. & Stotzka, R. (2015). Opt. Express, 23, 5368–5387.  CrossRef PubMed Google Scholar
First citationYu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791–2805.  Web of Science CrossRef PubMed Google Scholar
First citationZeng, G. L. (2010). Medical Image Reconstruction: A Conceptual Tutorial. Berlin: Springer.  Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds