The present invention relates to a method and a system for quantitative X-ray ra- diology using the absorption and the forward scattering information.

In X-ray imaging, the absorption contrast is well known to follow the Beer- Lambert law. For monochromatic beam, we have

where / and I ₀ are the intensity signal on the detector and the incident X-rays intensity, respectively. / is the direction along the path of the X-ray beam and μ is the linear attenuation coefficient of the underlying object. Practically, in medical imaging, the beam generated by X-ray tube is always rather polychromatic. In that case, the linear attenuation coefficient will depend on the energy of the photons. Eq. (1 ) becomes

where Ε(ε) is the spectrum of the X-rays, and 77(e) is the energy response of the detector at photon energy ε . Eq. (2) is no longer a simple linear integral, namely the Beer-Lambert law. However, in the range of certain applications, for instance mammography, it has been shown that Eq. (1 ) approximately holds for polychromatic beam by theory as well as experiments [1]. Therefore, we have

where μ _4ί is defined as the effective attenuation coefficient at the mean energy of ε . For single material, Eq. (3) can be simplified to

Τ = μ _4ί {ε ) - Σ , (4) where L is the thickness of the sample along the X-ray beam direction.

Further, any in-homogeneities in the micro- or nanometer range within the sample will cause X-ray photons to scatter. This forward scattering concentrates in very small angles, and therefore it is called small-angle scattering. Measurements of the local small-angle scattering power can deliver important structural information about the sample, and are thus widely used in the characterization of materials. It has been shown that the small-angle scattering signal, obtained for instance by a grating interferometer or diffraction enhanced imaging, follows a similar way as the absorption signal for monochromatic beam [2],

where a ² s the second moment of the scattering angle distribution of the outgoing beam and s(s,l) \s the generalized scattering parameter which is equiva- lent to the attenuation coefficient in the absorption case [3].

Usually the Gaussian scattering approximation is used to model the outgoing beam. The scattering angle distribution p(Q ) s modeled by

For polychromatic beam, the total scattering angular distribution ?(θ ) will be the incoherent addition of all the p(Q ) at different photon energies:

ρ(θ)

where ε ₂ ] is the integration interval of the spectrum Ε(ε) . If is(8) is assumed to be a continuous function and normalized on the integration interval, then the mean value theorem of integrals indicates that

where ε is a constant within the range [ε _χ ε ₂ ] . Eq. (8) shows that the total scattering angular distribution ?(0 ) for polychromatic beam is still a Gaussian distribution. For single material, the second moment of ?(θ ) is

= Γ E(s)s(s)Lds

(9)

= s _eff (s ) - L

Here the mean value theorem of integrals is used again as well as the fact that the scattering angular distributions of each single energy are independent. s _e (g ) is defined as the effective scattering parameter of the materials at the mean energy ε .

The linear relationship of Eq. (4) and Eq. (9) hold in a wider sense as long as the X-ray spectrum function is(8) is generally continuous and the beam hardening effect is not severe.

In conventional absorption-based radiographic methods, quantitative imaging (for instance, discriminating materials or determining the material composition) is not possible without pre-calibrations. The difficulty lies on the fact that Eq. (4) has two unknown parameters: the effective attenuation coefficient and the thickness of the sample. The same absorption value could represent arbitrary combination of the thickness and attenuation coefficient.

It is therefore the objective of the present invention to provide a system and a method for quantitative X-ray radiology eliminating the drawback given above. It would be desirable to improve the prior art in a way that the same absorption value might be analyzed in order to identify whether the absorption information has been caused rather by the present thickness of the probe or the attenuation coefficient of the probe which allows enables the evaluation whether one material or diverse materials have been penetrated by the X-ray beam.

This objective achieved according to the present invention by the features given in the independent claims 1 for the method and 17 for the system. The way according to the present invention to overcome this problem and get quantitative information of the sample is to involve another physical quantity in order to decouple the thickness parameter. Ideally, this quantity should be obtained with the absorption information simultaneously. This requirement will max- imally reduce the errors due to imaging the specimen in different circumstances. A multiple-modality imaging system is a potential solution, such as the grating interferometer. A suitable physical quantity in this sense is the forward scattering information used simultaneously with the absorption information of the x-ray analysis.

According to the present invention, in conditions where the beam hardening effect is not severe, for instance, in medical imaging applications like mammography, Eq. (4) and Eq. (9) are both linear with the thickness L , therefore their "ratio image", which is defined as the R image here, cancels out the unknown parame- terZ , leading to Eq. (10):

_R _ T _ μ _// (ε )∑ _ μ (ε ) _ _{( 1 Q)} S s _eff (s )L s _eff (s )

Since the l value doesn't depend on the thickness of the sample, naturally it can be considered as the signature of the material of the sample in the radiographic imaging. The same material is expected to show the same l value regardless from the thickness.

It is worth mentioning that the l value is not expected to be as unique as the atomic number or electron density for material discrimination purpose. For compressible samples, e.g. the glandular and fat tissue of breast in medical imaging, both their attenuation coefficient and the generalized scattering parameter might change when different compressions are applied. Therefore, the l value may also change accordingly. However, despite of these small limitations, this concept of the thickness irrelevance in the R value is useful in most medical applications. Preferred embodiments of the present invention are given in the dependent claims 2 to 16 with respect to the method and 18 to 24 with respect to the system. Preferred embodiments of the present invention are described hereinafter with reference to the attached drawings which depict in: Figure 1 results from validation experiments with the use of five different plastic materials;

Figure 2 results on a quantitative radiography; Figure 3 an R image for one breast sample in the craniocaudal (CC) view (a) and the corresponding histogram of the R image (b);

Figure 4 a schematic illustration of the dynamic process for determining

^ _g and R _F ;

Figure 5 the mean and standard deviation of VBD values estimated from the

CC view for each diagnostic ACR group; and

Figure 6 a two-gratings interferometer (top) and a three-gratings interferome- ter (bottom) for X-ray imaging.

With respect to the equations given above for the calculation of the R image, the justification of Eq. (4) and Eq. (9) by experiments on several known plastic materials is given in Fig. 1 . The samples were imaged using a grating interferometer with a broad spectrum like one of the grating interferometers shown in Figure 6. Figure 1 shows the results of the validation experiments on five plastic materials. Figure 1 (a) shows the phantom design; (b) is the absorption image and (c) is the scattering image: Figure 1 (d) shows the resulting R image. In Figure 1 (e) the profiles of the absorption image for the five materials are illustrated. Forty hori- zontal lines were used for each material (the ROI was indicated in Figure 1 (c)) and averaged along the vertical direction. Figure 1 (f) are the profiles of the scattering image for the five materials, forty horizontal lines were used for each material and averaged along the vertical direction. Figure 1 (g) are the profiles of the resulting R images for each material. Noticeable noises were observed at the places where the samples were too thin, especially for Nylon and Derlin. This was because the scattering signals were close to the noise level at those places. Figure 1 (g) is therefore truncated from the thickness of 1 .2 cm as indicated by the right rectangle in Figure 1 (d). In Figure 1 (h) the tabulated R value for the five given materials is listed.

It can be seen that the l values of the five materials are independent of the thickness within a good approximation. And the materials can be easily distin- guished by their l values which confirms exactly the general expectations on this techniques.

An application of using R image for quantitative radiology is given in Fig. 2. Conventional X-ray absorption-based radiography can only provide qualitative infor- mation (e.g. morphology) of the underlying sample, but not quantitative information, such as telling the composition of the sample. A 4cmx3cm wedge phantom and a 3cmx3cm one of the same material were overlapped with each other and imaged using a grating interferometer set up as shown in Figure 6. The results were shown in Fig.2. For the absorption image only (Fig. 2(a)), different ab- sorption properties were observed. However, there was not enough information to tell whether it was a single material with different thicknesses or just two different materials. By taking the R image of the same sample, it was clearly shown that the phantom was composed of one material (Fig. 2(b)). Using the calibrated data in Fig.1 (h), additionally, the unknown material could be identified, which was Teflon in the present case.

Using the R image, additional quantitative information can be obtained compared to conventional absorption-based radiographic methods. A particular application in medical imaging is given below regarding volumetric breast density (VBD) es- timation.

Breast density has gained increasing attentions in breast cancer screening and diagnosis because it is a strong indicator of the breast cancer risk [4, 5]. Conven- tional 2D breast density estimation methods have the problems of subjectivity and low inter-reader agreements [6]. Volumetric breast density (VBD) determines the breast density in a 3D manner, which provides more accurate result [7]. The R image can be used to achieve quantitative volumetric breast density (VBD) estimation in mammography.

On the pixel base, the VBD is defined by (volumetric percentage of the glandular tissue)

where L _g and L _f are the thicknesses of the glandular tissue and fat tissue along the X-ray beam, respectively. is the total compression thickness of the breast. The volumetric breast density of the whole breast is given by averaging m over the whole breast region

VBD = ^ ,

N

where N is the pixel numbers of the breast region.

With the two-material composition assumption of the breast, Eq. (10) is

T _ _g L _g + _f L _f _ μ τηΐ + μ _/ (1 - m)L

(12) S s L + s _f L _f s mL + s _f (l - m)L

Taking Eq. (1 1 ) into Eq. (12), the connection between the VBD m and the measured quantity R is setup by m R < R < R f (13)

where R =— and R _F =— are defined as the l values for 100% glandular tis- s* ^s f

sue and 100% fat tissue, respectively. The subscripts g and / represent the glandular tissue and fat tissue, respectively. R _G and R _f can be determined dynamically from the histogram of the R image as described in the following sec- tions. Note that the fat tissue has a bigger R value than the glandular tissue. μ

— is a constant for a certain energy and can be determined by calibrated data μ

[7]. For instance,—is around 1 .43 for a mean photon energy of 28 KeV.

A typical histogram of the R image for the breast sample is shown in Fig. 3. Fig- ure 3(a) shows the R image for one breast sample in the craniocaudal (CC) view, Figure 3(b) is the corresponding histogram of the R image. The histogram H(R) \S usually an asymmetric Gaussian distribution due to the fact that the fat tissue (the right part of the histogram) is dominant in the breast compared to the glandular tissue (the left part) [8], therefore has a larger variance. The VBD information of the breast is actually implied in the asymmetric distribution of the histogram. For instance, the skewness (the third moment) of the distribution, which is defined as can give a good estimation of how dominant the glandular or fat tissue is in the whole breast.

To get a quantitative value of the breast density from Eq. (12), we determine R _G and R _F in the following way with the consideration that noises always present in the R image and the dynamic process for determining R _G and R _F . is illustrated in Fig.4.:

(1 ) The peak position R _PEAK which gives the maximal H(R) is first decided.

(2) The histogram H(R) is divided into the left and right parts by R _PEAK . These two parts of the histogram are considered separately. For each half, the dis- tribution is mirrored around the peak value to form a symmetric Gaussian distribution. Then Gaussian fitting is applied to the two resulting Gaussian distributions and their standard deviations a _g and o _f are calculated respectively.

(3) R _g and R _f are decided using the Νσ criteria,

R _f = R _pmk + Na _f '

where N is a positive integer. N = 3 is considered to give the best approximation of the i? _g and R _f when the noises present. The 3σ criteria covers 99.7% of the area of the histogram. With the determined R _g and R _f , the 3D VBD of the whole breast is calculated by weighting the normalized histogram H(R) using Eq. (15),

VBD =∑m(R) - H(R) . (16)

The proposed method was evaluated by clinical mastectomy breast dataset. 27 patients were included in the evaluation with the diagnostic breast densities: ACR 1 (4 samples), ACR 2 (13 samples), ACR 3 (10 samples) and ACR 4 (0 sample). The VBD results for the CC view using the proposed method were given in Fig.5 which shows the mean and standard deviation of VBD values estimated from the CC view for each diagnostic ACR group. The calculated VBD values fell into the range of 10%~30% which is consistent with the previous studies using X-ray tomography and MRI [8]. A positive correlation was also clearly seen, which meant statistically the proposed method was consistent with the diagnostic results.

In this invention, the absorption and small-angle scattering signals obtained with multiple-modality systems are adopted to distinguish two types of microcalcifica- tions. Such multiple-modality systems have been developed in the last fifteen years, including techniques based on analyzer crystal [9], gratings [10, 1 1 , 12] and coded apertures [13, 14]. The described invention is in context with these techniques and uses the set-ups disclosed in the respective documents which form therefore part of this disclosure with respect to the set-ups. Practical applications in non-destructive testing and medical imaging demand techniques which can work well with conventional X-ray tubes. For this reason, grating-based methods are especially promising. Without loss of generality, we will discuss the practical aspects of the method using gratings-based interferome- try as an example.

Grating-based X-ray imaging setups can generate three different signals: the conventional absorption contrast (AC) signal, the differential phase contrast (DPC) signal caused by refraction due to phase shifts[1 1], and the small-angle scattering contrast (SC) signal (also named dark-field signal) caused by scattering from in-homogeneities in the sample [10].

Grating Interferometer set-ups with two gratings (G1 and G2) or three gratings (GO, G1 , and G2) as schematically shown in Figure 6 can be applied to record the deflection of the X-rays. In the case of a two-grating set-up, the source needs to fulfill certain requirements regarding its spatial coherence, wherein a GO grating is required, when the source size is bigger than p2 ^* l/d, where p2 is the period of G2, / is the distance between the source and G1 , and d is the distance between G1 and G2, while in a three grating setup no spatial coherence is required [15,16]. Therefore, the three grating set-up is suited for use with incoherent X-ray sources, in particular with standard X-ray tubes for medical investigations.

To separate the conventional attenuation contrast (AC) from the DPC and SC contrast, a phase-stepping approach is applied. One of the gratings is displaced transversely to the incident beam whilst acquiring multiple images. The intensity signal at each pixel in the detector plane oscillates as a function of the displacement. The average value of the oscillation represents the AC. The phase of the oscillation can be directly linked to the wave-front phase profile and thus to the DPC signal. The amplitude of the oscillation depends on the scattering of X-rays in the object and thus yields the SC signal.

For the (two or three) gratings, several approaches have been proposed and applied. The grating GO (if required) is the one closest to the source. It usually con- sists of a transmission grating of absorbing lines with the period pO. It can be replaced by a source that emits radiation only from lines with the same period [17]. The grating G1 is placed further downstream of the source. It consists of lines with a period p1. The grating G2 is the one most downstream of the setup. It usually consists of a transmission grating of absorbing lines with the period p2. It can be replaced by a detector system that has a grating-like sensitivity with the same period.

Two regimes of setups can be distinguished: in the so called "near field regime" and the "Talbot regime". In the "near field regime", the grating period p, grating distances d and the X-ray wavelength λ are chosen such, that diffraction effects are negligible. In this case, all gratings need to consist of absorbing lines. In the "Talbot regime", diffraction from the grating structures is significant. A sharp distinction between the two regimes is not easily given, as the exact criterion de- pends on the duty cycle of the grating structure, and whether the gratings are absorbing or phase shifting. E.g., for a grating with absorbing lines and a duty cycle of 0.5, the condition for the "near field regime" is d≥ ρ ² /2λ. Here G1 should consist of grating lines that are either absorbing or, preferentially, phase shifting. Several amounts of phase shift are possible, preferentially π/2 or multiples there- of. The grating periods must be matched to the relative distances between the gratings. In the case of setups in the "Talbot regime" the Talbot effect needs to be taken into account to obtain good contrast. The formulae for the grating periods and distances are described in [18]. The sample is mostly placed between GO and G1 (or upstream of G1 in the case of a two-grating set-up), however it can be advantageous to place it between G1 and G2 [19].

The presented invention is relevant in all of the abovementioned cases, i.e. in the two- and three-gratings case, in the case of the "nearfield regime" and the "Talbot regime", and for the sample placed upstream or downstream of G1. In addition, the invention presented here also works in combination with scanning-based systems as suggested in [20, 21] or for planar grating geometries, as suggested in [22]. Intensity curves (with and without sample) are usually obtained with "phase stepping" methods or alternative techniques. Defining for each pixel on the detector the mean, phase and visibility of the intensity curve with sample as / _s ,O _s ,F _s , and without sample as I _b , _b ,V _b , yields:

DPC = <& _s - <& _b (18)

SC = -log(¾ . (19)

For both the AC signal and SC signal, the valid data range is [0,+∞] , while for the DPC it is [-π, +π] . Images obtained by plotting such signals are all perfectly registered.

A similar way to generate these multiple information signals can be found in diffraction enhanced imaging where the equivalent of the intensity curve is named the rocking curve [9].

References

[I ] J. J. Heine and M. Behera, "Effective X-ray attenuation measurements with full field digital mammography," Med. Phys. 33, 4350-4366 (2006).

[2] Z-T Wang, K-J Kang, Z-F Huang, and Z-Q Chen, "Quantitative grating- based X-ray dark-field computed tomography," Appl. Phys. Lett. 95, 094105 (2009).

[3] G. Khelashvili, J. Brankov, D. Chapman, M. Anastasio, Y. Yang, Z. Zhong, and M. Wernick, "A physical model of multiple-image radiography," Phys. Med. Biol. 51 , 221 (2006).

[4] C. Byrne et al, "Mammographic features and breast cancer risk: effects with time, age, and menopause status," J. Nat. Cancer Inst. 87(21 ), 1622- 1629 (1995).

[5] C. H. van Gils et al, "Changes in mammographic breast density and concomitant changes in breast cancer risk," Eur. J. Cancer 8(6), 509

(1999).

[6] A. Redondo et al, "Inter- and intraradiologist variability in the BI-RADS assessment and breast density categories for screening mammograms, "

Br. J. Radiol. 85(1019) 1465, (2012).

[7] M. J. Yaffe, "Mammographic density. Measurement of mammographic density." Breast Cancer Res, 10(3), 209 (2008).

[8] M. J. Yaffe, et al. "The myth of the 50-50 breast," Medical Physics 36(12),

5437 (2009).

[9] D. Chapman et al, "Diffraction enhanced X-ray imaging", Physics in Medicine and Biology, 42, 2015 (1997).

[10] F. Pfeiffer, M. Bech, O. Bunk et al., "Hard-X-ray dark-field imaging using a grating interferometer," Nature Materials, 7(2), 134-137 (2008).

[I I ] T. Weitkamp, A. Diaz, C. David et al., "X-ray phase imaging with a grating interferometer," Optics Express, 13(16), 6296-6304 (2005).

[12] F. Pfeiffer, C. Kottler, O. Bunk et al., "Hard X-ray phase tomography with low-brilliance sources," Physical Review Letters, 98(10), - (2007). [13] A. Olivo and R. Speller, "A coded-aperture technique allowing X-ray phase contrast imaging with conventional sources", Appl. Phys. Lett. 91 , 074106 (2007).

[14] A. Olivo and R. Speller, "Image formation principles in coded-aperture based X-ray phase contrast imaging", Phys. Med. Biol. 53, 6461 (2008).

[15] F. Pfeiffer, T. Weitkamp, O. Bunk et al., "Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources," Nature Physics,

2(4), 258-261 (2006).

[16] C. David, T. Weitkamp, and F. Pfeiffer, [Interferometer for quantitative phase contrast imaging and tomography with an incoherent polychromatic

X-ray source].

[17] C. David, [Strukturierung der Anodenoberflache einer Rontgenrohre urn eine gitterformige Emission fur die Phasenkontrastbildgebung im

Rontgenbereich], (2008).

[18] T. Weitkamp, C. David, C. Kottler et al., "Tomography with grating interferometers at low-brilliance sources." Proc. of SPIE, 6318, 63180S

(2006).

[19] C. David, [Optimierte Anordnung von Gittern fur die

Phasenkontrastbildgebung im Rontgenbereich], Siemens

Erfindungsmeldung 2008E07939 DE (2008).

[20] C. Kottler, F. Pfeiffer, O. Bunk et al., "Grating interferometer based scanning setup for hard X-ray phase contrast imaging," Review of

Scientific Instruments, 78(4), 043710 (2007).

[21] C. David, and F. Pfeiffer, [An interferometer for X-rays for obtaining quantitative X-ray images from an object].

[22] C. David, and M. Stampanoni, [A method for X-ray phase contrast and dark-field imaging using an arrangement of gratings in planar geometry],

(2010).