research papers
A modified discrete tomography for improving the reconstruction of unknown multigraylevel material in the `missing wedge' situation
^{a}National 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 email: yongg@ustc.edu.cn
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 threedimensional 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, highquality reconstructions can be obtained from the simulations, and the reconstructed images based on the rmwDART algorithm can be obtained from soft Xray nanocomputed tomography experiments. The results showed that there is the potential for discrete tomography.
Keywords: discrete algebraic reconstruction technique (DART); missing wedge; oversegmentation; tomography; limitedangle.
1. Introduction
Combined with computed tomography (CT), highresolution microscopy tools such as the electron microscope (EM) and transmission Xray microscope (TXM) have been developed to obtain threedimensional 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 Xray microscope, mesh grids are used as the specimen carriers. It is difficult to obtain hightiltangle projection images because the samples are blocked by the mesh grids when imaging at a high tilt angle (Zhuge et al., 2017). In addition, fullrotation tilt series projections cannot be acquired because of radiation damage due to the longer exposure times for fullrotational tilt series projections (Dent et al., 2014). Thus, it is impossible to collect hightilt 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 anglelimited 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 highquality reconstructions. However, two types of prior knowledge are still required in order to obtain highquality 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 highquality 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 multigraylevel 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 variationbased simultaneous algebraic reconstruction technique (SARTTV), 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.
2.1. SARTTV 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 SARTTV method was used to conduct the initial reconstruction. Let x = (x_{j}) ∈ R^{n} denote a twodimensional image, where n is the number of the pixels. The line integral of the p (Natterer, 1986) from the measured projection data can be represented by
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 . The SART method can be expressed as follows,
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, TVmbased methods can be expressed as follows,
Here, TV is the l_{1}norm of the gradient image, and it can be expressed as
Thus, the SARTTV method can be expressed as
where TV(…) represents the TVm norm, and γ is the weight for controlling the tradeoff 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
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 anglelimited 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,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 peakpicking from the smoothed histogram h_{s}(x). w_{1} is a very important parameter in rmwDART because 2w_{1} represents the peakpicking resolution in the peakpicking process. c_{thr} is a nonnegative 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_{1} directly determines the degree of oversegmentation. By choosing a very small peakdetection parameter 2w_{1}, it can make the process of peakpicking 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 , 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 fourconnected 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).
of the region) of the image are calculated, and the image regions are redistributed. A flowchart of this step is shown in Fig. 2The 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 subiterations are used in this paper. Since the regional distribution and the regional values will probably change and become stable in the subiterations, in order to avoid misconnections between adjacent regions due to the relatively unstable regional values in the initial subiteration, T_{c,n} is set to be relatively small in the initial subiteration, and then it gradually increases in the subsequent subiterations. A discussion about this is detailed in §3.2.2. In this paper, the connectivity threshold is set to be the following,
and the fourconnected 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 subiterations 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.
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 nonartifact 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 onetime iteration for the operator {1}).
The SART method is then applied to update the values of the located artifact areas while the values of the nonartifact 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 nonartifact 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 multigray 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).
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 SARTTV reconstructions of these five images are shown in Fig. 6. The total number of iterations for the SARTTV reconstruction was 500 times.
Fig. 5(a) is used to explain rmwDART's procedures in one loop. First, the SARTTV 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_{1} [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_{1} in equation (8), and only the platforms with widths no smaller than w_{1} were kept and the width of the platforms decreased in y_{1}(x). The resulting smaller platforms of y_{1}(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_{1}(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_{1}(x) coincide at the same place. If the distance between two isolated peaks is smaller than the peakpicking resolution, the platform of the peak with the smaller count in y_{1}(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 peakpicking resolution, these two isolated peaks in h_{s}(x) will coincide with their own platforms in y_{1}(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 peakpicking resolution 2w_{1} 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.
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 higherquality 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.
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 segmentationinput 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).
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),
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.
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 i54570 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.

3.2. Parameters discussion
In rmwDART, w_{1} and T_{c,n} are the most important parameters since they directly determine the quality of the reconstruction. Therefore, the effect of different w_{1} 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 SARTTV reconstruction in Fig. 6(b) was used as the initialization.
3.2.1. Experiments on different w_{1}
2w_{1} represents the peakpicking resolution in the peakpicking process. Therefore, w_{1} directly determines the degree of oversegmentation. The smaller that 2w_{1} is set, the higher the degree of oversegmentation. In the following rmwDART reconstructions for phantom 2, 2w_{1} 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_{1} 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_{1} of 15% is relatively large, misconnections appear in the resulting segmentation.

In this study, 90 iterations were used for these rmwDART reconstructions. The convergence rates for the rwmDART reconstruction with different 2w_{1} are shown in Fig. 14(a), and the time consumption of each iteration for the rwmDART reconstruction with different 2w_{1} is shown in Fig. 14(b). The pixel errors K of the resulting reconstructions and the total time consumption of the 90iteration reconstruction are shown in Table 3. As shown in Fig. 14(a), when 2w_{1} is small enough, the degree of oversegmentation is high enough to segment the image, and then a highquality reconstruction can be obtained. For 2w_{1} = 0.5%, 2% and 5%, the results converged in ten iterations and obtained very high quality reconstructions. When 2w_{1} is not small enough, misconnections appear in the segmentation, and the reconstruction error can be very large. For 2w_{1} = 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_{1} = 0.5%, compared with 2w_{1} = 2% and 5%, the reconstruction needs a few more iterations to converge. This outcome occurs because, when 2w_{1} 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_{1} is recommended, which will ensure the high degree of oversegmentation, especially when the computer's processing capabilities are high enough.

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_{1} is the same as that in §3.1.
First, the number of subiterations for the redistribution of regions is discussed. The subiteration number is set to 8 and T_{c,n} is set to T_{c,n} = [2, 2, 2, 2, 2, 2, 2, 2] × 10^{−3}. After the oversegmentation step, the numbers of regions in each subiteration for the redistribution of regions are shown in Fig. 15. As shown in Fig. 15, after three subiterations, the change in the number of regions becomes very small, and the number of regions becomes stable. Thus, setting the number of subiterations for the redistribution of regions to five is rational.
Second, the same T_{c,n} in the subiterations is discussed. T_{c,n} is set to be the same value in the five subiterations, 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^{−3} and T_{c,n} = [8, 8, 8, 8, 8] × 10^{−3} 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 subiteration 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.

Thirdly, different T_{c,n} in different subiterations are discussed. As the misconnections may appear in the initial subiteration for the redistribution of regions, T_{c,n} is set to be relatively small in the initial subiteration, and it gradually increases in the subsequent subiterations. 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 subiteration 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^{−3} and T_{c,n} = [8, 8, 8, 8, 8] × 10^{−3}, which failed to converge and resulted in large errors, the resulting reconstructions with T_{c,n} = [1, 2, 3, 4, 5] × 10^{−3} and T_{c,n} = [2, 3, 4, 6, 8] × 10^{−3} successfully converged and had very high quality results. This shows that setting a relatively small T_{c,n} in the initial subiteration can improve rmwDART reconstructions. For the reconstructions T_{c,n} = [4, 8, 12, 16, 20] × 10^{−4} and T_{c,n} = [2, 3, 4, 6, 8] × 10^{−5}, 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.

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 signaltonoise ratio (SNR) were selected for the Poisson noise. The angular projection ranges of phantoms 1 and 3, the projection length, 2w_{1} 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 singlepixel regions, which will significantly increase the time consumption.
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 Kedge) Xray illumination. The pixel resolution of the image is 14.6 nm [for displaying the value of the resulting reconstruction in SI units (µm^{−1}), the projection data should be multiplied by 1000/14.6]. A sinogram of the projection is shown in Fig. 22(a). The SARTTV 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 fullangle reconstruction in Fig. 22(b) compared with the result in Fig. 22(d).
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 highfidelity multigraylevel 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 highquality 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üllerReichert 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., DuninBorkowski, 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.