Towards high-resolution synchrotron radiation imaging with statistical iterative reconstruction

Rashed, E.A.; Kudo, H.

doi:10.1107/S0909049512041301

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 20| Part 1| January 2013| Pages 116-124

https://doi.org/10.1107/S0909049512041301

Towards high-resolution synchrotron radiation imaging with statistical iterative reconstruction

Essam A. Rashed ^a,^b ^* and Hiroyuki Kudo ^a

^aDivision of Information Engineering, Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba 305-8573, Japan, and ^bDepartment of Mathematics, Faculty of Science, Suez Canal University, Ismailia 41522, Egypt
^*Correspondence e-mail: essam@imagelab.cs.tsukuba.ac.jp

(Received 2 March 2012; accepted 2 October 2012; online 10 November 2012)

Synchrotron radiation (SR) X-ray micro-computed tomography (CT) is an effective imaging modality for high-resolution investigation of small objects, with several applications in medicine, biology and industry. However, the limited size of the detector field of view (FOV) restricts the sample dimensions to only a few millimeters. When the sample size is larger than the FOV, images reconstructed using conventional methods suffer from DC-shift and low-frequency artifacts. This classical problem is known as the local tomography or the interior problem. In this paper, a statistical iterative reconstruction method is introduced to eliminate image artifacts resulting from the local tomography. The proposed method, which can be used in several SR imaging applications, enables high-resolution SR imaging with superior image quality compared with conventional methods. Real data obtained from different SR micro-CT applications are used to evaluate the proposed method. Results indicate a noteworthy quality improvement in the image reconstructed from the local tomography measurements.

Keywords: synchrotron radiation micro-computed tomography; local tomography; statistical iterative reconstruction; high resolution.

1. Introduction

X-ray computed tomography (CT) is an effective cross-sectional imaging modality with a wide range of applications in medicine, biology and industry. The synchrotron radiation (SR) micro-CT system provides high-resolution imaging for small objects with a spatial resolution of less than 1 µm. Nevertheless, a major disadvantage of SR micro-CT is the problem of a limited field of view (FOV) owing to current hardware limitations for the detector size (Li et al., 2007 ). This limits the imaged samples to those of size small enough to be located completely inside the FOV. When the sample size is extended outside the small FOV, the reconstructed image suffers from DC-shift and low-frequency artifacts. This is a classical problem in tomographic imaging known as the interior problem or local tomography (Natterer, 1986 ). Moreover, owing to the side effects generated from the ionizing radiation of X-ray CT, it is important to limit the dose as much as possible in several medical and biological imaging applications. A possible effective way to minimize the radiation dose is to reduce the region to be imaged to only the region of interest (ROI), which represents a local tomography problem when the ROI is located completely inside the imaged object. In the last few decades several efforts have been dedicated to obtaining exact and stable reconstruction for the local tomography of several imaging applications (Lewitt, 1979 ; Ogawa et al., 1984 ; Louis & Rieder, 1989 ; Manglos, 1992 ; Maass, 1992 ; Olson & DeStafano, 1994 ; Ohnesorge et al., 2000 ; Zhang & Zeng, 2007 ; Yu & Wang, 2009 ; Zeng & Gullberg, 2009 ; Wang et al., 2009 ; Yang et al., 2010 ; Ritschl et al., 2011 ; Lang et al., 2012 ).

The conventional image reconstruction methods, derived from an analytical inversion formula, such as filtered backprojection (FBP), are known for their non-local properties. In other words the reconstruction of the attenuation coefficient corresponding to a single image pixel is affected by all projection rays, including rays that do not intersect with the concerned pixel. This is why it was believed for a long time that the solution of the local tomography problem was not unique. However, the theories based on the concept of differentiated backprojection using the Hilbert transform (Noo et al., 2004 ) and backprojection-filtration (Pan et al., 2005 ) relax the data requirements for exact and stable ROI reconstruction to feasible ones. It became clear that the position of the ROI is an important factor in determining the possibility of accurate ROI reconstruction (Defrise et al., 2006 ). These results were extended later to solve the local tomography problem in X-ray CT imaging (Ye et al., 2007 ; Kudo et al., 2008 ). In our earlier work we proposed an iterative reconstruction algorithm for ROI reconstruction when the projection data are partially truncated (Rashed et al., 2009 ). In this paper we extend our previous work to fit with the interior problem in SR micro-CT. In the literature it has been demonstrated that the interior problem can be exactly and stably solved under the assumption that the object support and the internal region having known intensity value are known, as detailed later in §2. When either the object support or the region of known intensity is unknown, they should be estimated during the reconstruction. This estimation problem introduces some additional instability to the reconstruction problem. This appears in the form that accurate estimation of the object support and the region of known intensity can be achieved when the problem parameters are correctly selected. We have found an effective strategy for parameter selection which ensures that the object support and regions with air intensities are detected. Experimental results using real data obtained from high-resolution SR imaging indicate a significant reduction in the image artifacts when the proposed method is used.

This paper is organized as follows. In §2 the local tomography in SR micro-CT is introduced and the recent findings in the uniqueness results are discussed. The iterative reconstruction method is detailed in §3. Experimental studies and results are presented in §4 and the paper is concluded in §5.

2. Local tomography in SR micro-CT

The SR micro-CT imaging system is operated by radiating the target object with a specific X-ray power and measuring the unattenuated counts using a charge-coupled-device (CCD) camera located on the opposite side of the beam generator. The beam intensity loss represents the density of the scanned object. The projection data measurements are acquired along many angles by rotating the sample object. The acquired projection data are then used to reconstruct the attenuation coefficients that represent a cross-sectional image or volume of the scanned object. In many cases, when the sample size extends the tiny FOV of the CCD camera, the projection data are truncated and important information corresponding to the regions outside the FOV is lost, as shown in Fig. 1. To demonstrate the effect of the interior problem on the acquired projection data, two projection views acquired from an aluminium alloy sample are shown in Fig. 2. The same sample is measured in two different settings named global and local scans. A global scan is acquired such that the entire object is located completely inside the FOV, and a local scan is acquired by adding an additional surrounding shield of the same material to extend the sample size outside the FOV to almost double the size of the original object. Line integrals, computed from the measured raw data, indicate a significant change in the measured values, as shown in the profile plots in Fig. 2. Owing to the loss of complete data measurements corresponding to the regions located outside the FOV, the image reconstruction in local tomography is a challenging problem. In the following we discuss the recent advances in the uniqueness results of the local tomography problem. Most of the discussions in this section are introduced primarily for general-purpose X-ray CT imaging and associated problems in clinical medical investigations. The essential framework is also applicable for SR micro-CT imaging.

Figure 1
(a) Simple illustration of the SR micro-CT imaging system with cross-sectional views corresponding to (b) local tomography and (c) global tomography.

Figure 2
Single projection view representing line integrals of SR micro-CT data measured from a microsample of aluminium alloy (Al-5 wt% Cu) corresponding to local tomography (left) and global tomography (right). Images are shown with the same gray-level window. The corresponding central horizontal profile plots are shown below. Details of the data-acquisition parameters are discussed in §4

The recent theories derived from the Hilbert transform (Defrise et al., 2006; Ye et al., 2007; Kudo et al., 2008) and associated image reconstruction algorithms provide an interesting framework for exact ROI reconstruction from locally truncated data under a few constraints. Independent of the reconstruction algorithm, the exact and stable reconstruction from local data is achieved under the following constraints: (i) all projection rays passing through the ROI (Ω) are measured; (ii) for iterative reconstruction the image matrix should be large enough to include the entire extended object, not only the ROI (Ω); (iii) a compact support (μ) of the object, named the object support (OS), is available; and (iv) a small sub-region (ω) inside the internal ROI (Ω) is a priori known. A simple illustration that concludes the recent developments in exact ROI reconstruction results is shown in Fig. 3.

Figure 3
Different ROI imaging configurations. (a) Global scan where the entire object is located inside the FOV (standard configuration and the image reconstruction is exact). (b) The FOV contains two regions (shaded) located outside a roughly estimated OS (μ) where exact reconstruction is assured for the ROI (Ω) that is defined as a union of parallel vertical Hilbert lines that connect shaded regions without intersection with the remaining regions of the object (Noo et al., 2004

). (c) The FOV contains a single region (shaded) outside the OS (μ) and solution exactness for ROI (Ω) is assured (Defrise et al., 2006

). (d) The FOV is located completely inside the OS (μ) where the exact OS is estimated and reconstruction of ROI (Ω) is exact (Rashed et al., 2009

). (e) Local tomography is exact if a sub-region (ω) inside the ROI (Ω) is known and a rough OS (μ) is defined (Kudo et al., 2008

). (f) The local tomography problem to be investigated in this paper with ROI (Ω) is located completely inside the object and both OS (μ) and sub-region (ω) are estimated for exact ROI reconstruction.

The first constraint reduces the required data measurements for exact and stable ROI reconstruction. It was believed that for exact and stable ROI reconstruction all projection rays passing through the whole object should be measured and used for the reconstruction. This rule is directly indicated from the inversion formula of analytical reconstruction, where a single image pixel is computed not only from projection rays passing through this pixel but also from all rays passing through the whole object. However, this strong constraint is relaxed such that it is sufficient to include the projection rays passing through the ROI (Ω) only (Noo et al., 2004; Defrise et al., 2006). As for the second constraint, owing to the layout of iterative reconstruction schemes it is important to use an extended image grid that is large enough to contain the entire support of the object. This formulation is essential to enable the enforcement of the OS constraint presented in the third constraint. The first and second constraints can be easily satisfied in the practical implementation in the local tomography for SR micro-CT imaging.

The most challenging part concerns the third and fourth constraints. On one hand the third constraint is not always available, as the definition of the exact object contour requires additional effort to obtain the object measurements (e.g. using additional hardware devices). When the OS (μ) is known, the third constraint is simply enforced by setting the image pixels located outside the OS to the attenuation value of air, which is almost zero. The accuracy of the OS knowledge has a large influence on the exactness of the ROI reconstruction. The closer the OS is to the exact one, the higher the quality of the reconstructed ROI (Rashed et al., 2009). On the other hand it is not always possible to satisfy the fourth constraint as it requires both intensity and spatial information of a small sub-region (ω) inside the ROI (Ω).

The iterative algorithm for ROI reconstruction, detailed by Rashed et al. (2009), provides an automatic and efficient way to detect the exact OS during the iterative reconstruction. This method was used in ROI reconstruction when the projection data were partially truncated (i.e. the ROI includes some region outside the exact OS), as shown in Fig. 3(d). In this paper, not only the exact OS but also the interior sub-region is estimated during the iterative reconstruction. A good estimation of the OS and interior sub-region in the early iterations will provide the required constraints for theoretically exact and stable local tomography.

3. Materials and methods

In this paper we consider the interior problem where both the OS (μ) and the prior information of the internal sub-region (ω) are unavailable. The proposed method can estimate the OS with high accuracy during the iterative reconstruction. By considering imaging applications where the image includes a sub-region of air-attenuated value, the proposed method can also satisfy the fourth constraint defined above for exact reconstruction by estimating a sub-region inside the FOV with air-attenuation value. If the region with the air-attenuation value inside the FOV is not available, it is sufficient to consider the attenuation value of some other known materials of the scanned object (Kudo et al., 2008). The difference is only a simple modification in the implementation of the reconstruction algorithm. In this section we review the main concepts of statistical iterative reconstruction required to introduce the proposed algorithm.

3.1. Statistical iterative reconstruction

The SR micro-CT imaging system can be defined with the following transmission CT model,

$[y_{i}\simeq{\rm{Poisson}}\left[b_{i}\exp\left(-\textstyle\sum\limits_{j\,=\,1}^na_{ij}x_j\right)\right],\quad i=1,\ldots,m,\eqno(1)]$

where $[{\bf{x}}]$ = $[(x_{1},\ldots,x_{n})]$ is the image vector representing the attenuation coefficients of the imaged object, $[{\bf{y}}]$ = $[(\,y_{1},\ldots,y_{m})]$ is the measured photon counts with the CCD camera, $[{\bf{b}}]$ = $[(b_{1},\ldots,b_{m})]$ is the blank scan, which is acquired by data measurements without the object, and A = { a_ij} is the m×n system matrix that models the imaging system. In these statistical models we ignore the effect of scattered photons and other background events for simplicity; however, the model generalization is direct and easy.

The statistical model in (1) represents an inverse problem which can be solved using several approaches. The log-likelihood function for the observed photon counts (Lange & Carson, 1984 ) is defined as

$[l({\bf{x}}) = \textstyle\sum\limits_{{i\,=\,1}}^{m}\left[y_{i}\log(b_{i})-y_{i}\langle{\bf{a}}_{i},{\bf{x}}\rangle-\log(y_{i}!)-b_{i}\exp(-\langle{\bf{a}}_{i},{\bf{x}}\rangle)\right],\eqno(2)]$

where $[\langle{\bf{a}}_i,{\bf{x}}\rangle]$ = $[\textstyle\sum_{{j\,=\,1}}^{n}a_{{ij}}x_{j}]$ is the inner product of the ith row of matrix A and image vector $[{\bf{x}}]$ . The solution of the inverse problem is found through minimizing the negative log-likelihood function after ignoring the irrelevant terms and is given by

$[\hat{x}=\arg\min_{{{{\bf{x}}}\,\geq\,0}}L({\bf{x}}),\eqno(3)]$

$[L({\bf{x}})=\textstyle\sum\limits_{{i\,=\,1}}^{m}\left[b_{i}\exp(-\langle{\bf{a}}_{i},{\bf{x}}\rangle)+y_{i}\langle{\bf{a}}_{i},{\bf{x}}\rangle\right].\eqno(4)]$

When the projection data suffer from data inconsistency, such as the case of data limitation, the well known Bayesian methods, that include a priori information of the scanned object, can be used. The penalty term is mainly used to limit the number of feasible solutions to those that satisfy the penalty term. The solution is found by maximizing a posteriori (MAP) function

$[P({\bf{x}}|{\bf{y}})={{P({\bf{y}}|{\bf{x}})P({\bf{x}})} \big/ {P({\bf{y}})}},\eqno(5)]$

and the maximum a priori solution of the image reconstruction problem is found by

$[\hat{x}=\arg\min_{{{\bf{x}}\,\geq\,0}}f({\bf{x}}),\qquad f({\bf{x}})=L({\bf{x}})+\beta U({\bf{x}}),\eqno(6)]$

where $[U({\bf{x}})]$ is known as the penalty term and β is a parameter that controls the balance between the data fidelity term enforced by the log-likelihood function and the penalty term.

3.2. Majorization–minimization approach

The minimization of cost function in equation (6) is an essential task to be considered for the derivation of the image reconstruction algorithm. As shown in the next section, we are interested in using a penalty term $[U({\bf{x}})]$ based on the l₀ norm function to enforce the image sparsity. It is difficult to use ordinary gradient methods to minimize the cost function $[f({\bf{x}})]$ as the penalty term is neither convex nor differentiable. An alternative useful approach is the use of the majorization–minimization (MM) strategy (De Pierro, 1995 ; Daubechies et al., 2004 ). The MM approach is an iterative approach that operates as follows. At each iteration k the non-separable part of the cost function $[f({\bf{x}})]$ is approximated by a separable quadratic function $[\tilde{f}({\bf{x}}\semi{\bf{x}}^{k})]$ around $[{\bf{x}}]$ = $[{\bf{x}}^{k}]$ . This operation is called majorization. The resulting quadratic function is easier to be minimized. We may refer the reader to Hunter & Lange (2004 ) and Lange (2004 , ch. 6) for a more detailed explanation of the MM approach.

3.3. Reconstruction algorithm

To develop the cost function for the image reconstruction problem under study we consider the following formulation of the penalty term,

$[U({\bf{x}})=\textstyle\sum\limits_{{j\,=\,1}}^{{n}}s(x_{j}),\qquad s(t)=\lim\limits_{{\varepsilon\rightarrow+0}}|t|^{{\varepsilon}}\equiv\Big\{\matrix{ 1,\hfill & t\neq0,\hfill \cr 0,\hfill & t=0.\hfill}\eqno(7)]$

This design of the cost function for image reconstruction simply indicates that the reconstructed image is selected to be the sparsest image that fits with the measured projection data. This formulation is proven to be useful in automatic estimation of the OS during image reconstruction. The practical implementation yields a thresholding function that trims the image pixels with intensity values close to zero to the value of zero (Rashed et al., 2009). Moreover, it is expected that this approach would satisfy the fourth constraint for exact reconstruction of local tomography, stated in the previous section, when the local ROI includes a small region with zero attenuation coefficient such as air holes.

The main obstacle in minimizing the cost function in (6) is that the penalty term is based on the l₀ norm. The MM approach is used to handle this problem. At each iteration k the non-separable part of the cost function $[f({\bf{x}})]$ is approximated by a separable function $[\tilde{f}({\bf{x}}\semi{\bf{x}}^{k})]$ around $[{\bf{x}}]$ = $[{\bf{x}}^{k}]$ . Then $[\tilde{f}({\bf{x}}\semi{\bf{x}}^{k})]$ is analytically minimized to compute the image estimate for the next iterate k + 1. The iterative reconstruction algorithm can be summarized as follows:

(i) Initialization. Set initial image estimate x⁰ = ∊, with ∊ > 0 and initialize iteration number k = 0.

(ii) Majorization. The cost function $[f({\bf{x}})]$ is majorized to obtain the following separable cost function:

$[\tilde{f}({\bf{x}}\semi{\bf{x}}^{k})= \textstyle\sum\limits_{{j\,=\,1}}^{{n}}\beta\left[t_{j}(x_{j}-p_{j})^{2}+s(x_{j})\right]+C({\bf{x}}^{{k}}),\eqno(8)]$

$[t_{j} = {{1}\over{2\beta x_{j}^{k}}} \textstyle\sum\limits_{{i\,=\,1}}^{m}a_{{ij}}\langle{\bf{a}}_{i},{\bf{x}}^{k}\rangle b_{i}\exp\left(-\langle{\bf{a}}_{i},{\bf{x}}^{k}\rangle\right),\eqno(9)]$

$[p_{j} = x_{j}^{k}+x_{j}^{k}{{\textstyle\sum_{{i\,=\,1}}^{m}a_{{ij}}b_{i}\exp(-\langle{\bf{a}}_{i},{\bf{x}}^{k}\rangle)-y_{i}} \over {\textstyle\sum_{{i\,=\,1}}^{m}a_{{ij}}\langle{\bf{a}}_{i},{\bf{x}}^{k}\rangle b_{i}\exp(-\langle{\bf{a}}_{i},{\bf{x}}^{k}\rangle)}},\eqno(10)]$

where $[C({\bf{x}}^{k})]$ is a term independent of $[{\bf{x}}]$ . Equations (8)–(10) were derived by approximating the likelihood part with a quadratic function and keeping the penalty function part as it is. See Rashed & Kudo (2012 ) for the mathematical derivation.

(iii) Minimization. The separable cost function $[\tilde{f}({\bf{x}}\semi{\bf{x}}^{k})]$ in (8) is minimized over $[{\bf{x}}\geq]$ 0 to obtain the image estimate for the next iterate,

$[{\bf{x}}^{{k+1}}= \arg\min_{{{\bf{x}}\,\geq \,0}}\tilde{f}({\bf{x}}\semi{\bf{x}}^{k}).\eqno(11)]$

(iv) Enforce prior knowledge. When prior knowledge of the sub-region (ω) is available, assign image pixels located inside ω with the a priori known intensity value (μ),

$[x_{j}= \mu\quad \forall\quad j\in\omega.\eqno(12)]$

(v) Iteration condition. Increase k and go to step (ii) until reaching a stopping criterion.

An interesting observation in the practical implementation of the proposed algorithm is that the minimization in (11) can be implemented directly through following a simple hard-thresholding function,

$[x_{{j}}^{{k+1}}= \Bigg\{\matrix{ 0,\hfill & {\rm{for}}\quad p_{j}\leq(1/t_{j})^{1/2}\hfill \cr p_{j},\hfill & {\rm{otherwise.}}\hfill}\eqno(13)]$

Therefore, the reconstruction algorithm is performed as a combination of two main steps in an alternating manner. The first step is a single iteration of the conventional convex algorithm (Lange & Fessler, 1995 ) using (10). The second step is a pixel-by-pixel thresholding function using (13). To speed up the reconstruction, we employ the ordered subsets approach (Beekman & Kamphuis, 2001 ) in the experimental studies presented within this paper. We also note that in the experimental studies we have used a dynamic value of parameter β, which is proved to be a good approach to selecting β in practice. First, we initiate β with a relatively large value to be sure that the object support and internal regions with air are completely detected. The value of β is then gradually decreased to enforce the likelihood function for the recovery of pixels that incorrectly thresholded in the early iterations and retrieve missing image details. Details of this approach can be found by Rashed & Kudo (2012). The proposed algorithm has several advantages. First, the implementation is simple as it coincides with the conventional convex algorithm with additional thresholding operation applied after each iteration. Second, the computation cost is almost the same as the conventional iterative method since the cost of the thresholding operation can be ignored relative to the computation of a single iteration of the convex algorithm. Third, it is possible to speed up the reconstruction using image updates with the ordered subsets approach.

4. Experimental studies

In this section the developed algorithm is evaluated using real data obtained from different SR micro-CT applications. In all the experimental studies introduced in this paper there exists a sub-region with air-attenuated values so that the reconstruction algorithm becomes more simpler and step (iv) is not required.

4.1. Image quality measurements

For quantitative evaluation of image quality improvement a set of measurements are calculated. We use the FBP image obtained from global tomography as the golden truth. The relative root-mean-square error (RRME) is calculated as

$[{\rm{RRME}}({\bf{x}}) = \left[\textstyle\sum\limits_{{j\,=\,1}}^{{n}}(x_{j}-x_{j}^{{*}})^{2}/\textstyle\sum\limits_{{j\,=\,1}}^{{n}}(x_{j}^{*})^{2}\right]^{1/2},\eqno(14)]$

where $[{\bf{x}}^{*}]$ represents the FBP image obtained from global tomography, when available. The image contrast is also measured using the following formula,

$[{\rm{Contrast}}({\bf{x}}) = {{|\bar{x}_{\rm{s}}-\bar{x}_{\rm{b}}|} \over {\bar{x}_{\rm{s}}+\bar{x}_{\rm{b}}}},\eqno(15)]$

where $[\bar{x}_{\rm s}]$ and $[\bar{x}_{\rm b}]$ are the mean pixel values of selected resolution pixels and background pixels, respectively, and computed using the following equation,

$[\bar{x}_{\rm s}= {{1} \over {n_{\rm{s}}}}\,\,\textstyle\sum\limits_{{j\,=\,1}}^{{n_{\rm{s}}}}x_{j},\eqno(16)]$

where n_s is the number of pixels within the s-ROI. To evaluate the image contrast and noise properties, we also compute the contrast-to-noise ration (CNR) as follows,

$[{\rm{CNR}}({\bf{x}})= {{2|\bar{x}_{\rm{s}}-\bar{x}_{\rm{b}}|} \over {\delta_{\rm{s}}\sigma_{\rm{s}}+\delta_{\rm{b}}\sigma_{\rm{b}}}},\quad\delta_{\rm{s}} = {{n_{\rm{s}}} \over {n_{\rm{s}}+n_{\rm{b}}}},\quad\delta _{\rm{b}} = {{n_{\rm{b}}} \over {n_{\rm{s}}+n_{\rm{b}}}},\eqno(17)]$

where n_b is the number of pixels in b-ROI, and $[\sigma_{\rm{s}}]$ and $[\sigma_{\rm{b}}]$ are the standard deviation over s-ROI and b-ROI. The value of $[\sigma_{\rm{s}}]$ is computed as follows,

$[\sigma_{\rm{s}}= \left[ {{1} \over {n_{\rm{s}}-1}}\,\,\textstyle\sum\limits_{{j\,=\,1}}^{{n_{\rm{s}}}}(x_{j}-\bar{x}_{\rm{s}})^{2}\right]^{1/2}.\eqno(18)]$

The CNR value evaluates the ratio of the contrast of a target structure in the image and the standard deviation of statistical noise.

4.2. Aluminium alloy sample

The first experiment was carried out at beamline BL20XU of SPring-8 synchrotron radiation source located in Hyogo, Japan. In this experiment the object is a microsample of aluminium alloy (Al-5 wt% Cu). The microsample was imaged using two different imaging settings corresponding to global and local tomography for effective and reliable evaluation of the proposed reconstruction algorithm. First, the global tomography was obtained by using a microsample of size 1 mm that is located completely inside the FOV. Later, the local tomography was set by adding an additional surrounding shield of the same material to extend the sample size to 2 mm. This set-up provides the ground truth information for the 1 mm-size sample through global scanning for effective evaluation of the image reconstructed later from local scanning. The sample was scanned with an X-ray energy of 35 keV and projection data were obtained using a CCD camera of 4000 × 2624 pixels with 2 × 2 binning mode and 0.5 µm × 0.5 µm pixel size. The original parallel-beam projection data used in this experiment were acquired over 1500 uniformly spaced view angles over 180° with an exposure time of 300 ms per view. Two blank scans were obtained through pre- and post-data measurements for each data-acquisition process. The detector array was resampled into 500 × 328 pixels to reduce the computation cost for the image reconstruction.

We have used the proposed method to reconstruct images from both global and local projection data. The number of iterations was ten with five subsets and β = 0.1. The reconstructed images corresponding to the central two-dimensional slice using conventional FBP and the proposed method from global and local projection data are shown in Figs. 4 and 5, respectively. We note that the reconstruction using the proposed method from local data is assigned to an image grid of 1000 × 1000 pixels which is trimmed to the ROI size of 500 × 500 pixels for illustration in Fig. 5(b). Image quality measurements including RRME, contrast and CNR are shown in Fig. 6. In global data measurements it is clear that the reconstructed images are of equivalent quality in both methods. However, the FBP image reconstructed from local data suffer from sever DC-shift and low-frequency artifacts. These artifacts are significantly suppressed when using the proposed algorithm.

Figure 4
Images reconstructed from global projection data using (a) FBP and (b) proposed algorithm shown within the same gray-scale window, and the corresponding central horizontal (left) and vertical (right) profile plots. The yellow square region is magnified in the top right-hand corner of (a) showing the s-ROI in red and b-ROI in blue.

Figure 5
Images reconstructed from local projection data using (a) FBP and (b) proposed algorithm shown within the same gray-scale window, and the corresponding central horizontal (left) and vertical (right) profile plots. The profiles of the image shown in Fig. 4(a)

are also included for comparison.

Figure 6
Image quality measurements for aluminium alloy data reconstructed using different configurations.

4.3. Calcite in talc powder sample

In this experiment the sample is a calcite fluid inclusion embedded in talc powder. The data measurements were performed at beamline BL20B2 at SPring-8. The X-ray beam used was of power 20 keV. Two sets of data were acquired which we refer to as global and local data. The global data were measured using a CCD camera with 1560 × 680 pixels and a pixel size of 6.66 µm. The sample was rotated to acquire 2250 projections over 180°. The FOV of 10.4 mm was sufficient to enclose the entire sample. The global data were used as a true reference to evaluate the accuracy of the image reconstruction using local data. The local data measurements were obtained using a CCD camera with higher resolution and a size of 2000 × 1312 pixels with a pixel size of 3.33 µm that provides a FOV of 6.66 mm. The local data were obtained through 3000 projections over 180°. The local measurements provide a higher resolution of the sample. However, the measured data are fully truncated and present a local tomography. Blank scans were measured in sequence with sample measurements. A periodic blank scan was acquired after every 20 projections in both global and local measurements. Every periodic blank scan was used in association with the proceeding detector measurements.

First, a direct reconstruction using the FBP algorithm was employed using both local and global projection data and the results are shown in Figs. 7(a) and 7(c), respectively. Reconstruction from local projection data provides a higher resolution image with clearer image details than the global data acquisition. However, local image reconstructed using the FBP algorithm suffers from typical local tomography artifacts, especially in the peripheral regions. The proposed algorithm was used to improve the image quality with image reconstruction in an extended image grid of 3120 × 3120 pixels. Iterative reconstruction was implemented through ten iterations with ten subsets and β = 0.05. The reconstructed image was trimmed to the accurate FOV size and is shown in Fig. 7(b). For better visual evaluation of image quality, the reconstructed images from local projection data using FBP and proposed method were down-sampled using 2 × 2 pixels binning and embedded into the reference true image reconstructed using the FBP algorithm from global data. Results are shown in Figs. 7(d) and 7(e). It is clear that the image reconstructed using the conventional FBP method contains DC-shift artifacts which are clearly shown in the mismatch between the inner and outer regions of the image shown in Fig. 7(d), while this change is not observed in the image reconstructed using the proposed method in Fig. 7(e).

Figure 7
Images reconstructed from local projection data using (a) FBP and (b) proposed algorithm, and (c) image reconstructed from global projection data using the FBP algorithm. Images (a) and (b) are down-sampled (2 × 2) and then embedded in (c), and are shown in (d) and (e), respectively.

The proposed method can significantly suppress image artifacts known for local tomography. Moreover, as the proposed reconstruction algorithm is based on an accurate statistical model, results have better noise properties. This can be easily recognized from the image in Fig. 7(e) where the internal region reconstructed using the proposed algorithm is of less noise content compared with the external region reconstructed using the FBP algorithm. To demonstrate how exact is the estimation of the OS (μ) and the pixels with air attenuation value (ω) by using the proposed method, we track pixels with zero intensity value through the iterative reconstruction. These pixels are marked in red in the images shown in Fig. 8. We observe that, by iterating, the estimation of the OS as well as the internal pixels representing air are approaching the golden truth of the FBP global image. However, we should note that this observation is true after the optimization of the parameter β.

Figure 8
Estimation of OS and internal pixels with zero attenuation. Image pixels marked red have zero attenuation value. Images (a)–(e) are local reconstructions from one, two, four, six and ten iterations of the proposed method, respectively. (f) Global FBP image with pixel values ≤ 0 marked.

The RRME, image contrast and CNR values for the images in Figs. 7(a)–7(c) are shown in Fig. 9. The resolution of the images reconstructed from local and global data does not match; therefore the images obtained from local data are down-sampled using (2 × 2) binning to be able to compute the RRME. The image contrast is computed using selected ROIs that represent resolution and background pixels as shown in Fig. 10. Central horizontal and vertical profiles of the images in Figs. 7(a)–7(c) are also shown in Fig. 10. It is clear from these measurements that the proposed method provides a higher image quality in terms of RRME as well as image contrast.

Figure 9
Image quality measurements of the ROI reconstructions shown in Figs. 7(a)–7(c)

Figure 10
ROIs used to compute the image contrast are shown on the left. The red rectangle represents the resolution pixel region (s-ROI) while the blue rectangle is assigned to the background pixel regions (b-ROI). Central horizontal (top) and vertical (bottom) profile plots corresponding to the images shown in Figs. 7(a)–7(c)

are shown on the right.

4.4. Mouse lung

This experiment is implemented using projection data acquired from mouse lung. The target of the imaging application is to investigate lung alveoli. The sample was imaged at beamline BL20B2 at SPring-8 and projection data were acquired using a detector of 2000 × 300 pixels with a pixel size of 3.33 µm. Projection data were measured through 1500 projections over 180°. A single blank scan and silent data measurements were acquired during data acquisition. Unfortunately the ground truth information (i.e. global projection data) regarding this experiment is unavailable for quantitative evaluation. A single slice was reconstructed using standard FBP and the proposed method. The proposed method was used through ten iterations and five subsets with β = 0.05. The reconstructed images and the corresponding central horizontal and vertical profile plots are shown in Fig. 11 while the image quality measurements are plotted in Fig. 12.

Figure 11
Reconstructed images using (a) FBP and (b) proposed algorithm corresponding to the central slice of the mouse lung data. In the top right-hand corner of (a) is a magnification showing the s-ROI in red and the b-ROI in blue. Central horizontal (middle) and vertical (bottom) profiles are also shown.

Figure 12
Image quality measurements for mouse lung data.

In the FBP reconstruction the image suffers from DC-shift and low-frequency artifacts, which are effectively suppressed when using the proposed method. The profile plots indicate a large variation in the FBP image and possibly incorrect estimation, especially in peripheral regions. However, the image reconstructed using the proposed method provides consistent intensity values for pixels of the same anatomy structures, even for pixels located in the peripheral regions.

5. Conclusion

In this work we propose a new statistical iterative reconstruction algorithm for local tomography. The proposed algorithm is based on a simple approach that effectively estimates the required a priori known information for theoretically accurate and stable reconstruction from fully truncated projection data. The proposed method is tested using SR micro-CT real data, and an effective improvement in image quality is achieved compared with the conventional method. However, the computation cost of the proposed method is still the major challenge. Future work includes reducing the computation cost to a reasonable level for practical implementation in high-resolution micro-CT imaging.

Acknowledgements

The authors would like to thank Dr Hiroyuki Toda, Toyohashi University of Technology, Japan, Dr Toshihiro Sera and Dr Akira Tsuchiyama, Osaka University, Japan, Dr Tsukasa Nakano, AIST, Japan, and Dr Kentaro Uesugi, Japan Synchrotron Radiation Research Institute, Japan, for providing the real data used in this study. This work was partially supported by a Grant-in-Aid from the Japan Society for the Promotion of Science (JSPS) for foreign post-doctoral fellows (No. P10052).

References

Beekman, F. J. & Kamphuis, C. (2001). Phys. Med. Biol. 46, 1835–1844. Web of Science CrossRef PubMed Google Scholar
Daubechies, I., Defrise, M. & De Mol, C. (2004). Commun. Pure Appl. Math. 57, 1413–1457. Web of Science CrossRef Google Scholar
Defrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037–1053. Web of Science CrossRef Google Scholar
De Pierro, A. R. (1995). IEEE Trans. Med. Imaging, 14, 132–137. CrossRef PubMed CAS Web of Science Google Scholar
Hunter, D. R. & Lange, K. (2004). Am. Stat. 58, 30–37. Web of Science CrossRef Google Scholar
Kudo, H., Courdurier, M., Noo, F. & Defrise, M. (2008). Phys. Med. Biol. 53, 2207–2231. Web of Science CrossRef PubMed Google Scholar
Lang, S., Dominietto, M., Cattin, P., Ulmann-Schuler, A., Weitkamp, T. & Müller, B. (2012). J. Synchrotron Rad. 19, 114–125. Web of Science CrossRef IUCr Journals Google Scholar
Lange, K. (2004). Optimization. New York: Springer-Verlag. Google Scholar
Lange, K. & Carson, R. (1984). J. Comput. Assist. Tomogr. 8, 306–316. CAS PubMed Web of Science Google Scholar
Lange, K. & Fessler, J. A. (1995). IEEE Trans. Image Process. 4, 1430–1438. CrossRef PubMed CAS Web of Science Google Scholar
Lewitt, R. M. (1979). Med. Phys. 6, 412–417. CrossRef CAS PubMed Google Scholar
Li, L., Toda, H., Ohgaki, T., Kobayashi, M. & Kobayashi, T. (2007). J. App. Phys. 102, 114908. Web of Science CrossRef Google Scholar
Louis, A. K. & Rieder, A. (1989). Num. Math. 56, 371–383. CrossRef Web of Science Google Scholar
Maass, P. (1992). SIAM J. Appl. Math. 52, 710–724. CrossRef Web of Science Google Scholar
Manglos, S. H. (1992). Phys. Med. Biol. 37, 549–562. CrossRef PubMed CAS Web of Science Google Scholar
Natterer, F. (1986). The Mathematics of Computerized Tomography. Philadelphia: Wiley. Google Scholar
Noo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903–3923. Web of Science CrossRef PubMed Google Scholar
Ogawa, K., Nakajima, M. & Yuta, S. (1984). IEEE Trans. Med. Imaging, 3, 34–40. CrossRef PubMed CAS Google Scholar
Ohnesorge, B., Flohr, T., Schwarz, K., Heiken, J. P. & Bae, K. T. (2000). Med. Phys. 27, 39–46. Web of Science CrossRef PubMed CAS Google Scholar
Olson, T. & DeStafano, J. (1994). IEEE Trans. Acoust. Speech Signal Process. 42, 2055–2067. CrossRef Google Scholar
Pan, X., Zou, Y. & Xia, D. (2005). Med. Phys. 32, 673–684. Web of Science CrossRef PubMed Google Scholar
Rashed, E. A. & Kudo, H. (2012). Phys. Med. Biol. 57, 2039–2061. Web of Science CrossRef PubMed Google Scholar
Rashed, E. A., Kudo, H. & Noo, F. (2009). Med. Imaging Tech. 27, 321–331. Google Scholar
Ritschl, L., Bergner, F., Fleischmann, C. & Kachelriess, M. (2011). Phys. Med. Biol. 56, 1545–1561. Web of Science CrossRef PubMed Google Scholar
Wang, G., Yu, H. & Ye, Y. (2009). Med. Phys. 36, 3575–3581. Web of Science CrossRef PubMed Google Scholar
Yang, J., Yu, H., Jiang, M. & Wang, G. (2010). Inverse Probl. 26, 035013. Web of Science CrossRef Google Scholar
Ye, Y., Yu, H., Wei, Y. & Wang, G. (2007). Int. J. Biomed. Imaging, 2007, 63634. PubMed Google Scholar
Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791–2805. Web of Science CrossRef PubMed Google Scholar
Zeng, G. L. & Gullberg, G. T. (2009). Phys. Med. Biol. 54, 5805–5814. Web of Science CrossRef PubMed Google Scholar
Zhang, B. & Zeng, G. L. (2007). Med. Phys. 34, 935–944. Web of Science CrossRef PubMed CAS 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 OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 20| Part 1| January 2013| Pages 116-124

https://doi.org/10.1107/S0909049512041301

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

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

Towards high-resolution synchrotron radiation imaging with statistical iterative reconstruction

1. Introduction

2. Local tomography in SR micro-CT

3. Materials and methods

3.1. Statistical iterative reconstruction

3.2. Majorization–minimization approach

3.3. Reconstruction algorithm

4. Experimental studies

4.1. Image quality measurements

4.2. Aluminium alloy sample

4.3. Calcite in talc powder sample

4.4. Mouse lung

5. Conclusion

Acknowledgements

References

research papers