1. Introduction

Analyzer-based imaging (ABI) is a phase-contrast imaging technique that has attracted much interest in the scientific community in the recent years. Although the main physical principles of the technique have been known for many years [1], it is just with the development and spread-out of third generation synchrotron radiation facilities that this technique has been thoroughly theoretically and experimentally investigated.

In conventional absorption imaging, the contrast is generated by variations of the X-ray absorption coefficient that arise from density differences and from changes in thickness and composition of the sample.

Unlike absorption-based imaging, ABI derives contrast also from the phase modulations induced by the object onto the transmitted X-ray beam. These two effects can be described in terms of a complex index of refraction, which can be indicated as n = 1-δ-iβ. The real part δ corresponds to the phase shift due to refraction and the imaginary part β to the absorption. Since the δ term is much larger than β and has stronger variations in the hard X-ray regime, ABI is in principle able to provide improved contrast compared to absorption-based imaging, and therefore holds great promises in many application fields such as materials science and clinical imaging.

The typical ABI setup consists of a parallel monochromatic X-ray beam, which is used to irradiate the sample, and a perfect crystal, called the analyzer, placed between the sample and the detector. The analyzer acts as an angular filter of the radiation transmitted through the object, since only the X-rays travelling in a narrow angle window close to the Bragg condition are diffracted onto the detector [2]. Before being detected, the beam is modulated by the angular-dependent reflectivity of the crystal, its rocking curve (RC), which has a full-width half maximum (FWHM) typically of the order of a few microradians.

Besides absorption and refraction (the latter providing contrast thanks to the modulation given by the rocking curve (RC)) also X-ray scattering plays an important role in the generation of the image contrast. Scattering originates by very small internal structures that are not resolved by the detector. Diffuse scattering in the milliradian range, the so-called small-angle X-ray scattering, which is due to structures at the nanometer scale, is rejected by the analyzer and gives rise to what is called the extinction contrast. Ultra-small angle X-ray scattering (USAXS), instead, which is due to structures ranging in the order of magnitude of hundreds of nanometers up to micrometers and is characterized by angles in the microradian scale, partially falls within the acceptance of the crystal analyzer and has the effect of broadening the observed rocking curve.

The contrast in the recorded AB images is therefore given by a mixture of absorption, refraction, USAXS and small-angle X-ray scattering rejection (the latter three effects having the same physical nature as demonstrated by Davis [3]). The information contained in the images is thus very rich but the image interpretation can be in some cases ambiguous due to the signals superposition. In order to both effectively separate the different physical effects and accurately quantify them, several mathematical methods have been so far proposed. They allow calculating the corresponding physical parameters by means of combining two or more images acquired at different positions on the analyzer RC.

Most of the proposed methods are based on the geometrical optics (GO) approximation, which imposes some restrictions on the imaged object. This approximation, in fact, is strictly valid only if the phase of the wave incident on the crystal analyzer is a slowly varying function on the length scale of the extinction length of the crystal, which has been shown to be equivalent to the condition N_T >> 1, where N_T is the so-called Takagi number [4, 5].

Other methods based on different approximations have also been developed. The method proposed by Nesterets et al. [6] is based on the so-called weak object (WO) approximation, which requires that the phase shifts introduced by the object are weak. Another method introduced in literature [7] is, instead, based on the linear transfer function (LTF) approximation, which requires that the transfer function of the imaging system can be linearized in the Fourier space. The application of these last two methods, however, requires the knowledge of the imaging system PSF (point spread function), which is a complex function that needs specific measurements and calculations. Furthermore, additional terms to the expression of the imaging system PSF need to be considered if the incoming beam is not perfectly parallel and monochromatic [8], which is usually the case in common experimental conditions. Under the GO approximation, instead, just the knowledge of the imaging system RC, which can be easily measured, is needed.

The different GO algorithms have been widely applied to extract quantitative information from ABI images in different application fields, and particularly in biological tissues imaging. A comprehensive and systematic quantitative comparison of the most used GO extraction algorithms is not available in the literature to the best of our knowledge. In particular, no complete quantitative comparison of the different methods in terms of their accuracy in the presence of variable amounts of refraction and scattering has been performed. This is an aspect of crucial importance especially in biological and medical imaging. The choice of the algorithm to use may significantly influence the image quality and the accuracy of the extracted information, and therefore, ultimately, the image interpretation. Moreover, biological samples may differ greatly in the amount of absorption, refraction and scattering they produce. It is thus important to assess the performance of each algorithm under different experimental conditions simulating the cases that can be encountered with biological tissues. The final goal of the present study is to identify which algorithm could be the most suitable for a particular biological application.

In this paper, a theoretical and experimental comparison of the main GO algorithms is presented. Different experimental conditions characterized by various amounts of absorption, refraction and scattering are considered. The quantitative comparison is performed by using simple geometrical phantoms consisting of plastic materials (exhibiting absorption and refraction) and of paper (exhibiting absorption and scattering). The same algorithms are then applied to a human bone-cartilage sample, as a challenging example in biomedical imaging.

The paper is organized in the following way. In section 2 the various GO algorithms are introduced and discussed, with particular attention paid to their underlying assumptions and applicability range. In section 3, the experimental setup as well as the geometrical phantoms and the human bone-cartilage sample that have been used in the experiments are described. In sections 4 and 5, results obtained by applying the various algorithms to experimental images of the considered samples are shown and discussed.

2. Algorithms for quantitative analysis of AB images

2.1 Diffraction-enhanced imaging (DEI) algorithm

The diffraction-enhanced imaging (DEI) algorithm was developed by Chapman et al. [9] and is based on linearizing the rocking curve recorded without the sample at its two slopes. In this approach just two images of the sample, acquired with the analyzer positioned respectively at the “low” and at the “high” slopes of the RC, where the second derivative of the RC is approximately zero, are needed. The images of absorption and refraction angle are then obtained analytically [9]. The main hypothesis on which the algorithm is based is that the refraction angles are small compared to the FWHM of the RC, so that the RC can be well described by a first-order Taylor approximation. Additionally, this method does not take into account the USAXS produced by the object, which has the effect of broadening the observed RC compared to the reference one that is used in the calculations.

This algorithm requires just two input images, which is an asset in terms of dose deposited to the sample and in terms of overall duration of the image acquisitions. The applicability field is anyway limited because the assumptions concerning the sample are quite restrictive. Many biological tissues, in fact, create a non-negligible amount of USAXS. Moreover, if the refraction angles are of the order of the RC FWHM or bigger, the first-order Taylor approximation fails to reproduce accurately the RC and results are quantitatively (and also qualitatively) incorrect [10–12].

2.2 Extended DEI (E-DEI) algorithm

In order to overcome the intrinsic limitations of the DEI algorithm, other analytical algorithms have been proposed that allow for separation of the absorption and refraction contributions by using two input images without imposing a Taylor approximation to the analyzer RC [13, 14]. These methods are usually referred as extended-DEI (E-DEI) algorithms.

In this paper, the approach proposed by Hu et al. is followed [14]. In this case, the imaging system RC, measured without the sample, is fitted by a Gaussian function. If two images are taken at different angular positions of the analyzer, θ_1,2, the following expression for the intensity recorded after the analyzer can be written:

Since this method does not impose any restrictions on the values of the refraction angles (provided that the GO conditions are still fulfilled [4, 5]) the range of refraction angles that can be calculated is wider than in the case of the DEI algorithm.

Furthermore, the two images do not need to be acquired exactly at the two slopes of the RC but can be chosen anywhere along it.

2.3 Generalized DEI (G-DEI) algorithm

Since the amount of scattering introduced by the sample is in many cases not negligible (notably, for biological tissues), different algorithms have been developed in order to attempt to separate this effect from absorption and refraction.

The algorithm proposed by Rigon et al. [15] is based on a second-order Taylor approximation for the RC. In this case, the intensity recorded after the analyzer can be expressed as:

If three images are acquired at three arbitrary analyzer angular positions, the parametric images of refraction, absorption and scattering can be analytically calculated [15]. Since the second-order Taylor expansion is a good approximation of the RC just for small angular deviations, the accuracy of this algorithm is limited to values of σ _ΔθS and Δθ_R that are small compared to the FWHM of the RC.

2.4 Multiple image radiography (MIR) algorithm

Pagot et al. [11] and Wernick et al. [10] independently developed two statistical methods which allow for reconstructing the RC on a pixel-by-pixel basis, by conveniently combining several images at different positions along the RC.

In this paper, the approach of Pagot et al. [11] will be followed. The method consists in taking two series of N (N ≥ 3) images at different positions of the crystal analyzer with and without the sample, respectively. The ‘reference’ RC (acquired without the sample) and the ‘object’ RC (acquired with the sample) are then compared. If it is assumed that for each pixel the angular distribution of the diffracted intensity is the convolution of the object angular spectrum with the imaging system RC, then the refraction, integrated absorption and scattering images can be calculated respectively from the zeroth-, first- and second- moments of the ‘reference’ and ‘object’ RCs [11]. An additional maximum absorption image can be also calculated, which is defined as the ratio of the ‘object’ RC maximum with respect to the ‘reference’ one. Whereas the integrated absorption image computes the area under the measured RC, the maximum absorption computes the value of the measured RC peak.

Similarly to the G-DEI algorithm, the MIR method has the advantage that the USAXS contribution is explicitly taken into account. In addition, it is in principle very stable with respect to noise, since many images are combined in order to calculate the different parameters. A notable drawback of this method is that, since many images need to be acquired, both the dose to the sample and the acquisition time are increased.

2.5 Pixel-by-pixel Gaussian curve fitting (GCF) algorithm

As Huang et al. pointed out [12], high refraction angles can be underestimated in the MIR algorithm. Under this approach, in fact, the refraction angle is calculated as the shift of the ‘object’ RC centroid relative to the ‘reference’ RC centroid. Theoretically, the centre of the ‘object’ RC and its centroid are equal only if the sampling along the RC is continuous and made on an infinite range. Since the number of sampling points is limited in the experiment, the centroid and centre values may differ. It can actually be shown that, in particular for high refraction angles, the calculated values always underestimate the actual ones.

The situation is schematically presented in Fig. 1 . Let suppose that the ‘reference’ RC is a Gaussian function with a FWHM of 2.94 μrad (this corresponds to the FWHM measured in our experiment, as described in section 3) and that two different ‘object’ RCs are considered: a first one characterized by a shift of the centre of 1 μrad, simulating a refraction effect, and a second one characterized by the same centre shift plus a broadening, simulating USAXS. For the latter ‘object’ RC, the FWHM_obj is 1.5 times the FWHM_ref. Let further suppose that the ‘reference’ and ‘object’ RCs are sampled with 31 analyzer positions ranging from −3 to +3 μrad with a 0.2 μrad step.

 

Fig. 1 Example of the finite sampling of a ‘reference’ RC (FWHM_ref = 2.94 μrad), of a first ‘object’ RC characterized by refraction (Δθ_z = 1 μrad) and a second ‘object’ RC characterized by refraction and USAXS (Δθ_z = 1 μrad, FWHM_obj = 1.5 FWHM_ref). The calculated refraction angle is Δθ_est,1 = 0.88 μrad and Δθ_est,2 = 0.62 μrad respectively in the absence and in the presence of scattering.

Download Full Size | PDF

When applying the MIR method to retrieve the refraction angle in the case of the first ‘object’ RC (case of pure refraction), a Δθ_z=0.88 μrad is calculated instead of the expected value of 1 μrad. The calculated value is smaller than the actual one and this is due to the fact that the ‘object’ RC is no longer symmetrically sampled with respect to its centre (peak) since the curve is shifted.

What no publication has pointed out yet, to our knowledge, is that this smoothing effect is even enhanced when, besides a high refraction angle, we are in the presence of considerable scattering. In the case of the second ‘object’ RC, in fact, where high refraction and a high amount of USAXS mix up, the MIR method gives an even lower calculated refraction angle, Δθ_z=0.62, than in the previous case. This means that, under this approach, high amounts of USAXS can increase the inaccuracy of the refraction angle estimate.

Similarly, it can be shown that under the conditions of high refraction and/or high USAXS, also the extracted USAXS signal can be underestimated. Furthermore, if the sampling range and/or the number of sampling points are decreased, the inaccuracy is expected to become more important.

An alternative approach which allows for overcoming these problems, the so-called Gaussian curve fitting (GCF) algorithm, has been proposed by Nesterets et al. [5]. This method consists in fitting, pixel-by-pixel, a Gaussian function to the ‘reference’ and ‘object’ RCs.

Under this approach, for each pixel the sampled ‘reference’ and ‘object’ RCs are fitted with the following Gaussian expression:

The quantities obtained with the GCF algorithm correspond to those calculated using the MIR algorithm, but they are expected to provide more accurate estimates of the different parameters since they are exempt from the above-described limitations.

An important drawback of this method is that it is computationally very intensive, since the fitting procedure has to be repeated for each pixel. For example, for an image of 10³ × 10³ pixels², the fitting calculation is repeated 10⁶ times and this leads to much longer computation times compared to those needed with the other algorithms.

3. Experimental methods

3.1 Set-up

The experiment was performed at the biomedical beamline (ID17) of the European Synchrotron Radiation Facility (ESRF, Grenoble, France). The X-ray source is a 21-pole variable-field wiggler with B_max = 1.6 T. A fixed-exit Laue-Laue Si(111) double-crystal monochromator allows to select energies in the range 25-150 keV. In this experiment, 26 keV X-rays were used (energy resolution ΔE/E ~2·10⁻⁴). The long source-to-sample distance (~151.5 m) in combination with the very small dimensions of the source (56.3 (H) x 10.3 (V) µm², as standard deviation) leads to a highly coherent X-ray beam.

Images have been recorded with a FReLoN CCD-based X-ray detector system [16] with 2048 x 2048 pixel² coupled with a 8 μm optics producing images with a spatial resolution of about 16 μm [17].

A Si(333) additional crystal monochromator and a Si(333) crystal analyzer were used to monochromatize the incoming beam and to analyze the refracted beam, respectively. The measured analyzer RC width (as FWHM) was 2.94 μrad.

3.2 Plastics samples

In order to perform a quantitative comparison between the different extraction algorithms, simple geometrical phantoms giving rise to either refraction or scattering signals were used.

The first phantom consisted of two cylindrical nylon wires of diameters of 350 μm and 200 μm, respectively (Fig. 2a ). They can be considered as pure phase objects, since their absorption contrast is almost negligible at the used energy. Instead, they show a very high refraction signal especially at the edges, where the refraction angles become very large. In this paper, only results for the 350 μm diameter wire are reported. The parameters for this wire are, for 26 keV: δ = 3.94 × 10⁻⁷; β = 1.39 × 10⁻¹⁰; absorption contrast = 1.2%.

 

Fig. 2 Sketches of the phantoms used in the experiment: 200 μm and 350 μm diameter nylon wires with overlapped paper layers (a), and Lucite grooves phantom with overlapped paper layers (b).

Download Full Size | PDF

The second phantom is a Lucite parallelepiped of section 40 × 40 mm² and thickness 2.9 mm, in which five almost cylindrically-shaped grooves (holes) of different radius are made (Fig. 2b). The section of the holes can be well approximated by a part of circumference. Unlike the polymer wires, this phantom produces both refraction and absorption. The refraction angles at the edges of the grooves are lower than in the case of the wires: the biggest groove, for example, shows a theoretical maximum refraction angle of 0.94 μrad. The characteristics of the grooves phantom are summarized in Table 1 .

Table 1. Some physical parameters of the groove phantom. The absorption contrast is calculated at the groove centre with respect to the maximum thickness of the groove phantom. For Lucite at 26 keV, δ=3.94 × 10⁻⁷; β=1.39 × 10⁻¹⁰

View Table

The third phantom is made of eight partially overlapping layers of newspaper. Paper is known to produce a large quantity of scattering and has been used for this scope in various publications concerning ABI extraction algorithms [10, 18]. Paper thickness ranges from 40±5 μm, in the region where just one layer of paper is present, to 320±5 μm, where all the eight layers of paper overlap. The wires and the grooves phantoms were imaged superimposed to the newspaper layers, in order to produce a variety of combinations of absorption, refraction and ultra-small angle X-ray scattering on the same image, as schematically presented in Figs. 2a and 2b.

Images of the wires were acquired at 9 different positions along the RC, respectively at the peak (100%, “top” position), ±75%, ±50%, ±30% and ±15% (positions relative to the maximum RC intensity). Images of the grooves phantom were acquired at 7 positions along the RC, respectively at the peak, ±50%, ±30% and ±15%.

3.3 Bone-cartilage sample

As a biological sample, a cylinder-shaped healthy cartilage-on-bone sample with a diameter of 7 mm was used. The sample was extracted from the lateral facet of a human patella using a shell auger. The cylinder was trimmed to a total height of 12 mm including the complete cartilage tissue and about 6 mm of subchondral bone.

During the images acquisition the sample was placed in a cylindrical container and was dipped into a saline solution. Images of the sample were acquired at 5 different positions on the RC, ±50%, ±15% and 100%, respectively.

3.4 Computer implementation

MATLAB (The MathWorks, Inc.) was used for the image analysis using the different extraction algorithms. Calculations were performed on a Dual-Core AMD Opteron 2218 Processor (64-bit, 2.6 GHz).

Concerning the GCF method, a built-in MATLAB function was used for the fitting of Eq. (3) on a pixel-by-pixel basis. This iterative fitting routine required initial estimates of the parameters to be calculated: the RC peak value A, centre Δθ_z and standard deviation σ obtained from the MIR method were used for this purpose.

DEI, E-DEI, G-DEI and MIR algorithms run fast on the used PC: the calculation time for any of these three methods, applied for example to a 10³ × 10³ pixels² image, was a few seconds only.

The GCF method, instead, requires much longer computation time. In fact, about 10 hours were needed to process the same image size.

4. Results and discussion

4.1 350 μm diameter nylon wire

The absorption, refraction and USAXS angle distribution images obtained by applying the different ABI algorithms to the images of the 350 μm diameter wire are reported in Fig. 3 . For DEI and E-DEI algorithms, only the AB images at the positions ±50% on the RC were used. These two images and the additional “top” one were used for the G-DEI algorithm, while all images at the 9 different RC positions were used for the MIR and GCF algorithms. Concerning the DEI and E-DEI methods, only the absorption and the refraction angle images are presented since no USAXS image is calculated.

 

Fig. 3 Calculated absorption, refraction angle and ultra-small angle images of the 350 μm diameter nylon wire overlapped to the newspaper layers.

Download Full Size | PDF

Concerning the MIR and the GCF absorption images, only the integrated absorption is shown in Fig. 3.

The paper layers are clearly visible in the absorption images. Also the nylon wire can be clearly seen but mainly because of the artefacts at the edges. In these regions the refraction angles are very high and this fact leads to the violation of the GO assumptions. Strong artefacts are therefore visible on the absorption images, in which the signal intensity itself is quite low.

In the calculated USAXS images, the different paper layers, which give different amounts of scattering, and the nylon wire edges are again clearly visualized. The strong USAXS signal generated at the edges of the wire is due to the very strong variations of the refraction angles, which may considerably increase the range of refracted X-rays falling within a single pixel. This leads to a broadening of the observed RC.

In the refraction images, just the nylon wire is visible, while the paper does not produce any visible signal. Images of the wire phase have been calculated by simple integration of the refraction angle images [19]. The profiles of the phase calculated with the different algorithms are reported in Figs. 4a and 4b, respectively in the case in which scattering can be neglected (no paper) and in the case in which scattering is maximum (8 layers). The calculated profiles are obtained by taking the average value over 30 horizontal pixels in order to reduce the contribution of noise. In Figs. 4a and 4b also the theoretical profile for the phase is shown, which can be calculated using the real part δ of the complex refractive index for nylon and the known pixel-by-pixel cross-section thickness of the wire [20].

 

Fig. 4 Cross-sections of phase images of the 350 μm diameter wire, when (a) no overlapping paper layers are present and when (b) 8 overlapping paper layers are present. The results extracted by DEI, E-DEI, G-DEI, MIR and GCF algorithms are compared with the theoretical profile. The plotted values are obtained by taking the average over 30 horizontal pixels.

Download Full Size | PDF

In the absence of scattering (Fig. 4a), it can be seen that for high refraction angles (comparable with the FWHM of the RC), like those given by the nylon wire, the agreement between the theoretical and the calculated values strongly depends on the used extraction algorithm. The best agreement (13% difference in the centre of the wire) is obtained when the experimental RC is fitted pixel-by-pixel to a Gaussian function (GCF method). Calculated values are most inaccurate in the case of the DEI algorithm, which is based on the linear approximation of the RC and fails when the refraction angles become comparable with the RC FWHM. The G-DEI algorithm, based on a second-order Taylor approximation that is strictly valid only if the refraction angles are very small, also significantly differ with respect to the theoretical values. The discrepancies between the GCF and the MIR algorithms are due to the fact that, for high refraction angles and when the number of the sampling points on the RC is limited, the calculated centroid of the sampled ‘object’ RC underestimates the actual centre of the ‘object’ RC. This last result perfectly agrees with the theoretical considerations presented in Section 2.5.

The difference between the algorithms is even more pronounced when high scattering is superimposed to the refraction signal (Fig. 4b). Both the DEI and E-DEI methods, which do not take into account the scattering, show the highest discrepancy compared to the theoretical values. Refraction angles are highly underestimated also in the case of the G-DEI and MIR algorithms. In the first case the difference can be explained by the fact that the second-order Taylor approximation fails when the refraction angle or the scattering distribution (or both) are not small compared to the RC FWHM. In the latter case, the discrepancy is due to the difference between the calculated centroid of the sampled ‘object’ RC and the peak position. As discussed in Section 2.5, such a difference increases when, additionally to the high refraction angles, the large scattering contribution broadens the measured RC. The GCF algorithm maintains also in this case a good agreement with the theoretical values. The comparison of the various algorithms presented in Figs. 4a and 4b shows that, under conditions of high refraction angles (~RC FWHM) and especially when, besides refraction, also the amount of scattering is important, the GCF method can provide much more accurate results than the other methods considered here.

In Fig. 5a the values of the phase at the centre of the wire, calculated with the different ABI algorithms, are plotted versus the number of paper layers superimposed to the sample. The aim is to compare the accuracy of the refraction angle calculation when scattering progressively increases. The theoretical estimation of the phase value is therefore reported as well. The graph in Fig. 5a shows that the more the scattering increases the more the difference between the algorithms becomes pronounced. The values calculated with the GCF algorithm show no evident variations while the values calculated with the other methods present a decreasing magnitude as the layers of paper increase. A drop of 42% of the phase for the DEI algorithm, of 44% for the E-DEI, of 31% for the G-DEI algorithm and of 25% for MIR is calculated by comparing the case of no scattering (x=0) and the case of 8 layers of paper (x=8) superimposed to the wire. As expected, the percentage decrement is higher in the case of the algorithms that do not take the X-ray scattering into account (DEI and E-DEI) than in the case of algorithms that explicitly consider it (G-DEI and MIR).

 

Fig. 5 Extracted phase values (a) at the centre of the 350 μm diameter wire and (b) at the centre of groove 5, plotted with respect to the number of overlapping paper layers. Results obtained with the DEI, E-DEI, G-DEI, MIR and GCF algorithms are shown. The plotted values are obtained by taking the average over 30 horizontal pixels; the error bars represent the standard deviation of the mean. The theoretical value is also shown for comparison.

Download Full Size | PDF

The fluctuations of the plotted values of the phase corresponding to different layers of paper may be attributed not simply to the image noise but also to the inhomogeneity of the paper itself. Besides microscopic unresolved structures that mainly generate a scattering signal, the paper contains internal structures of dimensions bigger than the pixel size which determine a real measurable phase signal.

4.2 Grooves phantom

The ABI extraction algorithms were applied also to images of the grooves phantom superimposed to the newspaper layers. Results for the phase values calculated at the centre of the biggest groove (groove 5) versus the number of overlapping paper layers are shown and compared to the calculated theoretical values in Fig. 5b. In the case in which no scattering is superimposed to the phantom, all the algorithms show a good agreement with the theory. The refraction angles introduced by the groove are not very high (even at the edges) compared to the RC FWHM, and therefore the validity conditions of the different algorithms are never strongly violated. As soon as the fraction of scattered X-rays increases, the differences between the various algorithms become more marked. In cases of high scattering, the refraction and the scattering signals cannot be completely separated. The percentage decrement between the calculated phase values under conditions of no scattering (x = 0) and the values calculated under conditions of maximum scattering (x = 8) is more important for those algorithms that do not consider scattering (40% in DEI, 44% in E-DEI) than for the other ones (29% in G-DEI method, 20% in MIR), as already seen in the case of the wire phantom. Results obtained with the MIR method are particularly interesting in this case. They show that the broadening of the measured RC due to high levels of scattering induces an underestimation of the refraction angle even when the latter is not extremely high.

Unlike the nylon wires, the grooves phantom shows non-negligible absorption and may allow therefore for an analysis of efficiency of the different algorithms in separating the absorption component from the other contrast contributions. The calculated profiles of the absorption signal for the biggest groove (# 5) are compared in Fig. 6 with the theoretical profile. The experimental curves have been calculated by averaging 30 horizontal pixels in order to improve statistics. For the MIR and GCF algorithms, only the integrated absorption is considered. Theoretical absorption values have been calculated from the imaginary component β of the complex index of refraction [21]. The extracted profiles for absorption look more affected by noise and artefacts than the phase ones because of the very low contrast (5.4%) of the absorption signal. The image noise can be of statistical nature or due to spatial inhomogeneities and temporal instability of the X-ray beam.

 

Fig. 6 Cross-sections of absorption images of groove 5 (depth 1.44 mm). Values extracted by the DEI, E-DEI, G-DEI, MIR and GCF algorithms are compared with the theoretical profile.

Download Full Size | PDF

The DEI algorithm is the one that exhibits the highest artefacts at the two edges, since the validity conditions on which it is based (very low refraction angle) are violated in these regions. A similar behavior and an additional asymmetry between the two edges appear also in the plots of the E-DEI and MIR algorithms. The reason has to be probably attributed to the fact that the actual positions of the crystal analyzer during acquisition were slightly different from the nominal ones. As a consequence, some systematic errors can be present in the calculation of the absorption. Apart from noise and artefacts, the agreement between the theory and the extracted values in the points where the refraction signal is absent (or low) is very good.

In Fig. 7 the profile of the USAXS signal through the papers, as calculated with the G-DEI, MIR and GCF algorithms, is reported. The plotted profiles are the result of the average over 100 vertical pixels. No theoretical value is displayed since no modeling of the paper structure has been attempted. The results for the MIR and the G-DEI algorithms are very similar while the GCF, especially for high levels of scattering, leads to sensibly higher values. These results suggest that in both the MIR and G-DEI methods (but not in the GCF) there is some saturation effect when the amount of scattering is high.

 

Fig. 7 USAXS profiles from a horizontal line across images in Fig. 3. The results from G-DEI algorithm, MIR algorithm and GCF algorithm are shown.

Download Full Size | PDF

4.2 Bone-cartilage sample

Extraction algorithms were also applied to ABI images of a human bone-cartilage sample. In Fig. 8 the maximum absorption (a), the integrated absorption (b), the refraction angle (c) and the USAXS (d) images obtained by using the MIR algorithm are reported as an example. In all the images the bone and cartilage parts of the sample can be distinguished.

 

Fig. 8 Images calculated by applying MIR algorithm to the human bone-cartilage sample: (a) maximum absorption, (b) integrated absorption, (c) refraction angle and (d) USAXS images.

Download Full Size | PDF

Because of the particular and different elemental composition and structure of the two tissues (bone and cartilage), the bone presents very intense signal for all contrast contributions (absorption, USAXS and refraction) with respect to cartilage.

Nevertheless, on the refraction image (Fig. 8 (c)) it is noteworthy the depiction of details inside the cartilage tissues that are instead completely invisible on the other calculated images. The refraction from these fine cartilage structures, which are compatible with the cells (chondrocytes) and collagen fibers present in the tissues, produces sufficient contrast and well represent the different architecture of the tissue along its height (from the bone-cartilage interface up to the surface) [22].

Horizontal profiles in the cartilage and in the bone were obtained from the refraction angle images calculated using the DEI, MIR and GCF algorithms and are reported in Figs. 9 and 10 . Two important characteristics are clearly visible in the profile of the cartilage: a fine amplitude oscillation at high frequency and a low frequency decrease from the left to the right side. The fine amplitude oscillations likely represent the signal given by the chondrocytes cells embedded in the cartilage matrix. The approximate estimated dimensions of these structures from the profile, about 35-50 μm, are in fact compatible with the common dimensions of chondrocytes. The overall refraction angle variation visible in the plot is instead given by the cylindrical shape of the sample.

 

Fig. 9 Comparison of cartilage refraction angle values extracted by DEI, MIR and GCF algorithms. Profiles position corresponds to line 1 in Fig. 8.

Download Full Size | PDF

 

Fig. 10 Comparison of bone refraction angle values extracted by DEI, MIR and GCF algorithms. Profiles position corresponds to line 2 in Fig. 8.

Download Full Size | PDF

A careful comparison of the signal amplitude obtained with the different algorithms shows that the MIR and GCF lead to very similar values, while with DEI algorithm the refraction angles are underestimated, especially at the borders (see the left and right sides of the cartilage profiles). This is in agreement with the results obtained for the plastics phantoms, and is due to the fact that the linearization of the RC introduced in the DEI method is accurate only if the refraction angles are very small, as already discussed.

The horizontal profile in the bone shows a very different scenario. The refraction angles, here originating by the trabecular structure of the bone, are characterized by very high amplitude compared with angles produced in the cartilage. In particular, values extracted using the GCF algorithm are in some cases bigger than 7 μrad, which is considerably higher than the FWHM of the experimental RC (about 2.94 μrad). It is clear that in this situation the validity conditions of the DEI algorithm are strongly violated and the algorithm fails. Also the MIR algorithm underestimates the refraction signal with respect to the GCF algorithm. As seen theoretically in section 2.5 and experimentally on phantoms in section 4.1, this is due to the inaccuracy of the estimation of the actual RC centre when calculating the RC centroid. The error increases when the refraction angles are high compared to the angles sampling range.

5. Conclusions

Results show that in the presence of low refraction and low USAXS, the five algorithms studied here are all in good agreement with the theory. When the refraction angles are larger, the discrepancies among the methods also increase. For refraction angles of the order of the RC FWHM, only the GCF algorithm is in good agreement with the theoretical estimations. This effect is still more pronounced when, besides high refraction angles, considerable scattering is also present: the algorithms (DEI and E-DEI) that do not take USAXS into account show the worst agreement with the theory, while the GCF method gives the most accurate results.

The calculated absorption images show strong artefacts at the edges of samples, where the refraction angles are high and thus the validity conditions of the GO approximation are violated. The agreement with the theoretical values is instead good in the regions where the refraction angles are less important.

The USAXS signal from the newspaper layers has been extracted for the G-DEI, MIR and GCF methods. Both G-DEI and MIR methods underestimate the signal compared to the GCF algorithm, especially when the USAXS contribution is important. This fact suggests a sort of saturation effect similar to the one encountered in the refraction angle calculation.

It has been shown that the calculated images in the case of biological samples are able to provide important information concerning tissues structure. In the specific case of the sample used in this work, the bone and cartilage are clearly discriminated and the internal structure of cartilage is visible as well. This may represent an important result for an improved diagnosis of diseases affecting the cartilage, like osteoarthritis.

The refraction angle images show that the differences between the algorithms are reduced when considering the cartilage tissue, where the refraction angles are relatively small. The discrepancies are important on the bone tissue, where the refraction angles are very high compared to the RC FWHM. The signal provided by the GCF algorithm is here much higher than that given by the DEI and MIR algorithms. These results are in agreement with the theoretical considerations and with the results obtained for the plastics phantoms.

The analysis and the comparison of the different algorithms, under variable conditions of refraction and scattering, clearly identify the GCF method as the one having the best agreement with the theory. This result is confirmed also when refraction and scattering contributions are very important.

A drawback of this method is that it is computationally intensive compared to the other ABI extraction methods, since a fitting of Eq. (3) has to be performed for each image pixel. It can be estimated, anyway, that the large computational time needed by our routine could be sensibly reduced if parallel processing on multi-core CPUs, computer clusters or GPUs are used. This is an important point, especially in the view of the application of this method to tomographic images that will be discussed in a work in preparation.

The high accuracy of the GCF algorithm in calculating the refraction signal with respect to the other considered extraction methods (also in the presence of high scattering) makes this algorithm of great interest for different kinds of applications especially where the quantitative accuracy of the results is essential. The use of the GCF method for the separation of the different contrast contributions (refraction, absorption and scattering) from phase contrast AB images of clinical samples may help in the detection of those structures and pathologies that are otherwise difficult to visualize by improving the obtainable image contrast.