research papers
Fast iterative reconstruction of data in full interior tomography
aInstitute for Biomedical Engineering, ETH Zurich, 8092 Zurich, Switzerland, and bSwiss Light Source, Paul Scherrer Institute, 5232 Villigen, Switzerland
*Correspondence e-mail: filippo.arcadu@psi.ch
This paper introduces two novel strategies for iterative reconstruction of full interior tomography (FINT) data, i.e. when the field of view is entirely inside the object support and knowledge of the object support itself or the attenuation coefficients inside specific regions of interest are not available. The first approach is based on data edge-padding. The second technique creates an intermediate virtual sinogram, which is, then, reconstructed by a standard iterative algorithm. Both strategies are validated in the framework of the alternate direction method of multipliers plug-and-play with gridding projectors that provide a speed-up of three orders of magnitude with respect to standard operators implemented in real space. The proposed methods are benchmarked on synchrotron-based X-ray tomographic microscopy datasets of mouse lung alveoli. Compared with analytical techniques, the proposed methods substantially improve the reconstruction quality for FINT underconstrained datasets, facilitating subsequent post-processing steps.
Keywords: X-ray tomographic microscopy; interior tomography; reconstruction artifacts; iterative tomographic reconstruction; gridding method.
1. Introduction
The term interior tomography (INT) refers to the problem of reconstructing an object function, when its support (S) is not a subset of the field of view (FOV) (Natterer & Wübbeling, 2001). This imaging modality is broadly applied in medical screening, material science and biology, as it allows high-resolution investigations of small regions of interest (ROIs).
Since INT projections contain information FOV, reconstructions with filtered backprojection (FBP) suffer from a DC shift and low-frequency artifacts (Natterer & Wübbeling, 2001), that compromise the quantitativeness of the results and make further post-processing, visualization and rendering of the investigated object complicated (Xiao et al., 2007).
Two different analytical reconstruction methods can address this problem: the differentiated backprojection (DBP) and FBP of edge-padded projections (FBP-E). The idea of DBP (Noo et al., 2004) is to backproject the derivative of the data and, then, recover the original object by Hilbert transform techniques. This method provides quantitative reconstructions, when S is known and only specific geometries are involved: two opposite boundaries of the object are inside the FOV (Noo et al., 2004); the FOV exceeds the object only from one side (Defrise et al., 2006); the FOV and the attenuation coefficients are known in a subregion FOV (Ye et al., 2007; Kudo et al., 2008; Courdurier et al., 2008). FBP-E (Seger, 2002; Marone et al., 2010; Kyrieleis et al., 2011) corresponds to FBP of projections, that have been padded with the edge pixels. This strategy prevents the appearance of low-frequency artifacts within the reconstructed FOV, although results are not fully quantitative.
Iterative methods have also been proposed for INT datasets, characterized by low photon statistics and/or with a low number of views (e.g. low-dose scan). Usually, these underconstrained datasets cannot be reconstructed with sufficient accuracy with analytical methods. The separable paraboloidal surrogate (SPS) (Xu et al., 2011), the expectation-maximization (EM-ML) (Rashed & Kudo, 2013) and the penalized weighted least square (PWLS) method (Zhang et al., 2014) have been modified to deal with truncated projections. In Xu et al. (2011), the regularized SPS is initialized through the projection-onto-convex-sets (POCS) (Defrise et al., 2006); the knowledge of S, the attenuation coefficients in one or multiple subregions and the position of the background pixels are needed for this reconstruction technique. In Rashed & Kudo (2013), the standard EM-ML is combined with thresholding, acting on pixels of selected ROIs. In Zhang et al. (2014), the PWLS is regularized with a complete high-quality scan of the same object. All these methods work for very specific (mostly medical) applications and require a priori knowledge obtained from previous scans. Iterative reconstruction of INT datasets without knowledge of either the object support or the attenuation coefficients inside specific ROIs (case of `full' interior tomography, abbreviated as FINT) has been initiated by Yu & Wang (2009). In this work, the authors prove that an iterative scheme working with total variation (TV) regularization yields a unique reconstruction of a FINT dataset, in case the object under study is piecewise constant.
At synchrotron imaging beamlines, time-resolved high-resolution investigations of a large variety of dynamic systems spanning different fields [e.g. biology (Walker et al., 2014; Hesse et al., 2015), material sciences (Aagesen et al., 2010; Campi et al., 2015), energy research (Ebner et al., 2013; Lu et al., 2016), earth and environmental sciences (Berg et al., 2012; Ganne et al., 2012)] are becoming routine. Due to the large variability of the examined samples and to the large datasets (several tens of TB) often associated with these experiments, efficient iterative reconstruction algorithms, not bounded to a specific application, are highly demanded.
This paper introduces two fast strategies for iterative reconstruction of FINT data, i.e. the FOV S, S is unknown as well as the attenuation coefficients inside specific ROIs of the FOV. In this regard, the presented results can be considered a step further on the line of research launched by Yu & Wang (2009). Differently from previous studies, this work shows the importance of data pre-processing (either edge-padding or differentiation) for analytical and iterative reconstruction of FINT datasets. The proposed iterative strategies are fast and provide high-quality non-quantitative reconstructions, best suited for subsequent post-processing and analysis steps (e.g. segmentation and morphological studies), independently of the imaged system.
1.1. Contributions
The contributions of this manuscript are summarized as follows:
(i) Analyze the importance of FINT data pre-processing through edge-padding and differentiation for analytical and iterative tomographic reconstruction.
(ii) Introduction of two fast forward gridding projectors for iterative reconstruction of FINT data: one implements the derivative of the Radon transform, the other acts as Radon transform combined with edge-padding (Arcadu et al., 2016a,b).
(iii) Introduction of the virtual strategy, as a more efficient alternative to specific forward projectors for INT data. This strategy is independent of the chosen iterative scheme and regularization.
(iv) Validation of the proposed methods within the framework of the alternate direction method of multipliers plug-and-play (ADMP) (Venkatakrishnan et al., 2013), on simulated and real datasets of full interior X-ray tomographic microscopy.
The edge-padding and virtual strategies studied in this work are neither bounded to the ADMP nor applicable only to the case of piecewise constant objects, but they can be utilized within any kind of iterative reconstruction algorithm and regularization scheme.
2. Reconstruction artifacts in interior tomography
2.1. Preliminaries
This work focuses on the reconstruction of a two-dimensional slice from line projections acquired in parallel beam geometry. However, results can be generalized to more complex tomographic configurations (e.g. fan- and cone-beam).
The collection of line projections, measured at different angles , is called a sinogram. The object to be reconstructed is a finite integrable real function, whose support, S, is assumed to be compact.
In interior tomography, the object support lies outside the FOV, i.e. FOV and FOV . Different INT configurations, listed by Defrise et al. (2006) and displayed in Fig. 1, have been studied in the literature. The case of full interior tomography, often characterizing microtomographic scans (Fig. 1d, FOV ), is considered in this work.
2.2. Case of analytical reconstruction methods
FBP reconstruction of a FINT sinogram is affected by a DC shift and low-frequency artifacts (Seger, 2002; Marone et al., 2010; Kyrieleis et al., 2011), due to the effect of the ramp filter on truncated projections (Seger, 2002). To clarify this point, the example discussed by Seger (2002) is reproposed here.
Fig. 2 shows a complete and truncated projection of a homogeneous circle before and after ramp filtering. In standard tomography (Fig. 2, left), the filtered projection has a constant, positive profile in the middle with symmetric negative tails. Once all projections in are smeared back to the image grid, the circle is filled with constant pixel values and the negative tails of each projection are zeroed-out by the positive contributions from all other projections (Fig. 3b). In interior tomography, the filtered projection (Fig. 2, right) is instead everywhere positive, characterized by a bowl-shaped profile and tails with very high values. After backprojection, these positive tails are still present, because no negative compensation in the reconstruction process is possible, yielding the bowl artifact displayed in Fig. 3(c). The DC shift and bowl-shaped profile, shown in Fig. 3(d), can severely compromise the interpretation/morphological analysis of reconstructed FINT datasets.
2.3. Case of iterative reconstruction methods
Although iterative algorithms do not involve explicit ramp filtering, reconstructions of FINT data are affected by the same artifacts as for analytical methods. Fig. 4 shows an iterative reconstruction with SPS of the same sinogram used for Fig. 3(c) and the corresponding line profile along the central row. The image is once again affected by a DC shift and a bowl artifact (Fig. 4b). The reconstruction with any other iterative algorithm would present the same problems.
This observation can be explained by the fact that iterative methods `mimic' the effect of the ramp filter: the reconstruction after the first iteration corresponds to a blurred version of the object, similar to the result of backprojection when omitting the filtering step; after several iterations, the object becomes sharper, until a similar spatial resolution as for the FBP reconstruction is reached. This `mimicking' behavior leads to the same problems characterizing FBP, when dealing with FINT data.
3. Proposed approach
The proposed approach builds upon existing methods, extended and combined to address the FINT problem. Although here for illustration purposes, specific projectors, iterative scheme and regularization have been chosen, the proposed strategy is more general and not bounded to these choices.
3.1. Flexible iterative reconstruction scheme
The alternate direction method of multipliers plug-and-play (ADMP) (Venkatakrishnan et al., 2013) offers an optimization framework, where the reconstructive and the regularization tasks are neatly separated in two subproblems. This structure allows the direct use of any forward projector for the reconstructive subproblem and any denoising method for the regularization subproblem. The ADMP can be viewed as a particular case of the classical formulation of the alternate direction method of multipliers (Boyd et al., 2011).
In parallel beam geometry, the tomographic problem for a piecewise constant object has the following form,
where m is the number of views, nz is the number of slices, n is the number of pixels along a row or a column (assuming the image grid to be square), is the unknown three-dimensional object, is the unknown ith slice, is the forward projection operator ( is the adjoint operator or backprojector), is the sinogram and R is the functional enforcing the object to be piecewise constant. The dual variable is introduced to transform the constrained optimization problem (1) into
The constraint = is incorporated in the functional by augmenting F with an array of multipliers . In this way, the following Lagrangian is obtained,
The original problem (1) is, therefore, mapped into the minimization of the Lagrangian (3). The ADMP solves (3) by cyclically minimizing two subproblems with respect to the variable and and by updating the multipliers until convergence,
where the superscript (k) refers to the kth iteration of a selected variable. The first and third terms of (3) will contribute to the -subproblem, whereas the second and the third terms will contribute to the -subproblem. To explicitly derive the form of the -subproblem, the gradient of L with respect to is calculated,
Since is symmetric and positive-definite, the conjugate gradient (CG) technique (Hestenes & Stiefel, 1952) is adopted to iteratively find the solution of the system = . For the -subproblem, it results that
The last term in (6) corresponds exactly to a denoising problem, where represents the input noisy image and = is the regularization strength. The form of the -subproblem gives the `plug-and-play' qualification to the ADMP, as the type of denoising can be changed without altering the structure of the entire algorithm.
The iterative procedure is stopped when the L2-norm of the relative difference between reconstructions of subsequent iterations, and , is smaller than a certain threshold, i.e. = 0.01. One of the main advantages of ADMM-type methods is that satisfactory results in terms of image quality can be achieved after very few iterations.
For the experiments presented in §4, ADMP with split Bregman total variation (SBTV) (Goldstein & Osher, 2009) as plug-and-play regularization has been used. τ strongly determines the quality of the iterative reconstruction since it controls the trade-off between spatial resolution and noise removal.
3.2. Fast forward projectors for interior tomography
The DC shift compromises the quantitativity of FINT reconstructions, but, except for an offset constant throughout the tomographic slice, the information is preserved. The bowl artifact leads, instead, to non-quantitative reconstructions characterized by varying grey level values within homogeneous regions and prevents, therefore, any reliable morphological analysis, unless `decupping' algorithms or sophisticated segmentation approaches are utilized.
Two analytical methods can be adopted to avoid this bowl artifact: DBP and FBP-E. Reconstructions with these techniques are non-quantitative, i.e. exact up to an unknown constant, but can be safely used for image analysis if the knowledge of the attenuation coefficients in an absolute sense is not needed, as is often the case. Fig. 5 shows the reconstruction of a FINT sinogram created from a Shepp–Logan (SL) phantom using FBP, DBP and FBP-E. The original SL with 2048 × 2048 pixels (Fig. 5a) is forward projected over 800 angles in and only the central 512 pixels of the sinogram are retained. The bowl artifact is visible at the corners of the FBP reconstruction (Fig. 5c) and in a line profile (Fig. 5f, red). This artifact is instead not present in the DBP and FBP-E reconstructions (Figs. 5d, 5e; Fig. 5f, green, black). For the iterative reconstruction of FINT data, the proposed idea is to use tomographic forward operators based on the principles of DBP and FBP-E, i.e. a projector that implements the derivative of the Radon transform and a Radon transform acting on edge-padded projections, respectively.
To guarantee fast reconstructions, the forward gridding projector (FGP) (Arcadu et al., 2016a), that works in the Fourier space and has a complexity of O(N2log2N), has been selected for this work. Studies conducted by Arcadu et al. (2016a) showed that the FGP is significantly faster than standard space-based projectors [complexity of O(N3)] (Toft, 1996) and its accuracy results comparable with that of very sophisticated operators, when used in iterative schemes.
In standard tomography, the FGP works on an oversampled grid created with zero-padding. For FINT data, the following modifications of the original forward gridding projector are proposed: (i) the differentiated FGP (DFGP) (Arcadu et al., 2016b), that implements the derivative of the Radon transform (still works on an oversampled grid created with zero-padding) and computes a differential sinogram; (ii) the FGP combined with edge-padding (FGP-E), i.e. the oversampled grid is created by edge-padding of the object. If these two operators are used within the ADMP framework, two iterative schemes can be defined: ADMP-D implementing the DFGP and ADMP-E implementing FGP-E.
As a proof of principle, the experiment in Fig. 5 is repeated with the standard ADMP, ADMP-D and ADMP-E. The bowl artifact, appearing in the reconstruction with the standard ADMP (Fig. 6c), is not present in the images computed with ADMP-D (Fig. 6d) and ADMP-E (Fig. 6e).
Edge-padding has two main advantages with respect to the differentiated operator. First, the computation of the differential sinogram by finite differentiation is very likely to enhance the noise affecting the original sinogram; if a sophisticate technique is used instead [like the Savitzky–Golay filter (Savitzky & Golay, 1964)], additional parameters, ruling the trade-off between noise and spatial resolution, are added to the reconstruction problem. This argument is valid for both analytical and iterative reconstructions. Second, the number of sub-iterations required by the CG step in the ADMP is related to the conditioning number (CN) of the operator : for ADMP-D, corresponds to the derivative of the Radon transform; for ADMP-E, is the Radon transform itself, having a lower CN than its derivative. It has been experimentally observed that to reach convergence with the same number of iterations the CG loop of ADMP-D requires at least 15 sub-iterations, whereas four are enough for the CG loop of ADMP-E. For these reasons, only ADMP-E is utilized for the experiments discussed in §4.
3.3. Iterative virtual method
For the iterative reconstruction of FINT data without a priori knowledge of the object support or of the attenuation coefficients inside specific ROIs, an alternative four-step strategy is also proposed: the iterative virtual method (abbreviated as ADMP-V). The steps of ADMP-V are (Fig. 7): (i) data reconstruction by an analytical method, like DBP or FBP-E; (ii) zeroing of all pixels outside the reconstruction circle; (iii) forward projection of this newly computed image to obtain a `virtual' sinogram, simulating a non-FINT dataset; (iv) reconstruction with ADMP using projectors for standard tomography (FGP, in this case) and using physical constraints (e.g. zeroing all pixels outside the reconstruction circle at each iteration). After steps (i), (ii) and (iii), the initial FINT dataset is transformed into a standard tomographic dataset with known support, that can be reconstructed with any analytical or iterative algorithm.
4. Experiments
4.1. Analysis framework
Metrics to assess the reconstruction accuracy are computed over the square inscribed inside the reconstruction circle (SQRES). The contrast-to-noise ratio (CNR) is defined as
where Sri and are the mean and standard deviation of the ROI ri, assumed to be homogeneous. For the computation of the CNR, ri1 and ri2 must be neighbouring. The values of CNR, reported in the tables of §4, represent averages over multiple ROIs at different distances from the image centre.
In the case of simulated data, the peak signal-to-noise ratio (PSNR) (Huyn-Thu & Ghanbaru, 2008) and the mean structural similarity index (MSSIM) (Wang et al., 2004) are also used as figures of merit. Since results are not quantitative due to the DC shift, a linearly regressed version, , of the reconstruction, I, is used as input for the analysis. Calling O the phantom, is computed as
where . In this way, PSNR and MSSIM scores are not biased by the DC shift. Edge profiles and histograms are also displayed to provide additional information about the spatial resolution and the segmentability of the images. Reconstructions are mapped to the interval [0, 1] to facilitate the display and direct comparison of profiles and histograms.
4.2. Data and settings
The phantoms used to generate the simulated data are displayed in Fig. 8. The first dataset is the SL phantom (Shepp & Logan, 1974), which is often utilized to validate tomographic reconstruction algorithms. The second dataset is a segmented reconstruction of mouse lung tissue and it is labelled MLT. This phantom is evidently characterized by a high structural complexity and is related to the real data used in §4.7. The procedure to create FINT sinograms from a reference image has been described in §3.2.
The real datasets have been acquired at the TOMCAT beamline of the Swiss Light Source in the framework of the SNF project `in vivo study of lung physiology with fast X-ray tomographic microscopy' (Lovric et al., 2013). The three sinograms correspond to scans of mouse lung tissue at the micrometre scale with an effective detector pixel size of 2.9 µm (MOUSE-1, MOUSE-2) and 1.1 µm (MOUSE-3). The first two datasets consist of 901 Paganin phase-retrieved (Paganin et al., 2002) projections and 2016 pixels, whereas the third one has 501 projections × 2016 pixels.
The regularization strength τ (the only free parameter in ADMP) is tuned so that reconstructions with ADMP-E and ADMP-V look visually as similar as possible.
4.3. Optimal edge-padding length for analytical and iterative reconstructions
The edge-padding length rules the trade-off between the reconstruction accuracy and computational efficiency for both analytical and iterative reconstruction methods. An insufficient amount of padding fails at removing the bowl artifact, whereas an excessive padding substantially increases the computational cost with insignificant gain in the reconstruction quality.
In this work, analytical reconstructions of FINT data are performed by means of the ramp-filtered adjoint FGP-E. The method, labelled GRID-E, is equivalent to FBP-E. To improve the signal-to-noise ratio (SNR) of the retrieved image, projections are additionally windowed with a Hamming function superimposed on the ramp filter. GRID-E is also used in the first step of ADMP-V.
When performing analytical reconstructions with GRID-E, the edge-padding required for the INT dataset and the oversampling needed by the gridding backprojector coincide. The edge-padding or oversampling factor, i.e. the ratio between the number of pixels of edge-padded and original projections, used for GRID-E is indicated with (the superscript stands for `analytical'). Fig. 9 shows the PSNR and MSSIM scores as a function of for the reconstructions of a SL and MLT sinograms with 1500 views × 512 pixels created from initial phantoms with 4096 × 4096 pixels. The smallest corresponding to the highest values of both PSNR and MSSIM in Fig. 9 is 2.32 (marked with a dashed red line).
An iterative reconstruction method utilizing the FGP-E as forward projector (like ADMP-E) depends on two distinct edge-padding factors: (the superscript stands for `internal'), corresponding to the oversampling factor required by each call of the gridding implementations of and for both standard and interior tomography; (the superscript stands for `external'), which is the edge-padding factor required to address the reconstruction of FINT datasets, can be considered a simple data extrapolation approach and is used for the entire duration of the iterative procedure. A sinogram , reconstructed by ADMP-E, is first padded by a factor , becoming [ = (number of pixels b′)/(number of pixels b)]. Then, every time (analogously it works for ) is invoked inside the CG, is padded again by a factor , becoming . The second padding is only temporary (since it is a requirement of the gridding operators) and, once () has completed its own calculation, the resulting image (sinogram) is cropped to remove the additional -padding. The -padding, instead, remains for the entire run of ADMP-E. The double-padding mechanism required by ADMP-E is shown in Fig. 10.
Fig. 11 shows the PSNR and MSSIM maps as a function of and for the reconstruction of the same SL sinogram as in the experiment of Fig. 9. The maps clearly point out that the reconstruction accuracy strongly varies with and only in a weaker way with . We selected = 1.7 and = 1.87 as optimal edge-padding parameters, because they are the smallest values (highest computational efficiency) guaranteeing simultaneously maximum accuracy in terms of both PSNR and MSSIM.
The optimal edge-padding factors determined here are used for the following reconstructions with ADMP-V and ADMP-E. For ADMP-V, = 2.32 is utilized inside GRID-E and = 1.7 inside ADMP. ADMP-E utilizes = 1.87 as external and = 1.7 as internal edge-padding factors.
4.4. Validation for different zoom-in factors
The reconstruction accuracy of ADMP-E and ADMP-V has been tested for different zoom-in factors (ZIFs), defined as the ratio between the number of pixels of the sinogram in standard tomography and the corresponding FINT one. For example, given a sinogram in standard tomography with 4096 pixels, ZIF = 4 describes a FINT sinogram consisting of the central 1024 pixels of the original dataset.
Fig. 12 shows reconstructions computed by the two iterative methods for increasing ZIFs. The original MLT sinogram has 1500 views × 4096 pixels. At visual inspection the reconstructions are not affected by FINT artefacts and look almost identical, showing that both methods can satisfactorily tackle the reconstruction of FINT data for different ZIFs.
The reconstructions with ZIF = 8 [Figs. 12(c), 12(f)] are overall characterized by a sort of background pattern. For the simulated data used here, this problem, starting to be barely visible with ZIF = 8, can be neglected up to ZIF ≃ 32. This does not represent a problem for the real datasets used in §4.7, as the highest ZIF is roughly 9.
4.5. Validation for different conditions of asymmetry
The accuracy of ADMP-E and ADMP-V has been tested also for the reconstruction of FOVs placed at various distances from the phantom centre. In this way, the iterative methods can be validated for different conditions of asymmetry, i.e. not symmetrical distribution of the attenuation coefficients around the selected FOV.
An example of such a test, performed with the SL phantom, is shown in Fig. 13. The sinograms for FOV-1 and FOV-2 both have 1500 views × 512 pixels. Once again, the reconstructions are not affected by FINT artefacts and look practically identical, proving that both methods can handle FINT cases with pronounced feature asymmetry around the FOV.
4.6. Reconstruction of underconstrained sinograms
Iterative tomographic algorithms are mostly designed to address underconstrained datasets, providing insufficient direct information for accurate reconstruction with analytical methods. An underconstrained sinogram is either undersampled, noisy or a combination of both factors. A sinogram in standard tomography is considered undersampled if m < , where m is the number of views and n is the number of detector pixels (Kak & Slaney, 2001). An FINT sinogram is undersampled when m < with = : represents the number of detector cells required to `cover' the entire object in standard tomography and n the number of available detector cells from the FINT scan (Xiao et al., 2007). Since low-dose fast scans of FINT tomographic microscopy are usually characterized by , undersampling combined with local tomography artifacts result in a sort of `background texture', that can severely affect FINT reconstructions. Iterative reconstruction in standard tomography can greatly reduce the amount of this background texture through the usage of constraints and regularization. Since no constraints are available for the FINT datasets, the removal of the background texture relies entirely on the regularization: for the same amount of noise in the data, the regularization strength τ should be bigger, the higher the undersampling factor, defined as UF = . If τ is too low, the jump associated with the background texture can be regarded by the TV as a collection of faint edges to be preserved and the regularization will only remove the noise component.
ADMP-E and ADMP-V are here tested for the reconstruction of underconstrained datasets against GRID-E. The regularization strength, τ, is separately chosen for ADMP-E and ADMP-V to guarantee a higher quality reconstruction in terms of CNR (while maintaining a similar spatial resolution) compared with GRID-E.
In the first experiment, a FINT sinogram of the SL phantom with ZIF = 4, 200 views × 512 pixels, UF = 94% and corrupted by Poisson noise with = 2.5% of the sinogram mean is considered. Fig. 14 shows that the reconstructions computed by ADMP-E and ADMP-V look more similar to the phantom compared with the GRID-E reconstruction. They have higher MSSIM, PSNR and CNR (increased by a factor of 4.4) scores, as reported in Table 1. The edge profiles (Fig. 15a) indicate that the proposed iterative approaches lead to improved CNR without deterioration of the spatial resolution. Moreover, the peaks corresponding to the different phases of the SL phantom can be clearly identified in the histograms for the iterative reconstructions in Figs. 15(c) and 15(d), whereas only two peaks appear in Fig. 15(b) for the analytical reconstruction.
|
For the second experiment, a FINT sinogram of the MLT phantom with the same condition of undersampling as considered in the previous case (200 views × 512 pixels, ZIF = 4) and a larger amount of Poisson noise (3.5% of the sinogram mean) is used. ADMP-E and ADMP-V substantially reduce the noise in the reconstructions displayed in Figs. 16(c) and 16(d) compared with GRID-E (Fig. 16b) making the light structures clearer against the dark ones. The metric scores in Table 2 are largely higher for the iterative reconstructions, with an improvement of the CNR by a factor of nearly 4.4. The edge profiles (Fig. 17a) for the three different reconstructions coincide almost perfectly at the edge position, demonstrating the efficacy of the split Bregman total variation regularization in removing noise while preserving edges. The higher quality of the iterative reconstructions is also highlighted by the histograms in Figs. 17(c) and 17(d) showing clear peaks corresponding to the two phases of the MLT phantom.
|
For the third experiment, the dataset of MOUSE-1, characterized by ZIF ≃ 3.2, 901 views × 2016 pixels, UF = 91%, a pronounced asymmetry and highly absorbing structures (e.g. ribs) outside the FOV, is considered. Due to the feature size and the pattern complexity, the different quality of the reconstructions displayed in Fig. 18 can be best appreciated in the blow-ups below the images. Features in the iterative reconstructions can be more easily identified thanks to the reduced noise and to the CNR increased by a factor of 5.3 (Table 3). The line profiles at the edge position (Fig. 19a) for the ADMP-E and ADMP-V results match that for the GRID-E reconstruction. Moreover, the double peak in Figs. 19(c) and 19(d) shows that the two main phases of the lung tissue are better separated in terms of grey level in the iterative reconstructions, whereas a single peak characterizes the histogram for the analytical reconstruction in Fig. 19(b).
|
In the fourth experiment, the reconstructed dataset of MOUSE-2 (901 views × 2016 pixels, ZIF ≃ 3.2) has a very similar pattern complexity as MOUSE-1, whereas the morphology, e.g. the shape of the small structures, is quite different. Reconstructions with the ADMP methods in Fig. 20 clearly show higher quality compared with the analytical one, thanks to the decreased noise level and the higher CNR (increased by a factor of 4.6, Table 4). Small dark features are well identifiable in Figs. 20(b) and 20(c), whereas they are mainly covered by noise and undersampling/local tomography artefacts in Fig. 20(a). The edge profiles in Fig. 21 for the iterative reconstructions overlap almost exactly at the edge position with that for the analytical result, indicating also in this case that the total variation regularization operates with negligible loss in terms of spatial resolution. The histograms corresponding to the reconstructions with ADMP-E and ADMP-V in Figs. 21(c) and 21(d) show the presence of two phases. This is not the case in the histogram for GRID-E in Fig. 21(b).
|
For the fifth experiment, the dataset of MOUSE-3, characterized by 500 views × 2016 pixels and ZIF ≃ 9.0, is used. In this case, UF = 98%, but features are on average much larger compared with MOUSE-1 and MOUSE-2, since projections were acquired with a higher magnification. The iterative reconstructions in Fig. 22 offer a clearer vision of the object; nevertheless, all structures recognizable in Figs. 22(b) and 22(c) can be visually identified in Fig. 22(a) as well. The CNR is improved in this example by a factor of 3.7 (Table 5). In Fig. 23(a) an edge can hardly be recognized for the reconstruction with GRID-E, whereas those corresponding to ADMP-E and ADMP-V, practically identical, do show a flank. Histograms in Figs. 23(b) and 23(c) show a distinction between the two main phases of the lung tissue, whereas a single peak dominates in the histogram for the analytical reconstruction in Fig. 23(a).
|
These five experiments with simulated and real sinograms show that: (i) ADMP-E manages to substantially improve the image quality compared with an analytical method like GRID-E when tackling the reconstruction of underconstrained FINT datasets with different ZIF, noise level and feature complexity; (ii) the virtual strategy incorporated in ADMP-V can provide comparable results with those achievable with the ADMP-E.
4.7. Computational cost
To give an idea of the superior computational performance of ADMP-V compared with ADMP-E, the time required for a single iteration, , has been measured for two different datasets on an Intel(R) Core(TM) i7-3520M CPU 2.90 GHz. For a sinogram with 800 views × 504 pixels (convergence reached after eight iterations), (ADMP-E) = 46.7 s and (ADMP-V) = 1.45 s; for a sinogram with 1584 views × 1008 pixels (convergence reached after nine iterations), (ADMP-E) = 173.2 s and (ADMP-V) = 5.7 s. Thus, ADMP-V for small and medium datasets is approximately 30 times faster with respect to ADMP-E.
5. Conclusion
This work introduces two novel strategies for iterative reconstructions of datasets in full interior tomography (FINT), when the FOV is completely inside the object support and no a priori knowledge regarding the support itself and the distribution of the attenuation coefficients in certain ROIs is available. FINT represents a very general case, frequently encountered in many tomographic applications like synchrotron-based X-ray tomographic microscopy.
One strategy works with an edge-padding forward projector. The second is a four-step procedure, requiring the creation of an intermediate virtual sinogram, simulating a full tomography dataset; this sinogram is, then, reconstructed by a standard iterative algorithm, while enforcing a tight constraint outside the reconstruction circle. Both strategies, although not bounded to a specific iterative scheme, have been implemented in this work inside the alternate direction method of multipliers plug-and-play (ADMP), that offers a versatile optimization framework, where different forward projectors and regularizers can be easily plugged in, without altering the structure of the iterative solver. The two resulting iterative methods have been labelled ADMP-E (edge-padding strategy) and ADMP-V (virtual strategy).
The forward gridding projector with minimal oversampling (FGP) is used as standard and edge-padding forward operator for the ADMP. The FGP guarantees fast iterative reconstructions, while keeping the same accuracy of the results achieved with more sophisticated, but much slower, implementations of the Radon transform.
ADMP-E and ADMP-V have been, first, validated for the reconstruction of FINT datasets with different zoom-in factors and asymmetry conditions around the FOV. The methods have, then, been tested on underconstrained simulated and real FINT datasets. Results show that both iterative techniques yield reconstructions of higher quality compared with a standard analytical method: the CNR is greatly improved (on average four times higher), while preserving the spatial resolution, and small features can be more easily identified. The reconstruction quality achieved with the two proposed iterative strategies is comparable. ADMP-V provides, though, superior computational efficiency (about 30 times faster), since it requires a much smaller grid for the computations inside the iterative procedure.
References
Aagesen, L. K., Johnson, A. E., Fife, J. L., Voorhees, P. W., Miksis, M. J., Poulsen, S. O., Lauridsen, E. M., Marone, F. & Stampanoni, M. (2010). Nat. Phys. 6, 796–800. Web of Science CrossRef CAS Google Scholar
Arcadu, F., Nilchian, M., Studer, A., Stampanoni, M. & Marone, F. (2016a). IEEE Trans. Image Process. 25, 1207–1218. Web of Science CrossRef Google Scholar
Arcadu, F., Stampanoni, M. & Marone, F. (2016b). Opt. Express. Submitted. Google Scholar
Berg, S., Ott, H., Klapp, S. A., Schwing, A., Neiteler, R., Brussee, N., Makurat, A., Leu, L., Enzmann, F., Schwarz, J. O., Ersten, M. K., Irvine, S. & Stampanoni, M. (2012). Proc. Natl Acad. Sci. 110, 3755–3759. Web of Science CrossRef Google Scholar
Boyd, S., Parikh, N., Chu, E., Peleato, B. & Eckstein, J. (2010). Mach. Learn. 3, 1–122. CrossRef Google Scholar
Campi, G., Bianconi, A., Poccia, N., Bianconi, G., Barba, L., Arrighetti, G., Innocenti, D., Karpinski, J., Zhigadlo, N. D., Kazakov, S. M., Burghammer, M., Zimmermann, M., v, Sprung, M. & Ricci, A. (2015). Nature (London), 525, 359–362. Web of Science CrossRef CAS Google Scholar
Courdurier, M., Noo, F., Defrise, M. & Kudo, H. (2008). Inverse Probl. 24, 65001. 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
Ebner, M., Marone, F., Stampanoni, M. & Wood, V. (2013). Science, 342, 716–720. Web of Science CrossRef CAS PubMed Google Scholar
Ganne, J., De Andrade, V., Weinberg, R. F., Vidal, O., Dubacq, B., Kagambega, N., Naba, S., Baratoux, L., Jessell, M. & Allibon, J. (2012). Nat. Geosci. 5, 60–65. Web of Science CrossRef CAS Google Scholar
Goldstein, T. & Osher, S. (2009). SIAM J. Imaging Sci. 2, 323–343. Web of Science CrossRef Google Scholar
Hesse, B., Varga, P., Langer, M., Pacureanu, A., Schrof, S., Männicke, N., Suhonen, H., Maurer, P., Cloetens, P., Peyrin, F. & Raum, K. (2015). J. Bone Miner. Res. 30, 346–356. Web of Science CrossRef CAS PubMed Google Scholar
Hestenes, M. R. & Stiefel, E. (1952). J. Res. Natl. Bur. Stan. 49, 409–436. CrossRef Web of Science Google Scholar
Huyn-Thu, Q. & Ghanbaru, M. (2008). Electron. Lett. 44, 800–801. Google Scholar
Kak, A. C. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics. Google Scholar
Kudo, H., Courdurier, M., Noo, F. & Defrise, M. (2008). Phys. Med. Biol. 53, 2207–2231. Web of Science CrossRef PubMed Google Scholar
Kyrieleis, A. V., Titarenko, V., Ibison, M., Connolley, T. & Withers, P. J. (2011). J. Microsc. 241, 69–82. Web of Science CrossRef CAS Google Scholar
Lovric, G., Barré, S. F., Schittny, J. C., Roth-Kleiner, M., Stampanoni, M. & Mokso, R. (2013). J. Appl. Cryst. 46, 856–860. Web of Science CrossRef CAS IUCr Journals Google Scholar
Lu, J., Lee, Y. J., Luo, X., Lau, K. C., Asadi, M., Wang, H. H., Brombosz, S., Wen, J., Zhai, D., Chen, Z., Miller, D. J., Jeong, Y. S., Park, J. B., Fang, Z. Z., Kumar, B., Salehi-Khojin, A., Sun, Y. K., Curtiss, L. A. & Amine, K. (2016). Nature (London), 529, 377–382. Web of Science CrossRef CAS Google Scholar
Marone, F., Münch, B. & Stampanoni, M. (2010). Proc. SPIE, 7804, 780410. CrossRef Google Scholar
Natterer, F. & Wübbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia: Society for Industrial and Applied Mathematics. Google Scholar
Noo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903–3923. Web of Science CrossRef PubMed Google Scholar
Paganin, D., Mayo, S., Gureyev, T., Miller, P. & Wilkins, S. W. (2002). J. Microsc. 206, 33–40. Web of Science CrossRef PubMed CAS Google Scholar
Rashed, E. A. & Kudo, H. (2013). J. Synchrotron Rad. 20, 116–124. Web of Science CrossRef CAS IUCr Journals Google Scholar
Savitzky, A. & Golay, M. J. E. (1964). Anal. Chem. 36, 1627–1639. CrossRef CAS Web of Science Google Scholar
Seger, M. M. (2002). Proceeding of the Symposium on Image Analysis (SSAB02), Lund, Sweden, pp. 33–36. Google Scholar
Shepp, L. & Logan, B. F. (1974). IEEE Trans. Nucl. Sci. 21, 2143. Google Scholar
Toft, P. (1996). PhD thesis, Denmark Technical University, Denmark. Google Scholar
Venkatakrishnan, S. V., Bouman, C. A. & Wohlberg, B. (2013). Proceeding of the 2013 IEEE Global Conference on Signal and Information Processing (Globalsip 2013), Austin, TX, USA, 3–5 December 2013, pp. 945–948. Google Scholar
Walker, S. M., Schwyn, D. A., Mokso, R., Wicklein, M., Müller, T., Doube, M. & Stampanoni, M. (2014). PLoS Biol. 12, e1001823. Web of Science CrossRef Google Scholar
Wang, Z., Bovik, A. C., Sheikh, H. R. & Simoncelli, E. P. (2004). IEEE Trans. Image Process. 13, 600–612. Web of Science CrossRef PubMed Google Scholar
Xiao, X., De Carlo, F. & Stock, S. (2007). Rev. Sci. Instrum. 78, 063705. Web of Science CrossRef Google Scholar
Xu, Q., Mou, X., Wang, G., Sieren, J., Hoffman, E. A. & Yu, H. (2011). IEEE Trans. Med. Imaging, 30, 1116–1128. Web of Science CrossRef Google Scholar
Ye, Y., Yu, H., Wei, Y. & Wang, G. (2007). Intl J. Biomed. Imaging, 2007, 63634. Google Scholar
Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791–2805. Web of Science CrossRef PubMed Google Scholar
Zhang, H., Huang, J., Ma, J., Bian, Z., Feng, Q., Lu, H., Liang, Z. & Chen, W. (2014). IEEE Trans. Biomed. Eng. 61, 2367–2378. Web of Science CrossRef Google Scholar
This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.