On the use of flat-fields for tomographic reconstruction

Jailin, C.; Buffière, J.-Y.; Hild, F.; Poncelet, M.; Roux, S.

doi:10.1107/S1600577516015812

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 24| Part 1| January 2017| Pages 220-231

https://doi.org/10.1107/S1600577516015812

On the use of flat-fields for tomographic reconstruction

Clément Jailin,^a ^* Jean-Yves Buffière,^b François Hild,^a Martin Poncelet ^a and Stéphane Roux ^a

^aLMT (ENS Cachan/CNRS/Université Paris-Saclay), 61 avenue du Président Wilson, F-94235 Cachan, France, and ^bMATEIS (INSA-Lyon/CNRS/Université C. Bernard, Lyon-I), F-69621 Villeurbanne, France
^*Correspondence e-mail: clement.jailin@lmt.ens-cachan.fr

Edited by V. Favre-Nicolin, CEA and Université Joseph Fourier, France (Received 1 June 2016; accepted 7 October 2016)

Seeking for quantitative tomographic images, it is of utmost importance to limit reconstruction artifacts. Detector imperfections, inhomogeneity of the incident beam, as classically observed in synchrotron beamlines, and their variations in time are a major cause of reconstruction bias such as `ring artifacts'. The present study aims at proposing a faithful estimate of the incident beam local intensity for each acquired projection during a scan, without revisiting the process of data acquisition itself. Actual flat-fields (acquired without specimen in the beam) and sinogram borders (when the specimen is present), which are not masked during the scan, are exploited to construct a suited instantaneous detector-wide flat-field. The proposed treatment is fast and simple. Its performance is assessed on a real scan acquired at ESRF ID19 beamline. Different criteria are used including residuals, i.e. difference between projections of reconstruction and actual projections. All confirm the benefit of the proposed procedure.

Keywords: flat-field normalization; artifact corrections; computed tomography.

1. Introduction

Computed tomography (CT) is a three-dimensional imaging technique based on X-ray absorption of materials that provides a three-dimensional volume of a material from a collection of two-dimensional X-ray projections. To deal with the reconstruction procedure (i.e. an inverse Radon transform), many algorithms are used such as filtered back projection (FBP) or algebraic methods (Kak & Slaney, 1988 ).

The quality of the reconstructed volume in CT is affected (Vidal et al., 2005 ; Boas & Fleischmann, 2012 ), for example, by noise, phase contrast, beam-hardening (due to polychromatic X-ray beam), unsteady beam intensity or defective pixels of the detector. Because of the reconstruction technique, a systematic error on a detector site will produce a trace over half a circle or a complete one depending on the range of angles used (usually π-rotation for parallel beam tomography, whereas a full $[2\pi]$ angular range for cone or fan beams). Such marks are called `ring artifacts' (RAs). They are caused by different phenomena originating either from the detector (i.e. scintillator and camera), from spatial inhomogeneities of the beam itself or the optical elements along the line (e.g. monochromator) (Antoine et al., 2002 ; Rack et al., 2011 ). They are always present, although the increasing quality of the detectors, beamlines or sources tend to make them less apparent. However, in the quest for faster tomography or any other challenging specimen image acquisition (Baruchel et al., 2006 ) they tend to reappear. When quantitative image analysis is to be performed (Maire & Withers, 2014 ), such artifacts may introduce a spurious image texture that is difficult to deal with because of extended spatial correlations. Even if not perceived by bare eye, RAs may be revealed in residuals when digital volume correlation is used to follow the kinematics of a specimen in time (Hild et al., 2011 ; Limodin et al., 2011 ). Let us stress that ring artefacts are to be distinguished from `ringing artefacts'; the latter are the manifestation of the Gibbs phenomenon (Hazewinkel, 2001 ) at sharp boundaries, and come from the band-limited nature of the representation of a signal.

Various procedures have been proposed to correct or reduce RAs, which are readily available in most synchrotron facilities or laboratory tomographs (Stock, 2008 ). Corrections are applied to the radiographs before reconstruction (Boin & Haibel, 2006 ; Münch et al., 2009 ; Yousuf & Asaduzzaman, 2009 ; Sadi et al., 2010 ; Rashid et al., 2012 ) where it is assumed that RAs are represented as vertical or horizontal lines in polar coordinates. These lines are then to be filtered out. Corrections are also carried out in a post-processing step after reconstruction (Antoine et al., 2002; Sijbers & Postnov, 2004 ; Axelsson et al., 2006 ). Yet, filtering RAs without degradation of the reconstruction accuracy is still challenging. Another approach aiming at reducing RAs consists of revisiting the experimental acquisition procedure (Davis & Elliott, 1997 ) but it only concerns future acquisitions.

A good understanding of the origin of RAs is believed to offer strategies that should outperform heuristic post-processing techniques. In particular, a detailed study of the spatial and temporal inhomogeneity of the incident beam/detected signal (flat-field) appears promising (Titarenko et al., 2010a ,b ; Baek et al., 2015 ). However, when flat-field images are shot at different instants, they are observed to vary with time (Davis & Elliott, 2006 ). A simple rescaling in time is not sufficient and, as discussed below, even if time and space are naturally well decoupled, different modes are present. This observation was reported in particular by Van Nieuwenhove et al. (2015 ) who propounded an elegant way of extracting those modes from a principal component analysis (Abdi & Williams, 2010 ).

The purpose of the present paper is to propose a treatment of the flat-fields in a similar spirit, yet with some differences. Contrary to standard flat-field correction methods, it is proposed to describe the beam intensity variation in time by a combination of the different acquired flat-fields. The method uses the borders of the projections that are not masked by the specimen during the test to find the correspondence with the flat-field borders so that the beam intensity can be faithfully interpolated in the region masked by the specimen. The procedure requires a less costly determination of the time evolution, which does not rely on a specific assumption of the actual sample projections. Further, it can incorporate additional enrichment whenever necessary without difficulty. After having introduced the notations in §2, a methodology for handling flat-fields and estimating the incoming beam intensity at each instant of time is proposed in §3. The procedure is illustrated in §4 in a case study, namely, a cast-iron specimen that is CT scanned at the European Synchrotron Radiation Facility (ESRF) in Grenoble (France). Comparison with standard handling of flat-fields shows the benefit of the proposed incident beam normalization.

2. Statement of the problem and notations

Reconstruction is based on the relative beam intensity attenuation for each detector position $[{\bf{r}}]$ = (r,z) (where z is parallel to the specimen rotation axis, and r is perpendicular to it) and rotation angle ϕ. The actual collected intensity $[I({\bf{r}},\phi)]$ at time t when the rotation angle $[\phi(t)]$ is to be normalized with the intensity $[I_0({\bf{r}},t)]$ recorded without the sample on the beamline,

$[s({\bf{r}},\phi) = -\log\left[ I({\bf{r}},\phi(t))/I_0({\bf{r}},t) \right]. \eqno(1)]$

Note that, in the sequel, rotation angle ϕ and time t are used equivalently through the correspondence $[\phi(t)]$ as t mostly refers to the sequential ordering of data acquisitions.

The difficulty is that I and I₀ cannot be acquired simultaneously. The assumption that I₀ does not vary with time is only an approximation and hence I₀ fields have to be estimated at each instant t, although they are not available. What is generally accessible is the acquisition of so-called flat-fields, $[f_i({\bf{r}})]$ = $[f({\bf{r}},t_i)]$ , i.e. projections captured when the imaged sample has been moved away from the beamline. This is achieved before t = t₀ or after t = t_N (i.e. the scan is completed). In some cases, the scan is interrupted at specific times, say every 100 acquisitions, to capture intermediate flat-fields, but this interruption requires the sample to be moved away from the beam, which causes some delay and calls for a very accurate repositioning of the sample to be able to carry out the reconstruction. Hence the number of these manipulations is to be kept to a minimum.

These flat-fields f_i provide an estimate of $[I_0({\bf{r}},t)]$ for arbitrary times. Various interpolation procedures are considered using either a piecewise constant or linear temporal variation. The underlying assumption is that the beam intensity is well approximated by such an interpolation scheme. However, in a synchrotron beamline where electrons are to be re-injected into the ring (Willmott, 2011 ) such a procedure may be a crude approximation. Moreover, other artifacts take place, either in the beam shape itself, the X-ray optical elements along the path, on the scintillator or the end imaging device optics/camera.

Very recently, a major step forward has been proposed by Van Nieuwenhove et al. (2015). Based on the acquisition of a long series of flat-fields prior to a scan, the authors proposed to account for the time evolution of these flat-fields, or rather the deviation from their temporal mean $[\overline{f({\bf{r}})}]$ , through a principal component analysis (PCA) such that only a few spatial fields $[\varphi_n({\bf{r}})]$ are to be extracted and it is assumed that those fields provide a good estimate of $[I_0({\bf{r}},t)]$ = $[\overline{f({\bf{r}})}]$ + $[\textstyle\sum w_n(t) \,\varphi_n({\bf{r}})]$ where the weights w_n(t) are unknown. The authors further proposed to estimate those amplitudes in order to minimize the total variation of the attenuation image $[I_0[{\bf{r}},\phi(t)]/I({\bf{r}},t)]$ . Additional steps are proposed in their procedure (e.g. non-local means, intensity rescaling) and the reader is referred to their original publication for details. The philosophy for accounting for the flat-field variation, in particular through PCA, is extremely well suited to the problem, and the discussed examples of Nieuwenhove et al. (2015) showed a drastic reduction in the ring artifacts at the modest expense of the determination of the library of flat `modes' $[\varphi_n({\bf{r}})]$ , and their amplitude w_n(t) in time.

The purpose of the present paper is to propose a treatment of the flat-fields in a similar spirit. However, the resulting procedure will not make a direct usage of PCA, and no specific assumption regarding the texture of the scanned specimen is to be postulated [such as a minimal total variation in $[{\bf{r}}]$ of $[s({\bf{r}},\phi)]$ ].

3. Space-time representation

The exploited data are a series of flat-fields, $[f_j({\bf{r}})]$ , j = $[1,\ldots,N_f]$ , and the entire set of projections $[I({\bf{r}},t)]$ . The complete detector area is called Ξ in the following. When the specimen is in the beamline, more often than not its transverse dimensions are smaller than the detector itself. Hence, left and right edge regions, $[\Omega_{\rm{l}}]$ and $[\Omega_{\rm{r}}]$ , with $[\Omega]$ = $[\Omega_{\rm{l}} \cup \Omega_{\rm{r}}]$ , are available since they are never masked by the sample. These regions may be extended to the top of the surface detector if the specimen height is smaller than the detector (on ex situ tests for example). The restriction of pixels $[{\bf{r}} \in \Xi]$ to these regions is denoted as $[\boldrho \in \Omega]$ . Those edges are very appealing for our purpose as they coincide precisely with $[I_0(\boldrho,t)]$ and will be called control regions. Therefore, on the one hand, a complete coverage in time is available but only over part of the detector (i.e. spatially). On the other hand, the flat-fields $[f_j({\bf{r}})]$ cover the entire image, but only at a few instants of time, t_j, as schematically shown in Fig. 1.

Figure 1
(a) Flat-field captured between t = 300 and t = 301 (the spatial inhomogeneity of the flat-field is quite pronounced with bright spots and an overall intensity gradient). (b) Raw projection at t = 300, where the two regions $[\Omega_{\rm{l}}]$ and $[\Omega_{\rm{r}}]$ are the blue shaded rectangles, and the central darker region is caused by the specimen absorption. (c) In a space-time ( r– t) diagram (sliced at a given height, z = 256 pixels) full flat-fields, f, are shown as red lines to indicate that they cover the entire detector but are captured only at some interruptions of the scan and sinogram borders $[\Omega_{\rm{l}}]$ and $[\Omega_{\rm{r}}]$ shown in blue to highlight that they are available at all times but only over a limited portion of the detector.

It is assumed that the time modulation of the beam affects the flat-fields in proportion to the intensity, so that a correction $[c({\bf{r}},t)]$ of a reference flat-field $[f_{\rm{R}}]$ is sought under a multiplicative form,

$[I_0({\bf{r}},t) = f_{\rm{R}}({\bf{r}})\, c({\bf{r}},t). \eqno(2)]$

In the same spirit as Nieuwenhove et al. (2015), space and time are assumed to be decoupled, so that PCA can provide a faithful representation of the logarithm of the correction through a set of spatial modes $[\Phi_n({\bf{r}})]$ (i.e. independent of time) modulated in time by weights w_n(t),

$[\log[c({\bf{r}},t)] = \textstyle\sum_n w_n(t) \,\Phi_n({\bf{r}}). \eqno(3)]$

The strategy is to exploit $[f(\boldrho,t)]$ and $[I(\boldrho,t)]$ images to estimate $[\Phi_n({\bf{r}})]$ and w_n(t). [Let us note that Nieuwenhove et al. (2015) propose an additive correction to c rather than to log(c) as performed herein.] Likewise, the above additive decomposition is rephrased as a series of multiplicative terms,

$[c({\bf{r}},t) = \textstyle\prod_n \exp\left[w_n(t) \,\Phi_n({\bf{r}})\right]. \eqno(4)]$

3.1. Control regions

Because the sinogram borders Ω are never masked by the specimen during the acquisition, the incident beam intensity $[I_0(\boldrho,t)]$ is equal to the measured one $[I(\boldrho,t)]$ over Ω. At zeroth order, these intensities are close to any flat-field and hence one is chosen as a reference, $[f_{\rm{R}}(x)]$ . This allows us to work only with corrections, c, that are close to unity,

$[I(\boldrho,t) = I_0(\boldrho,t) = f_{\rm{R}}(\boldrho)\,c(\boldrho,t). \eqno(5)]$

It was checked that choosing another flat-field as a reference does not lead to any difference. PCA decomposition (Abdi & Williams, 2010) is classically obtained through the minimization of the quadratic difference between the logarithm of the full data $[I(\boldrho,t)]$ , normalized by the reference flat-field, $[G(\boldrho,t)]$ = $[\log[I(\boldrho,t)]]$ − $[\log[\,f_{\rm{R}}(\boldrho)]]$ , and its space and time mode representation

$[\Gamma^2 = \left\Vert G(\boldrho,t)-\textstyle\sum_n w_n(t) \,\Phi_n(\boldrho) \right\Vert^2_{\Omega}, \eqno(6)]$

where the norm $[\Vert\ldots\Vert^2_{\Omega}]$ implies a summation over the control regions Ω and time. Usually a Euclidian norm is chosen. However, it may not be the most appropriate norm to correctly represent medium frequencies. Because the Euclidian norm can equivalently be computed in Fourier space, all frequencies are given the same weight and the standard white noise will contribute dominantly to the norm and hide medium-frequency features. Filtering of the images with a Gaussian convolution permits low and medium frequencies to be emphasized. The characteristic width of the Gaussian will be of the order of 1 to 2 pixels. This convolution is a direct product in Fourier space, and hence this low-pass filter can also be seen as using a weighted norm in Fourier space. Hence, it is proposed to convolve G with a Gaussian as a pre-processing step and compute its separated space and time modal representation.

Using the conventional normalization $[\textstyle\int w_n(t)^2\,{\rm{d}}t]$ = 1, the minimization of Γ leads to

$[\left\{ \matrix{ w_n(t) = {{1}\over{{\textstyle\int}\Phi_n(\boldrho)^2\,{\rm{d}}\boldrho_{\vphantom{\Big|}}}} \int G(\boldrho,t)\,G(\boldrho,t^{\,\prime})\, w_n(t^{\,\prime})\,{\rm{d}}\boldrho\,{\rm{d}}t^{\,\prime},\hfill \cr \Phi_n(\boldrho) = \int G(\boldrho,t)\, w_n(t)\,{\rm{d}}t.\hfill }\right. \eqno(7)]$

This equation is written using the covariance, $[C_G(t,t^{\,\prime})]$ = $[\textstyle\int G(\boldrho,t)\,G(\boldrho,t^{\,\prime})\,{\rm{d}}\boldrho]$ , and its eigenvalues, $[\lambda_n]$ ,

$[\int C_G(t,t^{\,\prime}) \,w_n(t^{\,\prime})\,{\rm{d}}t^{\,\prime} = \lambda_n w_n(t). \eqno(8)]$

This decomposition is exact and complete when N temporal modes are considered, with N the number of instants t where a frame is available. However, the interest of PCA is to retain only the dominant N_m modes and truncate the remaining ones. The choice of N_m can be performed based on the ranking of the eigenvalues $[\lambda_n]$ , keeping the largest ones, and checking that the residual Γ when the sum is limited to those modes becomes smaller than a threshold. A qualitative but insightful criterion may also be given by spatial modes $[\Phi_{N_m+1}({\bf{r}})]$ , which should look like white noise when N_m is properly chosen.

3.2. Control regions expressed as combinations of (border) flat-fields

The decomposition obtained in the previous section represents the intensity evolution of control regions using the spatial and temporal modes $[\Phi_n(\boldrho)]$ and w_n(t). However, the spatial modes are limited to the Ω region whereas it would be desirable to extrapolate them in the entire domain. The driving idea is now to use the flat-fields in order to perform this extrapolation. For this goal, it is required to `translate' the eigenmodes $[\Phi_n(\boldrho)]$ into their flat-field expressions. Thus, the flat-fields are first clipped to the Ω region, and they constitute the second set of N_f fields $[f_j(\boldrho)]$ . The sought correspondence between both sets,

$[\Phi_n(\boldrho) = \sum_{j\,=\,1}^{N_f} \alpha_{nj}\, F_j(\boldrho), \eqno(9)]$

with $[F_j(\boldrho)]$ = $[\log[\,f_j(\boldrho)]-\log[\,f_{\rm{R}}(\boldrho)]]$ , is computed in the least-squares sense. Thus, the amplitudes are written with the pseudo inverse,

$[\alpha_{n} = \left[(F\otimes F)_{\Omega}\right]^{-1} \left\{(F\cdot\Phi_n)_{\Omega}\right\}, \eqno(10)]$

where

$[\left[(F\otimes F)_{\Omega}\right]_{ij} = \int F_i(\boldrho)\, F_j(\boldrho) \,{\rm{d}}\boldrho \eqno(11)]$

and

$[\left\{ (F\cdot \Phi_n)_{\Omega}\right\}_i = \int F_i(\boldrho)\,\Phi_n(\boldrho)\,{\rm{d}}\boldrho.\eqno(12)]$

Let us note that the method can easily be extended to incorporate any suitable enrichment of the flat-field library $[f_j(\boldrho)]$ by adding extra fields (e.g. gradient, constant term) before the minimization.

3.3. Estimates of I₀ over the entire image

Combining the previous results, the sinogram borders are approximated as

$[\eqalignno{ G(\boldrho,t)&\approx \sum_{n\,=\,1}^{N_m} w_n(t) \sum_{j\,=\,1}^{N_f} \alpha_{nj} F_j(\boldrho) \cr & \equiv \sum_{j\,=\,1}^{N_f} \beta_{j}(t) \,F_j(\boldrho). &(13)}]$

In order to ensure that the sinogram control regions have been accurately described, the residuals

$[\xi_n(\boldrho) = \Phi_n(\boldrho) -\sum_{j\,=\,1}^{N_f} \alpha_{nj}\,F_j(\boldrho) \eqno(14)]$

are analyzed. When specific features are observed in those N_m fields, the residuals indicate that the set of flat-fields is not sufficient for accounting for the relevant modes and may suggest relevant enrichments such as left–right or up–down brightness gradients.

The very same formula is proposed to be extended over the entire detector,

$[G({\bf{r}},t)\approx \sum_{j\,=\,1}^{N_f}\beta_{j}\, F_j({\bf{r}}). \eqno(15)]$

It can be noticed from equation (13) that the control region fields are directly decomposed onto a combination of flat-fields. This decomposition provides the same results. Thus, PCA indicates the `dimensionality' of the problem, i.e. the number of basic fields to be considered. Similarly, the spectrum of eigenvalues of the pseudo-inverse also gives the dimension of the effective space generated by the flat-fields. This shows whether or not one should enrich and, if needed, how enrichment is to be designed. Finally, once the basis is chosen, the β coefficients are computed without reference to the w and α parameters. In a similar fashion to the decomposition of spatial modes, the control region fields are decomposed, in a least-squares sense, onto the same filtered fields (in order to weight the medium frequencies as needed).

The sum over the entire domain of the sinograms $[\textstyle\int s({\bf{r}},t)\,{\rm{d}}{\bf{r}}]$ (for all t) should be constant in time and equal to the total absorption, A, of the specimen. A possible extension of the proposed method is to strictly account for this constraint by normalizing the sinogram

$[\int s({\bf{r}},t)\,{\rm{d}}{\bf{r}} = \int \sum_{j\,=\,1}^{N_f} \gamma_j(t)\,F_j({\bf{r}})\,{\rm{d}}{\bf{r}} - \int G({\bf{r}},t)\,{\rm{d}}{\bf{r}} = A. \eqno(16)]$

This is achieved from the minimization of the following functional to obtain the new γ coefficients (instead of the previous β parameters),

$[\eqalignno{ \Gamma^2 = {}& \left\Vert G(\boldrho,t)-\sum_{j\,=\,1}^{N_f} \gamma_j(t)\,F_j(\boldrho) \right\Vert^2_{\rho} \cr& + \kappa(t) \left(\,\,\sum_{j\,=\,1}^{N_f}\gamma_j(t)\int F_j({\bf{r}})\,{\rm{d}}{\bf{r}} - \int G({\bf{r}},t)\,{\rm{d}}{\bf{r}} - A \right), &(17)}]$

where $[\kappa(t)]$ is a Lagrange multiplier depending on time. Hence, when deriving the above functional, with $[B\equiv \int G({\bf{r}},t)\,{\rm{d}}{\bf{r}}]$ , $[\lbrace F_{\Xi}\rbrace_i]$ = $[\textstyle\int F_i({\bf{r}})\,{\rm{d}}{\bf{r}}]$ and $[\lbrace H_{\Omega}\rbrace_i]$ = $[\textstyle\int F_k(\boldrho)\,G_{ik}(\boldrho,t)\,{\rm{d}}\boldrho]$ ,

$[\left\{ \matrix{ \left\lbrace \gamma(t)\right\rbrace_j = [(F\otimes F)_{\Omega}]^{-1}_{ij} \left(\left\lbrace H_{\Omega}\right\rbrace_i + \kappa(t) \left\lbrace F_{\Xi} \right\rbrace_i \right)_{\vphantom{\big|}},\hfill \cr \left\lbrace\gamma(t)\right\rbrace_j \left\lbrace F_{\Xi} \right\rbrace_j = A+B_\Xi(t),\hfill } \right. \eqno(18)]$

so that

$[\kappa(t) = {{A+B_\Xi(t)-\left\lbrace H_{\Omega}\right\rbrace_i[(F\otimes F)_{\Omega}]^{-1}_{ij} \left\lbrace F_{\Xi} \right\rbrace_j} \over {\left\lbrace F_{\Xi} \right\rbrace_i [(F\otimes F)_{\Omega}]^{-1}_{ij} \left\lbrace F_{\Xi} \right\rbrace_j}}, \eqno(19)]$

that can be injected into the expression of $[\gamma(t)]$ ,

$[\eqalignno{ \left\lbrace \gamma(t)\right\rbrace_j = {}& \left\lbrace \beta(t)\right\rbrace_j &(20)\cr & +[(F\otimes F)_{\Omega}]^{-1}_{ij} {{B_\Xi(t)+A-\left\lbrace \beta(t)\right\rbrace_k\left\lbrace F_{\Xi} \right\rbrace_k }\over{ \left\lbrace F_{\Xi} \right\rbrace_l\, [(F\otimes F)_{\Omega}]^{-1}_{lk} \left\lbrace F_{\Xi} \right\rbrace_k }} \left\lbrace F_{\Xi} \right\rbrace_i.}]$

4. Test case

To illustrate the above treatment, a cast-iron sample scanned on a synchrotron beamline (i.e. ID19 at ESFR) inside an in situ testing machine (Buffière et al., 2010 ) is chosen. The specimen shape is a parallelepiped, 1.6 mm × 0.8 mm in cross-sectional area and 10 mm long. The voxel size is equal to 5.1 µm. The complete scan consists of 600 projections captured at equally spaced angles ranging over half a revolution. Seven flat-fields have been acquired, interrupting the scan every 100 projections and moving the testing machine, and labeled by t = $[[0,100,\ldots,600]]$ as shown in Fig. 1(a). One dark-field is also acquired before the experiment.

4.1. Flat-fields

Fig. 2(a) shows the first flat-field acquired at time t = 0. The beam intensity displays a significant vertical gradient (roughly a factor of two in intensity). A few bright spots are also observed. Fig. 2(b) is the pixel-to-pixel ratio between the first and second flat-fields. First, it appears that the average value is not 1, meaning that the overall beam intensity has changed by about 1%. For applications seeking the sample motion, it is to be observed that an inaccurate correction for the beam intensity could lead to a spurious vertical motion (exploiting the intensity gradient). More generally, any bad account for flat-fields may induce false contrast on radiographs that can be mis-interpreted in the reconstruction or the analysis of the sample kinematics. It is therefore very important to be as accurate as possible. Second, in Fig. 2(b) the bright spots seen in the initial flat-field can still be distinguished. Thus, they are not simply scaled by the overall beam intensity. This observation implies in particular that several modes should be introduced in PCA. Fig. 2(c) shows a dark-field composed of four quadrants (due to the detector technology) and with a low-intensity amplitude.

Figure 2
(a) Flat-field acquired at time t = 0. (b) Pixel-to-pixel ratio of two flat-fields ( f₁/f₃). The mean intensity is observed to vary with time and inhomogeneous patterns (bright spots) are observed and are not well described by the intensity rescaling. (c) Dark-field at t = 0.

4.2. Control regions modes

PCA has been performed based on the borders of the sinogram, where the left and right regions $[\Omega_{\rm l,r}]$ have been stacked together to create the $[I_0(\boldrho,t)]$ set of images. Fig. 3 shows the five main spatial modes $[\Phi_n(\boldrho)]$ describing most of the variations of the beam intensity normalized by the third flat-field (acquired in the middle of the experiment). The two areas $[\Omega_{\rm{l}}]$ and $[\Omega_{\rm{r}}]$ are clearly seen on those modes [as they are in $[I_0(\boldrho,t)]$ ]. The control region areas represent one-seventh of the total surface, Ξ (512 × 512 pixels). As can be seen in Fig. 1(c), the border regions are kept ten pixels away from the projected boundary of the specimen, in order to keep clear from the influence of X-ray diffraction or scattering.

Figure 3
From (a) to (e): spatial modes $[\Phi_n(\boldrho)]$ for n = 1 to 5. The first spatial mode (a) has the highest eigenvalue $[\lambda_1]$ and essentially represents the image intensity, (b) and (c)–(d) correspond to a vertical and horizontal gradient and (e) is a spot mode.

Figs. 4(a)–4(e) display the corresponding temporal amplitudes for the 600 time steps w_n(t). Modes with n > 5 are, spatially and temporally, close to white noise, and are therefore discarded in the following. The former fields can roughly be interpreted in words as a global intensity variation of the X-ray beam (mode 1), a vertical and horizontal intensity gradient (modes 2 and 3–4) and a few spots with different intensity changes (mode 5). The first ten eigenvalues are plotted in Fig. 5.

Figure 4
From (a) to (e): temporal modes w_n(t)/600 for n = 1 to 5 associated with the spatial modes shown in Fig. 3

Figure 5
The ten highest eigenvalues on a log scale. The first five are selected for the proposed procedure.

The number of modes is used to estimate the number of degrees of freedom of the intensity variations. However, consideration of eigenvalues only may not be sufficient for selecting the relevant modes, as some highly concentrated in space may have a low eigenvalue and yet a large effect, while others that are more uniformly distributed may have a large power, and yet a modest detrimental influence. In the present case, even if modes 2 to 5 have small eigenvalues, $[\lambda_n]$ [equation (8)], the temporal variations are significant and the maximum intensity variation could be compared (i.e. the temporal amplitude times the gray level amplitude for each mode).

4.3. Projection of modes onto flat-fields

The above five eigen-modes are now translated into their flat-field expression to allow for a faithful extrapolation to the entire domain. Because modes 2 and 3 present a clear uniform vertical and horizontal gradient, two such simple fields (i.e. homogeneous intensity gradients) are generated and added to the set of flat- and dark-fields as enrichments. Finally, to allow for an overall gray-level offset, a uniform field (valued 1) is also added. This enriched set now consists of ten different fields shown in Fig. 6.

Figure 6
Enriched flat-field library $[F_j(\boldrho)]$ normalized with the third flat-field. The first six fields are $[F_j(\boldrho)]$ , image 7 is the dark-field normalized with the third flat-field and the last three are added fields to help account for horizontal and vertical gradients and uniform intensity offset.

The decomposition onto the ten fields revealed not to be very accurate for high-order modes. To give a larger weight to medium and low spatial frequencies (spots) than on the very high spatial frequency where noise has a higher power, a Gaussian filter with a modest characteristic length of two pixels has been applied to the spatial modes $[\Phi_n(\boldrho)]$ and to the normalized flat-fields $[F_j(\boldrho)]$ . Last, the reconstruction has been performed with the initial noisy images.

Fig. 7 shows the obtained reconstruction of the previous modes (Fig. 3). The analysis of the residuals essentially shows uniform noise without spots, meaning that the modes have been described in a satisfactory manner.

Figure 7
First row: decomposition of the five filtered modes onto enriched flat-fields [equation (13)

]. Second row: residual [equation (14)

Finally, the extrapolation to the entire detector area is now straightforward, keeping the same weights. These modes are displayed in Fig. 8.

Figure 8
The five spatial modes reconstructed over Ξ from an enriched flat-field set [equation (15)

The participation of the dark-field is negligible compared with other flat-fields contributions (two orders of magnitude). This means that the usual dark-field correction can be omitted.

5. Comparison of different reconstructions

It is now possible to compute the ratio of the measured intensity $[I({\bf{r}},t)]$ during the scan with the above derived beam intensity at the corresponding time $[I_0({\bf{r}},t)]$ . From the corrected sinogram, $[s({\bf{r}},\phi)]$ , a tomographic reconstruction is performed to produce the three-dimensional image of the specimen. In order to assess the relevance of the proposed algorithm, different standard inferred choices for $[I_0({\bf{r}},t)]$ have been implemented to normalize $[I({\bf{r}},t)]$ .

5.1. Standard flat-field corrections

Two different standard flat-field normalization methods have been considered. The mostly used approach is to choose the average $[I_0^{\,\rm{Step}}]$ of the flat-fields acquired before and after the experiment (Chen et al., 2012 ; Weitkamp et al., 2011 ) (or the series of scans acquired in a row). In the test case, the six acquisition phases are corrected by the arithmetic average of the two flat-fields acquired just before, k(t), and after, k(t)+1, the considered group of 100 projections,

$[I_0^{\,\rm{Step}}({\bf{r}},t)\equiv {{\left[f_{k(t)}+f_{k(t)+1}\right]}\over{2}}. \eqno(21)]$

Another option consists of choosing a linear interpolation $[I_0^{\,\rm{Lin}}]$ of the flat-fields acquired before and after the experiment (Hammersley, 2001 ). In the test case, the six acquisition phases are corrected by the linear interpolation of the two flat-fields nearest in time, i.e. with a weight p proportional to the time delay via p(t) = [[t-t(k)]/[t(k+1)-t(k)]] ,

$[I_0^{\,\rm{Lin}}({\bf{r}},t) \equiv [1-p(t)]\,f_k+p(t)\,f_{k+1}. \eqno(22)]$

These two choices are challenged with the proposed estimate, $[I_0^{\,\rm{Prop}}]$ . The change of the mean intensity of these different methods is plotted in Fig. 9.

Figure 9
Average intensity for the three compared techniques: (a) piecewise constant average flat-field, (b) piecewise linear, (c) the proposed estimate and (d) proposed correction with a constant absorption constraint.

5.2. Comparison with the proposed procedure

In order to compare the different options for the normalization of the sinograms, three different procedures have been considered:

Procedure 1: comparison of sinograms normalized by the three previous estimates of I₀.

Procedure 2: comparison of reconstructed volumes from the three sinograms. The reconstructions are computed using the ASTRA tomography toolbox (Van Aarle et al., 2015 ) and a SIRT algorithm.

Procedure 3: comparison of projections of the reconstructed volumes with initial corrected radiographs. The projections at 0 and $[\pi/2]$ angles (i.e. nothing but the mere line-sum of all pixels) were first considered because they did not introduce any bias due to the interpretation of partly intercepted pixels by a ray.

5.2.1. Procedure 1: comparison of the sinograms

It is seen that the sinogram normalized with an average flat-field (Fig. 10a) has vertical stripes corresponding to intensity variations. For t = 300 and t = 500, two red arrows point at these global intensity changes that hold for the whole image. These stripes are also visible in Figs. 10(d) and 9, at the same abscissa where there is an important gap between the intensity of the image and the intensity of the proposed normalizing I₀. These variations are also perceived in Fig. 10(e) but their amplitude is much smaller. These global intensity variations are corrected with the first mode of the proposed method.

Figure 10
Step 1. First row: sinogram (x,t) slice at z = 114 pixel corrected with (a) $[I_0^{\,\rm{Step}}]$ , (b) $[I_0^{\,\rm{Lin}}]$ , (c) $[I_0^{\,\rm{Prop}}]$ . Second row: difference between the sinogram slice at z = 114 pixel: (d) proposed-step, (e) proposed-linear. The constant absorption constraint has not been plotted because it cannot be distinguished from (c).

Another feature to be noted is the presence of horizontal lines (blue arrow) in Figs. 10(d) and 10(e). The slice of the sinogram has been chosen to cut few spots. One such intersection is seen in slice z = 114 pixel at r = 300 pixel (Fig. 8). These spots give rise to horizontal lines and are not well corrected with the first two normalization methods $[I_0^{\,\rm{Step}}]$ and $[I_0^{\,\rm{Lin}}]$ .

Constant absorption criterion. As noted earlier, the sum over space of the sinogram $[\textstyle\int s({\bf{r}},\phi)\,{\rm{d}}{\bf{r}}]$ = A should be constant in time and equal to the total absorption, A, of the specimen. The change of this absorption is shown in Fig. 11. It is observed that the proposed reconstruction is much steadier (amplitude of 2%) than the standard correction procedure (amplitude of 7%) especially for high frequencies. The main intensity variation at t = 61 corresponds to a projection angle where the longitudinal axis of the rectangular cross-section sample is perfectly aligned with the beam. When the constant absorption constraint is introduced, the mean intensity is constant (the very small variation may originate from defective pixels and the positivity constraint of the radiographs).

Figure 11
Normalized sum over space of the radiographs. (a) Average correction, (b) linear correction, (c) proposed correction and (d) proposed correction with constant absorption constraint.

5.2.2. Procedure 2: comparison of reconstructed volumes

The SIRT reconstruction from the ASTRA tomography toolbox (Van Aarle et al., 2015) is used to compute the reconstructed volume $[f(\bf x)]$ . The three reconstructions and their differences are displayed in Fig. 12.

Figure 12
Step 2. First row: SIRT reconstruction slice at z = 114 pixel normalized with (a) $[I_0^{\,\rm{Step}}]$ , (b) $[I_0^{\,\rm{Lin}}]$ , (c) $[I_0^{\,\rm{Prop}}]$ . The red rectangle shows the region of interest. Second row: difference between the reconstructed slice at z = 114 pixel, (d) proposed-step, (e) proposed-linear.

In order to show the ring artifacts in the reconstruction, a zoom focusing on $[x\in[200,300]]$ pixels and $[y\in[200,300]]$ pixels is shown in Fig. 13. It is seen that ring artifacts are reduced from (a) to (c) although not completely suppressed with the proposed normalization.

Figure 13
Region of interest on the center of the reconstructed volume (red square in Fig. 12c

) normalized with (a) $[I_0^{\,\rm{Step}}]$ , (b) $[I_0^{\,\rm{Lin}}]$ , (c) $[I_0^{\,\rm{Prop}}]$ .

Entropy of reconstruction. A possible criterion to judge the relevance of the correction procedure is to study the Shannon entropy of the reconstructed volume. The gray level histogram p( f ) (i.e. probability of observing a gray level equal to f) of the reconstructed specimen is computed, and the Shannon entropy (i.e. relative to a uniform measure) S = $[-\textstyle\int p(\,f\,)\log[\,p(\,f\,)]\,{\rm{d}}f]$ is evaluated. Because cast iron consists of few phases, Shannon entropy is expected to be low when the reconstruction is accurate, and hence judging flat-field normalization quality can be carried out aiming at the lowest entropy. In the present case, entropy is based on an 8-bit digitization of the gray levels and is restricted to the specimen volume. The value of S normalized by that of $[S^{\,\rm{Step}}]$ , measured for normalization by the piecewise constant of I₀, is shown in Table 1. The proposed normalization has a markedly smaller value than the two standard procedures. Enforcing a constant absorption is most efficient with respect to the present criterion as it achieves the lowest entropy score.

Table 1
Variation of entropy for different normalizations: (a) piecewise constant, (b) piecewise linear, (c) proposed method, (d) proposed method enforcing constant absorption

Normalization	(a)	(b)	(c)	(d)
Normalized entropy S/S^Step	100%	99.8%	95.7%	92.3%

5.2.3. Procedure 3: comparison of projections of the reconstructed volumes with initial corrected radiographs

The `reconstruction residual', i.e. the difference between the projection of the reconstruction and the original sinogram, has been studied. Fig. 14 shows a cross section of the residual for t = 300 for the three reconstruction methods. The residual is homogeneously distributed over the image and does not display any region with significantly higher residuals.

Figure 14
Residual field of the projected reconstruction and the original sinogram at t = 300: (a) average correction, (b) linear correction, (c) proposed correction.

The norm of these residual fields is plotted for every angle in Fig. 15 for the four methods [i.e. piecewise constant (a), linear (b), proposed (c) and proposed with absorption constraint (d) estimates of the incoming beam intensity I₀]. Procedure (c) provides a good estimate and interestingly it seems to provide a smooth lower envelope of the other residuals (a) and (b). The constraint on the total absorption shows the best results.

Figure 15
Norm of the residual fields for every angle with different correction procedure: (a) average, (b) linear, (c) proposed correction, (d) proposed correction with constant absorption constraint.

Because the reconstructed volume does provide a constant absorption, the computed norm is expected to be strongly correlated with the single angle estimate of the absorption shown in Fig. 11 as observed here. The peaks that clearly appear at t = 61 and for some results at t = 361 correspond to the alignment of the specimen cross section principal axes (i.e. length and width, respectively) with the X-ray beam. Hence they may be attributed to phase contrast effects. When the constant absorption criterion is enforced in the reconstruction, it is observed that the norm of the projection residual reaches even lower values and is more evenly spread over angles.

In this experimental procedure, the in situ testing machine is moved away from the beam in order to acquire the flat-fields. The PMMA tube that constitutes the frame of this machine has a different projection on the flat-fields and on the scans during acquisition. This difference introduces a bias, particularly for standard correction techniques that only use flat-fields as reference. The final results (d) show a low-frequency residual that may originate from the off-centered position of the tube with respect to the rotation axis and an intermittent residual that may be attributed to scratches on the tube of the testing device (Buffière et al., 2010) and give temporally (angular-wise) incoherent residuals.

5.2.4. Perspectives

One challenge of the presented experiment is that the specimen has to be moved away from the beam at several instances in order to acquire the full flat-fields and has to be accurately repositioned (e.g. five times for the reported experiment). These steps may not be extremely precise and may degrade the quality of the reconstruction especially if high spatial resolutions are sought. With the enrichment of the flat-field library thanks to the PCA decomposition, only two flat-fields acquired before and after the experiment can be used and supplemented by additional uniform fields. Such a direction has been explored and it provides reconstructions of high quality. The Shannon entropy (again normalized by the standard piecewise constant procedure) for example gives 144 % if the normalization is performed with only two flat-fields and 94.4 % with the proposed method using the same two flat-fields and three enrichment fields. Thus the intermediate acquisition of flat-fields, at least for acquisitions comparable with the tested example, may be omitted, saving time without prejudice on the quality of reconstruction.

6. Conclusion

A flat-field intensity normalization procedure is proposed, which aims at enhanced tomographic reconstruction quality (through artifact reduction). A suited flat-field is estimated for every acquisition angle (e.g. time) using the borders of the sinogram called control regions. Because these areas are never masked by the scanned specimen, they provide a link between flat-fields acquired at specific scan interruptions and the current projection. Studying the different inhomogeneous beam variations with a space-time decomposition, the principal modes are extracted. The modes are further decomposed onto flat-fields, enriched by some additional modes, in order to extrapolate them over the entire detector area.

The method is applied to a synchrotron tomography scan acquired at ESRF (ID19 beamline). The proposed normalization by the reconstructed flat-field has been compared with alternative standard procedures at different stages considering (i) the sinogram, (ii) the reconstruction or (iii) the projection residual (difference between projection of reconstruction and acquired projection).

Acknowledgements

This work was made possible by an ESRF grant for the experiment MA-501 on beamline ID19. The authors would also like to thank Dr Jan Neggers for helpful comments.

References

Abdi, H. & Williams, L. J. (2010). WIREs Comput Stat, 2, 433–459. CrossRef Google Scholar
Antoine, C., Nygård, P., Gregersen, Ø. W., Holmstad, R., Weitkamp, T. & Rau, C. (2002). Nucl. Instrum. Methods Phys. Res. A, 490, 392–402. Web of Science CrossRef CAS Google Scholar
Axelsson, M., Svensson, S. & Borgefors, G. (2006). Proceedings of the 28th DAGM Symposium on Pattern Recognition, Berlin, Germany, 12–14 September 2006, pp. 61–70. Berlin, Heidelberg: Springer. Google Scholar
Baek, J., De Man, B., Harrison, D. & Pelc, N. J. (2015). Opt. Express, 23, 7514–7526. Web of Science CrossRef PubMed Google Scholar
Baruchel, J., Buffière, J.-Y., Cloetens, P., Di Michiel, M., Ferrié, E., Ludwig, W., Maire, E. & Salvo, L. (2006). Scr. Mater. 55, 41–46. Web of Science CrossRef CAS Google Scholar
Boas, F. E. & Fleischmann, D. (2012). Imaging Med. 4, 229–240. CrossRef Google Scholar
Boin, M. & Haibel, A. (2006). Opt. Express, 14, 12071–12075. Web of Science CrossRef PubMed Google Scholar
Buffière, J.-Y., Maire, E., Adrien, J., Masse, J.-P. & Boller, E. (2010). Exp. Mech. 50, 289–305. Google Scholar
Chen, R.-C., Dreossi, D., Mancini, L., Menk, R., Rigon, L., Xiao, T.-Q. & Longo, R. (2012). J. Synchrotron Rad. 19, 836–845. Web of Science CrossRef IUCr Journals Google Scholar
Davis, G. & Elliott, J. (1997). Nucl. Instrum. Methods Phys. Res. A, 394, 157–162. CrossRef CAS Web of Science Google Scholar
Davis, G. & Elliott, J. (2006). Mater. Sci. Technol. 22, 1011–1018. Web of Science CrossRef CAS Google Scholar
Hammersley, A. (2001). HST: High Speed Tomography Reference Manual V0.4. Technical report. ESRF, Grenoble, France. Google Scholar
Hazewinkel, M. (2001). Gibbs Phenomenon. Berlin: Springer. Google Scholar
Hild, F., Fanget, A., Adrien, J., Maire, E. & Roux, S. (2011). Arch. Mech. 63, 1–20. Google Scholar
Kak, A. C. & Slaney, M. (1988). Principles of Computerized Tomographic Imaging. Piscataway: IEEE Press. Google Scholar
Limodin, N., Réthoré, J., Adrien, J., Buffière, J.-Y., Hild, F. & Roux, S. (2011). Exp. Mech. 51, 959–970. Web of Science CrossRef Google Scholar
Maire, E. & Withers, P. (2014). Int. Mater. Rev. 59, 1–43. Web of Science CrossRef CAS Google Scholar
Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567–8591. Web of Science PubMed Google Scholar
Rack, A., Weitkamp, T., Zanette, I., Morawe, C., Vivo Rommeveaux, A., Tafforeau, P., Cloetens, P., Ziegler, E., Rack, T., Cecilia, A., Vagovič, P., Harmann, E., Dietsch, R. & Riesemeier, H. (2011). Nucl. Instrum. Methods Phys. Res. A, 649, 123–127. Web of Science CrossRef CAS Google Scholar
Rashid, S., Lee, S. Y. & Hasan, M. K. (2012). EURASIP J. Adv. Signal. Process. 2012, 1–18. Web of Science CrossRef Google Scholar
Sadi, F., Lee, S. Y. & Hasan, M. K. (2010). Comput. Biol. Med. 40, 109–118. Web of Science CrossRef PubMed Google Scholar
Sijbers, J. & Postnov, A. (2004). Phys. Med. Biol. 49, N247–N253. Web of Science CrossRef PubMed Google Scholar
Stock, S. R. (2008). Microcomputed Tomography: Methodology and Applications. Boca Raton: CRC Press. Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A. & Withers, P. J. (2010a). J. Synchrotron Rad. 17, 540–549. Web of Science CrossRef IUCr Journals Google Scholar
Titarenko, V., Titarenko, S., Withers, P. J., De Carlo, F. & Xiao, X. (2010b). J. Synchrotron Rad. 17, 689–699. Web of Science CrossRef CAS IUCr Journals Google Scholar
Van Aarle, W., Palenstijn, W. J., De Beenhouwer, J., Altantzis, T., Bals, S., Batenburg, K. J. & Sijbers, J. (2015). Ultramicroscopy, 157, 35–47. Web of Science CrossRef CAS Google Scholar
Van Nieuwenhove, V., De Beenhouwer, J., De Carlo, F., Mancini, L., Marone, F. & Sijbers, J. (2015). Opt. Express, 23, 27975–27989. Web of Science CrossRef CAS PubMed Google Scholar
Vidal, F., Létang, J. M., Peix, G. & Cloetens, P. (2005). Nucl. Instrum. Methods Phys. Res. B, 234, 333–348. Web of Science CrossRef CAS Google Scholar
Weitkamp, T., Haas, D., Wegrzynek, D. & Rack, A. (2011). J. Synchrotron Rad. 18, 617–629. Web of Science CrossRef CAS IUCr Journals Google Scholar
Willmott, P. (2011). An Introduction to Synchrotron Radiation: Techniques and Applications. new York: John Wiley and Sons. Google Scholar
Yousuf, M. & Asaduzzaman, M. (2009). J. Sci. Res. 2, 37–45. CrossRef 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 24| Part 1| January 2017| Pages 220-231

https://doi.org/10.1107/S1600577516015812