Information-based analysis of X-ray in-line phase tomography with application to the detection of iron oxide nanoparticles in the brain

: The study analyzes noise in X-ray in-line phase tomography in a biomedical context. The impact of noise on detection of iron oxide nanoparticles in mouse brain is assessed. The part of the noise due to the imaging system and the part due to biology are quantitatively expressed in a Neyman Pearson detection strategy with two models of noise. This represents a practical extension of previous work on noise in phase-contrast X-ray imaging which focused on the theoretical expression of the signal-to-noise ratio in mono-dimensional phantoms, taking account of the statistical noise of the imaging system only. We also report the impact of the phase retrieval step on detection performance. Taken together, this constitutes a general methodology of practical interest for quantitative extraction of information from X-ray in-line phase tomography, and is also relevant to assessment of contrast agents with a blob-like signature in high resolution imaging.


Introduction
X-ray phase contrast imaging is a field of coherent imaging which is receiving growing interest due to its greater sensitivity over conventional attenuation-based techniques (see [1] for an introduction to this field). In X-ray phase-contrast imaging, refraction of a partially coherent X-ray beam by the object of interest slightly modifies the original wavefront profile. These variations result in changes in the locally transmitted intensity of the wave which contains quantitative information on the phase shift induced by the object. X-ray phase contrast imaging has been shown to be particularly relevant for application in soft biological tissues composed of elements with low atomic numbers (see [2,3] for an overview on biomedical applications). X-ray phase contrast imaging in biomedical applications is currently mainly used for image rendering taking advantage of the edge-enhanced contrast at the interface between distinct tissues. Over the last two decades, a number of acquisition techniques coupled with phase retrieval algorithms have been developed. The performances of such acquisitions and reconstruction techniques are usually assessed for methodological purposes in terms of linear fidelity of 3D reconstruction with respect to calibrated phantoms [4].
Investigations of X-ray phase contrast imaging are progressively moving beyond simple linear reconstruction and visualization to consider this imaging technique from an informational point of view. Quantities such as contrast and signal-to-noise ratio (SNR) in various phase contrast imaging techniques have for instance only recently been investigated analytically [5][6][7][8] with the aim of achieving analytical expressions and experimental comparison with calibrated phantoms. In [6] and [5], theoretical analytical expressions were derived for the signal-to-noise ratio, resolution and contrast, as a function of the physical experimental parameters involved in X-ray in-line phase contrast imaging for a generic unidimensional phase edge. In [8] and [7], a theoretical and experimental study compare three different X-ray phase contrast imaging techniques in terms of signal-to-noise ratio also on a mono-dimensional phase edge. The noise model used in these works was assumed to be Gaussian or modeled as Poisson noise to account for the statistical nature of the photons. Other sources of perturbations, not included in [5][6][7][8], have been reported in phase contrast imaging. Because of the coherent nature of the X-ray source, speckle noise, controlling the second order statistics of the noise, can be recorded [12]. Such second order statistics of the noise for detection tasks with in-line X-ray phase-contrast imaging have been studied in [9][10][11] by simulations with homogeneous materials. In addition, when phase contrast imaging is coupled to computed tomography, it is well-known that spurious inhomogeneities in the sensitivity of the scintillator converting X-ray to photons produce so-called "ring artifacts" in the phase-contrast reconstruction (see for instance [13][14][15]). These sources of noise are combined with the electronic noise of the CCD camera placed at the end of the optical setup. It could therefore be possible to extend the realism of the noise model from [5][6][7][8][9][10][11]. Such considerations, common in other fields of coherent imaging [16], are very important in the perspective of nonlinear informational tasks such as detection.
For a specific pattern to be detected, however, signal-to-noise ratio has to be defined and detection performance has to be assessed. In [5][6][7][8] a simple mono-dimensional phase edge is considered. With such a mono-dimensional phantom, the detection scheme is only sensitive to the first order statistics of the noise. However, this does not take into account the possible correlation of the noise included in the second order statistics. In this article, we address this issue in a practical context with the detection problem of bidimensional heterogeneous objects with one phase contrast imaging technique, namely in-line phase tomography (IPT) in a biomedical application framework: the detection of iron oxide nanoparticles in excised mouse brain. The use of ultrasmall superparamagnetic particles of iron oxide as contrast agent in biomedical imaging in general [17] and in the study of stroke in particular [18] is a field of current interest, and was recently shown to be feasible and effective in a preclinical study of stroke [19]. The present study sought to go beyond simple visualization as in [19] and to undertake the automatic detection of these microstructures included in the lesion. We use this application of practical biomedical interest to study the influence of the noise and phase reconstruction parameter on the performance of the detection task. The paper is organized as follows. The first section, is a material and method section, which describes the experimental setup and the algorithmic pipeline for the phase reconstruction and the noise analysis and modeling. Then, in the results section, these noise models are used to assess the performance of detection of the nanoparticles in terms of probability at various signal-to-noise ratios and we explore the influence of the reconstruction parameter of the phase-contrast on the signal-to-noise ratios and on the detection performance of our algorithm.

Image acquisition
The experimental image acquisition setup is the one described in detail in [19]. Briefly, IPT acquisitions in the propagation based imaging mode are performed on the synchrotron radiation beamline ID19 of the European Synchrotron Radiation Facility (ESRF) in Grenoble, France. As shown in Fig. 1(a), the samples are 8 post-mortem mouse brains with permanent middle cerebral artery occlusion [20], which had received an intravenous injection of iron oxide nanoparticles. Iron oxide nanoparticles are engulfed by macrophages in vivo and thus used as imaging markers in inflammatory disorders associated with elevated phagocytic activity such as ischemic stroke. The cerebral artery occlusion causes a lesion in the cortex of the brain, colorized in purple in Fig. 1(c). As demonstrated in [19] and visible in Fig. 1(d) the nanoparticles accumulate in the lesion where they are detected as hyperintense small areas. These areas correspond to spherical microstructures of about 40 μm, compatible with macrophages which have engulfed iron oxide nanoparticles.
Image processing pipeline used in this study. Imaging is performed with 17.6 keV selected from undulator radiation using Al filters. The X-ray beam transmitted through the specimen is acquired on a detector using a LuAg scintillator screen, visible light optics, and a 2048 × 2048 CCD detector. The detector is positioned at 1 meter from the sample for in-line phase-contrast imaging. The 3D volume is a stack of 1,000 slices of 2048 × 2048 voxels and slice thickness is equal to pixel size (isotropic voxels) 8μm.

Image reconstruction
Phase retrieval is performed from a single phase-contrast image at each projection angle, using Paganin's method [21]. With this phase retrieval method, the complex refractive index n is used to describe attenuation and phase shift caused by the imaged object, where β is responsible for attenuation and δ for beam phase shift. The complex amplitude U of an attenuated plane wave of wavelength λ , propagating along direction x, in a homogeneous medium with complex refractive index n, can be described as The refractive index of a homogeneous sample of thickness d thus causes a phase shift Φ = −2πδ d/λ on the X-ray wave. The phase retrieval method chosen in this report assumes proportionality between the phase term δ and absorption term β of the complex refractive index. The proportionality coefficient is the δ /β coefficient. This δ /β coefficient is a tunable parameter either selected based on the prior knowledge of the type of tissue imaged by using tables like in [24,25] or adjusted empirically. As visible in Fig. 2, this influences the contrast of  [25] is given for linearity between grey level and phase, assuming the imaged object as homogeneous. However, it is not easy to determine the optimal choice of a single δ /β for biological samples which are highly heterogeneous. In [19], the δ /β was set at 321 to correspond to the highest nanoparticle concentration in stereotaxically injected animals. In our case, the final informational task is not to produce grey levels proportional to a physical quantity but to produce a nonlinear binary decision between pixels containing iron oxide nanoparticles and pixels in the background lesion. In the results section, we therefore propose to adjust the δ /β taking into account our detection task. After phase retrieval, the phase maps are used as input to a 3D parallel-beam tomographic reconstruction algorithm based on the filtered back projection algorithm described in [22]. For each scan, this procedure provides a reconstructed 3D volume. As depicted in Fig. 1(a), the reconstructed 3D volume suffers from spatially circular noise typical of reconstruction artifacts in tomography caused by differences in the individual pixel response of the detector. These ring artifacts are removed with a nonlinear median-based filter algorithm of [14] to produce Fig. 1(b). A 3D segmentation of the ischemic lesion is then performed. Since the mice in the stroke model systematically undergo right middle cerebral artery occlusion, lesion location is known, and could be used to create a lesion atlas, by aligning the 8 scans of mice brain, enabling the 3D average mask of the lesion presented in Fig. 1(c) to be created. The macrophages containing the iron oxide nanoparticles appear in the lesion in the form of blobs of sizes between 16 μm to 48 μm. For the detection of these microstructures, we use the multiscale blob detection algorithm introduced in [23] and described in Fig. 3. Using prior information concerning the search for spherical structures, the two-stage algorithm first convolves the reconstructed images with Laplacian of a Gaussian filters at multiple scales, r representing the standard deviation of the filters. Circular blobs of size equals to r produce local minima in the map resulting from this convolution. The second stage of the algorithm keeps local minima in these maps across the scales to produce the estimated positions (x c ,ỹ c ) and estimated radiusr of detected circular blobs.
The image processing pipeline presented in this section is an association of standard image processing tools and other choices of algorithm could be made. In what follows, the focus is on the impact of the noise and the reconstruction parameter δ /β on the performance of such a detection scheme. [r min = 16μm; r max = 48μm]. A blob is detected for each local minimum in the detection maps. Restriction is made on the amplitude down to a given percentage of the amplitude of the smallest minimum. This percentage is arbitrarily set here as 60%. The process is applied slice by slice.

Noise analysis
We propose to analyze the noise remaining after the tomographic reconstruction and the denoising of the ring artifacts as in Fig. 1(b). To this end, as depicted in Fig. 4, a region of interest is cropped in the background of the image, identified as background noise N 1 .
For statistical analysis, the region of interest in the background noise N 1 is 3685 pixels and is applied on all the slices of the 8 scans of the mouse brains. As illustrated in Fig. 5, the histograms for the background noise showed a wide variety of shape including mono-and bimodal distributions, symmetric and non-symmetric distributions. It is therefore not easy to propose a simple statistical model for this distribution. The Shapiro-Wilk test for normal distribution is negative with a p-value of the order of 10 −13 (whereas it would have to be larger than 0.05 for the test to be positive) for all the slices tested for the background noise N 1 . Figure 6 shows the empirical standard deviation σ N 1 and average M N 1 . The average value M N 1 is displaying nonmonotonic evolutions reproducible from one mouse brain to another. The standard deviation σ N 1 is restricted to a 1 to 8 interval.
From these observations, we propose the following "empirical" model to simulate noisy IPT images S 1 (x, y) where (x, y) are the spatial coordinates. In Eq. (3), B(x, y) are simulated binary blobs. M B is the average value recorded in the iron oxide nanoparticles areas. We assume in Eq. (3) an additive signal noise coupling and propose to simulate the background noise N 1 (x, y) by randomly picking centered and normalized instances of this noise from the 8 scans containing 1,000 slices each. However, the background noise N 1 does not take into account the biological spatial fluctuations due to the microstructures of the lesion distinct from the targeted iron oxide nanoparticles. To take account of this biological noise, we cropped a region of interest called lesion noise N 2 in Fig. 4 located in a part of the lesion that is visibly free of iron oxide nanoparticles. The distribution of the lesion noise N 2 are found similar to those presented for the background noise N 1 in Fig. 4, i.e. a large variety of distribution shapes and similar fluctuations for average M N2 and standard deviation σ N2 . The same kind of empirical model is therefore constructed for lesion noise N 2 using stacks of images of this noise and assuming additive coupling between the blobs with an equation similar to Eq. (3). The lesion noise N 2 shows some spatial correlation due to biological microstructures while background noise N 1 visually appears closer to a spatially uncorrelated white noise. The comparison of the performance of our detection scheme shown in Fig. 3 with background noise N 1 and lesion noise N 2 will therefore be useful to quantify the impact of the second order statistics introduced by biological tissues. For comparison with a statistical process of reference, we will also consider, in the results section, with the same additive coupling as in Eq. (3), a purely white noise N 3 with Gaussian distribution centered on the same average value M N 2 as the lesion noise N 2 and with standard deviation σ N 2 .

Detection performance
We propose to simulate the nanoparticles with B(x, y) in Eq. (3) taken as binary images with black background including white disks with center randomly positioned at (xc true , yc true ) in B(x, y) and radius r true uniformly drawn from the interval [16,48] μm, the expected scale range for iron oxide nanoparticles. The evaluation of the detection performance based on these simulated binary blobs are expressed in terms of accuracy of the estimated center position (xc,ỹc) and radiusr of the blobs as well as in terms of the probability of good detection and false alarms. The accuracy of the estimation is computed with the average error and where · stands for the empirical average. For the probability of good detection and false alarms, we call H 0 the hypothesis that there is no blob to be detected and H 1 the hypothesis that there is one blob to be detected. Let us call D 0 the decision to detect no blob and D 1 the decision to detect one blob. The probability of good detection is defined as Pr(D 1 | H 1 ) and the probability of false alarm is defined as Pr (D 1 | H 0 ). These probabilities are estimated empirically over some 10,000 trials.

Results
We are now ready to assess the detection scheme shown in Fig. 3 with the noise models of the previous section, i.e. background noise N 1 , lesion noise N 2 and white Gaussian noise N 3 . Figure 7 shows how the detection performances, Er radius , Er center , Pr(D 1 | H 1 ) and Pr(D 1 | H 0 ) degrade as the noise level is raised in the range of noise standard deviations observed experimentally. In the range of noise standard deviations σ ∈ [1,8], the detection scheme performs better with the Gaussian white noise model N 3 than with the background noise N 1 . This confirms that the background noise N 1 in IPT is not a white Gaussian noise. Also, in Fig. 7 we observe that the detection scheme performs better with the background noise than with the lesion noise N 2 for all the observed noise standard deviations larger than σ = 2. For small noise regimes the effect of the background noise N 1 due to the apparatus has the same impact as the lesion noise N 2 incorporating the biological microstructures distinct from nanoparticles. But for greater noise (σ > 2), in contrast, the microstructures in the lesion noise N 2 , when appearing in the spatial range [16,48]μm, bias the estimation of correctly detected blob or/and constitute false positive detection. The distance between the curve for the background noise N 1 and the lesion noise N 2 in Fig. 7 quantifies this effect.
A B C D In IPT, the signal-to-noise ratio controlling the detection performance is impacted not only by noise but also, as shown in Fig. 2, by the phase-retrieval algorithm. The impact of the choice of the δ /β coefficient on the detection scheme shown in Fig. 3 is quantitatively assessed using a the Fisher contrast-to-noise ratio (FR), defined as where μ is the mean and σ the standard deviation. The signal regions are manually segmented in areas considered as compatible with iron oxide nanoparticles. The noise regions are areas surrounding the signal region. The Fisher ratio FR of Eq. (6) therefore constitutes a local contrast to noise measure quantifying the detectability of iron oxide nanoparticles. Figure 8 presents the evolution of the local contrast to noise measure FR as a function of the δ /β coefficient: the smaller the δ /β the higher the detectability of the iron oxide nanoparticles. This detectability estimation from the deflection is confirmed in Fig. 9(a) in terms of receiver operator characteristic (ROC) curves which give the probability of good detection of nanoparticles Pr(D 1 | H 1 ) as a function of the probability of false detection of nanoparticles Pr(D 1 | H 0 ). ROC curves plotted for various δ /β coefficients were closer to the upper left corner for smaller δ /β values. The ROC curves obtained with the smallest δ /β values were significantly distant from the ROC curve obtained for δ /β = 321, which would be the recommended coefficient for linear measurement of the tissue. This enables quantification of the gain in detection of the information-based approach for tuning δ /β . In such an approach, it is important to understand that the tuning of δ /β is problem-dependent. To illustrate this, we have plotted in Fig. 8 the evolution of the deflection between the lesion and the healthy surrounding tissue. In this case, the larger the δ /β the better is the detectability. This is also confirmed in Fig. 9(b) in terms of ROC curve which are closer to the upper left corner for the largest δ /β . As visible in Fig. 8, the best contrast FR between the lesion and the healthy tissue (obtained for large δ /β ) is larger than the best contrast between nanoparticles and lesion (obtained for low δ /β ). This explains the fact that for high δ /β in Fig. 9(b), the ROC curves for large δ /β are very close to the upper left corner and almost superimposed. These findings demonstrate that the optimal value for the phase-retrieval parameter δ /β in IPT depends on the final informational task. Following this observation, various scenarios of informational-based reconstruction schemes can be imagined. If detection of multiple classes is targeted, it would be interesting to perform multiple reconstructions. For instance here, it would be interesting to perform one reconstruction at δ /β = 1250 for the detection of the lesion area and then another reconstruction with δ /β = 50 for the detection of nanoparticles within the lesion. This is a proposal similar to the one formulated in [26] in terms of signal-to-noise ratio and applied to multi-material phantoms. The study [26] proposes a numerically efficient phase retrieval algorithm to reconstruct the complex refractive index distribution of known homogeneous materials embedded in a second homogeneous medium from a single X-ray phase contrast image per projection. This algorithm makes use of a single propagation-based-imaging image per projection, separately and selectively reconstructing each interface between any given pair of distinct homogeneous materials. But this is not the case in the present study, where we are not working with known homogeneous materials. The biological tissues in the brain are heterogeneous and we therefore have no prior knowledge of the δ /β values corresponding to each feature (lesion, nanoparticles, etc.) and local background. The present study thus extends to heterogeneous material with unknown exact composition the idea introduced in [26], that, when multiple structures are to be detected, it is useful to optimize phase retrieval separately at each interface between the structures. Also, if only one reconstruction is computationally feasible, it could be interesting to optimize reconstruction in terms of detection for both regions. One criterion could be to choose a reconstruction δ /β value of 600, where the two contrast curves intersect and are equal in  Fig. 8. It is again important to note that this optimal reconstruction parameter value for detection is different from the one (δ /β = 321) that was optimal for visualization in the present case. Statistics computed on 18 lesions representing a total of 1,080 pixels for the signal and 272,700 pixels for the noise. The solid line stands for FR between nanoparticles and lesion and the dashed line represents FR between mouse brain lesion and surrounding healthy tissue.

Conclusion and discussion
The present study assessed the impact of noise in IPT applied to a detection problem in a biomedical context. The method quantified the part due to the statistical noise of the IPT system and the part due to the biological tissue. Tuning the phase-retrieval reconstruction parameter was shown to benefit to the detection task in comparison to the conventional tuning implemented for metrological or visualization purposes only. This used a very common phaseretrieval method in IPT [21]. The theoretical derivation of this Paganin method [21] is based on a number of assumptions: the X-ray beam is assumed to be partially coherent and monochromatic; the distance between the object and the detector should be close enough to fulfill the near-field condition; and the object is assumed to be composed of a single homogeneous material [27]. Knowing the energy of the beam and the composition of the sample enables the ideal reconstruction parameter δ /β to be selected [28]. However, in biomedical applications, the object is in general inhomogeneous and it is difficult to know its exact composition. Thus the question of how to select the δ /β value arises. This is generally done by selecting a range of δ /β values and heuristically choosing the best image from a visual point of view. One originality of the present study was to provide objective criteria for selecting the δ /β parameter based on the detection task.
In this informational approach, there is a range of open possibilities. From a biomedical point of view, the present detection scheme could be applied in semi-automatic segmentation of ironlabeled macrophages in the ischemic lesion area, a task that would be too time-consuming to perform manually. The detection scheme was applied in 2D images for purposes of illustration but could easily be extended to 3D by considering the 3D Laplacian of Gaussian filters. Also, the approach could easily be extended to multi-class detection by using multiple phaseretrieval reconstruction parameter values with a different δ /β value for each class. Assessment of nanoparticle detection performance in all the lesion tissue would require comparison with a global ground truth providing the total amount of nanoparticles fixed in the brain in comparison to the volume of nanoparticles injected. This would require having nanoparticles which very specifically target the inflammatory part of the lesion. Alternatively, further applications include the assessment of the ability of new contrast agents to specifically target a cerebral ischemic lesion, and the detection of iron oxide nanoparticles in other pathologies, such as atherosclerosis. From a statistical point of view, the present results showed that the first-order statistics of the background noise in IPT were not Gaussian. It would be interesting to further investigate this statistical characterization and seek to construct first-order statistical models of background noise. The noise in biological tissue and in the background appeared as non-white (i.e., correlated). This is in agreement with the simulation performed in [9][10][11] with homogeneous materials. It would therefore also be interesting to investigate the modeling of the second-order statistics and their dependence on the phase-retrieval parameters, but with heterogeneous material. However, X-ray Compton scattering caused by the observed samples may constitute a source of correlation, making the noise sample-dependent. Simulations of Xray in-line phase tomography including Compton scattering from a raw reconstruction of the sample could constitute another direction for noise modeling. The information-based approach for phase retrieval optimization was illustrated here with one specific technique: X-ray in-line phase tomography. Other X-ray phase-contrast imaging techniques, involving interference with a reference or the use of gratings, may also benefit from an informational approach where the phase reconstruction algorithm is optimized in terms of the final extracted information.