A modified discrete tomography for improving the reconstruction of unknown multi-gray-level material in the `missing wedge' situation

Liu, J.; Liang, Z.; Guan, Y.; Wei, W.; Bai, H.; Chen, L.; Liu, G.; Tian, Y.

doi:10.1107/S1600577518013681

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 25| Part 6| November 2018| Pages 1847-1859

https://doi.org/10.1107/S1600577518013681

A modified discrete tomography for improving the reconstruction of unknown multi-gray-level material in the `missing wedge' situation

Jianhong Liu,^a Zhiting Liang,^a Yong Guan,^a ^* Wenbin Wei,^a Haobo Bai,^a Liang Chen,^a Gang Liu ^a and Yangchao Tian ^a

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

Edited by S. M. Heald, Argonne National Laboratory, USA (Received 13 May 2018; accepted 25 September 2018; online 24 October 2018)

Full angular rotational projections cannot always be acquired in tomographic reconstructions because of the limited space in the experimental setup, leading to the `missing wedge' situation. In this paper, a recovering `missing wedge' discrete algebraic reconstruction technique algorithm (rmwDART) has been proposed to solve the `missing wedge' problem and improve the quality of the three-dimensional reconstruction without prior knowledge of the material component's number or the material's values. By using oversegmentation, boundary extraction and mathematical morphological operations, `missing wedge' artifact areas can be located. Then, in the iteration process, by updating the located areas and regions, high-quality reconstructions can be obtained from the simulations, and the reconstructed images based on the rmwDART algorithm can be obtained from soft X-ray nano-computed tomography experiments. The results showed that there is the potential for discrete tomography.

Keywords: discrete algebraic reconstruction technique (DART); missing wedge; over-segmentation; tomography; limited-angle.

1. Introduction

Combined with computed tomography (CT), high-resolution microscopy tools such as the electron microscope (EM) and transmission X-ray microscope (TXM) have been developed to obtain three-dimensional reconstructed structures. However, in many situations, the `missing wedge' problem (Zhuge et al., 2017 ) is inevitable. For instance, due to the morphological limitation of the specimen carrier, tilt series projection images cannot be acquired under full rotation. In the cases of the EM and soft X-ray microscope, mesh grids are used as the specimen carriers. It is difficult to obtain high-tilt-angle projection images because the samples are blocked by the mesh grids when imaging at a high tilt angle (Zhuge et al., 2017). In addition, full-rotation tilt series projections cannot be acquired because of radiation damage due to the longer exposure times for full-rotational tilt series projections (Dent et al., 2014 ). Thus, it is impossible to collect high-tilt projection images during CT data acquisition, resulting in `missing wedge' artifacts in the subsequent reconstruction process.

Many reconstruction methods have been developed to solve the artifacts caused by the `missing wedge' problem (Guan & Gordon, 1994 ; Panin et al., 1999 ; Liang et al., 2013 ). Iterative reconstruction algorithms are commonly used to address projection data in an angle-limited range, and they can slightly reduce the `missing wedge' problem. The total variation (TV) based method is able to provide relatively high fidelity reconstructions on the basis of an iterative algorithm (Liang et al., 2013; Panin et al., 1999). Recently, the Discrete Algebraic Reconstruction Technique (DART) was reported and shown to be capable of calculating very high quality reconstructions using very few projections (Batenburg & Sijbers, 2011 ; Batenburg et al., 2010 ). The DART method with total variation regularization (Zhuge et al., 2016 , 2017) and soft constraints (Bleichrodt et al., 2014 ; van Aarle et al., 2012 ) was proposed to achieve high-quality reconstructions. However, two types of prior knowledge are still required in order to obtain high-quality reconstructions with the DART method. One is that a sample consists of a limited number of materials, and the other is that sparse boundaries exist between materials (Batenburg et al., 2009 , 2010; Batenburg & Sijbers, 2007 , 2011; Nemeth, 2015 ). Liang et al. proposed a method (mDART) based on DART to improve the `missing wedge' situation without the need of prior knowledge and acquired high-quality reconstructions, but it requires manual segmentation (Liang et al., 2016 ).

The typical DART method contains a segmentation process. In general, to acquire an accurate segmentation using the DART method, prior knowledge of the material component's number and the material's values are necessary. In this paper, an advanced DART algorithm (rmwDART) is proposed to solve the `missing wedge' situation by computing the multi-gray-level tomographic reconstruction without the use of any prior knowledge. The automatic oversegmentation strategy was applied to the reconstructed image for each iteration because oversegmentation does not need prior knowledge of the sample for segmentation since it just oversegments the image. Generally, for image processing, oversegmentation will produce inappropriate segmentation because it will oversegment the artifact areas into too many small regions resulting in uncertainty. In this paper, this `disadvantageous' feature of oversegmentation becomes an exploitable advantage, and a method that combines the oversegmentation strategy, boundary extraction and mathematical morphological operations is proposed to locate the artifact areas. Thus, this method does not require knowledge of the values of different materials and the number of gray value levels. The recovery of these areas was then achieved by applying the simultaneous algebraic reconstruction technique (SART). This method helps to improve the accuracy of the artifact areas and to further improve the reconstruction quality. The simulation and sample experiment's details are described in the following sections.

2. Algorithm

The main reconstruction program of rmwDART contains four main steps in an iteration loop: the initialization of the total variation-based simultaneous algebraic reconstruction technique (SART-TV), the oversegmentation, the redistribution of regions and the calculation of the regional values, and the positioning and updating of artifact areas. A flowchart of rmwDART is shown in Fig. 1, where `i' is the iteration number.

Figure 1
Flowchart of rmwDART.

2.1. SART-TV reconstruction initialization

Typically, an initial reconstruction is required as input data before the main iterative reconstruction process of rmwDART proceeds. Usually, algebraic reconstruction methods can be used to conduct the initial reconstruction. In this paper, the SART-TV method was used to conduct the initial reconstruction. Let x = (x_j) ∈ Rⁿ denote a two-dimensional image, where n is the number of the pixels. The line integral of the linear absorption coefficient p (Natterer, 1986 ) from the measured projection data can be represented by

$[{\bf{p}} = {\bf{Wx}}\eqno(1)]$

where W = (w_ij) is the m × n projection matrix, and m is the total number of the measured detectors.

Basically, the CT reconstruction process can be interpreted by computing the solution x to $[\min \left\| {{\bf{Wx - p}}} \right\|]$ . The SART method can be expressed as follows,

$[x_j^{(k + 1)} = x_j^{(k)} + \lambda {{\sum\nolimits_i^p {\left [{{{{w_{ij}}\left({{p_i} - \sum\nolimits_h^n {{w_{ih}}\,x_h^{(k)}} } \right)} \Big/ {\sum\nolimits_h^n {{w_{ih}}} }}} \right]} } \over {\sum\nolimits_i^p {{w_{ij}}} }}, \eqno(2)]$

where x_j^(k) denotes the jth pixel of the current reconstruction, x_j^(k+1) denotes that of the next iteration, λ is the relaxation parameter, and p is the number of total detectors at each angle.

Generally, TVm-based methods can be expressed as follows,

$[\min {\left\| {\bf{x}} \right\|_{{\rm{TV}}}},\quad {\rm such\,\,that} \quad {\rm{ }}{\bf{p}} = {\bf{Wx}}. \eqno(3)]$

Here, TV is the l₁-norm of the gradient image, and it can be expressed as

$[{\left\| {\bf{x}} \right\|_{{\rm{TV}}}} = {\sum\limits_{i,j} {\left \{{{{\left[{{\bf{x}}(i + 1,j) - {\bf{x}}(i,j)} \right]}^2} + {{\left[{{\bf{x}}(i,j + 1) - {\bf{x}}(i,j)} \right]}^2}} \right\}} ^{1/2}}. \eqno(4)]$

Thus, the SART-TV method can be expressed as

$[{\hat {\bf x}} = \arg \mathop {\min }\limits_x \big\{ \left\| {{\bf{Wx - p}}} \right\|_2^2 \,\,+\,\, \gamma {\mathop{\rm TV}\nolimits} ({\bf{x}})\big\}, \eqno(5)]$

where TV(…) represents the TVm norm, and γ is the weight for controlling the trade-off between the SART and TVm norm.

2.2. Oversegmentation

Since the automatic oversegmentation strategy does not need prior knowledge of the material component number or the linear absorption coefficient values of the material, it was applied to the image reconstruction in each iteration. After the input image is set or the artifact areas are updated, oversegmentation will be applied to the image reconstruction. This means that many small regions will be segmented in the `missing wedge' areas. Because the values of the missing wedge areas are smooth instead of discrete in angle-limited reconstructions, the oversegmentation method would segment these areas into many small regions, which is necessary. This is because only these segmented small regions will be retained after subsequently applying the boundary extraction and mathematical morphological operations. Therefore, the oversegmentation method is very helpful in locating the missing wedge areas precisely. To reach the goal of oversegmentation, the thresholding segmentation method was chosen since it has fast processing and adequately meets the requirements. For thresholding segmentation, the thresholds must be selected. The procedures for computing the image's thresholds are as follows,

$[\eqalign{ h_s(x) & = {1\over3}\left[h(x)+\textstyle\sum\limits_{i\,=\,-1/2}^{1/2} h(x+i)\right]_{\vphantom{\Big|}}, \cr h(x+{1/2}) & = h(x)+\big[h(x+1)-h(x)\big]/2, } \eqno(6)]$

$[y(x) = \big\{ \max[h_{\rm{s}}(i)] \,|\, x-w_1 \le i \le x + w_1 \big\}, \eqno(7)]$

$[y_1(x) = \left\{ \matrix{ y(x),\hfill & \!\!\!{\rm{if}}\,\, y^{\,\prime}(i)=0, \forall \,\,i \in \left[x-w_1/2,x+w_1/2\right],\hfill \cr 0,\hfill & \!\!\!{\rm{otherwise}},\hfill } \right. \eqno(8)]$

$[\big\{x_{\rm{p}}|y_1\left(x_{\rm{p}}\right) > c_{\rm{thr}},\,\,y_1\left(x_{\rm{p}}\right) = h_{\rm{s}}\left(x_{\rm{p}}\right)\big\}, \eqno(9)]$

$[\big\{ x_{\rm{v}}|\,x_{\rm{v}}(i) = \mathop{\min}\limits_{x\,\in\,\left[x_{\rm{p}}(i),\,x_{\rm{p}}(i+1)\right]} \left[h_{\rm{s}}(x)\right] \big\}, \eqno(10)]$

h(x) is the histogram of an image, where x is the grayscale value and h(x) is the count number of the grayscale value x. Equation (6) smooths the histogram from h(x), and the smoothed histogram is denoted as h_s(x). Equations (7)–(9) represent the process of peak-picking from the smoothed histogram h_s(x). w₁ is a very important parameter in rmwDART because 2w₁ represents the peak-picking resolution in the peak-picking process. c_thr is a non-negative constant for filtering the small peaks, which could be considered the noise counts. In this paper, c_thr is set to 2.5. {x_p} is the point set of the peak locations in the histogram. According to the point set {x_p} and the smoothed histogram h_s(x), the locations of the minimal counts between two adjacent peak locations in {x_p} are considered the point set of the threshold values {x_v} by equation (6). Then, thresholding segmentation is applied on the image with the threshold values {x_v}. w₁ directly determines the degree of oversegmentation. By choosing a very small peak-detection parameter 2w₁, it can make the process of peak-picking more sensitive, thereby detecting more peaks and thresholds, which finally leads to the oversegmentation result.

2.3. Redistribution of regions and calculation of the regional values

After the image is segmented into a regionalized image by oversegmentation, the regional connectivity and the regional values (corresponding to the linear absorption coefficient of the region) of the image are calculated, and the image regions are redistributed. A flowchart of this step is shown in Fig. 2, and it contains three substeps: the evaluation of the regional connectivity, the redistribution of the image regions, and the calculation of the regional value. The regional connectivity is a state in which adjacent image regions are connected. It represents the connectivity of two adjacent regions. Let the connectivity threshold T_c,n represent the detection resolution of two regional values. If the difference between two adjacent regional values is smaller than T_c,n, then these two regions would still be regarded as the same substance and should be connected into one region; otherwise, these two regions are regarded as two different regions (corresponding to two different substances). In this paper, the four-connected neighborhood is used to determine whether two regions are spatially adjacent. Then, according to the regional connectivity, the image regions are redistributed. Therefore, T_c,n is another important parameter in rmwDART. After calculating the regional connectivity of the image, the image regions are redistributed. Let x_s represent the redistributed regional image, and the unknown regional values are then computed by solving the linear equations Wx_s = p using the LSQR (sparse equations and least squares) method (Paige & Saunders, 1982a ,b ).

Figure 2
Flowchart of the `Redistribution of regions and the calculation of the regional values' step.

The regional connectivity is based on the values of adjacent regions. Theoretically, after calculating the regional values, new regional values are generated. This generation means that the regional connectivity probably also changed, so the regional connectivity and the regional distribution should be renewed. However, after the renewal of the regional connectivity and the regional distribution, the image regions are changed, and so the regional values should also be renewed. Therefore, to acquire a stable distribution of the image regions and their regional values, this process of calculating the regional connectivity and regional values will be repeated more than once. However, after the first calculation, the changes in the regional connectivity and regional values are small because the oversegmentation will overly segment each region, which can avoid misconnections. Therefore, only a very few iterations are needed to reach a stable result, and five sub-iterations are used in this paper. Since the regional distribution and the regional values will probably change and become stable in the sub-iterations, in order to avoid misconnections between adjacent regions due to the relatively unstable regional values in the initial sub-iteration, T_c,n is set to be relatively small in the initial sub-iteration, and then it gradually increases in the subsequent sub-iterations. A discussion about this is detailed in §3.2.2. In this paper, the connectivity threshold is set to be the following,

$[\eqalignno{ T_{{\rm{c}},n} & = \big[ {T_{{\rm{c}},1}},{T_{{\rm{c}},2}},{T_{{\rm{c}},3}},{T_{{\rm{c}},4}},{T_{{\rm{c}},5}} \big] \cr & = \big[1,1.5,2,3,4\big]\times10^{-3}, \,\, n \in [1,5], &(11)}]$

and the four-connected neighborhood was used to judge the spatial pixel connectivity. Because T_c,n represents the detection resolution of the regional values, the parameters of the five sub-iterations and T_c,n can be adjusted by users. However, as T_c,n decreases, the time consumption for solving the linear equations will also increase. Therefore, in the simulation experiment of this paper, to balance the detection resolution and the time consumption, T_c,n was set to be 0.001, 0.0015, 0.002, 0.003 and 0.004.

2.4. Update on the located artifact areas

After the regional connectivity and regional values are calculated, if the stop criterion is not yet met, the next step is the positioning and updating of the artifact areas. A flowchart of this step is shown in Fig. 3, and contains five substeps: boundary extraction, erosion operation, dilation operation, SART updating and Gaussian smoothing.

Figure 3
Flowchart of the `Positioning and updating of artifact areas' step.

First, the differential gradient is used in the regionalized image to extract the boundaries of all regions. Then, mathematical morphology operations are applied on the extracted boundaries. The mathematical morphology operations used in this paper contain the erosion and dilation operations. Oversegmentation would segment many small regions that only occupy a few pixels or less in the `missing wedge' artifact areas. When operating the erosion operator on the boundaries, some of the boundaries in the `missing wedge' artifact areas would be retained, while the other boundaries in the non-artifact areas are eroded. This procedure helps to locate the positions of `missing wedge' artifact areas and subsequently focuses on recovering these areas. Furthermore, it does not affect the other boundaries during the image reconstruction process. The utilized erosion operator is operator {1} in the `Positioning and updating of artifact areas' in Fig. 3. Because the image boundaries are eroded, the dilation operator is then used on the eroded image boundaries to recover the size of the artifact areas.

In this paper, the dilation operators {1}, {2 3} and {4 0} were used in rmwDART, as shown in Fig. 4. Typically, using the dilation operator {1} can roughly recover the artifact areas. However, the dilation operators {2 3} and {4 0} can recover finer artifact areas in two opposite directions and improve the efficiency. Therefore, for the dilation operators {2 3} and {4 0}, two iterations are applied for each (after the one-time iteration for the operator {1}).

Figure 4
Three different dilation operators. The gray arrows represent the different directions of the dilation effect by different dilation operators.

The SART method is then applied to update the values of the located artifact areas while the values of the non-artifact areas remain unchanged. Subsequently, a Gaussian filter with a standard deviation of 0.5 pixel is applied on the updated artifact areas. The filtered updated artifact areas are then recombined with the non-artifact areas into an output image. Because the updated areas do not include the areas with clear boundaries, the reconstruction error will not affect the areas with clear boundaries. This means that the quality of the whole reconstruction can be improved by reducing the reconstruction error.

When the remainder of i divided by 5 equals {1 2 4}, the next iteration starts from the oversegmentation step, and, when it equals {3 0}, the next iteration starts from the connectivity calculation step. This process means that when a new iteration for the next dilation operator begins, an oversegmentation step must be applied. In the case with two iterations, the oversegmentation procedure is used only one time. When the stop criterion is met, the iteration process stops and the resulting reconstructed image is output. In this paper the maximum iteration number was set to be the stop criterion.

3. Result and discussion

3.1. Simulation results and analysis

Five simulated images were used to test rmwDART and are shown in Fig. 5. The size of all images is 256 × 256 pixels, and they are multi-gray phantoms. Phantom 1 is the classic Shepp–Logan phantom. Phantoms 2 and 3 contains 14 and 15 gray levels, respectively. Phantoms 4 and 5 are composed of many homogeneous circles, each of which has a different value in the image, and they contain 50 and 101 gray levels, respectively. Simulation experiments on phantoms 1–5 were also tested in the work of Liang et al. (2016).

Figure 5
Five different images used to test rmwDART. (a) Phantom 1 is the classic Shepp–Logan phantom containing seven gray levels. (b) Phantom 2 is the modified phantom with 14 gray levels. (c) Phantom 3 has 15 gray levels. (d) Phantom 4 is composed of many homogeneous circles with 50 gray levels. (e) Phantom 5 is composed of many homogeneous circles with 101 gray levels.

The projections of phantoms 1–5 were acquired. The projection length of the projections is 367 pixels, which is about the length of the diagonal line of the image. The angular ranges of the projections are (a) 0–138°, (b) 0–138°, (c) 0–120°, (d) 0–90° and (e) 0–90°. These angular ranges were chosen to demonstrate the maximum capability of rmwDART. The SART-TV reconstructions of these five images are shown in Fig. 6. The total number of iterations for the SART-TV reconstruction was 500 times.

Figure 6
SART-TV reconstruction images of the five test images. The number of iterations for the SART-TV reconstruction was 500. The angular ranges for the projections of phantoms 1–5 are as follows: (a) 0–138°, (b) 0–138°, (c) 0–120°, (d) 0–90° and (e) 0–90°, respectively.

Fig. 5(a) is used to explain rmwDART's procedures in one loop. First, the SART-TV reconstruction initialization is shown in Fig. 6(a). The histogram h(x) of Fig. 6(a) is shown in Fig. 7(a). Fig. 7(b) is the smoothed histogram h_s(x) from h(x). Then, h_s(x) is scanned using a scanning window with a width of 2w₁ [see equation (7)], and the maximum counts in the scanning window were recorded in y(x). The result of y(x) is shown in Fig. 7(c). As shown in Fig. 7(c), the count platforms were generated on the isolated count peaks after the scanning. Then, y(x) was scanned by a smaller scanning window with a smaller width of w₁ in equation (8), and only the platforms with widths no smaller than w₁ were kept and the width of the platforms decreased in y₁(x). The resulting smaller platforms of y₁(x) are shown in Fig. 7(d). In equation (9), after filtering the noise counts, the peak locations {x_p} where h_s(x) and y₁(x) coincide were acquired. As shown in Fig. 7(e), the red line represents the peak locations {x_p}, and h_s(x) and y₁(x) coincide at the same place. If the distance between two isolated peaks is smaller than the peak-picking resolution, the platform of the peak with the smaller count in y₁(x) will not coincide with the peak itself in h_s(x), which means that this smaller peak will not be picked. When the distance between two isolated peaks is larger than the peak-picking resolution, these two isolated peaks in h_s(x) will coincide with their own platforms in y₁(x), and thus these two isolated peaks are picked. Then, the locations of the minimal counts between two adjacent peak locations in {x_p} were acquired using equation (10). The threshold values {x_v} are these minimal count locations, and the result of {x_v} is shown in Fig. 7(f). In this segmentation procedure, the peak-picking resolution 2w₁ is set to be 0.5% of the grayscale value range in the histogram. The image is then segmented with these threshold values {x_v}. After the oversegmentation step, the segmented image is shown in Fig. 8(a). The image was segmented into 871 regions.

Figure 7
(a) Histogram h(x) of Fig. 6(a)

. (b) Smoothed histogram h_s(x) from h(x). (c) Maximum counts y(x) in the scanning window. (d) Smaller platforms y₁(x). (e) Peak locations {x_p}, where the red line represents the peak locations {x_p}, the purple line represents y₁(x) and the cyan curve represents h_s(x). (f) Threshold values {x_v}, where the red line represents the threshold values {x_v} and the cyan curve represents h_s(x).

Figure 8
Procedures of rmwDART in one loop. (a) Resulting segmentation of the input reconstruction of Fig. 6(a)

. (b) Resulting boundary extraction from (a). (c) Calculated image after the `Redistribution of regions and the calculation of the regional values' step. (d) Extracted boundaries image from (c). (e) Image of the areas that were extracted after applying the erosion and dilation operations on (d). (f) Resulting image of one loop iteration of rwmDART.

The boundaries extracted from the segmented regionalized image are shown in Fig. 8(b). It can be seen that many small regions assembled in the `missing wedge' artifact areas. Then, in the `Redistribution of regions and the calculation of the regional values' step, the connectivity thresholds were set by equation (7), and the calculated image after this step is shown in Fig. 8(c), which is much better than that in Fig. 6(a). However, the boundaries of some regions still need to be improved. In the procedure of calculating the regional values, 300 iterations of the LSQR method were applied. Then, in the next step, the outline of the extracted boundaries from the calculated image (Fig. 8c) was acquired, as shown in Fig. 8(d). After applying the erosion and dilation operations on the extracted outline of the boundaries, the areas that need to be updated are shown in Fig. 8(e), which were mostly coincident with the `missing wedge' artifact areas. Then, 15 iterations of the SART method and Gaussian smoothing were applied to update the values of these areas. Finally, the higher-quality reconstruction was obtained and is shown in Fig. 8(f). At this point, one loop iteration of rmwDART had been completed.

To demonstrate the effect of three different dilation operators, three resulting images after different dilation operations (by the dilation operators in Fig. 4) were used as examples, shown in Fig. 9. It can be seen that the dilated areas in Figs. 9(b) and 9(c) are smaller than that in Fig. 9(a), but they are two different parts of the dilated areas in Fig. 9(a) in two opposite directions. Using the dilation operators in Figs. 4(b) and 4(c) can help the convergence by restoring the artifacts in two opposite directions.

Figure 9
Explanation of the three different dilation operators. Panels (a), (b) and (c) are the resulting images from the dilation effect after the dilation operation by the operator in (a) Fig. 4(a)

, (b) Fig. 4(b)

and (c) Fig. 4(c)

, respectively.

After 30 iterations of rmwDART, the final reconstructions of the five test images are shown in Figs. 10(k)–10(o). In comparison, 100 iterations of mDART and DART were also applied on the projections of phantoms 1–5. The segmented images for the segmentation-input of mDART (Liang et al., 2016) were generated by the oversegmentation method proposed in this paper. The resulting reconstructions of DART and mDART are shown in Figs. 10(a)–10(e) and Figs. 10(f)–10(j).

Figure 10
Final reconstructions of five test images by different algorithms. (a)–(e) Resulting reconstructions of phantoms 1–5 by DART. (f)–(j) Resulting reconstructions of phantoms 1–5 by mDART. (k)–(o) Resulting reconstructions of phantoms 1–5 by rwmDART. The angular range for the projections of phantoms 1–5 are as follows: (a, f, k) 0–138°; (b, g, l) 0–138°; (c, h, m) 0–120°; (d, i, n) 0–90°; (e, j, o) 0–90°, respectively.

To quantitatively evaluate the quality of the reconstructed images in the simulation test, Liang et al. proposed using the following equation (Liang et al., 2016),

$[K = \left| {\left\{ {i\left| {\left [\,{\left| {x_i^{\,r} - {x_i}} \right| \gt \max (\omega \. \Delta \rho, 0.003)} \right]} \right.} \right\}} \right|, \eqno(12)]$

where K is the total number of error pixels, x_i^r is the pixel value of the reconstructed image, x_i is the pixel value of the original image, Δρ is the original images' minimum interval between two values, and ω is an indicator that ranges from 0 to 1. In this case, ω was set as 0.03, which is the same as in Liang et al.'s work (Liang et al., 2016). As shown in Fig. 11, the pixel errors K of Figs. 10(a)–10(e) are 3938, 4174, 7255, 10493 and 16268, respectively; the pixel errors K of Figs. 10(f)–10(j) are 4835, 6280, 17112, 4692 and 8546, respectively; and the pixel errors K of Figs. 10(k)–10(o) are 3, 1, 2, 5 and 1, respectively. As a result, in comparison with DART and mDART, rmwDART calculates better reconstructions of the simulated images of Fig. 3 in `missing wedge' situations. It should also be noted that mDART requires an appropriate manual segmentation for the input, and the segmentation result from the oversegmentation method is probably not suitable for the mDART method. This may be the reason why the mDART method did not perform well in the above simulation.

Figure 11
Pixel error K of the resulting reconstructions of phantoms 1–5 by rmwDART, mDART and DART. The angular range of the projections of phantoms 1–5 for the reconstructions were 0–138°, 0–138°, 0–120°, 0–90° and 0–90°, respectively.

As mDART may require an appropriate segmentation for the input image, experiments on mDART with the appropriate manual segmentation were assessed using the same phantoms 1–5 as the simulated images, as in Liang et al.'s work (Liang et al., 2016). The angular ranges for those mDART experiments on phantoms 1–5 were (a) 0–147°, (b) 0–145°, (c) 0–130°, (d) 0–135° and (e) 0–130°, respectively. In addition, the pixel errors K of these experiments on phantoms 1–5 were 212, 546, 161, 145 and 1580, respectively. It should be noted that the projection range used in Liang et al.'s work was larger than that of this paper. This shows that the above simulation results of rmwDART are better than those of mDART with the appropriate segmentation and DART.

The convergence rates of rwmDART for the resulting reconstructions in Figs. 10(k)–10(o) are shown in Fig. 12; in this figure the horizontal axis represents the iteration number and the vertical axis represents the pixel error K. As shown in Fig. 12, the five reconstruction simulations converged within 15 iterations. The time consumptions for the rwmDART reconstructions of phantoms 1–5 are shown in Table 1. In this study, a computer equipped with an Intel i5-4570 CPU and 16 GB RAM memory was used to compute the data. The time consumption for the reconstructions of phantoms 2 and 3 were relatively large, and the convergence rates also show that the reconstructions of phantoms 2 and 3 converged relatively late. It shows that, when the error of the iterated reconstruction is large, rmwDART consumes more time for its calculations. Compared with the time consumption of 15 iterations and 30 iterations, even if the reconstructions converge, rmwDART still needs significant time to calculate each iteration.

Table 1
Time consumption (in seconds) for the rwmDART reconstructions of phantoms 1–5

	Phantom
	1	2	3	4	5
15 iterations	383.05	686.58	574.75	284.27	341.40
30 iterations	722.22	1003.91	874.00	543.95	649.89

Figure 12
Convergence rates of rwmDART for the resulting reconstructions in Figs. 10(k)–10(o)

3.2. Parameters discussion

In rmwDART, w₁ and T_c,n are the most important parameters since they directly determine the quality of the reconstruction. Therefore, the effect of different w₁ and T_c,n on rmwDART are tested. Phantom 2 is used as an example to test the changed parameters, and the angular projection range and the projection length are the same as in §3.1. The SART-TV reconstruction in Fig. 6(b) was used as the initialization.

3.2.1. Experiments on different w₁

2w₁ represents the peak-picking resolution in the peak-picking process. Therefore, w₁ directly determines the degree of oversegmentation. The smaller that 2w₁ is set, the higher the degree of oversegmentation. In the following rmwDART reconstructions for phantom 2, 2w₁ was set to be 0.5%, 2%, 5% and 15%, and T_c,n is the same as that in §3.1. After applying the oversegmentation method on Fig. 6(b), the resulting segmentation images are shown in Fig. 13, and the numbers of regions of these resulting segmentation images are shown in Table 2. From the results of Fig. 13 and Table 2, as 2w₁ decreases, Fig. 6(b) is segmented into more regions, and the degree of oversegmentation becomes higher. As shown in Fig. 13(d), since the 2w₁ of 15% is relatively large, misconnections appear in the resulting segmentation.

Table 2
Number of regions of the resulting segmentations with different 2w₁

	2w₁
	0.5%	2%	5%	15%
Number of regions	1364	378	181	72

Figure 13
Different segmentation effects with different 2w₁. (a) 2w₁ is 0.5%. (b) 2w₁ is 2%. (c) 2w₁ is 5%. (d) 2w₁ is 15%.

In this study, 90 iterations were used for these rmwDART reconstructions. The convergence rates for the rwmDART reconstruction with different 2w₁ are shown in Fig. 14(a), and the time consumption of each iteration for the rwmDART reconstruction with different 2w₁ is shown in Fig. 14(b). The pixel errors K of the resulting reconstructions and the total time consumption of the 90-iteration reconstruction are shown in Table 3. As shown in Fig. 14(a), when 2w₁ is small enough, the degree of oversegmentation is high enough to segment the image, and then a high-quality reconstruction can be obtained. For 2w₁ = 0.5%, 2% and 5%, the results converged in ten iterations and obtained very high quality reconstructions. When 2w₁ is not small enough, misconnections appear in the segmentation, and the reconstruction error can be very large. For 2w₁ = 15%, misconnections appear in Fig. 13(d), the reconstruction has a large error until the 89th iteration, and it still does not have a small error at the 90th iteration. For 2w₁ = 0.5%, compared with 2w₁ = 2% and 5%, the reconstruction needs a few more iterations to converge. This outcome occurs because, when 2w₁ is very small, many regions are segmented, and the computer needs more time and iterations to redistribute the regions and calculate the regional values. For practical applications, a very small 2w₁ is recommended, which will ensure the high degree of oversegmentation, especially when the computer's processing capabilities are high enough.

Table 3
Pixel errors K and total time consumption of the resulting reconstructions with different 2w₁

	2w₁
	0.5%	2%	5%	15%
Pixel error K	1	1	1	152
Total time consumption (s)	2287.89	2801.25	2790.62	2215.33

Figure 14
(a) Convergence rates for the rwmDART reconstructions with different 2w₁. (b) Time consumption of each iteration for the rwmDART reconstructions with different 2w₁.

3.2.2. Experiments on different T_c,n

T_c,n is the connectivity threshold representing the detection resolution of the regional value. If T_c,n is too large, two adjacent regions that belong to two different substances may be misconnected. Conversely, if T_c,n is too small, it may take more iterations or a very long computational time to connect two adjacent regions that belong to the same substance. In the following rmwDART reconstruction for phantom 2, 2w₁ is the same as that in §3.1.

First, the number of sub-iterations for the redistribution of regions is discussed. The sub-iteration number is set to 8 and T_c,n is set to T_c,n = [2, 2, 2, 2, 2, 2, 2, 2] × 10⁻³. After the oversegmentation step, the numbers of regions in each sub-iteration for the redistribution of regions are shown in Fig. 15. As shown in Fig. 15, after three sub-iterations, the change in the number of regions becomes very small, and the number of regions becomes stable. Thus, setting the number of sub-iterations for the redistribution of regions to five is rational.

Figure 15
Number of regions in each sub-iteration for the redistribution of regions.

Second, the same T_c,n in the sub-iterations is discussed. T_c,n is set to be the same value in the five sub-iterations, and the settings for T_c,n are shown in Table 4. Forty iterations were used for the rmwDART reconstruction. The pixel errors K and the total time consumption of the resulting reconstructions with different T_c,n are also shown in Table 4. The convergence rates for the rmwDART reconstruction are shown in Fig. 16(a). The time consumption of each iteration for the rwmDART reconstruction is shown in Fig. 16(b). As shown in Fig. 16(a), the resulting reconstructions with T_c,n = [5, 5, 5, 5, 5] × 10⁻³ and T_c,n = [8, 8, 8, 8, 8] × 10⁻³ failed to converge, and the other conditions of T_c,n converged within 35 iterations. When T_c,n is relatively large, misconnections may appear in the computing process of the redistribution of regions. Especially in the initial sub-iteration for the redistribution of regions, the regional distribution is probably unstable, and this could lead to a large error in the resulting reconstruction. For the resulting reconstructions that converged, as T_c,n decreases, the resulting reconstruction needs more iterations to converge, and the time consumption of the computations also becomes greater.

Table 4
Pixel errors K and total time consumptions of the resulting reconstructions with different T_c,n

	Pixel error K	Total time consumption (s)
T_c,n = [2, 2, 2, 2, 2] × 10⁻³	1	1980.83
T_c,n = [2, 2, 2, 2, 2] × 10⁻⁴	1	6815.13
T_c,n = [2, 2, 2, 2, 2] × 10⁻⁵	1	10034.65
T_c,n = [5, 5, 5, 5, 5] × 10⁻³	65356	3245.41
T_c,n = [5, 5, 5, 5, 5] × 10⁻⁴	1	1755.67
T_c,n = [5, 5, 5, 5, 5] × 10⁻⁵	1	6313.53
T_c,n = [8, 8, 8, 8, 8] × 10⁻³	65444	2655.37
T_c,n = [8, 8, 8, 8, 8] × 10⁻⁴	1	1514.18
T_c,n = [8, 8, 8, 8, 8] × 10⁻⁵	1	5178.88

Figure 16
(a) Convergence rates for the rwmDART reconstruction with different T_c,n. (b) Time consumption of each iteration for the rwmDART reconstruction with different T_c,n.

Thirdly, different T_c,n in different sub-iterations are discussed. As the misconnections may appear in the initial sub-iteration for the redistribution of regions, T_c,n is set to be relatively small in the initial sub-iteration, and it gradually increases in the subsequent sub-iterations. The settings for T_c,n are shown in Table 5. In this study, 60 iterations were used for the rmwDART reconstructions. The pixel errors K and the total time consumptions of the resulting reconstructions with different T_c,n in different sub-iteration are also shown in Table 5. The convergence rates for the rmwDART reconstruction are shown in Fig. 17(a). The time consumption of each iteration for the rwmDART reconstruction is shown in Fig. 17(b). Compared with the resulting reconstructions with T_c,n = [5, 5, 5, 5, 5] × 10⁻³ and T_c,n = [8, 8, 8, 8, 8] × 10⁻³, which failed to converge and resulted in large errors, the resulting reconstructions with T_c,n = [1, 2, 3, 4, 5] × 10⁻³ and T_c,n = [2, 3, 4, 6, 8] × 10⁻³ successfully converged and had very high quality results. This shows that setting a relatively small T_c,n in the initial sub-iteration can improve rmwDART reconstructions. For the reconstructions T_c,n = [4, 8, 12, 16, 20] × 10⁻⁴ and T_c,n = [2, 3, 4, 6, 8] × 10⁻⁵, as T_c,n increases, the resulting reconstruction needs more iterations to converge, and the time consumption for the computations also increases. For practical applications, a very small T_c,n is recommended, which can avoid the misconnection during the redistribution of regions, especially when the computer's computational capabilities are high enough.

Table 5
Pixel errors K and total time consumptions of the resulting reconstructions with different T_c,n in different sub-iterations

	Pixel error K	Total time consumption (s)
T_c,n = [2, 3, 4, 6, 8] × 10⁻³	16	3625.39
T_c,n = [1, 2, 3, 4, 5] × 10⁻³	1	1715.66
T_c,n = [4, 8, 12, 16, 20] × 10⁻⁴	1	1705.33
T_c,n = [2, 3, 4, 6, 8] × 10⁻⁵	1	7435.25

Figure 17
(a) Convergence rates for the rwmDART reconstruction with different T_c,n in different sub-iterations. (b) Time consumption of each iteration for the rwmDART reconstruction with different T_c,n in different sub-iterations.

3.3. Experiments on conditions with noise

To investigate the noise resistance capability of rmwDART, phantoms 1 and 3 were used as examples to test the conditions with noise. In this study, Poisson noise was introduced in the projection data of phantoms 1 and 3, and 30 dB, 25 dB and 20 dB signal-to-noise ratio (SNR) were selected for the Poisson noise. The angular projection ranges of phantoms 1 and 3, the projection length, 2w₁ and T_c,n are the same as in §3.1. In comparison, mDART and DART were also applied on these noisy projections, the segmentation input for mDART (Liang et al., 2016) was generated by the oversegmentation method proposed in this paper, and the number of iterations for DART, mDART and rmwDART were 120. The resulting reconstructions of phantom 1 with different SNRs are shown in Fig. 18, and the resulting reconstructions of phantom 3 with different SNRs are shown in Fig. 19. The pixel errors K of phantoms 1 and 3 with different SNRs are shown in Figs. 20(a) and 20(b), respectively. The convergence rates and the time consumption of each iteration for these rmwDART reconstructions are shown in Figs. 21(a) and 21(b), respectively. As shown in Fig. 20, the pixel errors K in the resulting reconstructions with 20 dB, 25 dB and 30 dB SNR by rmwDART are smaller than those of mDART and DART. With respect to the mDART results, the resulting reconstruction in Fig. 19(d) have relatively serious distortions. It should also be noted that mDART requires an appropriate manual segmentation for inputs, and the oversegmentation method is probably not suitable for mDART reconstruction. With respect to the rmwDART and DART results, as the SNR of the projections decreases, the pixel error K in the resulting reconstructions also becomes larger because the noise may be amplified by the gradient operator. For the reconstruction of phantom 3 with a 20 dB SNR, although the pixel error K of the reconstruction by rwmDART is smaller than that by DART, the difference between them is not significant. As shown in Fig. 21(a), the fluctuation is still relatively large after convergence, especially when the SNR is relatively low. As shown in Fig. 21(b), as the SNR of the projections decreases, rmwDART needs more computational time, because the noise may create smaller single-pixel regions, which will significantly increase the time consumption.

Figure 18
Resulting reconstructions of phantom 1 with 20 dB, 25 dB and 30 dB SNR by different algorithms. (a)–(c) Resulting reconstructions using DART. (d)–(f) Resulting reconstructions using mDART. (g)–(i) Resulting reconstructions using rmwDART. The angular range for the projections of phantom 1 is 0–138°.

Figure 19
Resulting reconstructions of phantom 3 with 20 dB, 25 dB and 30 dB SNR by different algorithms. (a)–(c) Resulting reconstructions using DART. (d)–(f) Resulting reconstructions using mDART. (g)–(i) Resulting reconstructions using rmwDART. The angular range for the projections of phantom 3 is 0–120°.

Figure 20
Pixel error K of the resulting reconstructions of phantoms 1 and 3 with 20 dB, 25 dB and 30 dB SNR by different algorithms. (a) Pixel error K of the resulting reconstructions of phantom 1. (b) Pixel error K of the resulting reconstructions of phantom 3.

Figure 21
(a) Convergence rates for the rwmDART reconstruction of phantoms 1 and 3 with 20 dB, 25 dB and 30 dB SNR. (b) Time consumption of each iteration for the rwmDART reconstruction of phantoms 1 and 3 with 20 dB, 25 dB and 30 dB SNR.

4. Experiments on a solid oxide fuel cell (SOFC) anode

A SOFC anode was imaged at beamline 4W1A at the Beijing Synchrotron Radiation Facility (BSRF), China (Yuan et al., 2012 ). 180 projection images with an angular range from −89° to 90° were acquired by 8.4 keV (above the Ni K-edge) X-ray illumination. The pixel resolution of the image is 14.6 nm [for displaying the value of the resulting reconstruction in SI units (µm⁻¹), the projection data should be multiplied by 1000/14.6]. A sinogram of the projection is shown in Fig. 22(a). The SART-TV reconstruction of the SOFC anode with the 180 projections is shown in Fig. 22(b). 139 projections with angular range from −69 to 69° were selected to test the performances of rmwDART and mDART. The resulting reconstructions of rmwDART and mDART are shown in Figs. 22(c) and 22(d). The iteration number for the reconstructions was 140, and the parameters for rmwDART were the same as above. It can be seen that the specific detail pointed by the red arrow in Fig. 22(c) is closer to the full-angle reconstruction in Fig. 22(b) compared with the result in Fig. 22(d).

Figure 22
(a) Full-angle sinogram of the projections of the SOFC anode and the selected angular range. (b) Full-angle reconstruction by SART-TV. (c) Resulting reconstruction by rmwDART with the selected angular range from −69 to 69°. (d) Resulting reconstruction by mDART with the selected angular range from −69 to 69°.

5. Conclusion

In this paper, the rmwDART method that aimed to restore `missing wedge' artifacts has been presented. By using the oversegmentation strategy, prior knowledge of the sample for segmentation is unnecessary. Through the oversegmentation method, boundary extraction and mathematical morphological operations, positioning at the `missing wedge' artifact areas is realized. By updating the located areas, the reconstruction can be impressively improved.

The simulation results on five test images showed that rmwDART can efficiently help restore the `missing wedge' artifacts and provide high-fidelity multi-gray-level reconstructions. The important parameters of rmwDART were discussed. Simulations on noisy data showed the noise resistance capacity of rmwDART. Experimental results on an SOFC anode also proved the high-quality reconstruction of rmwDART.

Funding information

The following funding is acknowledged: National Key Research and Development Program of China (grant Nos 2016YFA0400902, 2017YFA0402904); National Natural Science Foundation of China (NSFC) (grant Nos 11475175, 11405175, 11275204, 11775224).

References

Aarle, W. van, Batenburg, K. J. & Sijbers, J. (2012). IEEE Trans. Image Process. 21, 4608–4621. Web of Science PubMed Google Scholar
Batenburg, K. J., Bals, S., Sijbers, J., Kübel, C., Midgley, P. A., Hernandez, J. C., Kaiser, U., Encina, E. R., Coronado, E. A. & Van Tendeloo, G. (2009). Ultramicroscopy, 109, 730–740. Web of Science CrossRef PubMed Google Scholar
Batenburg, K. J. & Sijbers, J. (2011). IEEE Trans. Image Process. 20, 2542–2553. Web of Science CrossRef PubMed Google Scholar
Batenburg, K. J., Sijbers, J. & IEEE (2007). 2007 IEEE International Conference on Image Processing (ICIP 2007), 16–19 September 2007, San Antonio, TX, USA, Vols. 1–7, pp. 1829–1832. New York: IEEE. Google Scholar
Batenburg, 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
Bleichrodt, F., Tabak, F. & Batenburg, K. J. (2014). Comput. Vis. Image Underst. 129, 63–74. Web of Science CrossRef Google Scholar
Dent, K. C., Hagen, C. & Grünewald, K. (2014). Methods in Cell Biology, Vol. 124, edited by T. Müller-Reichert and P. Verkade, pp. 179–216. New York: Academic Press. Google Scholar
Guan, H. Q. & Gordon, R. (1994). Phys. Med. Biol. 39, 2005–2022. CrossRef PubMed Web of Science Google Scholar
Liang, Z. T., Guan, Y., Liu, G., Bian, R., Zhang, X. B., Xiong, Y. & Tian, Y. C. (2013). Proc. SPIE, 8851, 885113. Google Scholar
Liang, Z., Guan, Y., Liu, G., Chen, X., Li, F., Guo, P. & Tian, Y. (2016). J. Synchrotron Rad. 23, 606–616. Web of Science CrossRef IUCr Journals Google Scholar
Natterer, F. (1986). The Mathematics of Computerized Tomography. Philadephia: Siam. Google Scholar
Nemeth, J. (2015). J. Math. Imaging Vis. 53, 314–331. Web of Science CrossRef Google Scholar
Paige, C. C. & Saunders, M. A. (1982a). ACM Trans. Math. Softw. 8, 43–71. CrossRef Web of Science Google Scholar
Paige, C. C. & Saunders, M. A. (1982b). ACM Trans. Math. Softw. 8, 195–209. CrossRef Web of Science Google Scholar
Panin, V. Y., Zeng, G. L. & Gullberg, G. T. (1999). IEEE Trans. Nucl. Sci. 46, 2202–2210. Web of Science CrossRef Google Scholar
Yuan, Q., Zhang, K., Hong, Y., Huang, W., Gao, K., Wang, Z., Zhu, P., Gelb, J., Tkachuk, A., Hornberger, B., Feser, M., Yun, W. & Wu, Z. (2012). J. Synchrotron Rad. 19, 1021–1028. Web of Science CrossRef CAS IUCr Journals Google Scholar
Zhuge, X. D., Jinnai, H., Dunin-Borkowski, R. E., Migunov, V., Bals, S., Cool, P., Bons, A. J. & Batenburg, K. J. (2017). Ultramicroscopy, 175, 87–96. Web of Science CrossRef PubMed Google Scholar
Zhuge, X., Palenstijn, W. J. & Batenburg, K. J. (2016). IEEE Trans. Image Process. 25, 455–468. Web of Science CrossRef PubMed 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.