Open Access
Add to BibSonomy
Abstract
Edge-illumination X-ray phase-contrast tomography (EIXPCT) is an emerging technique that enables practical phase-contrast imaging with laboratory-based X-ray sources. A joint reconstruction method was proposed for reconstructing EIXPCT images, enabling novel flexible data-acquisition designs. However, only limited efforts have been devoted to optimizing data-acquisition designs for use with the joint reconstruction method. In this study, several promising designs are introduced, such as the constant aperture position (CAP) strategy and the alternating aperture position (AAP) strategy covering different angular ranges. In computer-simulation studies, these designs are analyzed and compared. Experimental data are employed to test the designs in real-world applications. All candidate designs are also compared for their implementation complexity. The tradeoff between data-acquisition time and image quality is discussed.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
X-ray phase-contrast tomography (XPCT) seeks to estimate the complex-valued refractive index distribution $n=1-\delta +i\beta$ of a sample, where $i\equiv \sqrt {-1}$, $\delta$ denotes the sample’s refractive index decrement distribution, and $\beta$ denotes the absorption distribution. Compared with conventional computed tomography (CT) techniques, XPCT holds promise to provide better differentiation of soft tissues and generally low Z materials. In recent years, various types of XPCT methods have been developed, including grating-based X-ray phase-contrast tomography (GB-XPCT) [1–3], propagation-based X-ray phase-contrast tomography (PB-XPCT) [4–6] and edge illumination X-ray phase-contrast tomography (EIXPCT) [7–9]. XPCT can be implemented with microfocal sources, but it comes with the expense of greatly increased exposure times [4,10]. For example, a conventional X-ray source can be employed for XPCT imaging in GB-XPCT when a source grating [1] is employed to produce beamlets possessing the necessary spatial coherence. However, this greatly reduces the output of the source, resulting in prolonged exposure times. Compared with previous methods, EIXPCT has the advantage of being readily usable with laboratory-based X-ray sources without the need for high coherence or an extra grating that limits the source output, although some restrictions on maximum source size still exist mainly due to geometrical reasons, imposing a $\sim$ 100 $\mu$m limit for most commonly used setups. [9,11–13].
EIXPCT employs two carefully-positioned partially opaque masks to encode refraction and absorption information. Conventional reconstruction methods for EIXPCT employ a two-step approach [14,15]. The two steps are a phase retrieval step that separates the absorption and refractive properties, and a tomographic reconstruction step that reconstructs the distributions of both properties. The phase retrieval step requires at least two images to be taken at each tomographic view angle during the data-acquisition process. The two acquired images correspond to different mask positions. At these positions, different portions of the incident X-ray beams hit the sensitive areas of the detector pixels, thus creating different fringes in the projection images that correspond to the object’s boundaries. The differences in the two images taken at the same tomographic view angle contain refractive index information about the imaged sample, and enable the retrieval of both refractive and absorption information. In the second step, conventional tomographic image reconstruction methods are employed to estimate the $\delta$ and $\beta$ distributions, respectively.
An intrinsic limitation of this conventional reconstruction method is the need to acquire multiple images at each tomographic view angle, which slows the data-acquisition, hindering in vivo applications. To accelerate the data-acquisition process, single-shot imaging techniques have been proposed that require only one image for each tomographic view angle [16–19]. However, most existing single-shot technologies are restricted by certain assumptions – such as a parallel beam assumption [16] and a quasi-homogeneous object [17,20] – or require special hardware – such as energy sensitive detectors [18]. Recently, a joint reconstruction (JR) method was proposed for single-shot imaging [21,22]. Compared with other single-shot technologies, the JR method does not make restrictive assumptions that have been required in previous methods and does not rely on specialized hardware. Moreover, JR allows for flexible single-shot data-acquisition designs [23]. However, the advantages and disadvantages of different designs have not been comprehensively analyzed, and their impact on the stability and quality of the reconstruction process remains largely unknown. In this study, we propose a series of data-acquisition designs and assess them systematically.
The remainder of the article is organized as follows. The EIXPCT imaging system and the JR method are reviewed in Section 2. In Section 3., the simulation study and candidate mask displacement strategies are introduced. Then, in Section 4., the tomographic scanning requirements for stably implementing JR using the different mask displacement strategies are investigated in an idealized situation, and cross-talk is also discussed. Section 5. provides a deeper understanding of how the data-acquisition designs affect image quality in terms of contrast-to-noise ratio and resolution. Experimental results corresponding to different data-acquisition designs are presented in Section 6. Finally, a summary of the promising data-acquisition designs and their properties is provided in Section 7.
2. Background
2.1 Imaging model
An EIXPCT system is shown in Fig. 1 [14,24]. There are two parallel masks (also referred to as apertures): the sample mask, $A_1$, and the detector mask, $A_2$, which are placed immediately before the sample O and the detectors, respectively. The distance between the X-ray source S and $A_1$ is denoted by $l_{so}$, and that between $A_1$ and $A_2$ is denoted by $l_{od}$. For simplicity, the X-ray source is assumed to be monochromatic and point-like. The approximation is appropriate because the imaging model of EIXPCT has been shown to be tolerant of polychromatic beams and up to about $100\mu m$ focal spot size in most setups [11–13]. Define the stationary coordinate system $(x, y)$ so that the $x$-axis is perpendicular to the two masks, while the $y$-axis is parallel to them. During the data-acquisition process, the sample mask can be moved along the $y$-axis, and its displacement is denoted as $\Delta \xi$. The movement of the sample mask causes different portions of the incident X-ray beam to fall on the sensitive areas of the detector, encoding differential phase information in the measured X-ray intensities. A rotation stage is set below the sample O, allowing tomographic scanning. A rotated coordinate system $(x_r, y_r)$ is attached to the rotation stage, and the origin of both coordinate systems corresponds to the center of rotation for the tomographic scan. Denote the tomographic view angle measured from the positive $x$-axis as $\theta$. The two coordinate systems are related as $x_r=x\cos \theta -y\sin \theta$, and $y_r = y\cos \theta +x\sin \theta$. EIXPCT seeks to estimate the refractive index decrement distribution, $\delta$, and the absorption distribution, $\beta$, of the object. For a two-dimensional (2D) reconstruction problem, both are functions of the spatial location $(x, y)$ in the stationary coordinate system.
Fig. 1. Schematic of the experimental setup. $A_1$ is the sample mask and is placed before the sample. $A_2$ is the detector mask and is placed immediately before the detector, with a pixel size $P$. S is the X-ray source. The distance between the source and $A_1$ is denoted as $l_{so}$, and the distance between $A_1$ and the detector is denoted as $l_{od}$.
Consider a 2D object whose absorption and refractive index decrement distributions are denoted as $\beta (x, y)$ and $\delta (x, y)$ respectively. When $l_{so}$ is sufficiently large, the incident X-rays will travel along the direction of $x_r$-axis. Let $\mathcal {H}$ denote the 2D Radon transform operator that maps the object function into a sinogram, which is a function representing the line integrals through the object at all possible rotation angles and detector positions. The 2D Radon transform of $\beta (x, y)$ can be described as
where $L(y_r,\theta )$ is the line parallel to the $x_r$-axis that goes through $(0, y_r)$. DefineThe unitless $I_{\theta, y_r}(\beta (x, y),\delta (x, y))$ will denote the measured intensity normalized by the incident X-ray intensity at tomographic view angle $\theta$ and detector location $y_r$. For a thin biological sample, the canonical imaging model for EIXPCT is given by [24]
2.2 Discretized imaging model
The above section describes a continuous-to-continuous (C-C) imaging model. In practice, however, measurement data are described by vectors instead of continuous functions. In addition, reconstruction by digital computers requires the object functions to be modeled as vectors as well. Thus a discrete-to-discrete (D-D) model forms the basis for algebraic or iterative image reconstruction.
Consider that the absorption and refractive index decrement distributions of the object are sampled at N vertices, so the discretized object can be represented by the vectors $\boldsymbol \beta \in \mathbb {R}^N$ and $\boldsymbol \delta \in \mathbb {R}^N$. The object is assumed to be imaged at $P$ tomographic view angles. The collection of mask positions at all view angles can be denoted as $\boldsymbol {\Delta \xi } \in \mathbb {R}^P$, where $[\boldsymbol {\Delta \xi }]_i$ ($i=1,2,\ldots, P$) denotes the mask position for the $i$-th view angle. A collection of $Q$ detector pixels record $Q$ measured wavefield intensities for each view angle. Concatenating the measured intensities from all views and detectors creates a measurement vector $\mathbf {I}$, where $\mathbf {I}_{k}$ ($k=(i-1)Q+(j-1)$; $i=1,2,\ldots, P$; $j=1,2,\ldots, Q$) denotes the measured intensity at the $i$-th view angle and the $j$-th detector pixel. For a fixed data-acquisition design, the measurement vector can also be simulated. The simulated measurement vector is a function of $\boldsymbol \beta$ and $\boldsymbol \delta$. Thus it is denoted by $\mathbf {I}(\boldsymbol \beta ,\boldsymbol \delta )$. Denote the discrete representations of the 2D Radon transform and its first-order derivatives as $\mathbf {H}\in \mathbb {R}^{PQ\times N}$ and $\mathbf {D}\in \mathbb {R}^{PQ\times N}$, respectively. The D-D imaging model can then be expressed as [21]
2.3 Joint reconstruction (JR) method
The JR method seeks to estimate the distributions of absorption and refractive index decrement simultaneously by solving a penalized least square optimization problem [21] of the form:
To solve the optimization problem, the fast iterative soft-thresholding algorithm (FISTA) with backtracking is employed [25]. The algorithm is stopped when the step size determined by the backtracking is smaller than an empirically set threshold.
3. Methods
3.1 Single-shot mask displacement strategies and data-acquisition designs
Numerical studies were conducted to assess how data-acquisition designs affect image reconstruction quality when the JR method is employed. Different data-acquisition designs were considered that utilized distinct strategies to specify the mask displacement $\Delta \xi _\theta$ at each view angle. These strategies are referred to as mask displacement strategies. Previously, two single-shot mask displacement strategies have been proposed: the constant mask position (CAP) strategy, and the alternating mask position (AAP) strategy. (The masks were originally called apertures, and thus the “A” in the abbreviations stands for the word “aperture”. Now, we generally call them masks instead of apertures, but we still keep the same abbreviation to keep consistency.) This section reviews the two strategies and proposes a new single-shot mask displacement strategies: the piece-wise constant mask position (PCAP) strategy. The mask displacement strategy used for the conventional multiple-shot two-step reconstruction method will be referred to as the conventional (CVN) strategy.
There exist two mask positions for capturing highest contrast refractive index information, and they are located on opposite sides of the illumination curve. These positions correspond to highest absolute values of the derivatives of the illumination curve [24] and will be referred to as Position 1 and 2 below. In the CVN strategy, images are acquired at both of these positions at all tomographic view angles. The single-shot mask displacement strategies require acquisition of only one image at either of these two mask positions at each tomographic view angle. However, they vary in how the mask positions are chosen at each tomographic view angle. In the CAP strategy, the mask position is fixed and does not vary with tomographic view angle. In the AAP strategy, the mask position alternates between the two positions as a function of tomographic view angle. The PCAP strategy requires a more complicated data-acquisition process. The mask is placed at position 1 for $v$ consecutive view angles and $v$ images acquired. Then, it is switched to position 2 for the next $v$ view angles and another $v$ images taken. Then, it is switches back to position 1 for another $v$ view angles, and so on. In this study, two different $v$ values are tested, yielding two versions of PCAP strategies: PCAP60 for $v=60$ and PCAP180 for $v=180$. The mask displacement strategies are illustrated in Fig. 2.
Fig. 2. Different mask displacement strategies: (1) CVN strategy, (2) AAP strategy, (3) PCAP strategy, (4) CAP strategy. Each asterisk indicates one image acquired at position 1 or position 2 at the specified tomographic view angle.
Mask displacement strategies specify the general data-acquisition strategy. However, to uniquely describe a data-acquisition design, the angular step ($\Delta \theta$), number of tomographic view angles acquired ($N$), and total angular scanning range ($\theta _{a}$) need to be specified. For clarification, we denote a data-acquisition design as a combination of a mask displacement strategy together with the three acquisition parameters ($\Delta \theta , N, \theta _{a}$). For the single-shot designed considered here, the three parameters have the relationship:
Therefore, only two independent parameters are needed. In this study, the angular step is selected to be 0.5° in all test cases. Accordingly, a data-acquisition design can be fully determined by a mask displacement strategy and the angular scanning range. Data-acquisition designs will be denoted by the notation strategy($\theta _{a}$). For example, PCAP60($1\pi$) specifies the PCAP60 design with angular scanning range of $1\pi$. Given the fact that $\Delta \theta =0.5$, in total, 360 projection images are acquired during this data-acquisition process. The mask position is changed every $60$ views, or equivalently, every 30°.These data-acquisition designs would lead to varying data-acquisition times. In general, it is expected that, given a fixed angular step size and mask displacement strategy, smaller total angular scanning range will require less data-acquisition time; given a fixed angular step size and total angular scanning range, a strategy that moves the mask less-frequently would help shorten data-acquisition time. Which design would lead to the shortest total data-acquisition time is highly dependent on the hardware design and data acquisition implementation. The hardware selected determines the rotation speed, mask movement speed, and the required minimum exposure time per shot; even if the hardware is the same, we can determine the sequence of how the images are taken. The object keeps rotating for the total angular scanning range while the mask is moved according to the mask displacement strategy. Another option would be to let the object rotate for entire angular scanning range twice. The mask is fixed at position 1 during the first rotation session and then is fixed at position 2 during the second rotation session. During each rotation, exposures are only acquired at a subset of view angles required by the mask displacement strategy. For our experimental setup, the mask displacement and rotation times are much shorter than the exposure time. Thus given the same number of tomographic view angles, the increase in data-acquisition time of employing the AAP or PCAP strategy compared to the CAP strategy is minor, regardless of the sequence of image acquisition. On the other hand, increasing the number of tomographic view angles may significantly increase the total data-acquisition time.
3.2 Simulation study
3.2.1 Phantom designs
3.2.1.1 Digital and analytical phantoms
Both digital and analytic object phantoms were employed. Here, a digital phantom was modeled as a 2-dimensional (2D) matrix. The line integral data for a digital phantom was implemented as the weighted sum of the pixels’ values along the line trajectory, where the weight for any specific pixel was proportional to the line segment length inside the pixel. Analytical phantoms were formed by a collection of simple geometries where all parameters needed to determine the shapes and locations of the geometries were specified. The line integral data for an analytical phantom were calculated analytically.
3.2.1.2 First phantom: Analytical geometry phantom
The geometry phantom was an artificial phantom made up of simple shapes, such as ellipses and rectangles, as shown in Fig. 3(a) and Fig. 3(b). The largest circle was assigned $\delta$ and $\beta$ values representative of adipose tissue and all the small structures inside were assigned the same $\delta$ and $\beta$ values representative of soft tissues [26]. The bars at the top of the phantom were of dimension $1$ mm$\times 200\,\mu$m, and the bar at the bottom was $1$ mm$\times 300\,\mu$m. The ellipses-shaped phantom’s major axis was 2 cm, and the minor axis was 1.4 cm. The small ellipse’s major axis was 0.8 mm, and the minor axis was 0.6 mm.
Fig. 3. The phantoms for (a) (c) absorption distributions and (b) (d) refractive index decrement distributions. (a) and (b) correspond to the geometry phantom; (c) and (d) correspond to the digimouse phantom.
3.2.1.3 Second phantom: digital digimouse phantom
A digital digimouse phantom [27,28] employed in this study is shown in Fig. 3(c) and Fig. 3(d). This phantom had a pixel size of $100\,\mu$m$\times 100\,\mu$m.
3.2.2 Simulated measurement data
For the analytical phantom, the measurement data were simulated based on Eq. (3), where an analytical integration step was implemented to compute the Radon transformation. This simulation was closest to a real-world imaging process, where the objects are inherently continuous. For the discrete phantoms, the discrete integration replaced the analytical integration, and all other steps remained the same. To produce noisy measurement data, an uncorrelated Gaussian noise model was employed. The noise model had zero-mean and a standard deviation equal to 1% of the noise-less simulated measurement intensity was employed.
3.2.3 Data-acquisition design stability study
In a previous work [22], the optimization problem in Eq. (5) was shown to be non-convex. In this component of our study, the ability of the JR method to produce accurate image estimates when different data-acquisition designs are employed was investigated.
3.2.3.1 Effect of tomographic scanning range
As pointed out in Sec. 3.1, a data-acquisition design requires specifying the mask displacement strategy as well as the tomographic angular scanning range. The use of different mask displacement strategies may result in differing requirements over the tomographic angular scanning range. In this part of the study, the minimum required angular scanning range for each mask displacement strategy was investigated in simulation studies with idealized measurement data. Digital digimouse phantoms were employed as test objects. For each candidate mask displacement strategy (CAP, AAP, PCAP60, and PCAP180), five sets of simulated measurement data were produced, corresponding to angular range coverages of 180°, 225°, 270°, 315°, and 360°. Reconstructed images were obtained from noise-free measurement data by use of the JR method with no penalty employed. The algorithm was stopped when the stopping criterion was met or at a fixed maximum amount of iterations, whichever was earlier. The reconstructed images were analyzed based on their perceived visual quality as well as on quantitative measurements of their mean square error (MSE) and structural similarity (SSIM) [29]. Smaller MSE values and higher SSIM values indicated more accurate reconstructed images. Promising data-acquisition designs were selected for more detailed analysis.
3.2.3.2 Cross-talk for each design
The reconstructed $\boldsymbol \delta$ and $\boldsymbol \beta$ images can be interpreted as two reconstruction channels [30]. Cross-talk refers to image artifacts that occur when the structures in one channel erroneously appear in the other channel. Cross-talk is an important potential problem for any JR method and merits careful study [30]. For the data-acquisition designs identified as promising in the angular range study, cross-talk was assessed. Two sets of cross-talk studies were carried out, with both studies using analytical geometry phantoms.
In the first set of studies, we investigated how errors in one channel could affect the quality of the other channel. In the iterative reconstruction process of the JR method, all subsequent iterations are affected by intermediate reconstruction results. Thus, this study reveals how one channel might evolve when the other channel suffers from suboptimal quality from noise or system inconsistency. Noisy simulated measurement data (1% Gaussian noise) were employed. The noise level was selected to roughly match that in actual experiments.
Two single-channel reconstruction methods were defined:
In the second set of studies, the gradients computed during the reconstruction process were analyzed. If the gradient of $\boldsymbol \delta$ computed when solving Eq. (5) is almost independent of $\boldsymbol \beta$, and vise versa, the cross-talk effect is expected to be small. For different data-acquisition designs, we calculated gradients of the $\boldsymbol \delta$ channel when $\boldsymbol \beta$ was fixed, and analyzed how the value of $\boldsymbol \beta$ affects the $\boldsymbol \delta$ gradient. To be specific, $\boldsymbol \delta$ was selected to be an all-zero vector at the beginning, while $\boldsymbol \beta$ took different initial guesses:
- • All-zero vector, whose values were all 0s everywhere;
- • True phantom vector, which was the ground truth phantom;
- • Noisy vector as described above;
- • Blurred vector, which was a blurred version of the ground truth phantom. The blurring kernel was a Gaussian function with a standard deviation of $2\times 10^{-4}\,$m.
3.2.4 Data-acquisition design quality assessment
The candidate data-acquisition designs were compared based on the quality of the reconstructed images. The noise level and resolution of the reconstructed images were analyzed by the MSE analysis and reconstructed profiles for the bar region, respectively.
3.3 Experimental validation
Experimental data acquired in a previous study [20] were employed for comparing the data-acquisition designs selected in the computer-simulation studies. The source was a Rigaku MicroMax 007 HF X-ray tube with a molybdenum rotating target. The detector was a Hamamatsu C9732DK flat panel detector. The imaged object was a custom-built multiple-component phantom consisting of a hollow plastic cylinder, a plastic rod, a smaller hollow plastic cylinder, rolled plastic paraffin film, and chalk inside a plastic container. A more detailed description of the imaging study and parameters can be found in [20]. Unlike in the previous study, however, ten-step dithering was employed. Each exposure took 1 second to complete. Due to the dithering, 10s was spent on radiation exposure per tomographic view angle, and the pixel size in the reconstructed images was $7.9$ $\mu \,$m. The dataset was originally obtained from $720$ equally spaced tomographic view angles covering a $2\pi$ angular range, and two images were separately obtained at both of the optimal mask positions for each view angle. The original dataset was subsampled to mimic different data-acquisition designs.
4. Results: Data-acquisition stability study
4.1 Data-acquisition design tomographic scanning requirements in idealized situations
The reconstructed estimates of $\boldsymbol \delta$ for different mask displacement strategies and angular scanning ranges are shown in Fig. 4. The absorption distribution is of less interest for XPCT imaging, and is thus not shown here. However, it was found that the analysis of the reconstructed refractive index decrement distributions and the absorption distributions led to the same conclusions regarding data-acquisition designs.
Fig. 4. The reconstructed estimates of $\boldsymbol \delta$ for four different mask displacement strategies (columns) with different angular ranges (rows). From left to right, the mask position is changed less and less frequently, while the five rows correspond to 180°, 225°, 270°, 315°, and 360°, respectively. All images are displayed in the same gray-scale window. The panel at the bottom-right corner for each reconstructed image is a gray-scale-adjusted view of a uniform region in the reconstructed object, so that fluctuations in intensity can be better visualized. The image quality is higher when the angular range is larger or the masks are changed more frequently. Increasing angular range would increase dose, while moving the mask would not.
The MSE and SSIM values were calculated for the reconstructed $\boldsymbol \beta$ and $\boldsymbol \delta$, respectively. The relationship between the MSE/SSIM and the angular coverage is shown in Fig. 5.
Fig. 5. The MSE (top row) and SSIM (bottom row) for the reconstructed refractive index decrement distribution $\boldsymbol \delta$ and absorption distribution $\boldsymbol \beta$ for data-acquisition designs using different mask displacement strategies and tomographic angular scanning ranges.
For certain data-acquisition designs, image artifacts may appear. In general, the numerical studies suggested that the less frequently the mask is moved and the smaller the total angular coverage, the more artifacts appear in the final reconstructed images. This effect is most prominent with the CAP($1\pi$) design, as shown in Fig. 4. For some other designs, although the image is free from visually obvious artifacts, it can still contain artifacts. To observe this, a small panel is displayed for each data-acquisition design with a narrower gray-scale window, showing a small region that is supposed to be uniform; this region is marked with the red square in the first image of Fig. 4. Minor artifacts can be better visualized from the small panel. The quantitative metrics demonstrate a similar trend. Based on Figs. 4 and 5, when the angular coverage increases, all data-acquisition designs result in reduced MSE values and increased SSIM values, which is expected because more information is obtained at the cost of increased radiation time and doses.
The numerical results suggest that, the more frequently the mask position is changed, the more likely the data-acquisition design is to obtain an artifact-negligible reconstructed image with a smaller angular range coverage. On one extreme, when the AAP strategy is employed, an angular range of $1\pi$ was observed to be adequate for obtaining artifact-negligible images. On the other extreme, when the CAP strategy is employed, a $2\pi$ angular range was required for artifact-negligible images. The PCAP strategies stand in between. Thus, changing the mask position helps reduce the doses.
Since the reconstruction results were obtained from noise-free simulated measurement data utilizing a digital phantom, they reflect the upper-bound performance limit of all these data-acquisition designs. Although many potential designs work well in the idealized setting, only three representative designs were selected for the following in-depth study of the properties: CAP($2\pi$), AAP($2\pi$), and AAP($1\pi$). The differences between the reconstructed images utilizing AAP($2\pi$) and CAP($2\pi$) designs are representative of the effect of decreasing mask movement frequency, and the differences between those utilizing AAP($2\pi$) and AAP($1\pi$) show the effect of decreasing the total angular range. The performance of other data-acquisition designs can be deduced based on the trends revealed by the three representative designs.
4.2 Cross-talk effect
4.2.1 The effect of imperfection of one channel on the another channel
For the single-channel reconstruction problem in Eq. (8), the fixed estimates of $\boldsymbol \delta$ as detailed in Sec. 3.2.3 and the corresponding reconstructed $\boldsymbol \beta$ are shown in Fig. 6(a). For the single-channel reconstruction problem in Eq. (9), the fixed estimates of $\boldsymbol \beta$ and their corresponding reconstructed estimates of $\boldsymbol \delta$ are shown in Fig. 6(b).
Fig. 6. (a) Single-channel reconstructed estimates of $\boldsymbol \beta$ when imperfect values of $\boldsymbol \delta$ are assumed. The first image shows the imperfect assigned values for $\boldsymbol \delta$. The other three images show the reconstructed $\boldsymbol \beta$ images. Each corresponds to a data-acquisition design and a specific assigned $\boldsymbol \delta$. (b) Single-channel reconstructed estimates of $\boldsymbol \delta$ when imperfect values of $\boldsymbol \beta$ are assumed. The first image shows the imperfect assigned values for $\boldsymbol \beta$. The other three images shows the reconstructed $\boldsymbol \delta$ images. Each corresponds to a data-acquisition design and a specific assigned $\boldsymbol \beta$.
Based on the reconstructed images, although the assigned values for the other channel may not be ideal, the to-be-reconstructed channel was still well estimated. The only exception is for AAP($1\pi$) and a noisy vector assigned for either $\boldsymbol \beta$ or $\boldsymbol \delta$, where the reconstructed channel is still noisy. This problem, however, can be alleviated by increasing the value of regularization parameter. In conclusion, the quality of one channel does not exert a significant influence on the other channel in the CAP($2\pi$) and AAP($2\pi$) designs. For the AAP($1\pi$) design, the effect exists but can be mitigated by imposing higher regularization.
4.2.2 Impact of crosstalk on cost function gradients
To further understand the cross-talk, the gradients were visualized. For most cases, the gradients of one channel are very similar to each other when the other channel employs different estimates. The $\boldsymbol \delta$ gradient is only affected by a noisy $\boldsymbol \beta$ vector when the AAP($1\pi$) data-acquisition design is employed, which is best illustrated in the reduced conspicuity of small structures in the AAP($1\pi$) gradient image. This observation suggests that the AAP($1\pi$) design can be affected by cross-talk. In Fig. 7, the gradients of $\boldsymbol \delta$ with $\boldsymbol \beta$ being the noisy vector are compared with those with $\boldsymbol \beta$ being the all-zero vector. The other $\boldsymbol \delta$ gradient images when $\boldsymbol \beta$ is the true phantom vector or the blurred vector are not shown here, but provide results similar to those of the all-zero vector. This finding is consistent with the conclusion from the previous cross-talk effect test, suggesting that the AAP($1\pi$) is more easily affected by cross-talk, indicating a higher potential for getting trapped in local minima. Thus, special care needs to be taken when choosing the regularization parameters to stabilize the reconstruction.
Fig. 7. The gradient with respect to $\boldsymbol \delta$ when $\boldsymbol \delta =0$ and different values are assigned for $\boldsymbol \beta$: the first row shows when $\boldsymbol \beta =0$; the second row shows $\boldsymbol \beta$ being the noisy vector defined in Sec.3.2.3.
5. Results: Data-acquisition quality assessment study
5.1 Effect of regularization parameters on image quality
The regularization parameters can significantly affect image quality. Consequently, the effect of regularization parameters is considered first in this section. The trend is similar for all designs, but only the reconstructed images corresponding to the AAP($2\pi$) design are shown here.
The TV regularization parameters for the $\boldsymbol \beta$ and $\boldsymbol \delta$ distributions were swept on a 2D grid. Five representative regularization values were empirically selected for each of the parameters. The values employed for AAP($2\pi$) and CAP($2\pi$) designs are summarized in Table 1. For the AAP($1\pi$) design, the angular coverage range is only half that of the other two designs, thus the data fidelity term is also approximately half that of the other designs covering a $2\pi$ angular range. Thus, the candidate regularization parameters were also selected to be half of those for the other two designs.
Table 1. Candidate regularization parameters. In Fig. 8, the reconstructed images are only shown for underlined parameters.
A zoomed-in view of the bar region of the reconstructed $\boldsymbol \beta$ and $\boldsymbol \delta$ estimates corresponding 9 sets of representative regularization parameters are shown in Fig. 8(a) and Fig. 8(b), and a similarly zoomed-in view of the circle region is shown in Fig. 8(c) and Fig. 8(d).
Fig. 8. The bar region (a, b) and circle region (c, d) of the reconstructed $\boldsymbol \beta$ (a, c) and $\boldsymbol \delta$ (b, d) image of the analytical geometric phantom for the AAP($2\pi$) design with selected regularization parameters from Table 1. From left to right, the regularization parameter for $\boldsymbol \beta$ is increasing. From top to bottom, the regularization parameter for $\boldsymbol \delta$ is increasing.
The reconstructed images are intuitively consistent, in that higher regularization reduces noise, but setting the parameter too high can blur fine structures and reduce contrast. Moreover, the visual quality of the reconstructed absorption is mainly affected by the TV regularization parameter for $\boldsymbol \beta$, and that of the reconstructed refractive index decrement is mainly affected by the TV regularization parameter for $\boldsymbol \delta$. As a result, the parameters can be tuned separately.
Regularization parameters were chosen as a best compromise among their noise suppression effect, resolution, and contrast: $\alpha _1=1\times 10^6$ for $\boldsymbol \beta$ and $\alpha _2=1000$ for $\boldsymbol \delta$. The same procedure was repeated for the CAP($2\pi$) and AAP($1\pi$) data-acquisition designs, and the same set of regularization parameters were determined to be optimal for CAP($2\pi$). The optimal parameters for AAP($1\pi$) are $\alpha _1=5\times 10^5$ for $\boldsymbol \beta$ and $\alpha _2=500$ for $\boldsymbol \delta$.
5.2 Comparing data-acquisition designs based on image quality
The reconstructed images obtained by the identified regularization parameters in each candidate design were employed to compare the data-acquisition designs. A zoomed-in view of the bar region is given in Fig. 9(a), and that for the circle region is given in Fig. 9(b).
Fig. 9. The bar region (a) and the circle region (b) of the reconstructed absorption distribution (top) and refractive index decrement distribution (bottom) of the analytical geometric phantom for the CAP($2\pi$), AAP($2\pi$), and AAP($1\pi$) designs (from left to right), with an optimal regularization parameter. (c) The averaged profile of the bar region of the reconstructed absorption distribution (left) and refractive index decrement distribution (right) of the analytical geometric phantom for the CAP($2\pi$), AAP($2\pi$), and AAP($1\pi$) designs, with an optimal regularization parameter.
Judging from the visual appearance of the bar and circle regions, all of the designs resulted in images that revealed the bars and the circle clearly, with the reconstructed $\boldsymbol \delta$ image corresponding to the AAP($1\pi$) design having a higher level of noise than the other two designs. The MSE of the reconstructed images, as summarized in Table 2, also supports this finding. In addition, the shape of the boundary of the bar region is more distorted, indicating poorer reconstruction quality.
Table 2. MSEs for reconstructed images with the optimal regularization parameters
To better analyze the resolution, the profiles of the bar region along the vertical direction were averaged to produce the averaged profile shown in Fig. 9(c). The averaging operation minimizes the effect of noise on the analysis of the resolution. According to the profiles, the resolutions obtained by all three candidate designs are similar.
In summary, the image obtained by the AAP($1\pi$) design shows slightly inferior performance in terms of resistance to measurement data noise, but its resolution is mainly unaffected by the data acquisition designs.
6. Experimental verification
The JR method was applied to the experimentally acquired data and the reconstructed images are shown in Fig. 10. Consistent with the simulation studies, the AAP($1\pi$) design was more affected by noise, but the structure shapes can all be reliably reconstructed, which is clearly demonstrated especially by the fine details in the paraffin film region.
Fig. 10. The reconstructed absorption (first row) and refractive index decrement (second row) distributions using three data-acquisition designs. The third row shows a zoomed-in region of the reconstructed chalk section with varying data-acquisition designs.
Note that there are dark artifacts in the center of the reconstructed refractive index decrement image corresponding to the CAP($2\pi$) data-acquisition design, which are possibly due to the non-uniformity of the illumination curves across the whole field of view. Techniques such as local phase retrieval [31] may alleviate these artifacts. In comparison, the AAP designs benefit from the diversity in their data that is caused by frequent movement of the sample mask, and thus are less affected by the model mismatch problem in practice.
7. Conclusion
This work investigated a series of flexible data-acquisition designs for EIXPCT enabled by a JR method. The designs differ in terms of their data-acquisition simplicity and reconstruction stability and quality. The data-acquisition simplicity mainly refers to whether the mask position needs to be moved. Three representative designs were selected for the low artifact level of their reconstructed images from an idealized noise-free simulated measurement dataset. The three designs were further investigated for cross-talk, reconstructed image noise level, resolution, and ultimately tested on real experimental data. Table 3 summarizes the data-acquisition properties of the three single-shot methods.
Table 3. A comparison of the three single-shot methods (in the table, V is the number of views obtained for $\pi$ angular range)
Different data-acquisition designs have different advantages. The CAP($2\pi$) design has the simplest data-acquisition procedure because the masks do not need to be moved, but it may be more affected by model mismatch in practice. The AAP($2\pi$) design demonstrated the highest fidelity reconstructed images and is less affected by the model mismatch problem. However, it requires extra mask movement and therefore makes the data-acquisition process more complex. The AAP($1\pi$) design requires the minimum amount of images due to its small angular coverage range and single-shot design. However, images reconstructed from data acquired by the AAP($1\pi$) design can be susceptible to noise amplification. Ultimately, the selection of a data-acquisition design should be optimized based on the specific needs of the application. If we want to obtain highest quality and no constraint of time is available, such as imaging a detail-rich object, the conventional method would be ideal. If high quality is desired but there is still some time constraint, such as imaging the detail-rich object which may slowly dehydrate overtime, AAP($2\pi$) design would be desired. If the imaging time restriction is tighter, such as in an in vivo study, an AAP($1\pi$) would be the best fit. If system simplicity is most heavily appreciated, CAP($2\pi$) would help save some effort finding only one optimal mask location instead of two.
Funding
National Institutes of Health (EB020604, EB023045).
Acknowledgments
CKH is supported by the Royal Academy of Engineering under the RAEng Research Fellowship scheme.
References
1. C. Kottler, F. Pfeiffer, O. Bunk, C. Grünzweig, and C. David, “Grating interferometer based scanning setup for hard X-ray phase contrast imaging,” Rev. Sci. Instrum. 78(4), 043710 (2007). [CrossRef]
2. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2(4), 258–261 (2006). [CrossRef]
3. P. Zhu, K. Zhang, Z. Wang, Y. Liu, X. Liu, Z. Wu, S. A. McDonald, F. Marone, and M. Stampanoni, “Low-dose, simple, and fast grating-based X-ray phase-contrast imaging,” Proc. Natl. Acad. Sci. 107(31), 13576–13581 (2010). [CrossRef]
4. S. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature 384(6607), 335–338 (1996). [CrossRef]
5. P. Diemoz, A. Bravin, and P. Coan, “Theoretical comparison of three X-ray phase-contrast imaging techniques: propagation-based imaging, analyzer-based imaging and grating interferometry,” Opt. Express 20(3), 2789–2805 (2012). [CrossRef]
6. A. Ruhlandt, M. Krenkel, M. Bartels, and T. Salditt, “Three-dimensional phase retrieval in propagation-based phase-contrast imaging,” Phys. Rev. A 89(3), 033847 (2014). [CrossRef]
7. A. Olivo and R. Speller, “A coded-aperture technique allowing X-ray phase contrast imaging with conventional sources,” Appl. Phys. Lett. 91(7), 074106 (2007). [CrossRef]
8. A. Zamir, C. Hagen, P. C. Diemoz, M. Endrizzi, F. Vittoria, Y. Chen, M. A. Anastasio, and A. Olivo, “Recent advances in edge illumination X-ray phase-contrast tomography,” J. Med. Imag. 4(04), 1 (2017). [CrossRef]
9. A. Olivo and R. D. Speller, “A novel X-ray imaging technique based on coded apertures making phase contrast imaging feasible with conventional sources,” in 2008 IEEE Nuclear Science Symposium Conference Record, (IEEE, 2008), pp. 1447–1450.
10. T. Davis, D. Gao, T. Gureyev, A. Stevenson, and S. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature 373(6515), 595–598 (1995). [CrossRef]
11. P. R. Munro, K. Ignatyev, R. D. Speller, and A. Olivo, “The relationship between wave and geometrical optics models of coded aperture type X-ray phase contrast imaging systems,” Opt. Express 18(5), 4103–4117 (2010). [CrossRef]
12. P. C. Diemoz and A. Olivo, “On the origin of contrast in edge illumination X-ray phase-contrast imaging,” Opt. Express 22(23), 28199–28214 (2014). [CrossRef]
13. M. Endrizzi, P. C. Diemoz, T. P. Millard, J. Louise Jones, R. D. Speller, I. K. Robinson, and A. Olivo, “Hard X-ray dark-field imaging with incoherent sample illumination,” Appl. Phys. Lett. 104(2), 024106 (2014). [CrossRef]
14. C. K. Hagen, P. C. Diemoz, M. Endrizzi, L. Rigon, D. Dreoosi, F. Arfelli, F. C. Lopez, R. Longo, and A. Olivo, “Quantitative edge illumination X-ray phase contrast tomography,” in Developments in X-Ray Tomography IX, vol. 9212 (International Society for Optics and Photonics, 2014), p. 921205.
15. C. K. Hagen, P. C. Diemoz, and A. Olivo, “On the relative performance of edge illumination X-ray phase-contrast CT and conventional, attenuation-based CT,” Med. Phys. 44(5), 1876–1885 (2017). [CrossRef]
16. C. K. Hagen, M. Endrizzi, P. C. Diemoz, and A. Olivo, “Reverse projection retrieval in edge illumination X-ray phase contrast computed tomography,” J. Phys. D: Appl. Phys. 49(25), 255501 (2016). [CrossRef]
17. P. C. Diemoz, F. A. Vittoria, C. K. Hagen, M. Endrizzi, P. Coan, E. Brun, U. H. Wagner, C. Rau, I. K. Robinson, and A. Bravin, “Single-image phase retrieval using an edge illumination X-ray phase-contrast imaging setup,” J. Synchrotron Radiat. 22(4), 1072–1077 (2015). [CrossRef]
18. M. Das and Z. Liang, “Spectral X-ray phase contrast imaging for single-shot retrieval of absorption, phase, and differential-phase imagery,” Opt. Lett. 39(21), 6343–6346 (2014). [CrossRef]
19. M. Endrizzi, A. Astolfo, F. A. Vittoria, T. P. Millard, and A. Olivo, “Asymmetric masks for laboratory-based X-ray phase-contrast imaging with edge illumination,” Sci. Rep. 6(1), 25466 (2016). [CrossRef]
20. P. Diemoz, C. Hagen, M. Endrizzi, M. Minuti, R. Bellazzini, L. Urbani, P. De Coppi, and A. Olivo, “Single-shot X-ray phase-contrast computed tomography with nonmicrofocal laboratory sources,” Phys. Rev. Appl. 7(4), 044029 (2017). [CrossRef]
21. Y. Chen, H. Guan, C. K. Hagen, A. Olivo, and M. A. Anastasio, “Single-shot edge illumination X-ray phase-contrast tomography enabled by joint image reconstruction,” Opt. Lett. 42(3), 619–622 (2017). [CrossRef]
22. Y. Chen and M. A. Anastasio, “Properties of a joint reconstruction method for edge-illumination X-ray phase-contrast tomography,” Sens. Imaging 19(1), 7 (2018). [CrossRef]
23. Y. Chen, W. Zhou, and M. A. Anastasio, “Joint-reconstruction-enabled data acquisition design for single-shot edge-illumination X-ray phase-contrast tomography,” in Medical Imaging 2018: Physics of Medical Imaging, vol. 10573 (International Society for Optics and Photonics, 2018), p. 1057322.
24. P. R. Munro, C. K. Hagen, M. B. Szafraniec, and A. Olivo, “A simplified approach to quantitative coded aperture X-ray phase imaging,” Opt. Express 21(9), 11187–11201 (2013). [CrossRef]
25. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. on Imaging Sci. 2(1), 183–202 (2009). [CrossRef]
26. J. H. Hubbell and S. M. Seltzer, “Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients 1 keV to 20 MeV for elements Z= 1 to 92 and 48 additional substances of dosimetric interest,” Tech. rep., National Inst. of Standards and Technology-PL, Gaithersburg, MD (United States). Ionizing Radiation Div. (1995).
27. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007). [CrossRef]
28. D. Stout, “Creating a whole body digital mouse atlas with PET, CT and cryosection images,” Mol. Imaging Biol. 4(1), 27–33 (2002). [CrossRef]
29. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE T. Image Process. 13(4), 600–612 (2004). [CrossRef]
30. A. Sawatzky, Q. Xu, C. O. Schirra, and M. A. Anastasio, “Proximal ADMM for multi-channel image reconstruction in spectral X-ray CT,” IEEE T. Med. Imaging 33(8), 1657–1668 (2014). [CrossRef]
31. A. Zamir, M. Endrizzi, C. K. Hagen, F. A. Vittoria, L. Urbani, P. De Coppi, and A. Olivo, “Robust phase retrieval for high resolution edge illumination X-ray phase-contrast computed tomography in non-ideal environments,” Sci. Rep. 6(1), 31197 (2016). [CrossRef]
