Large-area single-photon-counting CdTe detector for synchrotron radiation computed tomography: a dedicated pre-processing procedure

Brombal, L.; Donato, S.; Brun, F.; Delogu, P.; Fanti, V.; Oliva, P.; Rigon, L.; Di Trapani, V.; Longo, R.; Golosio, B.

doi:10.1107/S1600577518006197

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 25| Part 4| July 2018| Pages 1068-1077

https://doi.org/10.1107/S1600577518006197

Large-area single-photon-counting CdTe detector for synchrotron radiation computed tomography: a dedicated pre-processing procedure

Luca Brombal,^a,^b ^* Sandro Donato,^a,^b Francesco Brun,^b Pasquale Delogu,^c,^d Viviana Fanti,^e,^f Piernicola Oliva,^g,^f Luigi Rigon,^a,^b Vittorio Di Trapani,^c,^d Renata Longo ^a,^b and Bruno Golosio ^e,^f

^aDepartment of Physics, University of Trieste, 34127 Trieste, Italy, ^bINFN Division of Trieste, 34127 Trieste, Italy, ^cDepartment of Physical Sciences, Earth and Environment, University of Siena, 53100 Siena, Italy, ^dINFN Division of Pisa, 34127 Pisa, Italy, ^eDepartment of Physics, University of Cagliari, 09042 Monserrato (CA), Italy, ^fINFN Division of Cagliari, 09042 Monserrato (CA), Italy, and ^gDepartment of Chemistry and Pharmacy, University of Sassari, 07100 Sassari, Italy
^*Correspondence e-mail: luca.brombal@ts.infn.it

Edited by P. A. Pianetta, SLAC National Accelerator Laboratory, USA (Received 28 February 2018; accepted 23 April 2018; online 1 June 2018)

Large-area CdTe single-photon-counting detectors are becoming more and more attractive in view of low-dose imaging applications due to their high efficiency, low intrinsic noise and absence of a scintillating screen which affects spatial resolution. At present, however, since the dimensions of a single sensor are small (typically a few cm²), multi-module architectures are needed to obtain a large field of view. This requires coping with inter-module gaps and with close-to-edge pixels, which generally show a non-optimal behavior. Moreover, high-Z detectors often show gain variations in time due to charge trapping: this effect is detrimental especially in computed tomography (CT) applications where a single tomographic image requires hundreds of projections continuously acquired in several seconds. This work has been carried out at the SYRMEP beamline of the Elettra synchrotron radiation facility (Trieste, Italy), in the framework of the SYRMA-3D project, which aims to perform the world's first breast-CT clinical study with synchrotron radiation. An ad hoc data pre-processing procedure has been developed for the PIXIRAD-8 CdTe single-photon-counting detector, comprising an array of eight 30.7 mm × 24.8 mm modules tiling a 246 mm × 25 mm sensitive area, which covers the full synchrotron radiation beam. The procedure consists of five building blocks, namely dynamic flat-fielding, gap seaming, dynamic ring removal, projection despeckling and around-gap equalization. Each block is discussed and compared, when existing, with conventional approaches. The effectiveness of the pre-processing is demonstrated for phase-contrast CT images of a human breast specimen. The dynamic nature of the proposed procedure, which provides corrections dependent upon the projection index, allows the effective removal of time-dependent artifacts, preserving the main image features including phase effects.

Keywords: single-photon counting; CdTe; pre-processing; dynamic correction; breast-CT.

1. Introduction

In recent years high-Z large-area single-photon-counting detectors have become appealing for imaging both in synchrotron and conventional sources experiments (Vedantham et al., 2016 ). These detectors offer remarkable advantages over conventional indirect detection and charge integration systems. Properly operated high-Z single-photon-counting detectors show minimum electronic noise (i.e. noise is Poisson dominated), energy discrimination of photons (i.e. spectral performances) and high detective efficiency (Ballabriga et al., 2016 ; Takahashi & Watanabe, 2001 ). Moreover, unlike scintillator-based detectors where an increase in the efficiency typically leads to a decrease in the spatial resolution due to the scintillating process regardless of the pixel dimension, in direct conversion devices the spatial resolution is mainly limited by the pixel size (Taguchi & Iwanczyk, 2013 ). The aforementioned features make these detectors suitable for low-dose phase-contrast imaging experiments, where both high efficiency for limiting the dose and high spatial resolution to detect phase-effects (e.g. edge enhancement) are needed.

At present, however, the data processing of large-area high-Z single-photon-counting detectors is still challenging. In fact, given the limited area of a single sensor (typically a few cm²) a large field of view is obtained via a multi-module architecture employing arrays or matrices of sensors (Delogu et al., 2017a ; Mozzanica et al., 2016 ). These arrangements lead to the presence of non-negligible gaps between the sensors and require the use of the close-to-edge pixels which often show worse efficiency, stability and gain constancy. Furthermore, these detectors usually suffer from local charge-trapping effects that depend, in general, on the polarization time and on the exposure (Astromskas et al., 2016 ; Park et al., 2014 ; Pennicard & Graafsma, 2011 ; Knoll, 2010 ): these effects cause a gain variation in time which is detrimental especially in CT, where the scan duration may be of the order of several seconds or more (Delogu et al., 2017b ). In the absence of a dedicated pre-processing procedure, all these effects cause artifacts which alter significantly the image quality, possibly impairing its scientific or diagnostic significance. In this work a pre-processing procedure tailored on the characteristics of a novel CdTe single-photon-counting detector (PIXIRAD-8; Bellazzini et al., 2013 ) used for tomographic applications is presented. The implementation and optimization of the pre-processing are performed within the framework of the SYRMA-3D project, whose aim is to perform the world's first in vivo synchrotron radiation breast-CT at the Elettra synchrotron (Trieste, Italy) (Longo et al., 2016 ; Sarno et al., 2016 ; Longo, 2016 ; Delogu et al., 2017a; Brombal et al., 2018 ). The beneficial effects of the pre-processing are shown directly on the reconstruction of a breast specimen which, due to the poor contrast given by adipose and glandular (or tumoral) tissues, represents a challenging sample. In the following, the procedure is described and the effects on the image reconstruction is reported step-by-step.

2. Material and methods

2.1. Beamline, acquisition mode and sample

The images shown in this work were acquired at the SYRMEP (SYnchrotron Radiation for MEdical Physics) beamline at Elettra (Abrami et al., 2005 ). The X-ray beam is produced by one storage ring bending magnet and monochromated by means of a Si(111) double-crystal monochromator allowing to tune the energy in the range 8.5–38 keV, with a resolution of 0.1%. The beam cross section in the patient's room is 220 mm (horizontal) × 3 mm (vertical, Gaussian shape, FWHM). The sample was imaged hanging from the patient support, composed of a rotating table with an ergonomically designed aperture at the rotation center. Thanks to the negligible divergence of the beam within the object (i.e. parallel-beam geometry), the projections were collected only over 180°, thus speeding up the acquisition. Each scan was performed in 40 s in continuous-rotation mode with an angular speed of 4.5° s⁻¹ (Delogu et al., 2017a). The object-to-detector distance of 1.6 m (the maximum available in the present configuration) allows detection of phase effects and, along with the laminar shape of the beam, working in a scatter-free geometry without the need of anti-scattering grids. The sample, fixed with formalin and sealed in a vacuum bag, was imaged at 32 keV delivering a mean glandular dose of 20 mGy. To calculate the dose, the air kerma at the breast position is multiplied by a conversion factor accounting for breast size and glandularity, obtained from an ad hoc developed Monte Carlo simulation based on the GEANT4 code optimized for breast dosimetry (Mettivier et al., 2016 ; Fedon et al., 2015 ). The Directive 2004/23/EC of the European Parliament and of the Council of 31 March 2004 on setting standards of quality and safety for the donation, procurement, testing, processing, preservation, storage and distribution of human tissues were followed. The images reported in this study were acquired to guide the pathologist in the localization of lesions for the histological examination, according to the standard procedures of the Pathology Unit of the Academic Hospital of Cattinara, Trieste University, accredited by JCI (Joint Commission International). The sample was derived from surgical material sent to the Pathology Unit according to local guidelines for histological examination.

2.2. Detector

The detector used in this project is a large-area high-efficiency direct-conversion CdTe photon-counting device (PIXIRAD-8) (Bellazzini et al., 2013). It is made up of eight modules and the pixels are arranged on a honeycomb matrix with 60 µm pitch. Each block has a hybrid architecture in which the 650 µm-thick CdTe sensor and the readout electronics are coupled by means of the flip-chip bump-bonding technique. The active area of each block is 30.7 mm × 24.8 mm, leading to a global active area of 246 mm × 25 mm, corresponding to 4096 × 476 pixels with a gap of 3 pixels between adjacent modules. Each pixel is associated with two independent 15-bit counters which can be used either in color mode (two different energy thresholds, useful for chromatic imaging) or in dead-time-free mode. When the latter is selected, the detector fills one counter while reading the other, thus providing a zero-dead-time acquisition. The detector shows a linear response up to 6.5 × 10⁷ counts mm⁻² s⁻¹ (30 keV photons, 5 keV threshold) while its absorption efficiency is almost 100% for photon energy below 50 keV (Bellazzini et al., 2013; Delogu et al., 2016 ; Vincenzi et al., 2015 ). The maximum photon flux used in this work corresponds to 1.5 × 10⁷ counts mm⁻² s⁻¹, while the threshold was fixed at 3 keV, which is sufficient to suppress the electronic noise. The projection images were collected at about 30 Hz (the maximum frame rate available) which, given the selected rotation speed, corresponds to 1200 projections per scan. Flat projections (i.e. without the sample) used for the flat-fielding procedure were obtained prior to the scan with the same acquisition parameters. All the images are first streamed to the control PC via a Gigabit Ethernet connection and then they undergo the pre-processing procedure described in the next section.

2.3. Pre-processing

The pre-processing procedure has a modular structure comprising five steps, specifically dynamic flat-fielding, gap seaming, dynamic ring removal, projection despeckling and around-gap equalization. All the steps will be described in the following sections. For the sake of portability and computational efficiency, the code is implemented the language C and is available upon request to the corresponding author. The complete pre-processing of 1200 16-bit raw projections, with dimensions of 2300 × 70 pixels each, requires about 4 min on an 8-core Intel Core i7-6700 CPU @ 3.40 GHz including loading and saving of the dataset.

2.3.1. Dynamic flat-fielding

The flat-fielding procedure is common to most of the X-ray imaging applications and it serves multiple purposes, namely to correct the beam shape and intensity inhomogeneity, to equalize the inhomogeneity in the pixel gain and to perform the normalization preparing the planar images to the CT reconstruction. The standard flat-fielding consists of a pixel-by-pixel division of each projection image with a constant flat image (i.e. acquired without the sample). Defining $[P(x, y\semi t)]$ as the projection image, with x,y the pixel coordinates and t the projection index proportional to the acquisition time, and $[\bar F_0(x,y)]$ as the constant flat image, the corrected image will be

$[f_{\rm{static}}(x,y\semi t) = {{P(x, y\semi t)} \over {\bar F_0(x,y)}}. \eqno(1)]$

Given a fixed detector's frame rate, the statistics of $[\bar F_0(x,y)]$ are increased computing the average of (2w+1) flat images, where w determines the width of the window,

$[\bar F_0(x,y) = {{1} \over {2w+1}}\sum\limits_{t\,=\,1}^{2w+1}F(x, y\semi t). \eqno(2)]$

The choice of an odd number as the window width has been made for the sake of notation coherence: in the following most of the presented filter windows are centered in a pixel of interest so that an odd filter dimension is required. With this procedure, hereinafter referred to as static flat-fielding, the presence of a detector's gain time-dependence in the projection images cannot be compensated since the flat image is not time-dependent. On the contrary, the dynamic flat-fielding approach requires as many flat-field images as the number of projections so that the denominator of equation (1) can be substituted with a moving average of 2w+1 flat images,

$[\bar F(x, y\semi t) = {{1} \over {2w+1}} \sum\limits_{t'\,=\,t-w}^{t+w}F(x, y\semi t'). \eqno(3)]$

In this way, assuming that the gain time-dependence is reproducible, each flat image has both high statistics and the same time-dependence as the projection images. The flat-fielded projections will be

$[f_{\rm{dynamic}}(x, y\semi t) = {{P(x, y\semi t)} \over {\bar F(x, y\semi t)}}. \eqno(4)]$

In order for this approach to be used, a slow time-dependence of gain is assumed so that, within the moving average window 2w+1, the flat images are considered to be constant. Namely, given a 30 Hz frame rate and a window of 2w+1 = 11 frames, the gain should not vary significantly for times of the order of 1 s. In addition, the fluctuations of the beam are assumed to be small in the time scale of the acquisition: this requirement is generally fulfilled at the Elettra synchrotron operated in top-up mode, where 1 mA of ring current is injected every 20 min, having a baseline of 140 mA at 2.4 GeV. A different approach to the dynamic flat-fielding based on principal component analysis is given by Van Nieuwenhove et al. (2015 ).

2.3.2. Gap seaming

The PIXIRAD-8 detector, as for most of the multi-module single-photon-counting devices, has a small gap (3 pixels wide) between each pair of modules which needs to be filled within the pre-processing procedure. The selected approach is a linear interpolation with a rectangular 9 × 8 pixels kernel. For each pixel within the gap, the interpolation window is chosen to be half in the left module and half in the right one (regions A and B in Fig. 1), then the mean value of each half is computed and the gap-pixel value is defined as the weighted average of the two mean values,

$[\eqalignno{ f_{\rm{gap}}(x,y\semi t) = {}& {{u(x)}\over{N\!_A}} \, \sum\limits_{(x'\!,\,y')\,\in\,A_{\vphantom{\big|}}} f\left(x',y'\semi t\right) \cr& + {{v(x)}\over{N\!_B}} \,\sum\limits_{(x'\!,\,y')\,\in\,B} f\left(x',y'\semi t\right), & (5) }]$

where N_A = N_B is the normalization factor while the weights u(x) and v(x) are the normalized distances between the pixel within the gap and the regions A and B. Despite its simplicity, this procedure represents a good compromise between image quality and computational load. Nevertheless, more sophisticated approaches, such as the inpainting technique described by Brun et al. (2017 ), may be considered if wider gaps or high-contrast dishomogeneities in the sample are present.

Figure 1
Illustration of the gaps seaming procedure. The gray region represents the gap while the rectangle is the interpolation window used for the pixel of interest (light blue). The figure is not to scale.

2.3.3. Dynamic ring removal

Ring artifacts, produced by gain inhomogeneities at the pixel level, are commonly encountered in tomographic reconstruction. In most of the cases the pixel (or group of pixels) producing the ring has a constant gain offset with respect to its neighbors, so that a single equalization is sufficient to remove or at least mitigate the artifact. In this case, despite the application of the dynamic flat-fielding, some pixels still show a time-dependent gain, resulting in rings with a non-constant intensity. To compensate for this artifact a dynamic (i.e. depending on the projection index) equalization factor has to be used. The implemented ring-removal algorithm makes use of the alpha-trimmed filter, which is a hybrid of the mean and median filters. For each pixel, this filter takes a window of nearest neighbors, sorts their values, excludes the largest and the smallest values and replaces the pixel with the average of the remaining ones. Let g(i) be a one-dimensional image, h and c two integers that represent, respectively, the filter window and the confidence window half-widths, with c < h. The alpha-trimmed filter algorithm can be described as follows.

(i) For each pixel i, consider the window of its 2h+1 neighbors,

$[w(\,j)=g(i+j),\qquad -h \leq j \leq h. \eqno(6)]$

(ii) Sort the values of w in ascending order,

$[w_{\rm{s}}={\rm{sort}}(w). \eqno(7)]$

(iii) Substitute the pixel i with the average of w_s within the confidence window of size 2c+1,

$[\bar{g}_{\rm{s}}(i)={{1}\over{2c+1}} \,\sum\limits_{j\,=\,-c}^c \, w_{\rm{s}}(\,j). \eqno(8)]$

Basically, in this average we are excluding the h-c smallest values and the h-c largest values. Note that if c = 0 the alpha-trimmed filter reduces to the median filter, while if c = h it reduces to the mean filter. In a two-dimensional or three-dimensional image, the alpha-trimmed filter can be applied along each dimension: we will call S_x[g], S_y[g] and S_t[g] the images filtered along the dimensions x, y and t, respectively. Furthermore, we define the filter applied along two or three dimensions as the composition of two or three one-dimensional alpha-trimmed filters, as for instance S_xy[g] = S_x[S_y[g]].

Given $[f(x,y\semi t)]$ , the three-dimensional image describing the whole set of projections, and $[G^\sigma_t[\,f]]$ , the convolution of the image f with a Gaussian function of standard deviation σ along the projection axis t, the ring-removal algorithm consists of the following steps:

(i) First apply the alpha-trimmed filter to the projections along the dimension t, then filter them with a Gaussian convolution along the same dimension,

$[f_1(x,y\semi t) = G^{\sigma}_t[S_t[\,f]](x,y\semi t), \eqno(9)]$

where σ should be a significant fraction of the number of projections.

(ii) Apply the alpha-trimmed filter to f₁ along the dimensions x and y,

$[f_2(x,y\semi t) = S_{xy}[\,f_1](x,y\semi t). \eqno(10)]$

(iii) f₁ is smooth along the dimension t by construction. It is also expected to be a smooth function along the dimensions x and y, therefore f₂ and f₁ should be close to each other, unless there is an equalization problem. Evaluate the equalization correction factor as

$[\alpha(x,y\semi t) = f_2(x,y\semi z)/f_1(x,y\semi t). \eqno(11)]$

(iv) Apply the correction factor to obtain the ring-corrected image,

$[f_{\rm{rc}}(x,y\semi t) = \alpha(x,y\semi t)\,f(x,y\semi t). \eqno(12)]$

In our implementation, we are using h_x = h_y = h_z = 10, c = h/2 for all dimensions and σ = N_p /10. Here we remark that the main advantage of this algorithm is that the equalization factor α varies with the projection index, allowing to cope with non-constant ring artifacts. The results of this approach will be compared with two of the most known filters which tackle the ring-removal problem from different perspectives, namely the one proposed by Rivers (Rivers, 1998 ; Boin & Haibel, 2006 ), based on a moving average filtering, and the one proposed by Münch, based on a combined wavelet-Fourier filtering (Münch et al., 2009 ).

2.3.4. Projection despeckling

In each projection image few (about 0.5%) pixels with an abnormal number of counts, either lower or higher than the neighboring pixels, are present. Their appearance is not reproducible neither in space nor in time and their content cannot be correlated with the actual number of impinging photons. To remove these speckles, which cause streaking artifacts in the reconstructed image, they first need to be recognized and then replaced. The procedure is based on a slightly different version of the alpha-trimmed filter described in the previous section, modified in order to filter only the bad pixels: for each pixel position i the average $[\bar{f}(i)]$ and standard deviation $[\sigma(i)]$ of the pixels within a confidence window are computed, then the pixel of interest is replaced only if its value differs from the mean value more than $[N\sigma(i)]$ , N being a parameter of the filter. In this way N acts as a filter sensitivity threshold, where if $[N\to 0]$ all the pixels will be filtered, as in the implementation reported in §2.3.3, while if $[N\to\infty]$ no pixels will be modified. Moreover, when calculating the average and standard deviation we are excluding the h-c smallest values and the h-c largest values, meaning that pixels with either abnormally high or low counts can be easily discarded. For the projection despeckling, the filter window is a 5×5 pixels square and the confidence window is a 3×3 pixels square, while the optimization of the parameter N is reported in the results section.

2.3.5. Around-gap equalization

The last step of the pre-processing is a dedicated procedure for equalizing the pixels around the gaps between modules. This further equalization is required since several adjacent columns of close-to-edge pixels show a non-optimal gain behavior. This effect involves a large number of pixel columns (30–40 columns across the gap), hence the action of the ring-removal filter, which operates with a 10 pixel window, is not sufficient. This procedure is based on a moving average along the projection axis and it is described as follows.

(i) Given a projection t, a volume C of width 2c = 40 pixels, height equal to the full height of the projection and depth N_p/3, where N_p is the number of projections, is selected across the gap between two modules. Two other volumes (A and B) with the same height, depth and a width of 2a = 10 pixels are selected adjacent to C (see Fig. 2).

Figure 2
Illustration of the equalization procedure: pixels of the projection t within the volume C are those to be equalized. See text for a complete description. The figure is not to scale.

(ii) The mean value along the x and t axes is computed for the volumes A and B,

$[\eqalign{ \bar{f}_A(y\semi t) & = {{1}\over{2aN_{\rm{p}}/3}} \,\,\sum\limits_{x\,=\,x_A-a}^{x_A+a} \,\, \sum\limits_{t'\,=\,t-N_{\rm{p}}/6}^{t+N_{\rm{p}}/6} f(x,y\semi t'), \cr \bar{f}_B(y\semi t) & = {{1}\over{2bN_{\rm{p}}/3}} \,\,\sum\limits_{x\,=\,x_B-b}^{x_B+b} \,\, \sum\limits_{t'\,=\,t-N_{\rm{p}}/6}^{t+N_{\rm{p}}/6} f(x,y\semi t'). } \eqno(13)]$

(iii) The mean value along t is computed for the volume C,

$[\bar{f}_C(x,y\semi t) = {{1}\over{N_{\rm{p}}/3}} \,\, \sum\limits_{t'\,=\,t-N_{\rm{p}}/6}^{t+N_{\rm{p}}/6} f(x,y\semi t'). \eqno(14)]$

(iv) The equalization factor is computed as

$[{\rm{Eq}}(x,y\semi t) = {{u(x)\,\bar{f}_A(y\semi t) + v(x)\,\bar{f}_B(y\semi t)} \over {\bar{f}_C(x,y\semi t)}}, \eqno(15)]$

where u(x) and v(x) are defined as in §2.3.2.

(v) The image is multiplied for the equalization factor,

$[f_{\rm{around}}(x,y\semi t) = f(x,y\semi t)\,{\rm{Eq}}(x,y\semi t). \eqno(16)]$

In order for this procedure to be effectively used, the pixels within the regions A and B must not show a non-optimal behavior. Moreover, as mentioned for the dynamic flat-fielding and ring-removal steps, the around-gap fixing equalization factor depends on the projection index, thus allowing to compensate for slow gain variations of close-to-gap pixels.

2.4. Image reconstruction

Prior to the tomographic reconstruction, the pre-processed projections can be optionally phase-retrieved. The benefit of a proper application of the phase-retrieval algorithm, based on the homogeneous transport of intensity equation approach (Paganin et al., 2002 ), is to increase the contrast-to-noise ratio while preserving spatial resolution, thus enhancing the visibility of low-contrast structures (Beltran et al., 2011 ). The projections, either with or without the phase retrieval, are reconstructed via a GPU-based filter back-projection (FBP) with a standard Shepp–Logan filtering. The reconstruction process, comprising the phase retrieval, requires less than 1 min using a Nvidia GeForce GTX 1080.

3. Results and discussion

In order to compare the flat-fielding procedures in the projection space, two sets of 1300 flat projections were acquired with different photon fluences: one is collected with a low photon fluence to simulate the sample's absorption, the other, acquired with four times higher statistics, is used for the flat-fielding. This choice is made to uncouple the effects of time and exposure on the detector's gain, thus having two datasets with the same acquisition time (i.e. acquired after the same time from the polarization of the CdTe sensor) but different exposures. In Figs. 3(a) and 3(b) details of the first projection normalized with the static and the dynamic flat-field approach are reported: at the center of both images a cluster of pixels with a gain lower than the neighboring ones is present. Observing the same region at a later time, it is evident that the cluster exhibits a gain variation which is more pronounced for the static flat-fielding, in Fig. 3(c), with respect to the dynamic flat-fielding, in Fig. 3(d). Focusing on the intensity plots as a function of time in Figs. 3(e) and 3(f) of a group of pixels within the cluster, it is clear that the gain variation of the statically flat-fielded (∼55%) dataset is significantly higher with respect to that (∼20%) of the dynamically flat-fielded projections. Moreover, as should be expected, the latter shows a smoother time-dependence which can be better compensated by the ring-removal procedure. The effects of each uncompensated crystal defect can be traced through the tomographic reconstruction process. In Fig. 4(a) a detail of the reconstructed image corresponding to a row through the defective pixel cluster obtained with the static flat-fielding is shown: a bright streak-like artifact embedded within a partial ring artifact, due to the uncompensated gain variation, is observed. Fig. 4(b) reports the same detail when the dynamic flat-field approach is used: in this case the streak is barely visible while the ring has been removed. In both images the whole pre-processing procedure has been applied in order to show only the effect of the flat-fielding in the final reconstruction.

Figure 3
Comparison between static and dynamic flat-fielding procedures in the projection space using two flat datasets with different statistics: in (a) and (c) the first and last projections when the static flat-field is applied, in (b) and (d) the first and last projections when the dynamic flat-field is applied. In (e) and (f) the average intensity of the bad pixel cluster as a function of time for the static and dynamic flat-field, respectively, is shown.

Figure 4
Details of a reconstruction obtained applying the static (a) and the dynamic (b) flat-fielding. The arrow indicates a streak artifact clearly visible in (a) while it is barely visible in (b).

Figs. 5(a) and 5(c) show, respectively, the sinogram and the tomographic reconstruction of the sample where only the flat-fielding has been applied. The sample was imaged using four modules of the detector, thus in the sinogram only three gaps are visible, producing marked ring artifacts in the reconstruction. The artifacts cover only half of the circumference because the projections are acquired over 180°. In Figs. 5(b) and 5(d) both the sinogram and the reconstruction are reported after the gap seaming: given the small size of the gaps (3 pixels wide) the interpolation does not introduce significant artifacts, thus preserving the anatomical information. Nevertheless, the resulting image is still affected from the presence of several artifacts which need to be corrected.

Figure 5
Sinograms and reconstructions obtained before (a, c) and after (b, d) the gap seaming.

Figs. 6(a) and 6(d) show the sinogram and the reconstruction where the Rivers ring-removal filter (Rivers, 1998; Boin & Haibel, 2006) has been applied with a window width of 11 pixels, while in Figs. 6(b) and 6(e) the Münch filter (Münch et al., 2009) has been applied with a decomposition level of 5 and a Gaussian bandpass function width of 3. From the sinograms it can be seen that neither the Rivers nor the Münch filter are optimal: in both cases most of the rings are only partially compensated resulting in arc (i.e. partial ring) artifacts. In particular, focusing on the Rivers approach where a constant equalization factor is used, the artifacts appear to be brighter at the top of the sinogram, well corrected in the central part and darker at the bottom. Again, this is due to the time gain variation which occurs to some pixels as previously described, e.g. the plot of Fig. 3(f). The Münch filter yields slightly better results on the rings but it introduces a low spatial frequency modulation strongly affecting the image quality. Comparing these results with Figs. 6(c) and 6(f), obtained using the procedure described in §2.3.4, it is clear that the latter yields the best results, substantially removing most of the ring artifacts. It is worth noting that the main advantage of this approach is the presence of an equalization factor varying with the projection index.

Figure 6
Sinograms and reconstructions obtained by applying the Rivers (a, d), Münch (b, e) and the dynamic (c, f) ring-removal filters. Sinograms are inverted and displayed on a logarithmic scale for better visualizing the action of the filters. The arrows in both the sinograms and the reconstructions indicate uncompensated ring artifacts.

As reported in §2.3.4, the parameter N of the despeckling filter should be optimized in order to remove only the bad pixels. For this purpose a dataset of 1300 flat projections has been acquired and subdivided into two datasets consisting of the even and the odd projection, respectively. Then the even projections were divided, pixel by pixel, by the odd projections. Thus, the gain dependence from time and exposure is matched and the distribution of the bad pixels alone can be studied. The gray-level histogram of the dataset obtained with this approach is plotted in Fig. 7(a) (black dashed line): if no bad pixels are present, the distribution should be a Gaussian centered around 1, whose width is only dependent on the photon statistics. On the contrary, the presence of bad pixels widens the distribution on both sides. The despeckling filter is expected to suppress the tails of the distribution without affecting the width of the Gaussian, i.e. the statistical noise. By varying continuously the filter parameter N, it is found that values of around 15 satisfy this request (blue solid line) while, for lower N (e.g. N = 3, red dashed line), the statistical noise is reduced, and thus the image is smoothed. The same overcorrection effect is observed when applying a common despeckling filter, such as the median filter, as reported in Fig. 7(b).

Figure 7
Histograms for the despeckling filter optimization. In (a) the non-filtered spectrum (black dashed line) is compared with the filtered ones (blue solid line for N = 15, red dashed line for N = 3); in (b) the median filtered spectrum (green dashed line) is also reported.

Once the parameter N has been optimized, the despeckling filter can be applied to the projections. Figs. 8(a) and 8(b) show details of the sinogram before and after the application of the filter respectively: the bad pixels have been removed without affecting the image noise and texture. The effect of the filter on the reconstruction is reported in Figs. 8(c) and 8(d), where in the unfiltered image several striking artifacts due to bad pixels are visible. Here, it has to be remarked that the optimization of the parameter N is crucial since an excessive smoothing of the projections may disrupt the edge-enhancement effect, which is one of the key features of the synchrotron radiation breast CT.

Figure 8
Sinograms and reconstructions before (a, c) and after (b, d) the application of the despeckling filter. The arrows indicate some of the speckles in the sinogram and some of the streaks in the reconstruction.

The last step of the pre-processing procedure is the around-gap equalization. In fact, referring to Fig. 9(a), two ring artifacts corresponding to the regions around the gaps between modules can still be observed. Once the equalization procedure is applied, the rings are removed and the final reconstructed image, reported in Fig. 9(b), is free from major artifacts.

Figure 9
Reconstructions before (a) and after (b) the around-gap equalization.

After the projections have been pre-processed, the phase-retrieval algorithm is applied [two materials approach, $[\delta/\beta]$ = 869 (Burvall et al., 2011 ; Brombal et al., 2018)]. Noticeably, the phase-retrieval algorithms produce a remarkable increase in the contrast-to-noise ratio, thus highlighting also uncompensated artifacts which may be barely visible in the phase-contrast images. In Figs. 10(a) and 10(b) details of the reconstruction processed only with the first two steps of the pre-processing procedure (namely flat-fielding and gap seaming) are reported, with and without phase retrieval: in both cases severe ring artifacts are observed but, when phase retrieval is applied, streaking artifacts arising from uncompensated speckles become evident, definitely impairing the image quality. Conversely, when the whole pre-processing is applied, both the absorption and phase-retrieved images in Figs. 10(c) and 10(d) do not report significant artifacts. In this context, it should be stressed that the optimization of the pre-processing procedure must account also for the subsequent image processing (e.g. phase retrieval) in order to yield a high image quality. In Fig. 11 the final result of the data processing, comprising the pre-processing and the phase-retrieval procedure, is shown: the extension, shape and boundaries of both the tumoral and glandular tissue (light gray) are clearly visible without artifacts.

Figure 10
Detail of a reconstruction without (a, c) and with (b, d) the phase retrieval. In (a) and (b) only the flat-fielding and gap seaming steps are applied; in (c) and (d) the whole pre-processing procedure is used.

Figure 11
Final reconstruction obtained subsequently applying the pre-processing procedure and the phase-retrieval.

4. Conclusions

In this work the effectiveness of the pre-processing procedure tailored for a multi-module single-photon-counting CdTe detector (PIXIRAD-8) was discussed step-by-step and demonstrated in the framework of a synchrotron radiation breast-CT experiment. The challenges concerning the presence of a dead space between adjacent modules, the time-dependent gain variation due to charge trapping and crystal impurities and the non-optimal behavior of close-to-edge pixels have been specifically addressed. In particular, the dynamic flat-fielding and ring-removal procedures yielded better results if compared with standard techniques. In fact, the main advantage of the implemented algorithms is their ability to cope with pixels showing a time gain variation during the tomographic scan, which is a feature common to most high-Z single-photon-counting detectors while it is rarely encountered in conventional CT applications. Moreover, great care is taken in the filters optimization in order to preserve phase effects, which are of paramount importance in a synchrotron-radiation-based experiment. The effects of the pre-processing on the image quality were demonstrated using a breast specimen, which represents a challenging sample to be imaged due to the poor contrast between different tissues. In addition to the presented images, the pre-processing procedure has shown its effectiveness for a number of different test objects and biological samples, within a wide range of beam energies, photon fluences and detector thresholds (Contillo et al., 2018 ; Brombal et al., 2018). Moreover, an effective artifact suppression has been obtained for both phase-contrast and phase-retrieved images. For the sake of portability, the whole code of the pre-processing procedure is implemented in the language C in order to be compiled and executed within all the operating systems and it is made available upon request to the corresponding author. It is the authors' belief that high-Z single-photon-counting detectors will be widely used in future CT applications, especially in medical imaging, due to their high-efficiency, low noise and spectral performances: in this context, the pre-processing procedure presented in this work may represent a useful approach to be used in other applications.

Acknowledgements

The authors thank Professor F. Zanconati MD and Dr Deborah Bonazza MD for providing and preparing the breast specimen.

Funding information

The SYRMA-3D project is supported by Istituto Nazionale di Fisica Nucleare (National Scientific Committee 5 for Technological and inter-disciplinary research) and Elettra-Sincrotrone Trieste SCpA. SD is partially supported by Consorzio per la Fisica Trieste

References

Abrami, A., Arfelli, F., Barroso, R., Bergamaschi, A., Bille, F., Bregant, P., Brizzi, F., Casarin, K., Castelli, E., Chenda, V., Dalla Palma, L., Dreossi, D., Fava, C., Longo, R., Mancini, L., Menk, R.-H., Montanari, F., Olivo, A., Pani, S., Pillon, A., Quai, E., Ren Kaiser, S., Rigon, L., Rokvic, T., Tonutti, M., Tromba, G., Vascotto, A., Venanzi, C., Zanconati, F., Zanetti, A. & Zanini, F. (2005). Nucl. Instrum. Methods Phys. Res. A, 548(1–2), 221–227. Web of Science CrossRef Google Scholar
Astromskas, V., Gimenez, E. N., Lohstroh, A. & Tartoni, N. (2016). IEEE Trans. Nucl. Sci. 63, 252–258. Web of Science CrossRef Google Scholar
Ballabriga, R., Alozy, J., Campbell, M., Frojdh, E., Heijne, E., Koenig, T., Llopart, X., Marchal, J., Pennicard, D., Poikela, T., Tlustos, L., Valerio, P., Wong, W. & Zuber, M. (2016). J. Instrum. 11, P01007. Google Scholar
Bellazzini, R., Spandre, G., Brez, A., Minuti, M., Pinchera, M. & Mozzo, P. (2013). J. Instrum. 8, C02028. Web of Science CrossRef Google Scholar
Beltran, M., Paganin, D., Siu, K., Fouras, A., Hooper, S., Reser, D. & Kitchen, M. (2011). Phys. Med. Biol. 56, 7353–7369. Web of Science CrossRef Google Scholar
Boin, M. & Haibel, A. (2006). Opt. Express, 14, 12071–12075. Web of Science CrossRef PubMed Google Scholar
Brombal, L., Golosio, B., Arfelli, F., Bonazza, D., Contillo, A., Delogu, P., Donato, S., Mettivier, G., Oliva, P., Rigon, L., Taibi, A., Tromba, G., Zanconati, F. & Longo, R. (2018). Proc. SPIE, 10573, 1057320. Google Scholar
Brun, F., Delogu, P., Longo, R., Dreossi, D. & Rigon, L. (2017). Meas. Sci. Technol. 29, 014001. Web of Science CrossRef Google Scholar
Burvall, A., Lundström, U., Takman, P. A., Larsson, D. H. & Hertz, H. M. (2011). Optics Express, 19, 10359–10376. Web of Science CrossRef Google Scholar
Contillo, A., Veronese, A., Brombal, L., Donato, S., Rigon, L., Taibi, A., Tromba, G., Longo, R. & Arfelli, F. (2018). Radiol. Oncol. doi:10.2478/raon-2018-0015. Google Scholar
Delogu, P., Brombal, L., Di Trapani, V., Donato, S., Bottigli, U., Dreossi, D., Golosio, B., Oliva, P., Rigon, L. & Longo, R. (2017b). J. Instrum. 12, C11014. Web of Science CrossRef Google Scholar
Delogu, P., Golosio, B., Fedon, C., Arfelli, F., Bellazzini, R., Brez, A., Brun, F., Lillo, F. D., Dreossi, D., Mettivier, G., Minuti, M., Oliva, P., Pichera, M., Rigon, L., Russo, P., Sarno, A., Spandre, G., Tromba, G. & Longo, R. (2017a). J. Instrum. 12, C01016. Web of Science CrossRef Google Scholar
Delogu, P., Oliva, P., Bellazzini, R., Brez, A., De Ruvo, P., Minuti, M., Pinchera, M., Spandre, G. & Vincenzi, A. (2016). J. Instrum. 11, P01015. Google Scholar
Fedon, C., Longo, F., Mettivier, G. & Longo, R. (2015). Phys. Med. Biol. 60, N311–N323. Web of Science CrossRef Google Scholar
Knoll, G. F. (2010). Radiation Detection and Measurement. New York: John Wiley and Sons. Google Scholar
Longo, R. (2016). Nucl. Instrum. Methods Phys. Res. A, 809, 13–22. Web of Science CrossRef Google Scholar
Longo, R., Arfelli, F., Bellazzini, R., Bottigli, U., Brez, A., Brun, F., Brunetti, A., Delogu, P., Di Lillo, F., Dreossi, D., Fanti, V., Fedon, C., Golosio, B., Lanconelli, N., Mettivier, G., Minuti, M., Oliva, P., Pinchera, M., Rigon, L., Russo, P., Sarno, A., Spandre, G., Tromba, G. & Zanconati, F. (2016). Phys. Med. Biol. 61, 1634–1649. Web of Science CrossRef CAS Google Scholar
Mettivier, G., Fedon, C., Di Lillo, F., Longo, R., Sarno, A., Tromba, G. & Russo, P. (2016). Phys. Med. Biol. 61, 569–587. Web of Science CrossRef CAS Google Scholar
Mozzanica, A., Bergamaschi, A., Brueckner, M., Cartier, S., Dinapoli, R., Greiffenberg, D., Jungmann-Smith, J., Maliakal, D., Mezza, D., Ramilli, M., Ruder, C., Schaedler, L., Schmitt, B., Shi, X. & Tinti, G. (2016). J. Instrum. 11, C02047. Web of Science CrossRef Google Scholar
Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567–8591. Web of Science PubMed Google Scholar
Paganin, D., Mayo, S., Gureyev, T. E., Miller, P. R. & Wilkins, S. W. (2002). J. Microsc. 206, 33–40. Web of Science CrossRef PubMed CAS Google Scholar
Park, S. E., Kim, J. G., Hegazy, M., Cho, M. H. & Lee, S. Y. (2014). Proc. SPIE, 9033, 90335. Google Scholar
Pennicard, D. & Graafsma, H. (2011). J. Instrum. 6, P06007. Google Scholar
Rivers, M. (1998). Tutorial introduction to X-ray computed microtomography data processing. https://www.mcs.anl.gov/research/projects/X-ray-cmt/rivers/tutorial.html. Google Scholar
Sarno, A., Mettivier, G., Golosio, B., Oliva, P., Spandre, G., Di Lillo, F., Fedon, C., Longo, R. & Russo, P. (2016). Phys. Med. 32, 681–690. Web of Science CrossRef Google Scholar
Taguchi, K. & Iwanczyk, J. S. (2013). Med. Phys. 40, 100901. Web of Science CrossRef PubMed Google Scholar
Takahashi, T. & Watanabe, S. (2001). IEEE Trans. Nucl. Sci. 48, 950–959. 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
Vedantham, S., Shrestha, S., Karellas, A., Shi, L., Gounis, M. J., Bellazzini, R., Spandre, G., Brez, A. & Minuti, M. (2016). Med. Phys. 43, 2118–2130. Web of Science CrossRef Google Scholar
Vincenzi, A., de Ruvo, P., Delogu, P., Bellazzini, R., Brez, A., Minuti, M., Pinchera, M. & Spandre, G. (2015). J. Instrum. 10, C04010. Web of Science 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.